Merge pull request #1023 from fvitt/waccmx_corrections

cam6_4_004: Corrections for WACCMX

doc/ChangeLog

@@ -1,5 +1,77 @@
 ===============================================================
 
+Tag name: cam6_4_004
+Originator(s): fvitt
+Date: 29 Jun 2024
+One-line Summary: Misc corrections for WACCM-X
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1023
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+ Implement corrections to:
+ - geometric height calculations (issue #987)
+ - thermosphere heating diagnostics (issue #1013)
+ - DTVKE vertical diffustion diagnostic
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: cacraigucar, nusbaume
+
+List all files eliminated: N/A
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+M       src/ionosphere/waccmx/ionosphere_interface.F90
+ - Hanli's formulation for geometric height
+
+M       src/physics/cam/vertical_diffusion.F90
+ - correction to DTVKE diagnostic
+
+M       src/physics/waccmx/ion_electron_temp.F90
+ - corrections to thermosphere heating diagnostics (issue #1013)
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam:
+  FAIL ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s
+  DIFF SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h
+  FAIL SMS_D_Ln9.T42_T42.FSCAM.derecho_intel.cam-outfrq9s
+  PEND SMS_D_Ln9_P1280x1.ne0ARCTICne30x4_ne0ARCTICne30x4_mt12.FHIST.derecho_intel.cam-outfrq9s
+  PEND SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s
+ - pre-existing failures
+
+  DIFF ERC_D_Ln9.f19_f19_mg17.QPX2000.derecho_intel.cam-outfrq3s
+  DIFF ERS_Ln9.f09_f09_mg17.FX2000.derecho_intel.cam-outfrq9s
+  DIFF ERS_Ln9.f19_f19_mg17.FXSD.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f19_f19_mg17.FXHIST.derecho_intel.cam-outfrq9s_amie
+  DIFF SMS_D_Ln9.ne16_ne16_mg17.QPX2000.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.ne16pg3_ne16pg3_mg17.FX2000.derecho_intel.cam-outfrq9s
+ - expected baseline failures in waccmx due to corrections in diagnostices
+
+izumi/nag/aux_cam:
+  FAIL DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae
+ - pre-existing failure
+
+  DIFF SMS_D_Ln3.ne5pg3_ne5pg3_mg37.QPX2000.izumi_nag.cam-outfrq3s
+ - expected baseline failure in waccmx due to corrections in diagnostices
+
+izumi/gnu/aux_cam: All PASS
+
+Summarize any changes to answers: bit-for-bit unchanged
+
+===============================================================
+===============================================================
+
 Tag name: cam6_4_003
 Originator(s): adamrher, jet
 Date: Thu Jun 28, 2024
@@ -206,7 +278,7 @@ Github PR URL: https://github.com/ESCOMP/CAM/pull/1028
 
 Purpose of changes (include the issue number and title text for each relevant GitHub issue):
 
-#813 - Introduce "-phys cam7" and remove "-phys cam_dev".  
+#813 - Introduce "-phys cam7" and remove "-phys cam_dev".
 https://github.com/ESCOMP/CAM/issues/813
 
 - The compset tokens CAM%DEV are replaced by CAM70
@@ -231,14 +303,14 @@ Issue #1039 - Change transient ne30np4 cam tests to ne30np4.pg3 #1039
 https://github.com/ESCOMP/CAM/issues/1039
 
 - ERP_Ln9.ne30_ne30_mg17.FCnudged.derecho_intel.cam-outfrq9s
-  changed to 
+  changed to
   ERP_Ln9.ne30pg3_ne30pg3_mg17.FCnudged.derecho_intel.cam-outfrq9s
 
 . resolves #813
 . resolves #1033
 . resolves #1038
 . resolves #1039
- 
+
 Describe any changes made to build system:
 
 . The physics package name 'cam_dev' is replaced by 'cam7'
@@ -316,7 +388,7 @@ bld/namelist_files/namelist_defaults_cam.xml
 
 bld/namelist_files/namelist_definition.xml
 . update description of do_clubb_sgs to indicate that it is not user
-  settable. 
+  settable.
 
 cime_config/SystemTests/mgp.py
 . change 'cam_dev' to 'cam7'
@@ -325,7 +397,7 @@ cime_config/config_component.xml
 cime_config/config_compsets.xml
 . change 'CAM%DEV' to 'CAM70'
 . modify compset matching so that %LT and %MT are only matched for CAM70
-  physics. 
+  physics.
 . remove %GHGMAM4 modifier (default chemistry set in configure)
 . F2000dev, FCLTHIST, FCMTHIST - change CLM50 to CLM51.  CLM no longer supports CLM50
   with CAM70 physics.

src/ionosphere/waccmx/ionosphere_interface.F90

@@ -20,6 +20,7 @@ module ionosphere_interface
    use pio,                 only: var_desc_t
    use perf_mod,            only: t_startf, t_stopf
    use epotential_params,   only: epot_active, epot_crit_colats
+   use shr_const_mod,  only: SHR_CONST_REARTH ! meters
 
    implicit none
 
@@ -95,6 +96,8 @@ module ionosphere_interface
 
    integer :: mag_nlon=0, mag_nlat=0, mag_nlev=0, mag_ngrid=0
 
+   real(r8), parameter :: rearth_inv = 1._r8/SHR_CONST_REARTH ! /meters
+
  contains
 
    !---------------------------------------------------------------------------
@@ -480,7 +483,6 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
       ! Gridded component call
       use edyn_grid_comp, only: edyn_grid_comp_run2
       use shr_assert_mod, only: shr_assert_in_domain
-      use shr_const_mod,  only: SHR_CONST_REARTH ! meters
 
       ! - pull some fields from pbuf and dyn_in
       ! - invoke ionosphere/electro-dynamics coupling
@@ -546,7 +548,6 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
       real(r8)          :: r8tmp
       real(r8), pointer :: tempm(:,:) => null() ! Temp midpoint field for outfld
       real(r8), pointer :: tempi(:,:) => null() ! Temp interface field for outfld
-      real(r8), parameter :: rearth_inv = 1._r8/SHR_CONST_REARTH ! /meters
       real(r8), parameter :: n2min = 1.e-6_r8  ! lower limit of N2 mixing ratios
 
       character(len=*), parameter :: subname = 'ionosphere_run2'
@@ -740,8 +741,8 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
                   ! Might need geometric height on midpoints for output
                   !------------------------------------------------------------
                   if (hist_fld_active('Z3GM')) then
-                     r8tmp = phys_state(lchnk)%zm(i, k) + phis(i)*rga
-                     tempm(i, k) = r8tmp * (1._r8 + (r8tmp * rearth_inv))
+                     ! geometric altitude (meters above sea level)
+                     tempm(i,k) = geometric_hgt(zgp=phys_state(lchnk)%zm(i,k), zsf=phis(i)*rga)
                   end if
                   ! physics state fields on interfaces (but only to pver)
                   zi_blck(k, j) = phys_state(lchnk)%zi(i, k) + phis(i)*rga
@@ -750,9 +751,9 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
                   !------------------------------------------------------------
                   ! Note: zht is pver instead of pverp because dynamo does not
                   !       use bottom interface
-                  hi_blck(k, j) = zi_blck(k, j) * (1._r8 + (zi_blck(k, j) * rearth_inv))
+                  hi_blck(k,j) = geometric_hgt(zgp=phys_state(lchnk)%zi(i,k), zsf=phis(i)*rga)
                   if (hist_fld_active('Z3GMI')) then
-                     tempi(i, k) = hi_blck(k, j)
+                     tempi(i,k) = hi_blck(k, j)
                   end if
                   omega_blck(k, j) = phys_state(lchnk)%omega(i, k)
                   tn_blck(k, j)    = phys_state(lchnk)%t(i, k)
@@ -1164,5 +1165,24 @@ end subroutine ionosphere_alloc
 
    !==========================================================================
 
+   ! calculates geometric height (meters above sea level)
+   pure function geometric_hgt( zgp, zsf ) result(zgm)
+
+     real(r8), intent(in) :: zgp ! geopotential height (m)
+     real(r8), intent(in) :: zsf ! surface height above sea level (m)
+     real(r8) :: zgm ! geometric height above sea level (m)
+
+     real(r8) :: tmp
+
+     ! Hanli's formulation:
+     ! Z_gm = 1/(1 - (1+Zs/r) * Z_gp/r) * (Zs + (1+Zs/r) * Z_gp)
+     ! Z_gm: geometric height
+     ! Zs: Surface height
+     ! Z_gp: model calculated geopotential height (zm and zi in the model)
+
+     tmp = 1._r8+zsf*rearth_inv
+     zgm = (zsf + tmp*zgp) / (1._r8 - tmp*zgp*rearth_inv)
+
+   end function geometric_hgt
 
 end module ionosphere_interface

src/physics/cam/vertical_diffusion.F90

@@ -1554,7 +1554,7 @@ subroutine vertical_diffusion_tend( &
   call outfld( 'KVT'          , kvt,                       pcols, lchnk )
   call outfld( 'KVM'          , kvm,                       pcols, lchnk )
   call outfld( 'CGS'          , cgs,                       pcols, lchnk )
-  dtk(:ncol,:) = dtk(:ncol,:) / cpair              ! Normalize heating for history
+  dtk(:ncol,:) = dtk(:ncol,:) / cpair / ztodt      ! Normalize heating for history
   call outfld( 'DTVKE'        , dtk,                       pcols, lchnk )
   dtk(:ncol,:) = ptend%s(:ncol,:) / cpair          ! Normalize heating for history using dtk
   call outfld( 'DTV'          , dtk,                       pcols, lchnk )

src/physics/waccmx/ion_electron_temp.F90

@@ -34,6 +34,7 @@ module ion_electron_temp
   use spmd_utils,       only : masterproc
   use cam_logfile,      only : iulog ! Output unit
   use ionos_state_mod,  only : ionos_state
+  use air_composition,only : cpairv
 
   implicit none
 
@@ -135,16 +136,16 @@ subroutine ion_electron_temp_init(pbuf2d)
     !-------------------------------------------------------------------------------
     !  Add history variables for ionosphere
     !-------------------------------------------------------------------------------
-    call addfld ('QIonElec' ,(/ 'lev' /), 'I', 'K/s', 'Electron Ion Thermal Heating Rate')
+    call addfld ('QIonElec' ,(/ 'lev' /), 'I', 'K sec-1', 'Electron Ion Thermal Heating Rate')
     call addfld ('TElec&IC'     ,(/ 'lev' /), 'I', 'K',      'Electron Temperature')
     call addfld ('TIon&IC'      ,(/ 'lev' /), 'I', 'K',      'Ion Temperature')
     call addfld ('TElec'        ,(/ 'lev' /), 'I', 'K',      'Electron Temperature')
     call addfld ('TIon'         ,(/ 'lev' /), 'I', 'K',      'Ion Temperature')
     call addfld ('ElecColDens'  ,horiz_only , 'I', 'TECU',   'Electron Column Density')
     if (.not.steady_state_ion_elec_temp) then
-       call addfld ('QIN'          ,(/ 'lev' /), 'I', 'J/kg/s', 'Ion-neutral Heating')
-       call addfld ('QEN'          ,(/ 'lev' /), 'I', ' ',      'Electron-neutral Heating')
-       call addfld ('QEI'          ,(/ 'lev' /), 'I', ' ',      'Electron-ion Heating')
+       call addfld ('QIN'          ,(/ 'lev' /), 'I', 'K sec-1','Ion-neutral Heating Rate')
+       call addfld ('QEN'          ,(/ 'lev' /), 'I', 'K sec-1','Electron-neutral Heating Rate')
+       call addfld ('QEI'          ,(/ 'lev' /), 'I', 'K sec-1','Electron-ion Heating Rate')
        call addfld ('LOSS_g3'      ,(/ 'lev' /), 'I', ' ',      'Loss Term g3')
        call addfld ('LOSS_EI'      ,(/ 'lev' /), 'I', ' ',      'Loss Term EI')
        call addfld ('LOSS_IN'      ,(/ 'lev' /), 'I', ' ',      'Loss Term IN')
@@ -334,7 +335,6 @@ end subroutine ion_electron_temp_inidat
 
   subroutine ion_electron_temp_tend(state, ptend, pbuf, ztodt)
 
-    use air_composition,     only: cpairv
     !-------------------------------------------------------------------------------------
     ! Calculate dry static energy and O+ tendency for extended ionosphere simulation
     !-------------------------------------------------------------------------------------
@@ -1037,9 +1037,9 @@ subroutine update_teti(state, dSETendIn, dSETendOut, ztodt, istate, tE, tI, teTi
     real(r8), dimension(pcols,pver)     :: delZ         ! Delta z: midpoints
 
     real(r8), dimension(pcols,pver)     :: qjoule       ! joule heating
-    real(r8), dimension(pcols,pver)     :: qen          ! electron-neutral heating
-    real(r8), dimension(pcols,pver)     :: qei          ! electron-ion Coulomb heating
-    real(r8), dimension(pcols,pver)     :: qin          ! ion-neutral heating
+    real(r8), dimension(pcols,pver)     :: qen          ! electron-neutral heating (units: ev/g/s)
+    real(r8), dimension(pcols,pver)     :: qei          ! electron-ion Coulomb heating (units: ev/g/s)
+    real(r8), dimension(pcols,pver)     :: qin          ! ion-neutral heating (units: ev/g/s)
     real(r8), dimension(pcols,pver)     :: rho          ! mass density
 
     real(r8), dimension(pcols,pver)     :: wrk2
@@ -1053,6 +1053,7 @@ subroutine update_teti(state, dSETendIn, dSETendOut, ztodt, istate, tE, tI, teTi
     logical, dimension(pcols)           :: colConv      ! flag for column converging
     logical                             :: converged    ! Flag for convergence in electron temperature
                                                         ! calculation iteration loop
+    real(r8) :: qrate(pcols,pver) ! heating rate diagnostic
 
     !---------------------------------------------------------------------------------------------------------
     !  Initialize arrays to zero and column convergence logical to .false.
@@ -1452,9 +1453,14 @@ subroutine update_teti(state, dSETendIn, dSETendOut, ztodt, istate, tE, tI, teTi
 
    dSETendOut(1:ncol,1:teTiBot) = (qei(1:ncol,1:teTiBot)+qen(1:ncol,1:teTiBot)) / sToQConv    ! J/kg/s
 
-   call outfld ('QEN', qen, pcols, lchnk)
-   call outfld ('QEI', qei, pcols, lchnk)
-   call outfld ('QIN', qin, pcols, lchnk)
+   qrate(:ncol,:) = qen(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
+   call outfld ('QEN', qrate, pcols, lchnk)
+
+   qrate(:ncol,:) = qei(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
+   call outfld ('QEI', qrate, pcols, lchnk)
+
+   qrate(:ncol,:) = qin(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
+   call outfld ('QIN', qrate, pcols, lchnk)
 
    return