Merge pull request NOAA-GFDL#29 from GEOS-ESM/feature/wmputman/zero-d…

…iff-levs-updates-20200624

updates in FV3 to support levels and optional remap codes

model/fv_control.F90

@@ -338,10 +338,6 @@ module fv_control_mod
    integer :: commID, max_refinement_of_global = 1.
    integer :: gid
 
-#ifdef MAPL_MODE
-   real   :: dyn_timer, comm_timer
-   public :: dyn_timer, comm_timer
-#endif
 !---- version number -----
    character(len=128) :: version = '$Id$'
    character(len=128) :: tagname = '$Name$'

model/fv_dynamics.F90

@@ -137,17 +137,9 @@ module fv_dynamics_mod
    use boundary_mod,        only: nested_grid_BC_apply_intT
    use fv_arrays_mod,       only: fv_grid_type, fv_flags_type, fv_atmos_type, fv_nest_type, fv_diag_type, fv_grid_bounds_type
    use fv_nwp_nudge_mod,    only: do_adiabatic_init
-#ifdef MAPL_MODE
-   use fv_control_mod,      only: dyn_timer, comm_timer
-#endif
 
 implicit none
 
-#ifdef MAPL_MODE
-  ! Include the MPI library definitons:
-  include 'mpif.h'
-#endif
-
    logical :: RF_initialized = .false.
    logical :: bad_range = .false.
    real, allocatable ::  rf(:)
@@ -283,8 +275,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
       integer :: is,  ie,  js,  je
       integer :: isd, ied, jsd, jed
       real :: dt2
-      real(kind=8) :: t1, t2
-      integer :: status
 
       is  = bd%is
       ie  = bd%ie
@@ -295,9 +285,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
       jsd = bd%jsd
       jed = bd%jed
 
-      dyn_timer = 0
-      comm_timer = 0
-
 ! Empty the accumulated mass flux and courant numbers
       mfx = 0.0
       mfy = 0.0
@@ -313,10 +300,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
       rdg = -rdgas * agrav
       allocate ( dp1(isd:ied, jsd:jed, 1:npz) )
 
-#ifdef MAPL_MODE
-! Begin Dynamics timer for GEOS history processing
-      t1 = MPI_Wtime(status)
-#endif
 #ifdef MOIST_CAPPA
       allocate ( cappa(isd:ied,jsd:jed,npz) )
       call init_ijk_mem(isd,ied, jsd,jed, npz, cappa, 0.)
@@ -925,10 +908,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
                          -50., 100., bad_range)
   endif
 
-#ifdef MAPL_MODE
-  t2 = MPI_Wtime(status)
-  dyn_timer = dyn_timer + (t2-t1)
-#endif
   end subroutine fv_dynamics
 
 #ifdef USE_RF_FAST

model/fv_mapz.F90

@@ -323,11 +323,11 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 ! Note: pt at this stage is Theta_v
              if ( hydrostatic ) then
 ! Transform virtual pt to virtual Temp
-             do k=1,km
+               do k=1,km
                    do i=is,ie
                       pt(i,j,k) = pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))/(akap*(peln(i,k+1,j)-peln(i,k,j)))
                    enddo
-             enddo
+               enddo
              else
 ! Transform "density pt" to "density temp"
                do k=1,km
@@ -542,6 +542,12 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
          do i=is,ie
             pkz(i,j,k) = (pk2(i,k+1)-pk2(i,k))/(akap*(peln(i,k+1,j)-peln(i,k,j)))
          enddo
+         if (remap_pt .and. last_step .and. (.not.adiabatic) ) then
+           ! Make pt T_v
+            do i=is,ie
+               pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
+            enddo
+         endif
       enddo
    else
     ! WMP: note that this is where TE remapping non-hydrostatic is invalid and cannot be run
@@ -573,7 +579,15 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
          do k=1,km
             do i=is,ie
                pkz(i,j,k) = exp( k1k*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)) )
+! Using dry pressure for the definition of the virtual potential temperature
+!              pkz(i,j,k) = exp(k1k*log(rrg*(1.-q(i,j,k,sphum))*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum))))
             enddo
+            if ( last_step .and. (.not.adiabatic) ) then
+              ! Make pt T_v
+               do i=is,ie
+                  pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
+               enddo
+            endif
          enddo
       endif
    endif
@@ -702,7 +716,16 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 
 !$OMP do
     do j=js,je
-      if (remap_t) then
+      if (remap_te) then
+          do i=is,ie
+             te_2d(i,j) = te(i,j,1)*delp(i,j,1)
+          enddo
+          do k=2,km
+             do i=is,ie
+                te_2d(i,j) = te_2d(i,j) + te(i,j,k)*delp(i,j,k)
+             enddo
+          enddo
+      else
        if ( hydrostatic ) then
             do i=is,ie
                gz(i) = hs(i,j)
@@ -755,57 +778,6 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 #endif
            enddo   ! k-loop
        endif  ! end non-hydro
-      elseif (remap_pt) then
-         if ( hydrostatic ) then
-            do i=is,ie
-               gz(i) = hs(i,j)
-               do k=1,km
-                  gz(i) = gz(i) + cp_air*pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))
-               enddo
-            enddo
-
-            do i=is,ie
-               te_2d(i,j) = pe(i,km+1,j)*hs(i,j) - pe(i,1,j)*gz(i)
-            enddo
-            do k=1,km
-               do i=is,ie
-                  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cp_air*pt(i,j,k)*pkz(i,j,k) +   &
-                               0.25*gridstruct%rsin2(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 +  &
-                                                v(i,j,k)**2+v(i+1,j,k)**2 -  &
-                            (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j)))
-               enddo
-            enddo
-         else
-!-----------------
-! Non-hydrostatic:
-!-----------------
-           do i=is,ie
-              phis(i,km+1) = hs(i,j)
-              do k=km,1,-1
-                 phis(i,k) = phis(i,k+1) - grav*delz(i,j,k)
-              enddo
-           enddo
-           do i=is,ie
-              te_2d(i,j) = 0.
-           enddo
-           do k=1,km
-              do i=is,ie
-                 te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cv_air*pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum)) + &
-                                 0.5*(phis(i,k)+phis(i,k+1) + w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*( &
-                                 u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -  &
-                                (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j))))
-              enddo
-           enddo
-         endif
-      elseif (remap_te) then
-          do i=is,ie
-             te_2d(i,j) = te(i,j,1)*delp(i,j,1)
-          enddo
-          do k=2,km
-             do i=is,ie
-                te_2d(i,j) = te_2d(i,j) + te(i,j,k)*delp(i,j,k)
-             enddo
-          enddo
       endif
 
          do i=is,ie
@@ -869,11 +841,11 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 endif        ! end last_step check
 
 ! Note: pt at this stage is T_v
-  if ( remap_t .and. (.not.do_adiabatic_init) .and. do_sat_adj ) then
-! if ( do_sat_adj ) then
+  if ( do_sat_adj ) then
                                            call timing_on('sat_adj2')
 
-#ifdef MAPL_MODE_FIX_SMALL_COND
+!#define MAPL_MODE_FIX_SMALL_COND
+#ifdef MAPL_MODE_FIX_SMALL_COND && USE_COND
    ! fix small cloud condensates
      ! Cloud
 !$OMP do
@@ -977,7 +949,25 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 
   if ( last_step ) then
        ! Output temperature if last_step
-      if (remap_t) then
+      if (remap_te) then
+!$OMP do
+         do j=js,je
+            do i=is,ie
+               gz(i) = hs(i,j)
+            enddo
+            do k=km,1,-1
+               do i=is,ie
+                  tpe = te(i,j,k) - gz(i) - 0.25*gridstruct%rsin2(i,j)*(    &
+                        u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -  &
+                       (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j) )
+                  dlnp = rg*(peln(i,k+1,j) - peln(i,k,j))
+                  tmp = tpe / ((cp - pe(i,k,j)*dlnp/delp(i,j,k))*(1.+r_vir*q(i,j,k,sphum)) )
+                  pt(i,j,k) =  tmp + dtmp*pkz(i,j,k) / (1.+r_vir*q(i,j,k,sphum))
+                  gz(i) = gz(i) + dlnp*tmp*(1.+r_vir*q(i,j,k,sphum))
+               enddo
+            enddo           ! end k-loop
+         enddo
+      else
 !$OMP do
         do k=1,km
            do j=js,je
@@ -1009,35 +999,9 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 #endif
            enddo   ! j-loop
         enddo  ! k-loop
-      elseif (remap_pt) then
-!$OMP do
-        do k=1,km
-          do j=js,je
-            do i=is,ie
-               pt(i,j,k) = (pt(i,j,k) + dtmp)*pkz(i,j,k)/(1.+r_vir*q(i,j,k,sphum))
-            enddo
-          enddo
-        enddo
-      elseif (remap_te) then
-!$OMP do
-         do j=js,je
-            do i=is,ie
-               gz(i) = hs(i,j)
-            enddo
-            do k=km,1,-1
-               do i=is,ie
-                  tpe = te(i,j,k) - gz(i) - 0.25*gridstruct%rsin2(i,j)*(    &
-                        u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -  &
-                       (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j) )
-                  dlnp = rg*(peln(i,k+1,j) - peln(i,k,j))
-                  tmp = tpe / ((cp - pe(i,k,j)*dlnp/delp(i,j,k))*(1.+r_vir*q(i,j,k,sphum)) )
-                  pt(i,j,k) =  tmp + dtmp*pkz(i,j,k) / (1.+r_vir*q(i,j,k,sphum))
-                  gz(i) = gz(i) + dlnp*tmp*(1.+r_vir*q(i,j,k,sphum))
-               enddo
-            enddo           ! end k-loop
-         enddo
       endif
     else  ! not last_step
+     ! Top of the loop expects PT to be Theta_V
       if (remap_t) then
 !$OMP do
          do k=1,km
@@ -1060,7 +1024,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
                        (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j) )
                   dlnp = rg*(peln(i,k+1,j) - peln(i,k,j))
                   tmp = tpe / (cp - pe(i,k,j)*dlnp/delp(i,j,k))
-                  pt(i,j,k) = (tmp/pkz(i,j,k) + dtmp)
+                  pt(i,j,k) = tmp/pkz(i,j,k)
                   gz(i) = gz(i) + dlnp*tmp
                enddo
             enddo           ! end k-loop
@@ -3593,7 +3557,7 @@ subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai
         cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
      enddo
   case default
-     call mpp_error (NOTE, 'fv_mapz::moist_cv - using default cv_air')
+    !call mpp_error (NOTE, 'fv_mapz::moist_cv - using default cv_air')
      do i=is,ie 
          qd(i) = 0.
         cvm(i) = cv_air

model/lin_cloud_microphys.F90

@@ -208,7 +208,8 @@ module gfdl_cloud_microphys_mod
     real :: tau_l2v = 300. !< cloud water to water vapor (evaporation)
     real :: tau_g2v = 900. !< graupel sublimation
     real :: tau_v2g = 21600. !< graupel deposition -- make it a slow process
-
+    real :: tau_revp = 150. !< rain re-evaporation
+
     ! horizontal subgrid variability
 
     real :: dw_land = 0.20 !< base value for subgrid deviation / variability over land
@@ -308,7 +309,7 @@ module gfdl_cloud_microphys_mod
         vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max,    &
         qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi,     &
         const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
-        tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
+        tau_g2v, tau_v2g, tau_revp, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
         tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs,           &
         z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice,        &
         rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof,         &
@@ -320,7 +321,7 @@ module gfdl_cloud_microphys_mod
         vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max,    &
         qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi,     &
         const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
-        tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
+        tau_g2v, tau_v2g, tau_revp, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
         tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs,           &
         z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice,        &
         rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof,         &
@@ -1279,7 +1280,8 @@ subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac,
 
     real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
     real :: qpz, dq, dqh, tin
-
+    real :: fac_revp
+
     integer :: k
 
     do k = ktop, kbot
@@ -1328,7 +1330,9 @@ subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac,
                 t2 = tin * tin
                 evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * &
                     exp (0.725 * log (qden))) / (crevp (4) * t2 + crevp (5) * qsat * den (k))
-                evap = min (qr (k), dt * evap, dqv / (1. + lcpk (k) * dqsdt))
+                fac_revp = 1. - exp (- dt / tau_revp)
+                evap = min (qr (k), dt * fac_revp * evap, dqv / (1. + lcpk (k) * dqsdt))
+              !!evap = min (qr (k), dt * evap, dqv / (1. + lcpk (k) * dqsdt))
                 ! -----------------------------------------------------------------------
                 ! alternative minimum evap in dry environmental air
                 ! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt))