Update master from dtc/develop 2020/03/17

physics/GFS_PBL_generic.F90

@@ -331,7 +331,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       character(len=*), intent(out) :: errmsg
       integer, intent(out) :: errflg
 
-      real(kind=kind_phys), parameter :: huge=1.0d30
+      real(kind=kind_phys), parameter :: huge=1.0d30, epsln = 1.0d-10
       integer :: i, k, kk, k1, n
       real(kind=kind_phys) :: tem, tem1, rho
 
@@ -499,32 +499,29 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       if (cplflx) then
         do i=1,im
           if (oceanfrac(i) > 0.0) then ! Ocean only, NO LAKES
-            if (fice(i) == oceanfrac(i)) then ! use results from CICE
+            if (fice(i) > 1.-epsln) then ! no open water, use results from CICE
               dusfci_cpl(i) = dusfc_cice(i)
               dvsfci_cpl(i) = dvsfc_cice(i)
               dtsfci_cpl(i) = dtsfc_cice(i)
               dqsfci_cpl(i) = dqsfc_cice(i)
-!           elseif (dry(i) .or. icy(i)) then   ! use stress_ocean from sfc_diff for opw component at mixed point
-            elseif (wet(i)) then                   ! use stress_ocean from sfc_diff for opw component at mixed point
-              if (icy(i) .or. dry(i)) then
-                tem1 = max(q1(i), 1.e-8)
-                rho = prsl(i,1) / (rd*t1(i)*(1.0+fvirt*tem1))
-                if (wind(i) > 0.0) then
-                  tem = - rho * stress_ocn(i) / wind(i)
-                  dusfci_cpl(i) = tem * ugrs1(i)   ! U-momentum flux
-                  dvsfci_cpl(i) = tem * vgrs1(i)   ! V-momentum flux
-                else
-                  dusfci_cpl(i) = 0.0
-                  dvsfci_cpl(i) = 0.0
-                endif
-                dtsfci_cpl(i) = cp   * rho * hflx_ocn(i) ! sensible heat flux over open ocean
-                dqsfci_cpl(i) = hvap * rho * evap_ocn(i) ! latent heat flux over open ocean
-              else                                       ! use results from PBL scheme for 100% open ocean
-                dusfci_cpl(i) = dusfc1(i)
-                dvsfci_cpl(i) = dvsfc1(i)
-                dtsfci_cpl(i) = dtsfc1(i)
-                dqsfci_cpl(i) = dqsfc1(i)
+            elseif (dry(i) .or. icy(i)) then   ! use stress_ocean from sfc_diff for opw component at mixed point
+              tem1 = max(q1(i), 1.e-8)
+              rho = prsl(i,1) / (rd*t1(i)*(1.0+fvirt*tem1))
+              if (wind(i) > 0.0) then
+                tem = - rho * stress_ocn(i) / wind(i)
+                dusfci_cpl(i) = tem * ugrs1(i)   ! U-momentum flux
+                dvsfci_cpl(i) = tem * vgrs1(i)   ! V-momentum flux
+              else
+                dusfci_cpl(i) = 0.0
+                dvsfci_cpl(i) = 0.0
               endif
+              dtsfci_cpl(i) = cp   * rho * hflx_ocn(i) ! sensible heat flux over open ocean
+              dqsfci_cpl(i) = hvap * rho * evap_ocn(i) ! latent heat flux over open ocean
+            else                                       ! use results from PBL scheme for 100% open ocean
+              dusfci_cpl(i) = dusfc1(i)
+              dvsfci_cpl(i) = dvsfc1(i)
+              dtsfci_cpl(i) = dtsfc1(i)
+              dqsfci_cpl(i) = dqsfc1(i)
             endif
 !
             dusfc_cpl (i) = dusfc_cpl(i) + dusfci_cpl(i) * dtf

physics/GFS_stochastics.F90

@@ -124,7 +124,7 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb,
                if (use_zmtnblck)then
                   sppt_wts(i,k)=(sppt_wts(i,k)-1)*sppt_vwt+1.0
                endif
-               sppt_wts_inv(i,km-k+1)=sppt_wts(i,k)
+               sppt_wts_inv(i,k)=sppt_wts(i,k)
 
                !if(isppt_deep)then
 
@@ -190,16 +190,16 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb,
          if (do_shum) then
            do k=1,km
              gq0(:,k) = gq0(:,k)*(1.0 + shum_wts(:,k))
-             shum_wts_inv(:,km-k+1) = shum_wts(:,k)
+             shum_wts_inv(:,k) = shum_wts(:,k)
            end do
          endif
 
          if (do_skeb) then
            do k=1,km
              gu0(:,k) = gu0(:,k)+skebu_wts(:,k)*(diss_est(:,k))
              gv0(:,k) = gv0(:,k)+skebv_wts(:,k)*(diss_est(:,k))
-             skebu_wts_inv(:,km-k+1) = skebu_wts(:,k)
-             skebv_wts_inv(:,km-k+1) = skebv_wts(:,k)
+             skebu_wts_inv(:,k) = skebu_wts(:,k)
+             skebv_wts_inv(:,k) = skebv_wts(:,k)
            enddo
          endif
 

physics/GFS_surface_composites.F90

@@ -70,37 +70,30 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cpl
       errmsg = ''
       errflg = 0
 
-      if (frac_grid) then  ! here cice is fraction of the whole grid that is ice
+      if (frac_grid) then  ! cice is ice fraction wrt water area
         do i=1,im
           frland(i) = landfrac(i)
           if (frland(i) > zero) dry(i) = .true.
-          tem = one - frland(i)
-          if (tem > epsln) then
+          if (frland(i) < one) then
             if (flag_cice(i)) then
-              if (cice(i) >= min_seaice*tem) then
+              if (cice(i) >= min_seaice) then
                 icy(i)  = .true.
               else
                 cice(i) = zero
               endif
             else
-              if (cice(i) >= min_lakeice*tem) then
+              if (cice(i) >= min_lakeice) then
                 icy(i) = .true.
-                cice(i) = cice(i)/tem  ! cice is fraction of ocean/lake
               else
                 cice(i) = zero
               endif
             endif
-!           if (icy(i)) tsfco(i) = max(tsfco(i), tisfc(i), tgice)
+            if (cice(i) < one ) then
+              wet(i)=.true. !there is some open ocean/lake water!
+              if (.not. cplflx) tsfco(i) = max(tsfco(i), tisfc(i), tgice)
+            end if
           else
             cice(i) = zero
-          endif
-
-                                       ! ocean/lake area that is not frozen
-
-          if (tem-cice(i) > epsln) then
-            wet(i) = .true.            ! there is some open water!
-!           if (icy(i)) tsfco(i) = max(tsfco(i), tgice)
-!           if (icy(i)) tsfco(i) = max(tisfc(i), tgice)
           endif
         enddo  
 
@@ -342,8 +335,8 @@ subroutine GFS_surface_composites_post_run (
 
           ! Three-way composites (fields from sfc_diff)
           txl = landfrac(i)
-          txi = cice(i)                 ! here cice is grid fraction that is ice
-          txo = one - txl - txi
+          txi = cice(i)*(one - txl) ! txi = ice fraction wrt whole cell
+          txo = max(zero, one - txl - txi)
 
           zorl(i)   = txl*zorl_lnd(i)   + txi*zorl_ice(i)   + txo*zorl_ocn(i)
           cd(i)     = txl*cd_lnd(i)     + txi*cd_ice(i)     + txo*cd_ocn(i)
@@ -397,12 +390,6 @@ subroutine GFS_surface_composites_post_run (
 
           if (.not. flag_cice(i)) then
             if (islmsk(i) == 2) then                           ! return updated lake ice thickness & concentration to global array
-              ! DH* NOT NEEDED? Sfcprop%hice(i)  = zice(i)
-! DH* is this correct? can we update cice in place or do we need separate variables as for IPD?
-!!             Sfcprop%fice(i)  = fice(i) * Sfcprop%lakefrac(i) ! fice is fraction of lake area that is frozen
-!              Sfcprop%fice(i)  = fice(i) * (one-Sfcprop%landfrac(i)) ! fice is fraction of wet area that is frozen
-              cice(i)  = cice(i) * (1.0-landfrac(i))           ! cice is fraction of wet area that is frozen
-! *DH
               tisfc(i) = tice(i)
             else                                               ! this would be over open ocean or land (no ice fraction)
               hice(i)  = zero

physics/cu_ntiedtke.F90

@@ -148,7 +148,7 @@ end subroutine cu_ntiedtke_finalize
 !-----------------------------------------------------------------------
 !          level 1 subroutine 'tiecnvn'
 !-----------------------------------------------------------------
-      subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
+      subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
                                  evap,hfx,zprecc,lmask,lq,ix,km,dt,dx,kbot,ktop,kcnv,&
                                  ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg)
 !-----------------------------------------------------------------
@@ -162,13 +162,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
       integer, dimension( lq ),   intent(in)  :: lmask
       real(kind=kind_phys), dimension( lq ),     intent(in ) :: evap, hfx, dx
       real(kind=kind_phys), dimension( ix , km ),     intent(inout) :: pu, pv, pt, pqv
-      real(kind=kind_phys), dimension( ix , km ),     intent(in )   :: poz, prsl, pomg, pqvf, ptf
+      real(kind=kind_phys), dimension( ix , km ),     intent(in )   :: tdi, qvdi, poz, prsl, pomg, pqvf, ptf
       real(kind=kind_phys), dimension( ix , km+1 ),   intent(in )   :: pzz, prsi
-      ! DH* TODO - check dimensions of clw, ktrac+2 seems to be smaller
-      ! than the actual dimensions (ok as long as only indices 1 and 2
-      ! are accessed here, and as long as these contain what is expected);
-      ! better to expand into the cloud-ice and cloud-water components *DH
-      real(kind=kind_phys), dimension( ix , km, ktrac+2 ),    intent(inout ) ::  clw
+      real(kind=kind_phys), dimension( ix , km, ktrac ),    intent(inout ) ::  clw
 
       integer, dimension( lq ),   intent(out)  :: kbot, ktop, kcnv
       real(kind=kind_phys), dimension( lq ),   intent(out)  :: zprecc
@@ -188,13 +184,13 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
       real(kind=kind_phys) ztp1(lq,km),    zqp1(lq,km),  ztu(lq,km),   zqu(lq,km),&
      &     zlu(lq,km),     zlude(lq,km), zmfu(lq,km),  zmfd(lq,km),  zmfude_rate(lq,km),&
      &     zqsat(lq,km),   zrain(lq)
-      real(kind=kind_phys) pcen(lq,km,ktrac),ptenc(lq,km,ktrac) 
+      real(kind=kind_phys),allocatable ::  pcen(:,:,:),ptenc(:,:,:) 
 
       integer icbot(lq),   ictop(lq),     ktype(lq),   lndj(lq)
       logical locum(lq)
 !
       real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt
-      integer i,j,k,k1,n,km1
+      integer i,j,k,k1,n,km1,ktracer
       real(kind=kind_phys) ztpp1
       real(kind=kind_phys) zew,zqs,zcor
 !
@@ -246,9 +242,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           zqs  = min(0.5,zqs)                                         
           zcor = 1./(1.-vtmpc1*zqs)                                    
           zqsat(j,k1)=zqs*zcor      
-          pqte(j,k1)=pqvf(j,k) 
+          pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst 
           zqq(j,k1) =pqte(j,k1)
-          ptte(j,k1)=ptf(j,k)
+          ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst
           ztt(j,k1) =ptte(j,k1)
           ud_mf(j,k1)=0.
           dd_mf(j,k1)=0.
@@ -258,16 +254,33 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
         end do
       end do
 
-      do n=1,ktrac
-        do k=1,km
-          k1=km-k+1
-          do j=1,lq
-            pcen(j,k1,n) = clw(j,k,n+2)
-            ptenc(j,k1,n)= 0.
+      if(ktrac > 2) then
+        ktracer = ktrac - 2
+        allocate(pcen(lq,km,ktracer))
+        allocate(ptenc(lq,km,ktracer))
+        do n=1,ktracer
+          do k=1,km
+            k1=km-k+1
+            do j=1,lq
+              pcen(j,k1,n) = clw(j,k,n+2)
+              ptenc(j,k1,n)= 0.
+            end do
           end do
         end do
-      end do
-
+      else
+        ktracer = 2
+        allocate(pcen(lq,km,ktracer))
+        allocate(ptenc(lq,km,ktracer))
+        do n=1,ktracer
+          do k=1,km
+            do j=1,lq
+              pcen(j,k,n) = 0.
+              ptenc(j,k,n)= 0.
+            end do
+          end do
+        end do
+      end if
+
 !      print *, "pgeo=",pgeo(1,:)
 !      print *, "pgeoh=",pgeoh(1,:)
 !      print *, "pap=",pap(1,:)
@@ -289,7 +302,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
      &     zqp1,     pum1,     pvm1,     pverv,   zqsat,&
      &     pqhfl,    ztmst,    pap,      paph,    pgeo, &
      &     ptte,     pqte,     pvom,     pvol,    prsfc,&
-     &     pssfc,    locum,    ktrac,    pcen,    ptenc,&
+     &     pssfc,    locum,    ktracer,  pcen,    ptenc,&
      &     ktype,    icbot,    ictop,    ztu,     zqu,  &
      &     zlu,      zlude,    zmfu,     zmfd,    zrain,&
      &     pcte,     phhfl,    lndj,     pgeoh,   zmfude_rate, dx)
@@ -314,7 +327,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst
           pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst
           ud_mf(j,k)= zmfu(j,k1)*ztmst
-          dd_mf(j,k)= zmfd(j,k1)*ztmst
+          dd_mf(j,k)= -zmfd(j,k1)*ztmst
           dt_mf(j,k)= zmfude_rate(j,k1)*ztmst
           cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1))
           cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k))
@@ -343,17 +356,21 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           end do
         end do
       endif
+
 !
-      if (ktrac > 0) then
-        do n=1,ktrac
-          do k=1,km
-            k1=km-k+1
-            do j=1,lq
-              clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
-            end do
-          end do
-        end do
-      end if
+! Currently, vertical mixing of tracers are turned off
+!      if(ktrac > 2) then
+!        do n=1,ktrac-2
+!          do k=1,km
+!            k1=km-k+1
+!            do j=1,lq
+!              clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
+!            end do
+!          end do
+!        end do
+!      end if
+      deallocate(pcen)
+      deallocate(ptenc)
 !
       return
       end subroutine cu_ntiedtke_run

physics/cu_ntiedtke.meta

@@ -80,6 +80,24 @@
   kind = kind_phys
   intent = inout
   optional = F
+[tdi]
+  standard_name = air_temperature
+  long_name = mid-layer temperature
+  units = K
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[qvdi]
+  standard_name = water_vapor_specific_humidity
+  long_name = water vapor specific humidity
+  units = kg kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
 [pqvf]
   standard_name = moisture_tendency_due_to_dynamics
   long_name = moisture tendency due to dynamics only
@@ -254,8 +272,8 @@
   intent = out
   optional = F
 [ktrac]
-  standard_name = number_of_total_tracers
-  long_name = number of total tracers
+  standard_name = number_of_tracers_for_convective_transport
+  long_name = number of tracers for convective transport
   units = count
   dimensions = ()
   type = integer

physics/dcyc2.meta

@@ -183,37 +183,37 @@
   intent = in
   optional = F
 [swh]
-  standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_time_step
+  standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_timestep
   long_name = total sky shortwave heating rate on radiation time step
   units = K s-1
-  dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation)
+  dimensions = (horizontal_dimension,vertical_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [swhc]
-  standard_name = tendency_of_air_temperature_due_to_shortwave_heating_assuming_clear_sky_on_radiation_time_step
+  standard_name = tendency_of_air_temperature_due_to_shortwave_heating_assuming_clear_sky_on_radiation_timestep
   long_name = clear sky shortwave heating rate on radiation time step
   units = K s-1
-  dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation)
+  dimensions = (horizontal_dimension,vertical_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [hlw]
-  standard_name = tendency_of_air_temperature_due_to_longwave_heating_on_radiation_time_step
+  standard_name = tendency_of_air_temperature_due_to_longwave_heating_on_radiation_timestep
   long_name = total sky longwave heating rate on radiation time step
   units = K s-1
-  dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation)
+  dimensions = (horizontal_dimension,vertical_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [hlwc]
-  standard_name = tendency_of_air_temperature_due_to_longwave_heating_assuming_clear_sky_on_radiation_time_step
+  standard_name = tendency_of_air_temperature_due_to_longwave_heating_assuming_clear_sky_on_radiation_timestep
   long_name = clear sky longwave heating rate on radiation time step
   units = K s-1
-  dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation)
+  dimensions = (horizontal_dimension,vertical_dimension)
   type = real
   kind = kind_phys
   intent = in