Add extended diagnostic output from Thompson MP, includes "Add optional scaling to RRTMGP flux adjustment" (#668), fix index bug in m_micro.F90

physics/dcyc2.f

@@ -178,7 +178,7 @@ subroutine dcyc2t3_run                                            &
      &       sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu,                   &
      &       sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd,                   &
      &       im, levs, deltim, fhswr,                                   &
-     &       dry, icy, wet,                                             &
+     &       dry, icy, wet, damp_LW_fluxadj, lfnc_k, lfnc_p0,           &
      &       minGPpres, use_LW_jacobian, sfculw, fluxlwUP_jac,          &
      &       t_lay, t_lev, p_lay, p_lev, flux2D_lwUP, flux2D_lwDOWN,    &
      &       pert_radtend, do_sppt,ca_global,                           &
@@ -213,10 +213,11 @@ subroutine dcyc2t3_run                                            &
 !     integer, intent(in) :: ipr
 !     logical lprnt
       logical, dimension(:), intent(in) :: dry, icy, wet
-      logical, intent(in) :: use_LW_jacobian, pert_radtend
+      logical, intent(in) :: use_LW_jacobian, damp_LW_fluxadj,          &
+     &     pert_radtend
       logical, intent(in) :: do_sppt,ca_global
       real(kind=kind_phys),   intent(in) :: solhr, slag, cdec, sdec,    &
-     &                                      deltim, fhswr, minGPpres
+     &     deltim, fhswr, minGPpres, lfnc_k, lfnc_p0
 
       real(kind=kind_phys), dimension(:), intent(in) ::                 &
      &      sinlat, coslat, xlon, coszen, tf, tsflw, sfcdlw,            &
@@ -253,11 +254,19 @@ subroutine dcyc2t3_run                                            &
       integer,          intent(out) :: errflg
 
 !  ---  locals:
-      integer :: i, k, nstp, nstl, it, istsun(im),iSFC
+      integer :: i, k, nstp, nstl, it, istsun(im),iSFC,iTOA
       real(kind=kind_phys) :: cns,  coszn, tem1, tem2, anginc,          &
      &                        rstl, solang, dT
       real(kind=kind_phys), dimension(im,levs+1) :: flxlwup_adj,        &
      &     flxlwdn_adj, t_lev2
+      real(kind=kind_phys) :: fluxlwnet_adj,fluxlwnet,dT_sfc,           &
+     &fluxlwDOWN_jac,lfnc,c1
+      ! Length scale for flux-adjustment scaling
+      real(kind=kind_phys), parameter ::                                &
+     &     L = 1.
+      ! Scaling factor for downwelling LW Jacobian profile.
+      real(kind=kind_phys), parameter ::                                &
+     &     gamma = 0.2
 !
 !===> ...  begin here
 !
@@ -267,9 +276,11 @@ subroutine dcyc2t3_run                                            &
 
 !     Vertical ordering?
       if (p_lev(1,1) .lt.  p_lev(1, levs)) then 
-         iSFC = levs
+         iSFC = levs + 1
+         iTOA = 1
       else
          iSFC = 1
+         iTOA = levs + 1
       endif
 
       tem1 = fhswr / deltim
@@ -375,32 +386,47 @@ subroutine dcyc2t3_run                                            &
          call cmp_tlev(im, levs, minGPpres, p_lay, t_lay, p_lev, tsfc,  &
      &        t_lev2)
 
+         ! Compute adjusted net LW flux foillowing Hogan and Bozzo 2015 (10.1002/2015MS000455)
+         ! Here we assume that the profile of the downwelling LW Jaconiam has the same shape
+         ! as the upwelling, but scaled and offset.
+         ! The scaling factor is 0.2
+         ! The profile of the downwelling Jacobian (J) is offset so that
+         !     J_dn_sfc / J_up_sfc = scaling_factor
+         !     J_dn_toa / J_up_sfc = 0
          !
-         ! Adjust up/downward fluxes (at layer interfaces).
-         !
-         do k = 1, levs+1
-            do i = 1, im
-               dT = t_lev2(i,k) - t_lev(i,k)
-               flxlwup_adj(i,k) = flux2D_lwUP(i,k) +                    &
-     &              fluxlwUP_jac(i,k)*dT
-            enddo
-         enddo
-         !
-         ! Compute new heating rate (within each layer).
-         !
-         do k = 1, levs
-            htrlw(1:im,k) =                                             &
-     &           (flxlwup_adj(1:im,k+1)  - flxlwup_adj(1:im,k) -        &
-     &           flux2D_lwDOWN(1:im,k+1) + flux2D_lwDOWN(1:im,k)) *     &
-     &           con_g / (con_cp * (p_lev(1:im,k+1) - p_lev(1:im,k)))
-         enddo
-
-         !
-         ! Add radiative heating rates to physics heating rate
-         !
-         do k = 1, levs
-            do i = 1, im
-               dtdt(i,k)  = dtdt(i,k)  + swh(i,k)*xmu(i)  + htrlw(i,k)
+         ! Optionally, the flux adjustment can be damped with height using a logistic function
+         ! fx ~ L / (1 + exp(-k*dp)), where dp = p - p0
+         ! L  = 1, fix scale between 0-1.      - Fixed
+         ! k  = 1 / pressure decay length (Pa) - Controlled by namelist
+         ! p0 = Transition pressure (Pa)       - Controlled by namelsit
+         do i = 1, im
+            c1 = fluxlwUP_jac(i,iTOA) / fluxlwUP_jac(i,iSFC)
+            dT_sfc = t_lev2(i,iSFC) - t_lev(i,iSFC)
+            do k = 1, levs
+               ! LW net flux
+               fluxlwnet = (flux2D_lwUP(i,  k+1) - flux2D_lwUP(i,  k) - &
+     &                      flux2D_lwDOWN(i,k+1) + flux2D_lwDOWN(i,k))
+               ! Downward LW Jacobian (Eq. 9)
+               fluxlwDOWN_jac = gamma *                                 &
+     &              (fluxlwUP_jac(i,k)/fluxlwUP_jac(i,iSFC) - c1) /     &
+     &              (1 - c1)
+               ! Adjusted LW net flux(Eq. 10)
+               fluxlwnet_adj = fluxlwnet + dT_sfc*                      &
+     &              (fluxlwUP_jac(i,k)/fluxlwUP_jac(i,iSFC) -           &
+     &              fluxlwDOWN_jac)
+               ! Adjusted LW heating rate
+               htrlw(i,k) = fluxlwnet_adj * con_g /                     &
+     &              (con_cp * (p_lev(i,k+1) - p_lev(i,k)))
+
+               ! Add radiative heating rates to physics heating rate. Optionally, scaled w/ height
+               ! using a logistic function
+               if (damp_LW_fluxadj) then
+                  lfnc = L / (1+exp(-(p_lev(i,k) - lfnc_p0)/lfnc_k))
+               else
+                  lfnc = 1.
+               endif
+               dtdt(i,k) = dtdt(i,k) + swh(i,k)*xmu(i) +                &
+     &              htrlw(i,k)*lfnc + (1.-lfnc)*hlw(i,k)
             enddo
          enddo
       else

physics/dcyc2.meta

@@ -370,6 +370,32 @@
   type = logical  
   intent = in
   optional = F
+[damp_LW_fluxadj]
+  standard_name = flag_to_damp_RRTMGP_LW_jacobian_flux_adjustment
+  long_name = logical flag to control RRTMGP LW calculation
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[lfnc_k]
+  standard_name = transition_pressure_length_scale_for_flux_damping
+  long_name = depth of transition layer in logistic function for LW flux adjustment damping
+  units = Pa
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[lfnc_p0]
+  standard_name = transition_pressure_for_flux_damping
+  long_name = transition pressure for LW flux adjustment damping
+  units = Pa
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
 [sfculw]
   standard_name = surface_upwelling_longwave_flux_on_radiation_time_step
   long_name = total sky sfc upward lw flux

physics/m_micro.F90

@@ -180,10 +180,7 @@ subroutine m_micro_run(   im,       lm,     flipv, dt_i           &
        real (kind=kind_phys), dimension(:,:),intent(in)  ::             &
      &                prsl_i,u_i,v_i,phil,   omega_i, QLLS_i,QILS_i,    &
      &                                       lwheat_i,swheat_i
-       real (kind=kind_phys), dimension(:,:),intent(in):: prsi_i, phii
-! GJF* These variables are conditionally allocated depending on whether the
-!     Morrison-Gettelman microphysics is used, so they must be declared 
-!     using assumed shape.
+       real (kind=kind_phys), dimension(:,0:),intent(in):: prsi_i, phii
        real (kind=kind_phys), dimension(:,:),  intent(in)  ::           &
      &       CNV_DQLDT_i, CLCN_i,     QLCN_i, QICN_i,                   &
      &       CNV_MFD_i,               cf_upi, CNV_FICE_i, CNV_NDROP_i,  &
@@ -210,9 +207,6 @@ subroutine m_micro_run(   im,       lm,     flipv, dt_i           &
        integer, dimension(:), intent(inout):: KCBL
        real (kind=kind_phys),dimension(:,:),intent(inout):: q_io, t_io,   &
      &                                             ncpl_io,ncpi_io,CLLS_io
-! GJF* These variables are conditionally allocated depending on whether the
-!     Morrison-Gettelman microphysics is used, so they must be declared 
-!     using assumed shape.
        real (kind=kind_phys),dimension(:,:),intent(inout):: rnw_io,snw_io,&
      &                                             ncpr_io, ncps_io,        &
      &                                             qgl_io,  ncgl_io
@@ -442,7 +436,7 @@ subroutine m_micro_run(   im,       lm,     flipv, dt_i           &
            END DO
          END DO
          DO K=0, LM
-           ll = lm-k+1
+           ll = lm-k
            DO I = 1,IM
              PLE(i,k)   = prsi_i(i,ll) * 0.01_kp      ! interface pressure in hPa
              zet(i,k+1) = phii(i,ll) * onebg