Merge branch 'oksanaguba/clubb-shflx' (PR #5009)

Fix surface heat flux conversion and use more accurate theta_l formulation in CLUBB

Two fixes in this PR:

1. change how SHFLX is used in clubb to make it consistent with clubb thermodynamic variable.
2. Redefine thetal to not to use an approximation https://glossary.ametsoc.org/wiki/Liquid_water_potential_temperature

A series of runs was performed first with fix 1 and then with both 1 and 2.
Using clubb energy fixer as a measure, fix 1 and fixes 1+2 reduce the fixer.
Also, 1+2 is the best option among the two.

[non-BFB] for any active atmosphere with clubb

components/eam/src/physics/cam/cam_diagnostics.F90

@@ -320,6 +320,10 @@ subroutine diag_init()
    call addfld ('PSL',horiz_only,    'A','Pa','Sea level pressure', &
       standard_name='air_pressure_at_mean_sea_level')
 
+
+   call addfld ('fixerCLUBB',horiz_only,    'A','J/m2','dTE fixed by CLUBB')
+
+
    call addfld ('T850',horiz_only,    'A','K','Temperature at 850 mbar pressure surface')
    call addfld ('T500',horiz_only,    'A','K','Temperature at 500 mbar pressure surface')
    call addfld ('T300',horiz_only,    'A','K','Temperature at 300 mbar pressure surface')

components/eam/src/physics/cam/clubb_intr.F90

@@ -1520,6 +1520,10 @@ subroutine clubb_tend_cam( &
    real(r8) :: sfc_v_diff_tau(pcols) ! Response to tau perturbation, m/s
    real(r8), parameter :: pert_tau = 0.1_r8 ! tau perturbation, Pa
 
+
+   real(r8) :: inv_exner_clubb_surf
+
+
 ! ZM gustiness equation below from Redelsperger et al. (2000)
 ! numbers are coefficients of the empirical equation
 
@@ -1840,7 +1844,20 @@ subroutine clubb_tend_cam( &
        rvm(i,k)     = state1%q(i,k,ixq)
        um(i,k)      = state1%u(i,k)
        vm(i,k)      = state1%v(i,k)
+
+#define NEWTHETAL
+#ifndef NEWTHETAL
        thlm(i,k)    = state1%t(i,k)*exner_clubb(i,k)-(latvap/cpair)*state1%q(i,k,ixcldliq)
+#else
+!NCAR
+!       thlm(i,k) = ( state1%t(i,k) &
+!                     - (latvap/cpairv(i,k,lchnk))*state1%q(i,k,ixcldliq) ) &
+!                   * inv_exner_clubb(i,k)
+
+       thlm(i,k) = ( state1%t(i,k) &
+                     - (latvap/cpair)*state1%q(i,k,ixcldliq) ) &
+                   * exner_clubb(i,k)
+#endif
 
        if (clubb_do_adv) then
           if (macmic_it .eq. 1) then
@@ -2028,6 +2045,15 @@ subroutine clubb_tend_cam( &
 
       !  Surface fluxes provided by host model
       wpthlp_sfc = real(cam_in%shf(i), kind = core_rknd)/(real(cpair, kind = core_rknd)*rho_ds_zm(1)) ! Sensible heat flux
+#if 1
+      inv_exner_clubb_surf = 1._r8/((state1%pmid(i,pver)/p0_clubb)**(rair/cpair)) !phl Option 2
+      wpthlp_sfc = wpthlp_sfc*inv_exner_clubb_surf
+#endif
+#if 0
+      inv_exner_clubb_surf = 1._r8/((state1%pint(i,pverp)/p0_clubb)**(rair/cpair)) !Peter B option
+      wpthlp_sfc = wpthlp_sfc*inv_exner_clubb_surf
+#endif
+
       wprtp_sfc  = real(cam_in%cflx(i,1), kind = core_rknd)/rho_ds_zm(1)                              ! Latent heat flux
       upwp_sfc   = real(cam_in%wsx(i), kind = core_rknd)/rho_ds_zm(1)                                 ! Surface meridional momentum flux
       vpwp_sfc   = real(cam_in%wsy(i), kind = core_rknd)/rho_ds_zm(1)                                 ! Surface zonal momentum flux
@@ -2436,7 +2462,11 @@ subroutine clubb_tend_cam( &
       wv_a = 0._r8
       wl_a = 0._r8
       do k=1,pver
+#ifdef NEWTHETAL
+         enthalpy = cpair*thlm(i,k)/exner_clubb(i,k) + latvap*rcm(i,k)
+#else
          enthalpy = cpair*((thlm(i,k)+(latvap/cpair)*rcm(i,k))/exner_clubb(i,k))
+#endif
          clubb_s(k) = enthalpy + gravit*state1%zm(i,k)+state1%phis(i)
 !         se_a(i) = se_a(i) + clubb_s(k)*state1%pdel(i,k)*invrs_gravit
          se_a(i) = se_a(i) + enthalpy * state1%pdel(i,k)*invrs_gravit
@@ -2540,6 +2570,10 @@ subroutine clubb_tend_cam( &
    enddo  ! end column loop
    call t_stopf('adv_clubb_core_col_loop')
 
+
+   call outfld('fixerCLUBB', te_a(:)-te_b(:), pcols, lchnk )
+
+
    ! Add constant to ghost point so that output is not corrupted
    if (clubb_do_adv) then
       if (macmic_it .eq. cld_macmic_num_steps) then