master: regain bit-for-bit identical results between IPD and CCPP for coupled model runs

physics/GFS_MP_generic.F90

@@ -191,11 +191,11 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
       end if
 
       if (lsm==lsm_ruc .or. lsm==lsm_noahmp) then
-            raincprv(:)   = rainc(:)
-            rainncprv(:)  = frain * rain1(:)
-            iceprv(:)     = ice(:)
-            snowprv(:)    = snow(:)
-            graupelprv(:) = graupel(:)
+        raincprv(:)   = rainc(:)
+        rainncprv(:)  = frain * rain1(:)
+        iceprv(:)     = ice(:)
+        snowprv(:)    = snow(:)
+        graupelprv(:) = graupel(:)
         !for NoahMP, calculate precipitation rates from liquid water equivalent thickness for use in next time step
         !Note (GJF): Precipitation LWE thicknesses are multiplied by the frain factor, and are thus on the dynamics time step, but the conversion as written
         !            (with dtp in the denominator) assumes the rate is calculated on the physics time step. This only works as expected when dtf=dtp (i.e. when frain=1).
@@ -341,8 +341,10 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
 
       if (cplflx .or. cplchm) then
         do i = 1, im
-          rain_cpl(i) = rain_cpl(i) + rain(i) * (one-srflag(i))
-          snow_cpl(i) = snow_cpl(i) + rain(i) * srflag(i)
+          drain_cpl(i) = rain(i) * (one-srflag(i))
+          dsnow_cpl(i) = rain(i) * srflag(i)
+          rain_cpl(i) = rain_cpl(i) + drain_cpl(i)
+          snow_cpl(i) = snow_cpl(i) + dsnow_cpl(i)
         enddo
       endif
 
@@ -376,15 +378,6 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
       if (do_sppt) then
 !--- radiation heating rate
         dtdtr(1:im,:) = dtdtr(1:im,:) + dtdtc(1:im,:)*dtf
-        do i = 1, im
-          if (t850(i) > 273.16) then
-!--- change in change in rain precip
-             drain_cpl(i) = rain(i) - drain_cpl(i)
-          else
-!--- change in change in snow precip
-             dsnow_cpl(i) = rain(i) - dsnow_cpl(i)
-          endif
-        enddo
       endif
 
     end subroutine GFS_MP_generic_post_run

physics/GFS_PBL_generic.F90

@@ -331,7 +331,10 @@ 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, epsln = 1.0d-10
+      real(kind=kind_phys), parameter :: zero  = 0.0d0
+      real(kind=kind_phys), parameter :: one   = 1.0d0
+      real(kind=kind_phys), parameter :: huge  = 9.9692099683868690E36 ! NetCDF float FillValue, same as in GFS_typedefs.F90
+      real(kind=kind_phys), parameter :: epsln = 1.0d-10 ! same as in GFS_physics_driver.F90
       integer :: i, k, kk, k1, n
       real(kind=kind_phys) :: tem, tem1, rho
 
@@ -486,7 +489,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       if (cplchm) then
         do i = 1, im
           tem1 = max(q1(i), 1.e-8)
-          tem  = prsl(i,1) / (rd*t1(i)*(1.0+fvirt*tem1))
+          tem  = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1))
           ushfsfci(i) = -cp * tem * hflx(i) ! upward sensible heat flux
         enddo
         ! dkt_cpl has dimensions (1:im,1:levs), but dkt has (1:im,1:levs-1)
@@ -498,22 +501,22 @@ 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) > 1.-epsln) then ! no open water, use results from CICE
+          if (oceanfrac(i) > zero) then ! Ocean only, NO LAKES
+            if (fice(i) > one - 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 (icy(i) .or. dry(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
+              rho = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1))
+              if (wind(i) > zero) 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
+                dusfci_cpl(i) = zero
+                dvsfci_cpl(i) = zero
               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

physics/GFS_PBL_generic.meta

@@ -1089,35 +1089,35 @@
   intent = in
   optional = F
 [dusfc_cice]
-  standard_name = surface_x_momentum_flux_for_coupling_interstitial
-  long_name = sfc x momentum flux for coupling interstitial
+  standard_name = surface_x_momentum_flux_for_coupling
+  long_name = sfc x momentum flux for coupling
   units = Pa
   dimensions = (horizontal_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [dvsfc_cice]
-  standard_name = surface_y_momentum_flux_for_coupling_interstitial
-  long_name = sfc y momentum flux for coupling interstitial
+  standard_name = surface_y_momentum_flux_for_coupling
+  long_name = sfc y momentum flux for coupling
   units = Pa
   dimensions = (horizontal_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [dtsfc_cice]
-  standard_name = surface_upward_sensible_heat_flux_for_coupling_interstitial
-  long_name = sfc sensible heat flux for coupling interstitial
+  standard_name = surface_upward_sensible_heat_flux_for_coupling
+  long_name = sfc sensible heat flux for coupling
   units = W m-2
   dimensions = (horizontal_dimension)
   type = real
   kind = kind_phys
   intent = in
   optional = F
 [dqsfc_cice]
-  standard_name = surface_upward_latent_heat_flux_for_coupling_interstitial
-  long_name = sfc latent heat flux for coupling interstitial
+  standard_name = surface_upward_latent_heat_flux_for_coupling
+  long_name = sfc latent heat flux for coupling
   units = W m-2
   dimensions = (horizontal_dimension)
   type = real

physics/GFS_debug.F90

@@ -402,7 +402,12 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                         call print_var(mpirank,omprank, blkno, 'Coupling%rain_cpl', Coupling%rain_cpl)
                         call print_var(mpirank,omprank, blkno, 'Coupling%snow_cpl', Coupling%snow_cpl)
                      end if
+                     if (Model%cplwav2atm) then
+                        call print_var(mpirank,omprank, blkno, 'Coupling%zorlwav_cpl' , Coupling%zorlwav_cpl  )
+                     end if
                      if (Model%cplflx) then
+                        call print_var(mpirank,omprank, blkno, 'Coupling%oro_cpl'     , Coupling%oro_cpl      )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%slmsk_cpl'   , Coupling%slmsk_cpl    )
                         call print_var(mpirank,omprank, blkno, 'Coupling%slimskin_cpl', Coupling%slimskin_cpl )
                         call print_var(mpirank,omprank, blkno, 'Coupling%dusfcin_cpl ', Coupling%dusfcin_cpl  )
                         call print_var(mpirank,omprank, blkno, 'Coupling%dvsfcin_cpl ', Coupling%dvsfcin_cpl  )
@@ -466,11 +471,24 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                         call print_var(mpirank,omprank, blkno, 'Coupling%shum_wts', Coupling%shum_wts)
                      end if
                      if (Model%do_skeb) then
-                        call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts)
-                        call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts)
+                        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%do_sfcperts) then
-                        call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts', Coupling%sfc_wts)
+                        call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts'  , Coupling%sfc_wts   )
+                     end if
+                     if (Model%do_ca) then
+                        call print_var(mpirank,omprank, blkno, 'Coupling%tconvtend', Coupling%tconvtend )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%qconvtend', Coupling%qconvtend )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%uconvtend', Coupling%uconvtend )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%vconvtend', Coupling%vconvtend )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_out   ', Coupling%ca_out    )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_deep  ', Coupling%ca_deep   )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_turb  ', Coupling%ca_turb   )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_shal  ', Coupling%ca_shal   )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_rad   ', Coupling%ca_rad    )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%ca_micro ', Coupling%ca_micro  )
+                        call print_var(mpirank,omprank, blkno, 'Coupling%cape     ', Coupling%cape      )
                      end if
                      if(Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then
                         call print_var(mpirank,omprank, blkno, 'Coupling%nwfa2d', Coupling%nwfa2d)

physics/GFS_stochastics.F90

@@ -79,10 +79,10 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb,
          real(kind_phys), dimension(1:im),      intent(inout) :: totprcpb
          real(kind_phys), dimension(1:im),      intent(inout) :: cnvprcpb
          logical,                               intent(in)    :: cplflx
-         ! rain_cpl, snow_cpl only allocated if cplflx == .true. or do_sppt == .true.
+         ! rain_cpl, snow_cpl only allocated if cplflx == .true. or cplchm == .true.
          real(kind_phys), dimension(:),         intent(inout) :: rain_cpl
          real(kind_phys), dimension(:),         intent(inout) :: snow_cpl
-         ! drain_cpl, dsnow_cpl only allocated if do_sppt == .true.
+         ! drain_cpl, dsnow_cpl only allocated if cplflx == .true. or cplchm == .true.
          real(kind_phys), dimension(:),         intent(in)    :: drain_cpl
          real(kind_phys), dimension(:),         intent(in)    :: dsnow_cpl
          ! tconvtend ... vconvtend only allocated if isppt_deep == .true.

physics/GFS_suite_interstitial.F90

@@ -228,15 +228,15 @@ subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplfl
 
         if (frac_grid) then
           do i=1,im
-            tem = one - cice(i) - frland(i)
+            tem = (one - frland(i)) * cice(i) ! tem = ice fraction wrt whole cell
             if (flag_cice(i)) then
-              adjsfculw(i) = adjsfculw_lnd(i) * frland(i) &
-                           + ulwsfc_cice(i)   * cice(i)   &
-                           + adjsfculw_ocn(i) * tem
+              adjsfculw(i) = adjsfculw_lnd(i) * frland(i)               &
+                           + ulwsfc_cice(i)   * tem                     &
+                           + adjsfculw_ocn(i) * (one - frland(i) - tem)
             else
-              adjsfculw(i) = adjsfculw_lnd(i) * frland(i) &
-                           + adjsfculw_ice(i) * cice(i)   &
-                           + adjsfculw_ocn(i) * tem
+              adjsfculw(i) = adjsfculw_lnd(i) * frland(i)               &
+                           + adjsfculw_ice(i) * tem                     &
+                           + adjsfculw_ocn(i) * (one - frland(i) - tem)
             endif
           enddo
         else

physics/GFS_surface_composites.F90

@@ -89,7 +89,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cpl
               endif
             endif
             if (cice(i) < one ) then
-              wet(i)=.true. !there is some open ocean/lake water!
+              wet(i)=.true. ! some open ocean/lake water exists
               if (.not. cplflx) tsfco(i) = max(tsfco(i), tisfc(i), tgice)
             end if
           else
@@ -414,7 +414,7 @@ subroutine GFS_surface_composites_post_run (
             fm10(i)   = fm10_lnd(i)
             fh2(i)    = fh2_lnd(i)
            !tsurf(i)  = tsurf_lnd(i)
-            tsfcl(i)  = tsfc_lnd(i)
+            tsfcl(i)  = tsfc_lnd(i) ! over land
             cmm(i)    = cmm_lnd(i)
             chh(i)    = chh_lnd(i)
             gflx(i)   = gflx_lnd(i)
@@ -426,9 +426,9 @@ subroutine GFS_surface_composites_post_run (
             hflx(i)   = hflx_lnd(i)
             qss(i)    = qss_lnd(i)
             tsfc(i)   = tsfc_lnd(i)
-            hice(i)   = zero
-            cice(i)   = zero
-            tisfc(i)  = tsfc(i)
+            !hice(i)   = zero
+            !cice(i)   = zero
+            !tisfc(i)  = tsfc(i)
           elseif (islmsk(i) == 0) then
             zorl(i)   = zorl_ocn(i)
             cd(i)     = cd_ocn(i)
@@ -441,7 +441,7 @@ subroutine GFS_surface_composites_post_run (
             fm10(i)   = fm10_ocn(i)
             fh2(i)    = fh2_ocn(i)
            !tsurf(i)  = tsurf_ocn(i)
-            tsfco(i)  = tsfc_ocn(i)
+            tsfco(i)  = tsfc_ocn(i) ! over lake (and ocean when uncoupled)
             cmm(i)    = cmm_ocn(i)
             chh(i)    = chh_ocn(i)
             gflx(i)   = gflx_ocn(i)
@@ -453,10 +453,10 @@ subroutine GFS_surface_composites_post_run (
             hflx(i)   = hflx_ocn(i)
             qss(i)    = qss_ocn(i)
             tsfc(i)   = tsfc_ocn(i)
-            hice(i)   = zero
-            cice(i)   = zero
-            tisfc(i)  = tsfc(i)
-          else
+            !hice(i)   = zero
+            !cice(i)   = zero
+            !tisfc(i)  = tsfc(i)
+          else ! islmsk(i) == 2
             zorl(i)   = zorl_ice(i)
             cd(i)     = cd_ice(i)
             cdq(i)    = cdq_ice(i)
@@ -468,30 +468,42 @@ subroutine GFS_surface_composites_post_run (
             fm10(i)   = fm10_ice(i)
             fh2(i)    = fh2_ice(i)
            !tsurf(i)  = tsurf_ice(i)
+            if (.not. flag_cice(i)) then
+              tisfc(i) = tice(i) ! over lake ice (and sea ice when uncoupled)
+            endif
             cmm(i)    = cmm_ice(i)
             chh(i)    = chh_ice(i)
             gflx(i)   = gflx_ice(i)
             ep1d(i)   = ep1d_ice(i)
             weasd(i)  = weasd_ice(i)
             snowd(i)  = snowd_ice(i)
+           !tprcp(i)  = cice(i)*tprcp_ice(i) + (one-cice(i))*tprcp_ocn(i)
             qss(i)    = qss_ice(i)
-            if (flag_cice(i)) then    ! this was already done for lake ice in sfc_sice
-              txi = cice(i)
-              txo = one - txi
-              evap(i) = txi * evap_ice(i) + txo * evap_ocn(i)
-              hflx(i) = txi * hflx_ice(i) + txo * hflx_ocn(i)
-              tsfc(i) = txi * tsfc_ice(i) + txo * tsfc_ocn(i)
-            else
-              evap(i)  = evap_ice(i)
-              hflx(i)  = hflx_ice(i)
-              tsfc(i)  = tsfc_ice(i)
-              tisfc(i) = tice(i)
-            endif
+            evap(i)   = evap_ice(i)
+            hflx(i)   = hflx_ice(i)
+            qss(i)    = qss_ice(i)
+            tsfc(i)   = tsfc_ice(i)
           endif
 
           zorll(i) = zorl_lnd(i)
           zorlo(i) = zorl_ocn(i)
 
+          if (flag_cice(i) .and. wet(i)) then ! this was already done for lake ice in sfc_sice
+            txi = cice(i)
+            txo = one - txi
+            evap(i) = txi * evap_ice(i) + txo * evap_ocn(i)
+            hflx(i) = txi * hflx_ice(i) + txo * hflx_ocn(i)
+            tsfc(i) = txi * tsfc_ice(i) + txo * tsfc_ocn(i)
+          else
+            if (islmsk(i) == 2) then
+              tisfc(i) = tice(i)
+            else ! over open ocean or land (no ice fraction)
+              hice(i)  = zero
+              cice(i)  = zero
+              tisfc(i) = tsfc(i)
+            endif
+          endif
+
         enddo
 
       endif ! if (frac_grid)