Skip to content

Commit

Permalink
add med_enthaply_mod
Browse files Browse the repository at this point in the history
  • Loading branch information
jedwards4b committed Aug 29, 2023
1 parent abaef5f commit 1943c2e
Show file tree
Hide file tree
Showing 6 changed files with 259 additions and 69 deletions.
14 changes: 14 additions & 0 deletions cime_config/namelist_definition_drv.xml
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
3 changes: 2 additions & 1 deletion mediator/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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}")
Expand Down
21 changes: 20 additions & 1 deletion mediator/med.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down
171 changes: 171 additions & 0 deletions mediator/med_enthalpy_mod.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
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_internalstate_mod, only : compocn, compatm, InternalState
use perf_mod, only : t_startf, t_stopf


implicit none
public :: med_compute_enthalpy

real(r8) :: global_htot_corr(1)
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(:)
real(r8), pointer :: snowl(:), snowc(:)
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
nmax = size(tocn)

call FB_GetFldPtr(is_local%wrap%FBImp(compocn,compocn), 'So_t', tocn, 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
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
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

do n = 1,nmax
! Need max to ensure that will not have an enthalpy contribution if the water is below 0C
hrain(n) = max((tocn(n) - tkfrz), 0._r8) * rain(n) * cpfw
hsnow(n) = min((tocn(n) - tkfrz), 0._r8) * snow(n) * cpice
hevap(n) = (tocn(n) - tkfrz) * min(evap(n), 0._r8) * cpwv
hcond(n) = (tocn(n) - tkfrz) * max(evap(n), 0._r8) * cpwv
hrofl(n) = max((tocn(n) - tkfrz), 0._r8) * rofl(n) * cpsw
hrofi(n) = min((tocn(n) - tkfrz), 0._r8) * rofi(n) * cpice
! GMM - note change in hcond
end do

! 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

deallocate(hcorr)
endif
call t_stopf(subname)

end subroutine med_compute_enthalpy

end module med_enthalpy_mod
56 changes: 8 additions & 48 deletions mediator/med_phases_prep_atm_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

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__
Expand Down Expand Up @@ -239,8 +237,14 @@ subroutine med_phases_prep_atm(gcomp, rc)
if (FB_FldChk(is_local%wrap%FBExp(compatm), 'Faxx_sen', rc=rc)) 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

! call med_compute_enthalpy(is_local, rc)
! if (ChkErr(rc,__LINE__,u_FILE_u)) return

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

Expand All @@ -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
Loading

0 comments on commit 1943c2e

Please sign in to comment.