Merge pull request #954 from grantfirl/misc_NRL_changes_20220824

Miscellaneous Bugfixes from NRL

physics/GFS_rad_time_vary.fv3.F90

@@ -16,7 +16,7 @@ module GFS_rad_time_vary
 !> \section arg_table_GFS_rad_time_vary_timestep_init Argument Table
 !! \htmlinclude GFS_rad_time_vary_timestep_init.html
 !!
-      subroutine GFS_rad_time_vary_timestep_init (                                     &
+      subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds,                     &
               lslwr, lsswr, isubc_lw, isubc_sw, icsdsw, icsdlw, cnx, cny, isc, jsc,    &
               imap, jmap, sec, kdt, imp_physics, imp_physics_zhao_carr, ps_2delt,      &
               ps_1delt, t_2delt, t_1delt, qv_2delt, qv_1delt, t, qv, ps, errmsg, errflg)
@@ -29,6 +29,8 @@ subroutine GFS_rad_time_vary_timestep_init (
          implicit none
 
          ! Interface variables
+         logical,                intent(in)    :: lrseeds
+         integer,                intent(in)    :: rseeds(:,:)
          integer,                intent(in)    :: isubc_lw, isubc_sw, cnx, cny, isc, jsc, kdt
          integer,                intent(in)    :: imp_physics, imp_physics_zhao_carr
          logical,                intent(in)    :: lslwr, lsswr
@@ -47,7 +49,7 @@ subroutine GFS_rad_time_vary_timestep_init (
 
          ! Local variables
          type (random_stat) :: stat
-         integer :: ix, j, i, nblks, ipseed
+         integer :: ix, j, i, ipseed
          integer :: numrdm(cnx*cny*2)
 
          ! Initialize CCPP error handling variables
@@ -60,18 +62,25 @@ subroutine GFS_rad_time_vary_timestep_init (
 
            !--- set up random seed index in a reproducible way for entire cubed-sphere face (lat-lon grid)
            if ((isubc_lw==2) .or. (isubc_sw==2)) then
-             ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0
-             call random_setseed (ipseed, stat)
-             call random_index (ipsdlim, numrdm, stat)
-
-             do ix=1,size(jmap)
-               j = jmap(ix)
-               i = imap(ix)
-               !--- for testing purposes, replace numrdm with '100'
-               icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx)
-               icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny)
-             enddo
-
+             !NRL If random seeds supplied by NEPTUNE
+             if(lrseeds) then
+               do ix=1,size(jmap)
+                 icsdsw(ix) = rseeds(ix,1)
+                 icsdlw(ix) = rseeds(ix,2)
+               enddo
+             else
+               ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0
+               call random_setseed (ipseed, stat)
+               call random_index (ipsdlim, numrdm, stat)
+
+               do ix=1,size(jmap)
+                 j = jmap(ix)
+                 i = imap(ix)
+                 !--- for testing purposes, replace numrdm with '100'
+                 icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx)
+                 icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny)
+               enddo
+             end if !lrseeds
            endif  ! isubc_lw and isubc_sw
 
            if (imp_physics == imp_physics_zhao_carr) then

physics/GFS_rad_time_vary.fv3.meta

@@ -7,6 +7,20 @@
 [ccpp-arg-table]
   name = GFS_rad_time_vary_timestep_init
   type = scheme
+[lrseeds]
+  standard_name = do_host_provided_random_seeds
+  long_name = flag to use host-provided random seeds
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[rseeds]
+  standard_name = random_number_seeds_from_host
+  long_name = random number seeds from host
+  units = none
+  dimensions = (horizontal_dimension, number_of_host_provided_random_number_streams)
+  type = integer
+  intent = in
 [lslwr]
   standard_name = flag_for_calling_longwave_radiation
   long_name = logical flags for lw radiation calls

physics/GFS_rad_time_vary.scm.F90

@@ -11,11 +11,12 @@ module GFS_rad_time_vary
       contains
 
 !>\defgroup mod_GFS_rad_time_vary GFS Radiation Time Update
+!! This module contains code related to GFS radiation setup.
 !> @{
 !> \section arg_table_GFS_rad_time_vary_timestep_init Argument Table
 !! \htmlinclude GFS_rad_time_vary_timestep_init.html
 !!
-      subroutine GFS_rad_time_vary_timestep_init (                                     &
+      subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds,                     &
               lslwr, lsswr, isubc_lw, isubc_sw, icsdsw, icsdlw, cnx, cny, isc, jsc,    &
               imap, jmap, sec, kdt, imp_physics, imp_physics_zhao_carr, ps_2delt,      &
               ps_1delt, t_2delt, t_1delt, qv_2delt, qv_1delt, t, qv, ps, errmsg, errflg)
@@ -28,6 +29,8 @@ subroutine GFS_rad_time_vary_timestep_init (
          implicit none
 
          ! Interface variables
+         logical,                intent(in)    :: lrseeds
+         integer,                intent(in)    :: rseeds(:,:)
          integer,                intent(in)    :: isubc_lw, isubc_sw, cnx, cny, isc, jsc, kdt
          integer,                intent(in)    :: imp_physics, imp_physics_zhao_carr
          logical,                intent(in)    :: lslwr, lsswr
@@ -59,26 +62,33 @@ subroutine GFS_rad_time_vary_timestep_init (
 
            !--- set up random seed index in a reproducible way for entire cubed-sphere face (lat-lon grid)
            if ((isubc_lw==2) .or. (isubc_sw==2)) then
-             ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0
-             call random_setseed (ipseed, stat)
-             call random_index (ipsdlim, numrdm, stat)
-
-             do ix=1,size(jmap)
-               j = jmap(ix)
-               i = imap(ix)
-               !--- for testing purposes, replace numrdm with '100'
-               icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx)
-               icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny)
-             enddo
-
+             !NRL If random seeds supplied by NEPTUNE
+             if(lrseeds) then
+               do ix=1,size(jmap)
+                 icsdsw(ix) = rseeds(ix,1)
+                 icsdlw(ix) = rseeds(ix,2)
+               enddo
+             else
+               ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0
+               call random_setseed (ipseed, stat)
+               call random_index (ipsdlim, numrdm, stat)
+
+               do ix=1,size(jmap)
+                 j = jmap(ix)
+                 i = imap(ix)
+                 !--- for testing purposes, replace numrdm with '100'
+                 icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx)
+                 icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny)
+               enddo
+             end if !lrseeds
            endif  ! isubc_lw and isubc_sw
 
            if (imp_physics == imp_physics_zhao_carr) then
              if (kdt == 1) then
                t_2delt  = t
                t_1delt  = t
-               qv_2delt = max(qmin,qv)
-               qv_1delt = max(qmin,qv)
+               qv_2delt = qv
+               qv_1delt = qv
                ps_2delt = ps
                ps_1delt = ps
              endif

physics/GFS_rad_time_vary.scm.meta

@@ -7,6 +7,20 @@
 [ccpp-arg-table]
   name = GFS_rad_time_vary_timestep_init
   type = scheme
+[lrseeds]
+  standard_name = do_host_provided_random_seeds
+  long_name = flag to use host-provided random seeds
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[rseeds]
+  standard_name = random_number_seeds_from_host
+  long_name = random number seeds from host
+  units = none
+  dimensions = (horizontal_dimension, number_of_host_provided_random_number_streams)
+  type = integer
+  intent = in
 [lslwr]
   standard_name = flag_for_calling_longwave_radiation
   long_name = logical flags for lw radiation calls

physics/GFS_rrtmg_post.meta

@@ -38,7 +38,7 @@
 [ltp]
   standard_name = extra_top_layer
   long_name = extra top layers
-  units = none
+  units = count
   dimensions = ()
   type = integer
   intent = in

physics/GFS_rrtmg_pre.F90

@@ -17,9 +17,9 @@ module GFS_rrtmg_pre
 !! \htmlinclude GFS_rrtmg_pre_run.html
 !!    
 !>\section rrtmg_pre_gen General Algorithm
-      subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
-        imfdeepcnv, imfdeepcnv_gf, me, ncnd, ntrac, num_p3d, npdf3d, ncnvcld3d,&
-        ntqv, ntcw,ntiw, ntlnc, ntinc, ntrnc, ntsnc, ntccn,                    &
+      subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
+        n_var_lndp, imfdeepcnv, imfdeepcnv_gf, me, ncnd, ntrac, num_p3d,       &
+        npdf3d, ncnvcld3d, ntqv, ntcw,ntiw, ntlnc, ntinc, ntrnc, ntsnc, ntccn, &
         ntrw, ntsw, ntgl, nthl, ntwa, ntoz,                                    &
         ntclamt, nleffr, nieffr, nseffr, lndp_type, kdt,                       &
         ntdu1, ntdu2, ntdu3, ntdu4, ntdu5, ntss1, ntss2,                       &
@@ -50,7 +50,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
 
       use physparam
 
-      use radcons,                   only: itsfc,ltp, lextop, qmin,  &
+      use radcons,                   only: itsfc, qmin,  &
                                            qme5, qme6, epsq, prsmin
       use funcphys,                  only: fpvs
 
@@ -84,8 +84,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
       use physparam,              only : iaermdl
       implicit none
 
-      integer,              intent(in)  :: im, levs, lm, lmk, lmp, n_var_lndp, &
-                                           imfdeepcnv,                         &
+      integer,              intent(in)  :: im, levs, lm, lmk, lmp, ltp,        &
+                                           n_var_lndp, imfdeepcnv,             &
                                            imfdeepcnv_gf, me, ncnd, ntrac,     &
                                            num_p3d, npdf3d, ncnvcld3d, ntqv,   &
                                            ntcw, ntiw, ntlnc, ntinc,           &
@@ -120,8 +120,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
 
       character(len=3), dimension(:), intent(in) :: lndp_var_list
 
-      logical,              intent(in) :: lsswr, lslwr, ltaerosol, lgfdlmprad, &
-                                          uni_cld, effr_in, do_mynnedmf,       &
+      logical,              intent(in) :: lextop, lsswr, lslwr, ltaerosol, lgfdlmprad, &
+                                          uni_cld, effr_in, do_mynnedmf,               &
                                           lmfshal, lmfdeep2, pert_clds
       logical,              intent(in) :: aero_dir_fdb
       real(kind=kind_phys), dimension(:,:), intent(in) :: smoke_ext, dust_ext
@@ -346,7 +346,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
           plyr(i,k1)    = prsl(i,k2)    * 0.01   ! pa to mb (hpa)
           tlyr(i,k1)    = tgrs(i,k2)
           prslk1(i,k1)  = prslk(i,k2)
-
+          rho(i,k1)     = prsl(i,k2)/(con_rd*tlyr(i,k1))
+          orho(i,k1)    = 1.0/rho(i,k1)
+
 !> - Compute relative humidity.
           es  = min( prsl(i,k2),  fpvs( tgrs(i,k2) ) )  ! fpvs and prsl in pa
           qs  = max( QMIN, con_eps * es / (prsl(i,k2) + epsm1*es) )
@@ -402,6 +404,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
           plyr(i,lyb)   = 0.5 * plvl(i,lla)
           tlyr(i,lyb)   = tlyr(i,lya)
           prslk1(i,lyb) = (plyr(i,lyb)*0.001) ** rocp ! plyr in hPa
+          rho(i,lyb)    = plyr(i,lyb) *100.0/(con_rd*tlyr(i,lyb))
+          orho(i,lyb)   = 1.0/rho(i,lyb)
           rhly(i,lyb)   = rhly(i,lya)
           qstl(i,lyb)   = qstl(i,lya)
         enddo

physics/GFS_rrtmg_pre.meta

@@ -44,6 +44,20 @@
   dimensions = ()
   type = integer
   intent = in
+[lextop]
+  standard_name = do_extra_top_layer_for_radiation
+  long_name = use an extra top layer for radiation
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[ltp]
+  standard_name = extra_top_layer
+  long_name = extra top layer for radiation
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [n_var_lndp]
   standard_name = number_of_perturbed_land_surface_variables
   long_name = number of land surface variables perturbed

physics/GFS_rrtmg_setup.F90

@@ -13,8 +13,6 @@ module GFS_rrtmg_setup
    &             iswcliq,                                             &
    &             kind_phys
 
-   use radcons, only: ltp, lextop
-
    implicit none
 
    public GFS_rrtmg_setup_init, GFS_rrtmg_setup_timestep_init, GFS_rrtmg_setup_finalize
@@ -50,7 +48,7 @@ subroutine GFS_rrtmg_setup_init (                          &
           icliq_sw, crick_proof, ccnorm,                      &
           imp_physics,                                        &
           norad_precip, idate, iflip,                         &
-          do_RRTMGP, me, errmsg, errflg)
+          do_RRTMGP, me, ltp, lextop, errmsg, errflg)
 ! =================   subprogram documentation block   ================ !
 !                                                                       !
 ! subprogram:   GFS_rrtmg_setup_init - a subprogram to initialize radiation !
@@ -141,6 +139,8 @@ subroutine GFS_rrtmg_setup_init (                          &
 !                     =0: index from toa to surface                     !
 !                     =1: index from surface to toa                     !
 !   me               : print control flag                               !
+!   ltp              : number of radiation extra top layers             !
+!   lextop           : control flag to denote extra top layers are used !
 !                                                                       !
 !  subroutines called: radinit                                          !
 !                                                                       !
@@ -171,6 +171,8 @@ subroutine GFS_rrtmg_setup_init (                          &
       integer, intent(in) :: iflip
       logical, intent(in) :: do_RRTMGP
       integer, intent(in) :: me
+      integer, intent(in) :: ltp
+      logical, intent(in) :: lextop
       character(len=*), intent(out) :: errmsg
       integer,          intent(out) :: errflg
 
@@ -243,7 +245,7 @@ subroutine GFS_rrtmg_setup_init (                          &
 
       call radinit                                                      &
 !  ---  inputs:
-     &     ( si, levr, imp_physics, me )
+     &     ( si, levr, imp_physics, me, ltp, lextop )
 !  ---  outputs:
 !          ( none )
 
@@ -324,7 +326,7 @@ end subroutine GFS_rrtmg_setup_finalize
 ! Private functions
 
 
-   subroutine radinit( si, NLAY, imp_physics, me )
+   subroutine radinit( si, NLAY, imp_physics, me, ltp, lextop )
 !...................................
 
 !  ---  inputs:
@@ -437,7 +439,8 @@ subroutine radinit( si, NLAY, imp_physics, me )
       implicit none
 
 !  ---  inputs:
-      integer, intent(in) :: NLAY, me, imp_physics 
+      integer, intent(in) :: NLAY, me, imp_physics, ltp
+      logical, intent(in) :: lextop 
 
       real (kind=kind_phys), intent(in) :: si(:)
 

physics/GFS_rrtmg_setup.meta

@@ -163,6 +163,20 @@
   dimensions = ()
   type = integer
   intent = in
+[ltp]
+  standard_name = extra_top_layer
+  long_name = extra top layer for radiation
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[lextop]
+  standard_name = do_extra_top_layer_for_radiation
+  long_name = use an extra top layer for radiation
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/cu_gf_deep.F90

@@ -2908,12 +2908,12 @@ subroutine cup_env(z,qes,he,hes,t,q,p,z1,                &
 !$acc declare copyin(p,t,q)
      real(kind=kind_phys),    dimension (its:ite,kts:kte)                &
         ,intent (out  )                   ::             &
-        he,hes,qes
-!$acc declare copyout(he,hes,qes)
+        hes,qes
+!$acc declare copyout(hes,qes)
      real(kind=kind_phys),    dimension (its:ite,kts:kte)                &
         ,intent (inout)                   ::             &
-        z
-!$acc declare copy(z)
+        he,z
+!$acc declare copy(he,z)
      real(kind=kind_phys),    dimension (its:ite)                        &
         ,intent (in   )                   ::             &
         psur,z1

physics/h2ophys.f

@@ -122,7 +122,7 @@ subroutine h2ophys_run(im, levs, kh2o, dt, h2o, ph2o, prsl,       &
           enddo
         endif
         do i=1,im
-          if (prsl(i,l) < prsmax) then
+          if (prsl(i,l) < prsmax .and. pltc(i,2) /= 0.0) then
             h2oib(i)  = h2o(i,l)            ! no filling
             tem       = 1.0 / pltc(i,2)     ! 1/teff
             h2o(i,l)  = (h2oib(i) + (pltc(i,1)+pltc(i,3)*tem)*dt)