MOSRS Trunk r8917

components/buoyancy/src/buoyancy.F90

@@ -29,8 +29,6 @@ module buoyancy_mod
   real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: tend_pr_tot_w
   logical :: l_tend_pr_tot_w
 
-  integer :: diagnostic_generation_frequency
-
   public buoyancy_get_descriptor
 
 contains
@@ -170,9 +168,6 @@ subroutine initialisation_callback(current_state)
       allocate( tend_pr_tot_w(current_state%local_grid%size(Z_INDEX)) )
     endif
 
-    ! Save the sampling_frequency to force diagnostic calculation on select time steps
-    diagnostic_generation_frequency=options_get_integer(current_state%options_database, "sampling_frequency")
-
   end subroutine initialisation_callback  
 
 
@@ -197,6 +192,10 @@ subroutine timestep_callback(current_state)
 
     integer :: k, n
     integer :: current_x_index, current_y_index, target_x_index, target_y_index
+    logical :: calculate_diagnostics
+
+    calculate_diagnostics = current_state%diagnostic_sample_timestep &
+                            .and. .not. current_state%halo_column
 
     current_x_index=current_state%column_local_x
     current_y_index=current_state%column_local_y
@@ -210,10 +209,8 @@ subroutine timestep_callback(current_state)
       endif
     endif  ! zero totals
 
-    if (mod(current_state%timestep, diagnostic_generation_frequency) == 0 .and. .not. current_state%halo_column) then
-      call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
-    end if
-
+    if (calculate_diagnostics) &
+        call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
 
 #ifdef W_ACTIVE
     if (.not. current_state%passive_th .and. current_state%th%active) then
@@ -252,9 +249,8 @@ subroutine timestep_callback(current_state)
     end if
 #endif
 
-    if (mod(current_state%timestep, diagnostic_generation_frequency) == 0 .and. .not. current_state%halo_column) then
-      call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
-    end if
+    if (calculate_diagnostics) &
+        call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
 
   end subroutine timestep_callback
 

components/casim/src/casim.F90

@@ -15,6 +15,7 @@ module casim_mod
   ! casim modules...
   use variable_precision, ONLY: wp
   use initialize, only: mphys_init
+  use mphys_parameters, only: cloud_params
   use mphys_switches, only: set_mphys_switches,  &
      l_warm,                                     &
      nq_l, nq_r, nq_i, nq_s, nq_g,               &
@@ -74,7 +75,17 @@ module casim_mod
      , l_pssub & ! sublimation of snow
      , l_pgsub & ! sublimation of graupel
      , l_pisub & ! sublimation of ice
-     , l_pimlt   ! ice melting
+     , l_pimlt & ! ice melting
+     ! New switches for sedimentation (these are sort-of temporary)
+     , l_gamma_online & ! when true use standard vn0.3.3 sed, when false use precalced gamma
+     , l_subseds_maxv & ! Use a CFL criteria based on max terminal velocity
+                    ! and sed_1M_2M 
+     , l_sed_eulexp & ! switch for eulexp sed based on UM. Default is false
+                      ! so standard casim sed used
+     , cfl_vt_max & ! cfl limit for sedimentation (default = 1.0)
+     , l_kfsm
+  use mphys_constants, only: fixed_cloud_number
+
 
   use micro_main, only: shipway_microphysics
   use generic_diagnostic_variables, ONLY: casdiags, allocate_diagnostic_space, &
@@ -93,7 +104,7 @@ module casim_mod
      , nc(:,:,:), qr(:,:,:), nr(:,:,:), m3r(:,:,:),rho(:,:,:) &
      , exner(:,:,:), w(:,:,:), tke(:,:,:)                               &
      , qi(:,:,:), ni(:,:,:), qs(:,:,:), ns(:,:,:), m3s(:,:,:) &
-     , qg(:,:,:), ng(:,:,:), m3g(:,:,:) 
+     , qg(:,:,:), ng(:,:,:), m3g(:,:,:), cfliq(:,:,:), cfice(:,:,:)
 
   REAL(wp), allocatable :: AccumSolMass(:,:,:), AccumSolNumber(:,:,:) ! Accumulation mode aerosol
   REAL(wp), allocatable :: ActiveSolLiquid(:,:,:)                      ! Activated aerosol
@@ -277,6 +288,8 @@ subroutine initialisation_callback(current_state)
     allocate(qg(kte,1,1))
     allocate(ng(kte,1,1))
     allocate(m3g(kte,1,1))
+    allocate(cfliq(kte,1,1))
+    allocate(cfice(kte,1,1))
 
     allocate(AccumSolMass(kte,1,1))
     allocate(AccumSolNumber(kte,1,1))
@@ -581,6 +594,7 @@ subroutine timestep_callback(current_state)
       iqx = iql
       qc(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
       dqc(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
+      cfliq(:,1,1) = 1.0
     end IF
     IF (nq_r > 0)then
       iqx = iqr
@@ -608,6 +622,7 @@ subroutine timestep_callback(current_state)
       iqx = iqi
       qi(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
       dqi(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
+      cfice(:,1,1) = 1.0
     end IF
     IF (nq_s > 0)then
       iqx = iqs
@@ -712,7 +727,7 @@ subroutine timestep_callback(current_state)
        pressure, rho,                              &
        w, tke,                                     &
        z_half, z_centre,                           &
-       dz,                                         &
+       dz, cfliq, cfice,                           &
                                 ! in/out
        dqv, dqc, dqr, dnc, dnr, dm3r,              &
        dqi, dqs, dqg, dni, dns, dng, dm3s, dm3g,   &
@@ -838,7 +853,7 @@ subroutine timestep_callback(current_state)
     ! and surface
     ! snow rate (precip_s), which is the sum of ice, snow and graupel (See micromain.F90 in casim for
     ! calculation).
-    if (l_warm) then
+    if (l_warm .or. .not. casdiags % l_surface_snow ) then
        surface_precip(target_y_index,target_x_index) = &
             casdiags % SurfaceRainR(1,1)
     else
@@ -877,7 +892,7 @@ subroutine read_configuration(current_state)
     sp2             = options_get_real(current_state%options_database, 'sp2')
     sp3             = options_get_real(current_state%options_database, 'sp3')
     max_mu          = options_get_real(current_state%options_database, 'max_mu')
-    fix_mu          = options_get_real(current_state%options_database, 'fix_mu')
+    cloud_params%fix_mu = options_get_real(current_state%options_database, 'fix_mu')
 
     l_aaut          = options_get_logical(current_state%options_database, 'l_aaut')
     l_aacc          = options_get_logical(current_state%options_database, 'l_aacc')
@@ -947,6 +962,12 @@ subroutine read_configuration(current_state)
     l_pgsub         = options_get_logical(current_state%options_database, 'l_pgsub')
     l_pisub         = options_get_logical(current_state%options_database, 'l_pisub')
     l_pimlt         = options_get_logical(current_state%options_database, 'l_pimlt')
+    l_gamma_online  = options_get_logical(current_state%options_database, 'l_gamma_online')
+    l_subseds_maxv  = options_get_logical(current_state%options_database, 'l_subseds_maxv')
+    l_sed_eulexp  = options_get_logical(current_state%options_database, 'l_sed_eulexp')
+    cfl_vt_max      = options_get_real(current_state%options_database, 'cfl_vt_max')
+    l_kfsm          = options_get_logical(current_state%options_database, 'l_kfsm')
+    fixed_cloud_number = options_get_real(current_state%options_database, 'fixed_cloud_number')
 
   end subroutine read_configuration
 

components/casim/src/casim_monc_diagnostics/casim_monc_dgs_space.F90

@@ -22,6 +22,7 @@ module casim_monc_dgs_space
      !--------------------------------
      ! Surface Precipitation rates
      REAL, ALLOCATABLE :: precip(:,:)
+     REAL, ALLOCATABLE :: SurfaceCloudR(:,:)
      REAL, ALLOCATABLE :: SurfaceRainR(:,:)
      REAL, ALLOCATABLE :: SurfaceSnowR(:,:)
      REAL, ALLOCATABLE :: SurfaceGraupR(:,:)
@@ -98,6 +99,11 @@ subroutine allocate_casim_monc_dgs_space(current_state, casdiags)
        casim_monc_dgs % precip(:,:) = 0.0_DEFAULT_PRECISION
     endif
 
+    if ( casdiags % l_surface_cloud ) then
+       allocate ( casim_monc_dgs % SurfaceCloudR(y_size_local, x_size_local) )
+       casim_monc_dgs % SurfaceCloudR(:, :) = 0.0_DEFAULT_PRECISION
+    endif
+
     if ( casdiags % l_surface_rain ) then
        allocate ( casim_monc_dgs % SurfaceRainR(y_size_local, x_size_local) )
        casim_monc_dgs % SurfaceRainR(:, :) = 0.0_DEFAULT_PRECISION
@@ -314,6 +320,10 @@ subroutine populate_casim_monc_dg(current_state, casdiags )
          casim_monc_dgs % SurfaceRainR(target_y_index,target_x_index) = &
          casdiags % SurfaceRainR(1,1)
 
+    if ( casdiags % l_surface_cloud ) & 
+         casim_monc_dgs % SurfaceCloudR(target_y_index,target_x_index) = &
+         casdiags % SurfaceCloudR(1,1)
+
     if ( casdiags % l_pcond ) & 
          casim_monc_dgs % pcond(:,target_y_index,target_x_index) = &
          casdiags % pcond(1,1,:)
@@ -354,9 +364,12 @@ subroutine populate_casim_monc_dg(current_state, casdiags )
          casdiags % dqr(1,1,:)
 
     if (.not. l_warm) then
-       if ( casdiags % l_precip ) & 
-            casim_monc_dgs % precip(target_y_index,target_x_index) = &
+       if ( casdiags % l_precip .and. casdiags % l_surface_snow ) & 
+         casim_monc_dgs % precip(target_y_index,target_x_index) = &
             casdiags % SurfaceRainR(1,1) + casdiags % SurfaceSnowR(1,1)
+      if ( casdiags % l_precip .and. .not. casdiags % l_surface_snow ) & 
+         casim_monc_dgs % precip(target_y_index,target_x_index) = &
+            casdiags % SurfaceRainR(1,1)
        if ( casdiags % l_surface_snow ) & 
             casim_monc_dgs % SurfaceSnowR(target_y_index,target_x_index) = &
             casdiags % SurfaceSnowR(1,1)
@@ -365,7 +378,7 @@ subroutine populate_casim_monc_dg(current_state, casdiags )
             casdiags % SurfaceGraupR(1,1)
        if ( casdiags % l_phomc ) & 
             casim_monc_dgs % phomc(:,target_y_index,target_x_index) = &
-            casdiags % psedr(1,1,:)
+            casdiags % phomc(1,1,:)
        if ( casdiags % l_pinuc ) & 
             casim_monc_dgs % pinuc(:,target_y_index,target_x_index) = &
             casdiags % pinuc(1,1,:)

components/casim/src/um_modules_core/lsp_sedim_eulexp.F90

@@ -0,0 +1,96 @@
+! *****************************COPYRIGHT*******************************
+! (C) Crown copyright Met Office. All rights reserved.
+! For further details please refer to the file COPYRIGHT.txt
+! which you should have received as part of this distribution.
+! *****************************COPYRIGHT*******************************
+!  Microphysics hydrometeor Eulerian sedimentation scheme
+MODULE lsp_sedim_eulexp_mod
+
+IMPLICIT NONE
+
+CHARACTER(LEN=*), PARAMETER, PRIVATE :: ModuleName='LSP_SEDIM_EULEXP_MOD'
+
+CONTAINS
+
+SUBROUTINE lsp_sedim_eulexp(                                                   &
+  points,m0,dhi,dhir,rho,rhor,                                                 &
+  flux_fromabove, fallspeed_thislayer,                                         &
+  mixratio_thislayer, fallspeed_fromabove,                                     &
+  total_flux_out)
+
+!USE lsprec_mod, ONLY: zero, one
+
+! Use in KIND for large scale precip, used for compressed variables passed down
+! from here
+!USE um_types,             ONLY: real_lsprec
+
+use variable_precision, only: wp, iwp, defp
+
+USE yomhook, ONLY: lhook, dr_hook
+USE parkind1, ONLY: jprb, jpim
+IMPLICIT NONE
+
+! Description:
+!  Dummy replacement for lsp_sedim_eulexp, to permit compile of 
+!  CASIM sedimentation in MONC.
+
+! Method:
+!   Based on method described in Rotstayn (1997)(QJRMS, 123, 1227-1282)
+!
+! Code Owner: Please refer to the UM file CodeOwners.txt
+
+! Subroutine arguments
+
+      ! Intent (In)
+INTEGER :: points       ! number of points to process
+
+REAL (KIND=wp) ::                                                     &
+  m0,                                                                          &
+                         ! Small mass (kg/kg) defined in c_lspmic
+  dhi(points),                                                                 &
+                         ! CFL limit (s m-1)
+  dhir(points),                                                                &
+                         ! 1.0/DHI (m s-1)
+  rho(points),                                                                 &
+                         ! Air density (kg m-3)
+  rhor(points),                                                                &
+                         ! 1.0/Rho
+  flux_fromabove(points),                                                      &
+  fallspeed_thislayer(points)
+
+    ! Intent (InOut)
+REAL (KIND=wp) ::                                                     &
+  mixratio_thislayer(points),                                                  &
+  fallspeed_fromabove(points)
+
+    ! Intent (Out)
+REAL (KIND=wp) ::                                                     &
+  total_flux_out(points)
+
+! Local variables
+
+REAL (KIND=wp) ::                                                     &
+  mixratio_fromabove,                                                          &
+                                 ! Mixing Ratio from above
+  flux_out,                                                                    &
+                                 ! Temporary flux out of layer
+  expfactor                  ! Exponential Factor
+
+INTEGER :: i                    ! Loop counter
+
+INTEGER(KIND=jpim), PARAMETER :: zhook_in  = 0
+INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
+REAL(KIND=jprb)               :: zhook_handle
+
+CHARACTER(LEN=*), PARAMETER :: RoutineName='LSP_SEDIM_EULEXP'
+
+
+!-----------------------------------------------------------------------
+
+IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle)
+
+
+IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle)
+RETURN
+END SUBROUTINE lsp_sedim_eulexp
+END MODULE lsp_sedim_eulexp_mod

components/casim/src/um_modules_core/um_types.F90

@@ -0,0 +1,25 @@
+! *****************************COPYRIGHT*******************************
+! (C) Crown copyright Met Office. All rights reserved.
+! For further details please refer to the file COPYRIGHT.txt
+! which you should have received as part of this distribution.
+! *****************************COPYRIGHT*******************************
+
+! Dummy version of um_types so CASIM will build in both MONC and UM
+! lsp_sed_eulexp (from the UM) requires the precision to be set using 
+! real_lsprec. In the UM and rose-stem, this can be either single 
+! or double. The following code ensures CASIM will build in MONC with 
+! lsp_sed_eulexp call and single precision microphysics 
+! 
+
+MODULE um_types
+
+use variable_precision, only: wp
+
+IMPLICIT NONE
+
+!Large scale precipitation scheme
+INTEGER, PARAMETER :: real_lsprec = wp
+
+
+END MODULE um_types
+