Merge pull request NCAR#378 from SMoorthi-emc/SM_Jan102020

Add RAS to CCPP, update SHOC and MG2/3 - this PR will break b4b reproducibility for several physics options between IPD and CCPP. This has been discussed and agreed upon.

physics/GFS_DCNV_generic.F90

@@ -104,7 +104,7 @@ end subroutine GFS_DCNV_generic_post_finalize
 !!
     subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, ras, cscnv, do_ca,               &
       isppt_deep, frain, rain1, dtf, cld1d, save_u, save_v, save_t, save_qv, gu0, gv0, gt0,          &
-      gq0_water_vapor, ud_mf, dd_mf, dt_mf, con_g, clw_ice, clw_liquid, npdf3d, num_p3d, ncnvcld3d,  &
+      gq0_water_vapor, ud_mf, dd_mf, dt_mf, con_g,  npdf3d, num_p3d, ncnvcld3d,                      &
       rainc, cldwrk, dt3dt, dq3dt, du3dt, dv3dt, upd_mf, dwn_mf, det_mf,                             &
       cnvw, cnvc, cnvw_phy_f3d, cnvc_phy_f3d,                                                        &
       cape, tconvtend, qconvtend, uconvtend, vconvtend, errmsg, errflg)
@@ -122,7 +122,6 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, ras, cscnv, do_c
       real(kind=kind_phys), dimension(im,levs), intent(in) :: gu0, gv0, gt0, gq0_water_vapor
       real(kind=kind_phys), dimension(im,levs), intent(in) :: ud_mf, dd_mf, dt_mf
       real(kind=kind_phys), intent(in) :: con_g
-      real(kind=kind_phys), dimension(im,levs), intent(in) :: clw_ice, clw_liquid
       integer, intent(in) :: npdf3d, num_p3d, ncnvcld3d
 
       real(kind=kind_phys), dimension(im), intent(inout) :: rainc, cldwrk
@@ -151,7 +150,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, ras, cscnv, do_c
       if (.not. ras .and. .not. cscnv) then
         if(do_ca) then
           do i=1,im
-            cape(i)=cld1d(i)
+            cape(i) = cld1d(i)
           enddo
         endif
         if (npdf3d == 3 .and. num_p3d == 4) then
@@ -186,13 +185,13 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, ras, cscnv, do_c
           do k=1,levs
             do i=1,im
               dt3dt(i,k) = dt3dt(i,k) + (gt0(i,k)-save_t(i,k)) * frain
-!              dq3dt(i,k) = dq3dt(i,k) + (gq0_water_vapor(i,k)-save_qv(i,k)) * frain
+!             dq3dt(i,k) = dq3dt(i,k) + (gq0_water_vapor(i,k)-save_qv(i,k)) * frain
               du3dt(i,k) = du3dt(i,k) + (gu0(i,k)-save_u(i,k)) * frain
               dv3dt(i,k) = dv3dt(i,k) + (gv0(i,k)-save_v(i,k)) * frain
 
-!              upd_mf(i,k)  = upd_mf(i,k)  + ud_mf(i,k) * (con_g*frain)
-!              dwn_mf(i,k)  = dwn_mf(i,k)  + dd_mf(i,k) * (con_g*frain)
-!              det_mf(i,k)  = det_mf(i,k)  + dt_mf(i,k) * (con_g*frain)
+!             upd_mf(i,k)  = upd_mf(i,k)  + ud_mf(i,k) * (con_g*frain)
+!             dwn_mf(i,k)  = dwn_mf(i,k)  + dd_mf(i,k) * (con_g*frain)
+!             det_mf(i,k)  = det_mf(i,k)  + dt_mf(i,k) * (con_g*frain)
             enddo
           enddo
         endif ! if (ldiag3d)

physics/GFS_DCNV_generic.meta

@@ -377,24 +377,6 @@
   kind = kind_phys
   intent = in
   optional = F
-[clw_ice]
-  standard_name = ice_water_mixing_ratio_convective_transport_tracer
-  long_name = moist (dry+vapor, no condensates) mixing ratio of ice water in the convectively transported tracer array
-  units = kg kg-1
-  dimensions = (horizontal_dimension,vertical_dimension)
-  type = real
-  kind = kind_phys
-  intent = in
-  optional = F
-[clw_liquid]
-  standard_name = cloud_condensed_water_mixing_ratio_convective_transport_tracer
-  long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) in the convectively transported tracer array
-  units = kg kg-1
-  dimensions = (horizontal_dimension,vertical_dimension)
-  type = real
-  kind = kind_phys
-  intent = in
-  optional = F
 [npdf3d]
   standard_name = number_of_3d_arrays_associated_with_pdf_based_clouds
   long_name = number of 3d arrays associated with pdf based clouds/mp

physics/GFS_MP_generic.F90

@@ -159,7 +159,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
       onebg = one/con_g
 
       do i = 1, im
-          rain(i) = rainc(i) + frain * rain1(i) ! time-step convective plus explicit
+        rain(i) = rainc(i) + frain * rain1(i) ! time-step convective plus explicit
       enddo
 
 !> - If requested (e.g. Zhao-Carr MP scheme), call calpreciptype() to calculate dominant 
@@ -211,20 +211,12 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
 
       if (cal_pre) then       ! hchuang: add dominant precipitation type algorithm
 !
-        call calpreciptype (kdt, nrcm, im, ix, levs, levs+1,           &
-                            rann, xlat, xlon, gt0,    &
-                            gq0(:,:,1), prsl, prsi,        &
-                            rain, phii, tsfc,           &  !input
-                            domr, domzr, domip, doms)                           ! output
+        call calpreciptype (kdt, nrcm, im, ix, levs, levs+1, &
+                            rann, xlat, xlon, gt0,           &
+                            gq0(:,:,1), prsl, prsi,          &
+                            rain, phii, tsfc,                &  ! input
+                            domr, domzr, domip, doms)           ! output
 !
-!        if (lprnt) print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS '
-!     &,DOMR(ipr),DOMZR(ipr),DOMIP(ipr),DOMS(ipr)
-!        do i=1,im
-!         if (abs(xlon(i)*57.29578-114.0) .lt. 0.2  .and.
-!     &    abs(xlat(i)*57.29578-40.0) .lt. 0.2)
-!     &    print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS ',
-!     &    DOMR(i),DOMZR(i),DOMIP(i),DOMS(i)
-!       end do
 !       HCHUANG: use new precipitation type to decide snow flag for LSM snow accumulation
 
         if (imp_physics /= imp_physics_gfdl .and. imp_physics /= imp_physics_thompson) then
@@ -270,7 +262,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
           do k=1,levs
             do i=1,im
               dt3dt(i,k) = dt3dt(i,k) + (gt0(i,k)-save_t(i,k)) * frain
-!              dq3dt(i,k) = dq3dt(i,k) + (gq0(i,k,1)-save_qv(i,k)) * frain
+!             dq3dt(i,k) = dq3dt(i,k) + (gq0(i,k,1)-save_qv(i,k)) * frain
             enddo
           enddo
         endif
@@ -281,7 +273,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
       do k = 1, levs-1
         do i = 1, im
           if (prsl(i,k) > p850 .and. prsl(i,k+1) <= p850) then
-            t850(i) = gt0(i,k) - (prsl(i,k)-p850) / &
+            t850(i) = gt0(i,k) - (prsl(i,k)-p850) /  &
                       (prsl(i,k)-prsl(i,k+1)) *      &
                       (gt0(i,k)-gt0(i,k+1))
           endif
@@ -299,7 +291,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
 ! determine convective rain/snow by surface temperature
 ! determine large-scale rain/snow by rain/snow coming out directly from MP
 
-        if (lsm/=lsm_ruc) then
+        if (lsm /= lsm_ruc) then
           do i = 1, im
             !tprcp(i)  = max(0.0, rain(i) )! clu: rain -> tprcp ! DH now lines 245-250
             srflag(i) = 0.                     ! clu: default srflag as 'rain' (i.e. 0)
@@ -326,7 +318,8 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
           enddo
         endif ! lsm==lsm_ruc
       elseif( .not. cal_pre) then
-        if (imp_physics == imp_physics_mg) then              ! MG microphysics
+        if (imp_physics == imp_physics_mg) then          ! MG microphysics
+          tem = con_day / (dtp * con_p001)               ! mm / day
           do i=1,im
             tprcp(i)  = max(0.0, rain(i) )     ! clu: rain -> tprcp
             if (rain(i)*tem > rainmin) then
@@ -355,7 +348,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
 
       if (cplchm) then
         do i = 1, im
-             rainc_cpl(i) = rainc_cpl(i) + rainc(i)
+          rainc_cpl(i) = rainc_cpl(i) + rainc(i)
         enddo
       endif
 
@@ -374,8 +367,6 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
         do i=1,im
           pwat(i) = pwat(i) + del(i,k)*(gq0(i,k,1)+work1(i))
         enddo
-!     if (lprnt .and. i == ipr) write(0,*)' gq0=',
-!    &gq0(i,k,1),' qgrs=',qgrs(i,k,1),' work2=',work2(i),' k=',k
       enddo
       do i=1,im
         pwat(i) = pwat(i) * onebg

physics/GFS_PBL_generic.F90

@@ -81,10 +81,10 @@ end subroutine GFS_PBL_generic_pre_finalize
 !!
       subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac,                       &
         ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc,                 &
-        ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm,                   &
+        ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm,          &
         imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,           &
-        imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, cplchm, ltaerosol, hybedmf, do_shoc,      &
-        satmedmf, qgrs, vdftra, errmsg, errflg)
+        imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, cplchm, ltaerosol, &
+        hybedmf, do_shoc, satmedmf, qgrs, vdftra, errmsg, errflg)
 
       use machine,                only : kind_phys
       use GFS_PBL_generic_common, only : set_aerosol_tracer_index
@@ -99,11 +99,11 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac,
       integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires
       logical, intent(in) :: cplchm, ltaerosol, hybedmf, do_shoc, satmedmf
 
-      real(kind=kind_phys), dimension(im, levs, ntrac), intent(in) :: qgrs
+      real(kind=kind_phys), dimension(im, levs, ntrac),  intent(in)    :: qgrs
       real(kind=kind_phys), dimension(im, levs, nvdiff), intent(inout) :: vdftra
 
       character(len=*), intent(out) :: errmsg
-      integer, intent(out) :: errflg
+      integer,          intent(out) :: errflg
 
       !local variables
       integer :: i, k, kk, k1, n
@@ -331,6 +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
       integer :: i, k, kk, k1, n
       real(kind=kind_phys) :: tem, tem1, rho
 
@@ -498,13 +499,13 @@ 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) == ceanfrac(i)) then ! 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
-            if (wet(i)) then                   ! use stress_ocean from sfc_diff for opw component at mixed point
+            if (fice(i) == oceanfrac(i)) then ! 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))
@@ -518,7 +519,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
                 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
+              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)
@@ -530,6 +531,12 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
             dvsfc_cpl (i) = dvsfc_cpl(i) + dvsfci_cpl(i) * dtf
             dtsfc_cpl (i) = dtsfc_cpl(i) + dtsfci_cpl(i) * dtf
             dqsfc_cpl (i) = dqsfc_cpl(i) + dqsfci_cpl(i) * dtf
+!
+          else
+            dusfc_cpl(i) = huge
+            dvsfc_cpl(i) = huge
+            dtsfc_cpl(i) = huge
+            dqsfc_cpl(i) = huge
 !!
           endif ! Ocean only, NO LAKES
         enddo
@@ -547,10 +554,6 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
           dtsfci_diag(i) = dtsfc1(i)
           dqsfci_diag(i) = dqsfc1(i)
         enddo
-  !       if (lprnt) then
-  !         write(0,*)' dusfc=',dusfc(ipr),' dusfc1=',dusfc1(ipr),' dtf=',
-  !    &     dtf,' kdt=',kdt,' lat=',lat
-  !       endif
 
         if (ldiag3d) then
           if (lsidea) then
@@ -565,9 +568,9 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
           endif
           do k=1,levs
             do i=1,im
-              du3dt_PBL(i,k) = du3dt_PBL(i,k) + dudt(i,k) * dtf
+              du3dt_PBL(i,k)  = du3dt_PBL(i,k)  + dudt(i,k) * dtf
               du3dt_OGWD(i,k) = du3dt_OGWD(i,k) - dudt(i,k) * dtf
-              dv3dt_PBL(i,k) = dv3dt_PBL(i,k) + dvdt(i,k) * dtf
+              dv3dt_PBL(i,k)  = dv3dt_PBL(i,k)  + dvdt(i,k) * dtf
               dv3dt_OGWD(i,k) = dv3dt_OGWD(i,k) - dvdt(i,k) * dtf
             enddo
           enddo

physics/GFS_SCNV_generic.F90

@@ -52,6 +52,7 @@ subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, gt0, gq0_water_vapor, &
 
     end subroutine GFS_SCNV_generic_pre_run
 
+
     end module GFS_SCNV_generic_pre
 
     module GFS_SCNV_generic_post

physics/GFS_phys_time_vary.fv3.F90

@@ -260,7 +260,7 @@ subroutine GFS_phys_time_vary_init (Data, Model, Interstitial, nthrds, errmsg, e
          do j = 1,Model%ny
            do i = 1,Model%nx
              ix = ix + 1
-             if (ix .gt. Model%blksz(nb)) then
+             if (ix > Model%blksz(nb)) then
                ix = 1
                nb = nb + 1
              endif