add med_enthaply_mod

cime_config/namelist_definition_drv.xml

@@ -727,6 +727,20 @@
       <value>.true.</value>
     </values>
   </entry>
+
+  <entry id="compute_enthalpy">
+    <type>logical</type>
+    <category>control</category>
+    <group>MED_attributes</group>
+    <desc>
+      Compute energy of enthalpy
+    </desc>
+    <values>
+      <value>.false.</value>
+      <value comp_ocn="mom">.true.</value>
+    </values>
+  </entry>
+
   <entry id="atm_nx" modify_via_xml="ATM_NX">
     <type>integer</type>
     <category>control</category>

mediator/CMakeLists.txt

@@ -16,7 +16,8 @@ set(SRCFILES esmFldsExchange_cesm_mod.F90 med_fraction_mod.F90
               med_phases_post_ocn_mod.F90  med_phases_ocnalb_mod.F90
               med_phases_post_atm_mod.F90 med_phases_post_ice_mod.F90
               med_phases_post_lnd_mod.F90 med_phases_post_glc_mod.F90
-              med_phases_post_rof_mod.F90 med_phases_post_wav_mod.F90)
+              med_phases_post_rof_mod.F90 med_phases_post_wav_mod.F90
+              med_enthalpy_mod.F90)
 
 foreach(FILE ${SRCFILES})
   if(EXISTS "${CASEROOT}/SourceMods/src.cmeps/${FILE}")

mediator/med.F90

@@ -662,7 +662,8 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc)
     use esmFlds, only : med_fldlist_init1, med_fld_GetFldInfo, med_fldList_entry_type
     use med_phases_history_mod, only : med_phases_history_init
     use med_methods_mod       , only : mediator_checkfornans
-
+    use med_enthalpy_mod      , only : mediator_compute_enthalpy
+
     ! input/output variables
     type(ESMF_GridComp)  :: gcomp
     type(ESMF_State)     :: importState, exportState
@@ -935,6 +936,24 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc)
     endif
 
 
+    ! Should enthalpy be calculated
+    call NUOPC_CompAttributeGet(gcomp, name="compute_enthalpy", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    if(isPresent .and. isSet) then
+       read(cvalue, *) mediator_compute_enthalpy
+    else
+       mediator_compute_enthalpy = .false.
+    endif
+    if(maintask) then
+       write(logunit,*) ' compute_enthalpy is ',mediator_compute_enthalpy
+       if(mediator_compute_enthalpy) then
+          write(logunit,*) ' Enthalpy calculation is ON'
+       else
+          write(logunit,*) ' Enthalpy calculation is OFF'
+       endif
+    endif
+
+
     if (profile_memory) call ESMF_VMLogMemInfo("Leaving "//trim(subname))
     call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO)
 

mediator/med_enthalpy_mod.F90

@@ -0,0 +1,196 @@
+module med_enthalpy_mod
+  use ESMF, only : ESMF_SUCCESS, ESMF_GridComp, ESMF_VMAllreduce, ESMF_REDUCE_SUM
+  use shr_kind_mod, only : R8=>shr_kind_r8
+  use shr_const_mod, only : tkfrz=>shr_const_tkfrz, cpfw=>shr_const_cpfw, cpice=>shr_const_cpice,&
+       cpwv=>shr_const_cpwv, cpsw=>shr_const_cpsw, pi=>shr_const_pi
+  use med_utils_mod         , only : chkerr        => med_utils_ChkErr
+  use med_methods_mod       , only : FB_fldchk     => med_methods_FB_FldChk
+  use med_methods_mod       , only : FB_GetFldPtr  => med_methods_FB_GetFldPtr
+  use med_methods_mod       , only : fldbun_getdata1d => med_methods_FB_getdata1d
+  use med_internalstate_mod, only : compocn, compatm, comprof, InternalState
+  use med_internalstate_mod , only : logunit, maintask
+  use perf_mod, only : t_startf, t_stopf
+
+
+  implicit none
+  public :: med_compute_enthalpy
+  logical, public :: mediator_compute_enthalpy = .false.
+
+  real(r8) :: global_htot_corr(1) = 0._r8
+  character(*), parameter :: u_FILE_u  = &
+       __FILE__
+contains
+  real(r8) function med_enthalpy_get_global_htot_corr()
+    ! Just return the latest global_htot_corr
+    med_enthalpy_get_global_htot_corr = global_htot_corr(1)
+  end function med_enthalpy_get_global_htot_corr
+
+  subroutine med_compute_enthalpy(is_local, rc)
+    type(InternalState), intent(in) :: is_local
+    integer, intent(out) :: rc
+
+    ! local variables
+
+    real(r8), pointer :: tocn(:), rain(:), snow(:), rofl(:), rofi(:), evap(:)
+    real(r8), pointer :: rainl(:), rainc(:), tbot(:)
+    real(r8), pointer :: snowl(:), snowc(:), ofrac(:)
+    real(r8), pointer :: hrain(:), hsnow(:), hevap(:), hcond(:), hrofl(:), hrofi(:)
+    real(r8), allocatable :: hcorr(:)
+    real(r8), pointer :: areas(:)
+    real(r8), parameter :: glob_area_inv = 1._r8 / (4._r8 * pi)
+    real(r8) :: local_htot_corr(1)
+
+    integer :: n, nmax
+    character(len=*), parameter:: subname = "med_compute_enthalpy"
+
+    call t_startf(subname)
+    rc = ESMF_SUCCESS
+
+    call FB_GetFldPtr(is_local%wrap%FBImp(compocn,compocn), 'So_t', tocn, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    nmax = size(tocn)
+
+    call FB_GetFldPtr(is_local%wrap%FBImp(compatm, compocn), 'Sa_tbot', tbot, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Faxa_rain', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Faxa_rain' , rain, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       call FB_GetFldPtr(is_local%wrap%FBImp(compatm, compatm), 'Faxa_rainl', rainl, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       call FB_GetFldPtr(is_local%wrap%FBImp(compatm, compatm), 'Faxa_rainl', rainc, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       allocate(rain(nmax))
+       rain = rainl + rainc
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrain', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrain', hrain, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hrain(nmax))
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_evap', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_evap' , evap, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       call FB_GetFldPtr(is_local%wrap%FBExp(compatm), 'Faxx_evap' , evap, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hevap', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hevap', hevap, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hevap(nmax))
+    endif
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hcond', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hcond', hcond, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hcond(nmax))
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Faxa_snow', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Faxa_snow' , snow, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       call FB_GetFldPtr(is_local%wrap%FBImp(compatm,compatm), 'Faxa_snowc' , snowc, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       call FB_GetFldPtr(is_local%wrap%FBImp(compatm,compatm), 'Faxa_snowl' , snowl, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       allocate(snow(nmax))
+       snow = snowc + snowl
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hsnow', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hsnow', hsnow, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hsnow(nmax))
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_rofl', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_rofl' , rofl, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       call FB_GetFldPtr(is_local%wrap%FBImp(comprof, comprof), 'Forr_rofl', rofl, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    endif
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofl', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrofl', hrofl, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hrofl(nmax))
+    endif
+
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_rofi', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_rofi' , rofi, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       call FB_GetFldPtr(is_local%wrap%FBImp(comprof, comprof), 'Forr_rofi', rofi, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    endif
+    if(FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofi', rc)) then
+       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrofi', hrofi, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       allocate(hrofi(nmax))
+    endif
+
+    call fldbun_getdata1d(is_local%wrap%FBImp(compocn,compocn), 'So_omask', ofrac, rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    do n=1,nmax
+       ! for F cases (docn) tocn is non-zero over land and so ofrac must be included
+       ! so that only ocean points are included in calculation
+       ! Need max to ensure that will not have an enthalpy contribution if the water is below 0C
+       hrain(n)  = max((tbot(n) - tkfrz), 0._r8) * rain(n)  * cpfw  * ofrac(n)
+       hsnow(n)  = min((tbot(n) - tkfrz), 0._r8) * snow(n)  * cpice * ofrac(n)
+       hevap(n)  = (tocn(n) - tkfrz) * min(evap(n), 0._r8)  * cpwv  * ofrac(n)
+       hcond(n)  = (tocn(n) - tkfrz) * max(evap(n), 0._r8)  * cpwv  * ofrac(n)
+       hrofl(n)  = max((tocn(n) - tkfrz), 0._r8) * rofl(n)  * cpfw  * ofrac(n)
+       hrofi(n)  = min((tocn(n) - tkfrz), 0._r8) * rofi(n)  * cpice * ofrac(n)
+       ! GMM - note change in hcond
+    end do
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Faxa_rain', rc)) deallocate(rain)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Faxa_snow', rc)) deallocate(snow)
+
+
+    ! Determine enthalpy correction factor that will be added to the sensible heat flux sent to the atm
+    ! Areas here in radians**2 - this is an instantaneous snapshot that will be sent to the atm - only
+    ! need to calculate this if data is sent back to the atm
+
+    if (FB_fldchk(is_local%wrap%FBExp(compatm), 'Faxx_sen', rc=rc)) then
+       allocate(hcorr(nmax))
+       areas => is_local%wrap%mesh_info(compocn)%areas
+       do n = 1,nmax
+          hcorr(n) = (hrain(n) + hsnow(n) + hcond(n) + hevap(n) + hrofl(n) + hrofi(n)) * &
+               areas(n) * glob_area_inv
+       end do
+
+       ! Determine sum of enthalpy correction for each hcorr index locally
+       local_htot_corr(1) = 0._r8
+       do n = 1,size(hcorr)
+          local_htot_corr(1) = local_htot_corr(1) + hcorr(n)
+       end do
+       call ESMF_VMAllreduce(is_local%wrap%vm, senddata=local_htot_corr, recvdata=global_htot_corr, count=1, &
+            reduceflag=ESMF_REDUCE_SUM, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       if (maintask) write(logunit, '(a,a,f21.13)') trim(subname),' global enthalpy correction: ',global_htot_corr(1)
+       deallocate(hcorr)
+    endif
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hsnow', rc)) deallocate(hsnow)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofl', rc)) deallocate(hrofl)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofi', rc)) deallocate(hrofi)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrain', rc)) deallocate(hrain)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hevap', rc)) deallocate(hevap)
+    if(.not. FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hcond', rc)) deallocate(hcond)
+
+    call t_stopf(subname)
+
+  end subroutine med_compute_enthalpy
+
+end module med_enthalpy_mod

mediator/med_internalstate_mod.F90

@@ -125,7 +125,8 @@ module med_internalstate_mod
     logical          :: ocn2glc_coupling = .false. ! obtained from attribute
     logical          :: lnd2glc_coupling = .false.
     logical          :: accum_lnd2glc = .false.
-
+    logical          :: docn_present               ! aoflux calc requires med_coupling_active true even for docn
+                                                   ! so we need an additional flag
     ! Mediator vm
     type(ESMF_VM) :: vm
 
@@ -282,8 +283,10 @@ subroutine med_internalstate_init(gcomp, rc)
     end if
     call NUOPC_CompAttributeGet(gcomp, name='OCN_model', value=ocn_name, isPresent=isPresent, isSet=isSet, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    is_local%wrap%docn_present = .false.
     if (isPresent .and. isSet) then
        if (trim(ocn_name) /= 'socn') is_local%wrap%comp_present(compocn) = .true.
+       if (trim(ocn_name) == 'docn') is_local%wrap%docn_present = .true.
     end if
     call NUOPC_CompAttributeGet(gcomp, name='ICE_model', value=ice_name, isPresent=isPresent, isSet=isSet, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return

mediator/med_phases_prep_atm_mod.F90

@@ -7,7 +7,7 @@ module med_phases_prep_atm_mod
   use med_kind_mod          , only : CX=>SHR_KIND_CX, CS=>SHR_KIND_CS, CL=>SHR_KIND_CL, R8=>SHR_KIND_R8
   use ESMF                  , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS
   use ESMF                  , only : ESMF_Field, ESMF_FieldGet, ESMF_FieldBundleGet
-  use ESMF                  , only : ESMF_GridComp, ESMF_GridCompGet
+  use ESMF                  , only : ESMF_GridComp, ESMF_GridCompGet, ESMF_FieldBundleIsCreated
   use med_constants_mod     , only : dbug_flag   => med_constants_dbug_flag
   use med_utils_mod         , only : memcheck    => med_memcheck
   use med_utils_mod         , only : chkerr      => med_utils_ChkErr
@@ -23,14 +23,12 @@ module med_phases_prep_atm_mod
   use perf_mod              , only : t_startf, t_stopf
   use med_phases_aofluxes_mod, only : med_aofluxes_map_xgrid2agrid_output
   use med_phases_aofluxes_mod, only : med_aofluxes_map_ogrid2agrid_output
+  use med_enthalpy_mod,       only : med_enthalpy_get_global_htot_corr, med_compute_enthalpy, mediator_compute_enthalpy
 
   implicit none
   private
 
   public :: med_phases_prep_atm
-  public :: med_phases_prep_atm_enthalpy_correction
-
-  real(r8), public :: global_htot_corr(1) = 0._r8  ! enthalpy correction from med_phases_prep_ocn
 
   character(*), parameter :: u_FILE_u  = &
        __FILE__
@@ -236,11 +234,17 @@ subroutine med_phases_prep_atm(gcomp, rc)
     end if
 
     ! Add enthalpy correction to sensible heat if appropriate
-    if (FB_FldChk(is_local%wrap%FBExp(compatm), 'Faxx_sen', rc=rc)) then
+    if (Mediator_compute_enthalpy) then
        call FB_getfldptr(is_local%wrap%FBExp(compatm), 'Faxx_sen', dataptr1, rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+       ! IF data ocn case compute first, otherwise computed in prep_ocn_mod
+       if(is_local%wrap%docn_present) then
+          call med_compute_enthalpy(is_local, rc)
+          if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       endif
        do n = 1,size(dataptr1)
-          dataptr1(n) = dataptr1(n) + global_htot_corr(1)
+          dataptr1(n) = dataptr1(n) + med_enthalpy_get_global_htot_corr()
        end do
     end if
 
@@ -255,48 +259,4 @@ subroutine med_phases_prep_atm(gcomp, rc)
 
   end subroutine med_phases_prep_atm
 
-  !-----------------------------------------------------------------------------
-  subroutine med_phases_prep_atm_enthalpy_correction (gcomp, hcorr, rc)
-
-    ! Enthalpy correction term calculation called by med_phases_prep_ocn_accum in
-    ! med_phases_prep_ocn_mod
-    ! Note that this is only called if the following fields are in FBExp(compocn)
-    ! 'Faxa_rain','Foxx_hrain','Faxa_snow' ,'Foxx_hsnow',
-    ! 'Foxx_evap','Foxx_hevap','Foxx_hcond','Foxx_rofl',
-    ! 'Foxx_hrofl','Foxx_rofi','Foxx_hrofi'
-
-    use ESMF            , only : ESMF_VMAllreduce, ESMF_GridCompGet, ESMF_REDUCE_SUM
-    use ESMF            , only : ESMF_VM
-
-    ! input/output variables
-    type(ESMF_GridComp) , intent(in)  :: gcomp
-    real(r8)            , intent(in)  :: hcorr(:)
-    integer             , intent(out) :: rc
-
-    ! local variables
-    type(InternalState) :: is_local
-    integer             :: n
-    real(r8)            :: local_htot_corr(1)
-    type(ESMF_VM)       :: vm
-    !---------------------------------------
-
-    rc = ESMF_SUCCESS
-
-    nullify(is_local%wrap)
-    call ESMF_GridCompGetInternalState(gcomp, is_local, rc)
-    if (chkErr(rc,__LINE__,u_FILE_u)) return
-
-    ! Determine sum of enthalpy correction for each hcorr index locally
-    local_htot_corr(1) = 0._r8
-    do n = 1,size(hcorr)
-       local_htot_corr(1) = local_htot_corr(1) + hcorr(n)
-    end do
-    call ESMF_GridCompGet(gcomp, vm=vm, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    call ESMF_VMAllreduce(vm, senddata=local_htot_corr, recvdata=global_htot_corr, count=1, &
-         reduceflag=ESMF_REDUCE_SUM, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
-
-  end subroutine med_phases_prep_atm_enthalpy_correction
-
 end module med_phases_prep_atm_mod