Land stochastic perturbations (wrapper PR for #65)

atmos_model.F90

@@ -627,9 +627,9 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
                         IPD_Interstitial, commglobal, mpp_npes(), Init_parm)
 
 !--- Initialize stochastic physics pattern generation / cellular automata for first time step
-   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr)
-   if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
-
+!   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr)
+!   if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
+!
 #else
    call IPD_initialize (IPD_Control, IPD_Data, IPD_Diag, IPD_Restart, Init_parm)
 #endif
@@ -684,6 +684,9 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
    ! Initialize the CCPP physics
    call CCPP_step (step="physics_init", nblks=Atm_block%nblks, ierr=ierr)
    if (ierr/=0)  call mpp_error(FATAL, 'Call to CCPP physics_init step failed')
+!--- Initialize stochastic physics pattern generation / cellular automata for first time step
+   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr)
+   if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
 #endif
 
    !--- set the initial diagnostic timestamp

ccpp/config/ccpp_prebuild_config.py

@@ -159,7 +159,6 @@
     'ccpp/physics/physics/ozphys_2015.f',
     'ccpp/physics/physics/precpd.f',
     'ccpp/physics/physics/phys_tend.F90',
-    'ccpp/physics/physics/tracer_sanitizer.F90',
     'ccpp/physics/physics/radlw_main.F90',
     'ccpp/physics/physics/radsw_main.F90',
     'ccpp/physics/physics/rascnv.F90',

ccpp/suites/suite_FV3_GSD_noah.xml

@@ -79,7 +79,6 @@
       <scheme>mp_thompson_post</scheme>
       <scheme>GFS_MP_generic_post</scheme>
       <scheme>cu_gf_driver_post</scheme>
-      <!-- <scheme>tracer_sanitizer</scheme> -->
       <scheme>maximum_hourly_diagnostics</scheme>
     </subcycle>
   </group>

ccpp/suites/suite_FV3_GSD_v0.xml

@@ -78,7 +78,6 @@
       <scheme>mp_thompson_post</scheme>
       <scheme>GFS_MP_generic_post</scheme>
       <scheme>cu_gf_driver_post</scheme>
-      <!-- <scheme>tracer_sanitizer</scheme> -->
       <scheme>maximum_hourly_diagnostics</scheme>
     </subcycle>
   </group>

gfsphysics/GFS_layer/GFS_typedefs.F90

@@ -811,7 +811,9 @@ module GFS_typedefs
     integer              :: lsoil_lsm       !< number of soil layers internal to land surface model
     integer              :: lsnow_lsm       !< maximum number of snow layers internal to land surface model
     integer              :: lsnow_lsm_lbound!< lower bound for snow arrays, depending on lsnow_lsm
-    real(kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model
+    real(kind=kind_phys), pointer :: zs(:)    => null() !< depth of soil levels for land surface model
+    real(kind=kind_phys), pointer :: pores(:) => null() !< max soil moisture for a given soil type for land surface model
+    real(kind=kind_phys), pointer :: resid(:) => null() !< min soil moisture for a given soil type for land surface model
     logical              :: rdlai           !< read LAI from input file (for RUC LSM or NOAH LSM WRFv4)
     logical              :: ua_phys         !< flag for using University of Arizona? extension to NOAH LSM WRFv4
     logical              :: usemonalb       !< flag to read surface diffused shortwave albedo from input file for NOAH LSM WRFv4
@@ -1110,6 +1112,8 @@ module GFS_typedefs
     integer              :: skeb_npass
     integer              :: lndp_type
     integer              :: n_var_lndp
+    logical              :: lndp_each_step    ! flag to indicate that land perturbations are applied at every time step,
+                                              ! otherwise they are applied only after gcycle is run
     character(len=3)     :: lndp_var_list(6)  ! dimension here must match  n_var_max_lndp in  stochy_nml_def
     real(kind=kind_phys) :: lndp_prt_list(6)  ! dimension here must match  n_var_max_lndp in  stochy_nml_def 
                                               ! also previous code had dimension 5 for each pert, to allow 
@@ -3469,8 +3473,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     logical :: do_shum      = .false.
     logical :: do_skeb      = .false.
     integer :: skeb_npass   = 11
-    integer :: lndp_type    = 0 
-    integer :: n_var_lndp   =  0 
+    integer :: lndp_type      = 0
+    integer :: n_var_lndp     = 0
+    logical :: lndp_each_step = .false.
 
 !--- aerosol scavenging factors
     character(len=20) :: fscav_aero(20) = 'default'
@@ -3555,7 +3560,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                                do_deep, jcap,                                               &
                                cs_parm, flgmin, cgwf, ccwf, cdmbgwd, sup, ctei_rm, crtrh,   &
                                dlqf, rbcr, shoc_parm, psauras, prauras, wminras,            &
-                               do_sppt, do_shum, do_skeb, lndp_type,  n_var_lndp,           & 
+                               do_sppt, do_shum, do_skeb,                                   &
+                               lndp_type,  n_var_lndp, lndp_each_step,                      &
                           !--- Rayleigh friction
                                prslrd0, ral_ts,  ldiag_ugwp, do_ugwp, do_tofd,              &
                           ! --- Ferrier-Aligo
@@ -3981,6 +3987,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     end if
     ! Set number of ice model layers
     Model%kice      = kice
+
+    ! Allocate variable for min/max soil moisture for a given soil type
+    allocate (Model%pores(30))
+    allocate (Model%resid(30))
+    Model%pores    = clear_val
+    Model%resid    = clear_val
     !
     if (lsnow_lsm /= 3) then
       write(0,*) 'Logic error: NoahMP expects the maximum number of snow layers to be exactly 3 (see sfc_noahmp_drv.f)'
@@ -4216,6 +4228,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%do_skeb          = do_skeb
     Model%lndp_type        = lndp_type
     Model%n_var_lndp       = n_var_lndp
+    Model%lndp_each_step   = lndp_each_step
 
     !--- cellular automata options
     Model%nca              = nca
@@ -5159,6 +5172,8 @@ subroutine control_print(Model)
       print *, ' aoasis            : ', Model%aoasis
       print *, ' fasdas            : ', Model%fasdas
       print *, ' kice              : ', Model%kice
+      print *, ' shape(pores)      : ', shape(Model%pores)
+      print *, ' shape(resid)      : ', shape(Model%resid)
 #endif
       print *, ' ivegsrc           : ', Model%ivegsrc
       print *, ' isot              : ', Model%isot
@@ -5316,7 +5331,8 @@ subroutine control_print(Model)
       print *, ' do_shum           : ', Model%do_shum
       print *, ' do_skeb           : ', Model%do_skeb
       print *, ' lndp_type         : ', Model%lndp_type
-      print *, ' n_var_lndp         : ', Model%n_var_lndp
+      print *, ' n_var_lndp        : ', Model%n_var_lndp
+      print *, ' lndp_each_step    : ', Model%lndp_each_step
       print *, ' '
       print *, 'cellular automata'
       print *, ' nca               : ', Model%nca

gfsphysics/GFS_layer/GFS_typedefs.meta

@@ -3285,6 +3285,20 @@
   type = real
   kind = kind_phys
   active = (flag_for_land_surface_scheme == flag_for_ruc_land_surface_scheme) 
+[pores]
+  standard_name = maximum_soil_moisture_content_for_land_surface_model
+  long_name = maximum soil moisture for a given soil type for land surface model
+  units = m
+  dimensions = (30)
+  type = real
+  kind = kind_phys
+[resid]
+  standard_name = minimum_soil_moisture_content_for_land_surface_model
+  long_name = minimum soil moisture for a given soil type for land surface model
+  units = m
+  dimensions = (30)
+  type = real
+  kind = kind_phys
 [rdlai]
   standard_name = flag_for_reading_leaf_area_index_from_input
   long_name = flag for reading leaf area index from initial conditions

stochastic_physics/stochastic_physics_wrapper.F90

@@ -13,10 +13,23 @@ module stochastic_physics_wrapper_mod
   real(kind=kind_phys), dimension(:,:,:), allocatable, save :: skebv_wts
   real(kind=kind_phys), dimension(:,:,:), allocatable, save :: sfc_wts
 
+  integer, save :: lsoil = -999
   real(kind=kind_phys), dimension(:,:,:), allocatable, save :: smc
   real(kind=kind_phys), dimension(:,:,:), allocatable, save :: stc
   real(kind=kind_phys), dimension(:,:,:), allocatable, save :: slc
+  !
   real(kind=kind_phys), dimension(:,:), allocatable, save :: vfrac
+  !albedo
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: snoalb
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: alvsf
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: alnsf
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: alvwf
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: alnwf
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: facsf
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: facwf
+  !emissivity
+  real(kind=kind_phys), dimension(:,:), allocatable, save :: semis
+
   real(kind=kind_phys), dimension(:,:), allocatable, save :: stype
 
   ! For cellular automata
@@ -58,7 +71,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
     use cellular_automata_global_mod, only: cellular_automata_global
     use cellular_automata_sgs_mod,    only: cellular_automata_sgs
     use lndp_apply_perts_mod, only: lndp_apply_perts
-    use namelist_soilveg, only: maxsmc
 
     implicit none
 
@@ -108,11 +120,24 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
          allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%n_var_lndp))
       end if
       if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
-          allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
-          allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
-          allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
+          if (GFS_Control%lsm == GFS_Control%lsm_noah) then
+            lsoil = GFS_Control%lsoil
+          elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
+            lsoil = GFS_Control%lsoil_lsm
+          endif
+          allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
+          allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
+          allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
           allocate(stype(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
           allocate(vfrac(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(snoalb(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(alvsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(alnsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(alvwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(alnwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(facsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(facwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
+          allocate(semis(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
       endif
 
       do nb=1,Atm_block%nblks
@@ -171,31 +196,76 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
 
              do nb=1,Atm_block%nblks
                 stype(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%stype(:)
-                smc(nb,1:GFS_Control%blksz(nb),:)  = GFS_Data(nb)%Sfcprop%smc(:,:)
-                slc(nb,1:GFS_Control%blksz(nb),:)  = GFS_Data(nb)%Sfcprop%slc(:,:)
-                stc(nb,1:GFS_Control%blksz(nb),:)  = GFS_Data(nb)%Sfcprop%stc(:,:)
                 vfrac(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%vfrac(:)
+                snoalb(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%snoalb(:)
+                alvsf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%alvsf(:)
+                alnsf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%alnsf(:)
+                alvwf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%alvwf(:)
+                alnwf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%alnwf(:)
+                facsf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%facsf(:)
+                facwf(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Sfcprop%facwf(:)
+                semis(nb,1:GFS_Control%blksz(nb))  = GFS_Data(nb)%Radtend%semis(:)
              end do
 
-             ! determine whether land paramaters have been over-written
-             if (mod(GFS_Control%kdt,GFS_Control%nscyc) == 1)  then ! logic copied from GFS_driver
-                    param_update_flag = .true.
+             if (GFS_Control%lsm == GFS_Control%lsm_noah) then
+               do nb=1,Atm_block%nblks
+                 smc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%smc(:,:)
+                 slc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%slc(:,:)
+                 stc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%stc(:,:)
+               end do
+             elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
+               do nb=1,Atm_block%nblks
+                 smc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%smois(:,:)
+                 slc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%sh2o(:,:)
+                 stc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%tslb(:,:)
+               end do
+             endif
+
+             ! determine whether land paramaters have been over-written to
+             ! trigger applying perturbations (logic copied from GFS_driver),
+             ! or if perturbations should be applied at every time step
+             if (mod(GFS_Control%kdt,GFS_Control%nscyc) == 1 ) then
+               param_update_flag = .true.
              else
-                    param_update_flag = .false.
+               param_update_flag = .false.
              endif
-             call lndp_apply_perts( GFS_Control%blksz, GFS_Control%lsm,  GFS_Control%lsoil, GFS_Control%dtf, &
-                             GFS_Control%n_var_lndp, GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, &
-                             sfc_wts, xlon, xlat, stype, maxsmc,param_update_flag, smc, slc,stc, vfrac, ierr)
+
+             call lndp_apply_perts(GFS_Control%blksz, GFS_Control%lsm, GFS_Control%lsm_noah, GFS_Control%lsm_ruc, lsoil,      &
+                               GFS_Control%dtf, GFS_Control%kdt, GFS_Control%lndp_each_step,                                  &
+                               GFS_Control%n_var_lndp, GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list,                  &
+                               sfc_wts, xlon, xlat, stype, GFS_Control%pores, GFS_Control%resid,param_update_flag,            &
+                               smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, snoalb, semis, ierr)
              if (ierr/=0)  then
                     write(6,*) 'call to GFS_apply_lndp failed'
                     return
              endif
+
              do nb=1,Atm_block%nblks
-                 GFS_Data(nb)%Sfcprop%smc(:,:) =  smc(nb,1:GFS_Control%blksz(nb),:)
-                 GFS_Data(nb)%Sfcprop%slc(:,:) =  slc(nb,1:GFS_Control%blksz(nb),:)
-                 GFS_Data(nb)%Sfcprop%stc(:,:) =  stc(nb,1:GFS_Control%blksz(nb),:)
-                 GFS_Data(nb)%Sfcprop%vfrac(:) =  vfrac(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%vfrac(:)  = vfrac(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%snoalb(:) = snoalb(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%alvsf(:)  = alvsf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%alnsf(:)  = alnsf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%alvwf(:)  = alvwf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%alnwf(:)  = alnwf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%facsf(:)  = facsf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Sfcprop%facwf(:)  = facwf(nb,1:GFS_Control%blksz(nb))
+               GFS_Data(nb)%Radtend%semis(:)  = semis(nb,1:GFS_Control%blksz(nb))
              enddo
+
+             if (GFS_Control%lsm == GFS_Control%lsm_noah) then
+               do nb=1,Atm_block%nblks
+                   GFS_Data(nb)%Sfcprop%smc(:,:) = smc(nb,1:GFS_Control%blksz(nb),:)
+                   GFS_Data(nb)%Sfcprop%slc(:,:) = slc(nb,1:GFS_Control%blksz(nb),:)
+                   GFS_Data(nb)%Sfcprop%stc(:,:) = stc(nb,1:GFS_Control%blksz(nb),:)
+               enddo
+             elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
+               do nb=1,Atm_block%nblks
+                   GFS_Data(nb)%Sfcprop%smois(:,:) = smc(nb,1:GFS_Control%blksz(nb),:)
+                   GFS_Data(nb)%Sfcprop%sh2o(:,:)  = slc(nb,1:GFS_Control%blksz(nb),:)
+                   GFS_Data(nb)%Sfcprop%tslb(:,:)  = stc(nb,1:GFS_Control%blksz(nb),:)
+               enddo
+             endif
+
          endif ! lndp block
       end if
 
@@ -313,6 +383,7 @@ subroutine stochastic_physics_wrapper_end (GFS_Control)
          if (allocated(skebv_wts)) deallocate(skebv_wts)
       end if
       if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast
+         lsoil = -999
          if (allocated(sfc_wts)) deallocate(sfc_wts)
       end if
       if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
@@ -321,6 +392,14 @@ subroutine stochastic_physics_wrapper_end (GFS_Control)
           if (allocated(stc)) deallocate(stc)
           if (allocated(stype)) deallocate(stype)
           if (allocated(vfrac)) deallocate(vfrac)
+          if (allocated(snoalb)) deallocate(snoalb)
+          if (allocated(alvsf)) deallocate(alvsf)
+          if (allocated(alnsf)) deallocate(alnsf)
+          if (allocated(alvwf)) deallocate(alvwf)
+          if (allocated(alnwf)) deallocate(alnwf)
+          if (allocated(facsf)) deallocate(facsf)
+          if (allocated(facwf)) deallocate(facwf)
+          if (allocated(semis)) deallocate(semis)
       endif
       call finalize_stochastic_physics()
    endif