diff --git a/atmos_model.F90 b/atmos_model.F90
index 051f5918d..c010789cd 100644
--- a/atmos_model.F90
+++ b/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
diff --git a/ccpp/config/ccpp_prebuild_config.py b/ccpp/config/ccpp_prebuild_config.py
index f649535ac..cfa0b5eb6 100755
--- a/ccpp/config/ccpp_prebuild_config.py
+++ b/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',
diff --git a/ccpp/physics b/ccpp/physics
index 0ac8068d5..acf281a01 160000
--- a/ccpp/physics
+++ b/ccpp/physics
@@ -1 +1 @@
-Subproject commit 0ac8068d53a19d342fc176dd62d4267f1b43008e
+Subproject commit acf281a01e19b840a1f0e0fb947d9672c6d10c05
diff --git a/ccpp/suites/suite_FV3_GSD_noah.xml b/ccpp/suites/suite_FV3_GSD_noah.xml
index 4d7d4e00f..42d55a5e4 100644
--- a/ccpp/suites/suite_FV3_GSD_noah.xml
+++ b/ccpp/suites/suite_FV3_GSD_noah.xml
@@ -79,7 +79,6 @@
mp_thompson_post
GFS_MP_generic_post
cu_gf_driver_post
-
maximum_hourly_diagnostics
diff --git a/ccpp/suites/suite_FV3_GSD_v0.xml b/ccpp/suites/suite_FV3_GSD_v0.xml
index 7838db77b..0d6531d19 100644
--- a/ccpp/suites/suite_FV3_GSD_v0.xml
+++ b/ccpp/suites/suite_FV3_GSD_v0.xml
@@ -78,7 +78,6 @@
mp_thompson_post
GFS_MP_generic_post
cu_gf_driver_post
-
maximum_hourly_diagnostics
diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90
index 2fef16ee4..1a63d5bc8 100644
--- a/gfsphysics/GFS_layer/GFS_typedefs.F90
+++ b/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
diff --git a/gfsphysics/GFS_layer/GFS_typedefs.meta b/gfsphysics/GFS_layer/GFS_typedefs.meta
index 4dfb5046e..dcacc8644 100644
--- a/gfsphysics/GFS_layer/GFS_typedefs.meta
+++ b/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
diff --git a/stochastic_physics/stochastic_physics_wrapper.F90 b/stochastic_physics/stochastic_physics_wrapper.F90
index 5a3701aa8..44c82ecbd 100644
--- a/stochastic_physics/stochastic_physics_wrapper.F90
+++ b/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