Updated land perturbation scheme

.gitmodules

@@ -8,5 +8,5 @@
 	branch = master
 [submodule "ccpp/physics"]
 	path = ccpp/physics
-	url =  https://github.com/NCAR/ccpp-physics.git 
-	branch = master 
+	url = https://github.com/NCAR/ccpp-physics
+	branch = master

atmos_model.F90

@@ -297,7 +297,9 @@ subroutine update_atmos_radiation_physics (Atmos)
 #endif
 
 !--- call stochastic physics pattern generation / cellular automata
-    call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block)
+    call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr) 
+    if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
+
 
 !--- if coupled, assign coupled fields
 
@@ -628,7 +630,8 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
 #endif
 
 !--- Initialize stochastic physics pattern generation / cellular automata for first time step
-   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block)
+   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr)
+   if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
 
    Atmos%Diag => IPD_Diag
 

gfsphysics/CMakeLists.txt

@@ -17,7 +17,6 @@ endif()
 set(CCPP_SOURCES
     physics/mersenne_twister.f
     physics/namelist_soilveg.f
-    physics/physparam.f
     physics/set_soilveg.f
 
     physics/noahmp_tables.f90

gfsphysics/GFS_layer/GFS_driver.F90

@@ -633,6 +633,7 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, &
       if (mod(Model%kdt,Model%nscyc) == 1) THEN
         call gcycle (nblks, Model, Grid(:), Sfcprop(:), Cldprop(:))
       endif
+        ! if not updating surface params through fcast, perturb params once at start of fcast
     endif
 
     !--- determine if diagnostics buckets need to be cleared

gfsphysics/GFS_layer/GFS_physics_driver.F90

@@ -660,6 +660,7 @@ subroutine GFS_physics_driver                                 &
       real    :: pshltr,QCQ,rh02
       real(kind=kind_phys), allocatable, dimension(:,:) :: den
 
+      real(kind=kind_phys) :: lndp_vgf  
       !! Initialize local variables (for debugging purposes only,
       !! because the corresponding variables Interstitial(nt)%...
       !! are reset to zero every time).
@@ -928,34 +929,28 @@ subroutine GFS_physics_driver                                 &
 !        alb1d(i)  = zero
          vegf1d(i) = zero
       enddo
-      if (Model%do_sfcperts) then
-        if (Model%pertz0(1) > zero) then
-          z01d(:) = Model%pertz0(1) * Coupling%sfc_wts(:,1)
-!          if (me == 0) print*,'Coupling%sfc_wts(:,1) min and max',minval(Coupling%sfc_wts(:,1)),maxval(Coupling%sfc_wts(:,1))
-!          if (me == 0) print*,'z01d min and max ',minval(z01d),maxval(z01d)
-        endif
-        if (Model%pertzt(1) > zero) then
-          zt1d(:) = Model%pertzt(1) * Coupling%sfc_wts(:,2)
-        endif
-        if (Model%pertshc(1) > zero) then
-          bexp1d(:) = Model%pertshc(1) * Coupling%sfc_wts(:,3)
-        endif
-        if (Model%pertlai(1) > zero) then
-          xlai1d(:) = Model%pertlai(1) * Coupling%sfc_wts(:,4)
-        endif
-! --- do the albedo percentile calculation in GFS_radiation_driver instead --- !
-!        if (Model%pertalb(1) > zero) then
-!          do i=1,im
-!            call cdfnor(Coupling%sfc_wts(i,5),cdfz)
-!            alb1d(i) = cdfz
-!          enddo
-!        endif
-        if (Model%pertvegf(1) > zero) then
-          do i=1,im
-            call cdfnor(Coupling%sfc_wts(i,6),cdfz)
-            vegf1d(i) = cdfz
-          enddo
-        endif
+      lndp_vgf=-999.
+
+      if (Model%lndp_type==1) then
+        do k =1,Model%n_var_lndp 
+           select case(Model%lndp_var_list(k)) 
+           case ('rz0') 
+                z01d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) 
+           case ('rzt') 
+                zt1d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) 
+           case ('shc') 
+                bexp1d(:) = Model%lndp_prt_list(k) * Coupling%sfc_wts(:,k) 
+           case ('lai') 
+                xlai1d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) 
+           case ('vgf') 
+! note that the pertrubed vegfrac is being used in sfc_drv, but not sfc_diff
+              do i=1,im
+                call cdfnor(Coupling%sfc_wts(i,k),cdfz)
+                vegf1d(i) = cdfz
+              enddo
+              lndp_vgf = Model%lndp_prt_list(k)
+           end select
+        enddo
       endif
 !*## CCPP ##
 !
@@ -1856,6 +1851,7 @@ subroutine GFS_physics_driver                                 &
 !    &,' pgr=',pgr(ipr),' sfcemis=',sfcemis(ipr)
 
 !## CCPP ##* sfc_drv.f/lsm_noah_run
+
           call sfc_drv                                                   &
 !  ---  inputs:
            (im, lsoil, Statein%pgr,                                      &
@@ -1867,7 +1863,8 @@ subroutine GFS_physics_driver                                 &
             Sfcprop%shdmin, Sfcprop%shdmax, Sfcprop%snoalb,              &
             Radtend%sfalb, flag_iter, flag_guess, Model%lheatstrg,       &
             Model%isot, Model%ivegsrc,                                   &
-            bexp1d, xlai1d, vegf1d, Model%pertvegf,                      &
+            bexp1d, xlai1d, vegf1d,lndp_vgf,                              &
+
 !  ---  input/output:
             weasd3(:,1), snowd3(:,1), tsfc3(:,1), tprcp3(:,1),           &
             Sfcprop%srflag, smsoil, stsoil, slsoil, Sfcprop%canopy,      &

gfsphysics/GFS_layer/GFS_radiation_driver.F90

@@ -1221,6 +1221,7 @@ subroutine GFS_radiation_driver                             &
 
       !  mg, sfc perts
       real(kind=kind_phys), dimension(size(Grid%xlon,1)) :: alb1d
+      real(kind=kind_phys) :: lndp_alb
 
       real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+ltp) :: cldtausw
       real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+ltp) :: cldtaulw
@@ -1846,14 +1847,18 @@ subroutine GFS_radiation_driver                             &
 !  ---  scale random patterns for surface perturbations with
 !  perturbation size
 !  ---  turn vegetation fraction pattern into percentile pattern
-      alb1d(:) = zero
-      if (Model%do_sfcperts) then
-        if (Model%pertalb(1) > zero) then
-          do i=1,im
-            call cdfnor(Coupling%sfc_wts(i,5),alb1d(i))
-          enddo
+      alb1d(:) = 0.
+      lndp_alb = -999.
+        if (Model%lndp_type ==1) then
+          do k =1,Model%n_var_lndp
+            if (Model%lndp_var_list(k) == 'alb') then 
+              do i=1,im
+                call cdfnor(Coupling%sfc_wts(i,k),alb1d(i)) 
+                lndp_alb = Model%lndp_prt_list(k) 
+              enddo
+            endif
+          enddo 
         endif
-      endif
 ! mg, sfc-perts
 !*## CCPP ##
 
@@ -1870,7 +1875,7 @@ subroutine GFS_radiation_driver                             &
                      Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf,    &
                      Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice,     &
                      Sfcprop%tisfc, IM,                              &
-                     alb1d, Model%pertalb,                           & !  mg, sfc-perts
+                     alb1d, lndp_alb,                                & !  mg, sfc-perts
                      sfcalb)                                           !  ---  outputs
 
 !> -# Approximate mean surface albedo from vis- and nir-  diffuse values.

gfsphysics/GFS_layer/GFS_typedefs.F90

@@ -517,7 +517,6 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: skebu_wts (:,:) => null()  !
     real (kind=kind_phys), pointer :: skebv_wts (:,:) => null()  !
     real (kind=kind_phys), pointer :: sfc_wts   (:,:) => null()  ! mg, sfc-perts
-    integer              :: nsfcpert=6                             !< number of sfc perturbations
 
     !--- aerosol surface emissions for Thompson microphysics
     real (kind=kind_phys), pointer :: nwfa2d  (:)     => null()  !< instantaneous water-friendly sfc aerosol source
@@ -1045,14 +1044,14 @@ module GFS_typedefs
     logical              :: do_shum
     logical              :: do_skeb
     integer              :: skeb_npass
-    logical              :: do_sfcperts
-    integer              :: nsfcpert=6
-    real(kind=kind_phys) :: pertz0(5)          ! mg, sfc-perts
-    real(kind=kind_phys) :: pertzt(5)          ! mg, sfc-perts
-    real(kind=kind_phys) :: pertshc(5)         ! mg, sfc-perts
-    real(kind=kind_phys) :: pertlai(5)         ! mg, sfc-perts
-    real(kind=kind_phys) :: pertalb(5)         ! mg, sfc-perts
-    real(kind=kind_phys) :: pertvegf(5)        ! mg, sfc-perts
+    integer              :: lndp_type
+    integer              :: n_var_lndp
+    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 
+                                              ! multiple patterns. It wasn't fully coded (and wouldn't have worked 
+                                              ! with nlndp>1, so I just dropped it). If we want to code it properly, 
+                                              ! we'd need to make this dim(6,5).
 !--- tracer handling
     character(len=32), pointer :: tracer_names(:) !< array of initialized tracers from dynamic core
     integer              :: ntrac           !< number of tracers
@@ -1948,6 +1947,8 @@ module GFS_typedefs
     real (kind=kind_phys), pointer      :: uustar_ocean(:)    => null()  !<
     real (kind=kind_phys), pointer      :: vdftra(:,:,:)      => null()  !<
     real (kind=kind_phys), pointer      :: vegf1d(:)          => null()  !<
+    real (kind=kind_phys)               :: lndp_vgf                      !<
+
     integer, pointer                    :: vegtype(:)         => null()  !<
     real (kind=kind_phys), pointer      :: w_upi(:,:)         => null()  !<
     real (kind=kind_phys), pointer      :: wcbmax(:)          => null()  !<
@@ -2800,9 +2801,9 @@ subroutine coupling_create (Coupling, IM, Model)
       Coupling%skebv_wts = clear_val
     endif
 
-    !--- stochastic physics option
-    if (Model%do_sfcperts) then
-      allocate (Coupling%sfc_wts  (IM,Model%nsfcpert))
+    !--- stochastic land perturbation option
+    if (Model%lndp_type .NE. 0) then
+      allocate (Coupling%sfc_wts  (IM,Model%n_var_lndp))
       Coupling%sfc_wts = clear_val
     endif
 
@@ -3314,15 +3315,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     logical :: use_zmtnblck = .false.
     logical :: do_shum      = .false.
     logical :: do_skeb      = .false.
-    integer :: skeb_npass = 11
-    logical :: do_sfcperts = .false.   ! mg, sfc-perts
-    integer :: nsfcpert    =  6        ! mg, sfc-perts
-    real(kind=kind_phys) :: pertz0   = -999.
-    real(kind=kind_phys) :: pertzt   = -999.
-    real(kind=kind_phys) :: pertshc  = -999.
-    real(kind=kind_phys) :: pertlai  = -999.
-    real(kind=kind_phys) :: pertalb  = -999.
-    real(kind=kind_phys) :: pertvegf = -999.
+    integer :: skeb_npass   = 11
+    integer :: lndp_type    = 0 
+    integer :: n_var_lndp   =  0 
 
 !--- aerosol scavenging factors
     character(len=20) :: fscav_aero(20) = 'default'
@@ -3399,7 +3394,7 @@ 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, do_sfcperts,                      &
+                               do_sppt, do_shum, do_skeb, lndp_type,  n_var_lndp,           & 
                           !--- Rayleigh friction
                                prslrd0, ral_ts,  ldiag_ugwp, do_ugwp, do_tofd,              &
                           ! --- Ferrier-Aligo
@@ -4007,21 +4002,15 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%e0fac            = e0fac
 
 !--- stochastic physics options
-    ! do_sppt, do_shum, do_skeb and do_sfcperts are namelist variables in group
+    ! do_sppt, do_shum, do_skeb and lndp_type are namelist variables in group
     ! physics that are parsed here and then compared in init_stochastic_physics
     ! to the stochastic physics namelist parametersto ensure consistency.
     Model%do_sppt          = do_sppt
     Model%use_zmtnblck     = use_zmtnblck
     Model%do_shum          = do_shum
     Model%do_skeb          = do_skeb
-    Model%do_sfcperts      = do_sfcperts ! mg, sfc-perts
-    Model%nsfcpert         = nsfcpert    ! mg, sfc-perts
-    Model%pertz0           = pertz0
-    Model%pertzt           = pertzt
-    Model%pertshc          = pertshc
-    Model%pertlai          = pertlai
-    Model%pertalb          = pertalb
-    Model%pertvegf         = pertvegf
+    Model%lndp_type        = lndp_type
+    Model%n_var_lndp       = n_var_lndp
 
     !--- cellular automata options
     Model%nca              = nca
@@ -5075,7 +5064,8 @@ subroutine control_print(Model)
       print *, ' do_sppt           : ', Model%do_sppt
       print *, ' do_shum           : ', Model%do_shum
       print *, ' do_skeb           : ', Model%do_skeb
-      print *, ' do_sfcperts       : ', Model%do_sfcperts
+      print *, ' lndp_type         : ', Model%lndp_type
+      print *, ' n_var_lndp         : ', Model%n_var_lndp
       print *, ' '
       print *, 'cellular automata'
       print *, ' nca               : ', Model%nca
@@ -7015,6 +7005,7 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     Interstitial%uustar_ocean    = huge
     Interstitial%vdftra          = clear_val
     Interstitial%vegf1d          = clear_val
+    Interstitial%lndp_vgf        = clear_val
     Interstitial%vegtype         = 0
     Interstitial%wcbmax          = clear_val
     Interstitial%weasd_ice       = huge