Merge branch 'gfsmeta' of https://github.com/XiaSun-NOAA/fv3atm into …

…gfsmeta

.gitmodules

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

atmos_model.F90

@@ -297,9 +297,7 @@ subroutine update_atmos_radiation_physics (Atmos)
 #endif
 
 !--- call stochastic physics pattern generation / cellular automata
-    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)
 
 !--- if coupled, assign coupled fields
 
@@ -630,8 +628,7 @@ 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, ierr)
-   if (ierr/=0)  call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
+   call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block)
 
    Atmos%Diag => IPD_Diag
 

gfsphysics/CMakeLists.txt

@@ -17,6 +17,7 @@ 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,7 +633,6 @@ 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,7 +660,6 @@ 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).
@@ -929,28 +928,34 @@ subroutine GFS_physics_driver                                 &
 !        alb1d(i)  = zero
          vegf1d(i) = zero
       enddo
-      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
+      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
       endif
 !*## CCPP ##
 !
@@ -1851,7 +1856,6 @@ 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,                                      &
@@ -1863,8 +1867,7 @@ 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,lndp_vgf,                              &
-
+            bexp1d, xlai1d, vegf1d, Model%pertvegf,                      &
 !  ---  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,7 +1221,6 @@ 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
@@ -1847,18 +1846,14 @@ subroutine GFS_radiation_driver                             &
 !  ---  scale random patterns for surface perturbations with
 !  perturbation size
 !  ---  turn vegetation fraction pattern into percentile pattern
-      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 
+      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
         endif
+      endif
 ! mg, sfc-perts
 !*## CCPP ##
 
@@ -1875,7 +1870,7 @@ subroutine GFS_radiation_driver                             &
                      Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf,    &
                      Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice,     &
                      Sfcprop%tisfc, IM,                              &
-                     alb1d, lndp_alb,                                & !  mg, sfc-perts
+                     alb1d, Model%pertalb,                           & !  mg, sfc-perts
                      sfcalb)                                           !  ---  outputs
 
 !> -# Approximate mean surface albedo from vis- and nir-  diffuse values.

gfsphysics/GFS_layer/GFS_typedefs.F90

@@ -517,6 +517,7 @@ 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
@@ -1044,14 +1045,14 @@ module GFS_typedefs
     logical              :: do_shum
     logical              :: do_skeb
     integer              :: skeb_npass
-    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).
+    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
 !--- tracer handling
     character(len=32), pointer :: tracer_names(:) !< array of initialized tracers from dynamic core
     integer              :: ntrac           !< number of tracers
@@ -1947,8 +1948,6 @@ 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()  !<
@@ -2801,9 +2800,9 @@ subroutine coupling_create (Coupling, IM, Model)
       Coupling%skebv_wts = clear_val
     endif
 
-    !--- stochastic land perturbation option
-    if (Model%lndp_type .NE. 0) then
-      allocate (Coupling%sfc_wts  (IM,Model%n_var_lndp))
+    !--- stochastic physics option
+    if (Model%do_sfcperts) then
+      allocate (Coupling%sfc_wts  (IM,Model%nsfcpert))
       Coupling%sfc_wts = clear_val
     endif
 
@@ -3315,9 +3314,15 @@ 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
-    integer :: lndp_type    = 0 
-    integer :: n_var_lndp   =  0 
+    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.
 
 !--- aerosol scavenging factors
     character(len=20) :: fscav_aero(20) = 'default'
@@ -3394,7 +3399,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, lndp_type,  n_var_lndp,           & 
+                               do_sppt, do_shum, do_skeb, do_sfcperts,                      &
                           !--- Rayleigh friction
                                prslrd0, ral_ts,  ldiag_ugwp, do_ugwp, do_tofd,              &
                           ! --- Ferrier-Aligo
@@ -4002,15 +4007,21 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%e0fac            = e0fac
 
 !--- stochastic physics options
-    ! do_sppt, do_shum, do_skeb and lndp_type are namelist variables in group
+    ! do_sppt, do_shum, do_skeb and do_sfcperts 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%lndp_type        = lndp_type
-    Model%n_var_lndp       = n_var_lndp
+    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
 
     !--- cellular automata options
     Model%nca              = nca
@@ -5064,8 +5075,7 @@ subroutine control_print(Model)
       print *, ' do_sppt           : ', Model%do_sppt
       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 *, ' do_sfcperts       : ', Model%do_sfcperts
       print *, ' '
       print *, 'cellular automata'
       print *, ' nca               : ', Model%nca
@@ -7005,7 +7015,6 @@ 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

gfsphysics/physics/GFS_debug.F90

@@ -468,7 +468,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                         call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts )
                         call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts )
                      end if
-                     if (Model%lndp_type .NE. 0) then
+                     if (Model%do_sfcperts) then
                         call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts'  , Coupling%sfc_wts   )
                      end if
                      if (Model%do_ca) then

gfsphysics/physics/radiation_surface.f

@@ -382,7 +382,7 @@ subroutine setalb                                                 &
      &       slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif,           &
      &       alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc,     &
      &       sncovr, snoalb, albPpert                                     ! sfc-perts, mgehne
-      real (kind=kind_phys),  intent(in) :: pertalb                     ! sfc-perts, mgehne
+      real (kind=kind_phys), dimension(5), intent(in) :: pertalb          ! sfc-perts, mgehne
 
 !  ---  outputs
       real (kind=kind_phys), dimension(IMAX,NF_ALBD), intent(out) ::    &
@@ -620,12 +620,12 @@ subroutine setalb                                                 &
 
 ! sfc-perts, mgehne ***
 ! perturb all 4 kinds of surface albedo, sfcalb(:,1:4)
-      if (pertalb>0.0) then
+      if (pertalb(1) > 0.0) then
         do i = 1, imax
           do kk=1, 4
             ! compute beta distribution parameters for all 4 albedos
             m = sfcalb(i,kk)
-            s = pertalb*m*(1.-m)
+            s = pertalb(1)*m*(1.-m)
             alpha = m*m*(1.-m)/(s*s)-m
             beta  = alpha*(1.-m)/m
             ! compute beta distribution value corresponding

gfsphysics/physics/sfc_drv.f

@@ -182,7 +182,7 @@ subroutine sfc_drv                                                &
 
 !  ---  input:
       integer, intent(in) :: im, km, isot, ivegsrc
-      real (kind=kind_phys), intent(in) :: pertvegf
+      real (kind=kind_phys), dimension(5), intent(in) :: pertvegf
 
       integer, dimension(im), intent(in) :: soiltyp, vegtype, slopetyp
 
@@ -383,10 +383,10 @@ subroutine sfc_drv                                                &
 !  perturb vegetation fraction that goes into sflx, use the same
 !  perturbation strategy as for albedo (percentile matching)
           vegfp  = vegfpert(i)                    ! sfc-perts, mgehne
-          if (pertvegf > zero) then
+          if (pertvegf(1) > zero) then
                 ! compute beta distribution parameters for vegetation fraction
                 mv = shdfac
-                sv = pertvegf*mv*(one-mv)
+                sv = pertvegf(1)*mv*(1.-mv)
                 alphav = mv*mv*(one-mv)/(sv*sv)-mv
                 betav  = alphav*(one-mv)/mv
 ! compute beta distribution value corresponding

stochastic_physics/makefile

@@ -18,7 +18,7 @@ $(info $$ESMF_INC is [${ESMF_INC}])
 
 LIBRARY  = libstochastic_physics_wrapper.a
 
-FFLAGS   += -I$(FMS_DIR) -I ../../stochastic_physics -I../ccpp/physics -I../ccpp/build/physics -I../atmos_cubed_sphere
+FFLAGS   += -I$(FMS_DIR) -I ../../stochastic_physics -I../ccpp/physics -I../atmos_cubed_sphere
 
 SRCS_f   =