Skip to content

Helfand Monin-Obukhov derivatives for GEOSldas offline runs #918

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
biljanaorescanin wants to merge 17 commits into from
Closed

Helfand Monin-Obukhov derivatives for GEOSldas offline runs #918

Changes from 7 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
aefa0da
helfand derivatives
biljanaorescanin Jan 18, 2024
c38db94
cleanup
biljanaorescanin Jan 23, 2024
c891499
move calls to run1
biljanaorescanin Jan 29, 2024
a250940
remove prints
biljanaorescanin Jan 29, 2024
9084ad6
Merge branch 'develop' into feature/boresca_helfand
biljanaorescanin Mar 22, 2024
8f8c9ae
code clenup
biljanaorescanin Mar 28, 2024
8e45c15
more comments clenup
biljanaorescanin Mar 28, 2024
b5bd6fe
change small_TC and small_QC
weiyuan-jiang Apr 4, 2024
b066a23
change pert value
biljanaorescanin Apr 15, 2024
534f271
add timer to louis and helfand
weiyuan-jiang Apr 19, 2024
8999b44
Merge branch 'develop' into feature/boresca_helfand
biljanaorescanin Apr 25, 2024
f23520f
Merge branch 'develop' into feature/boresca_helfand
biljanaorescanin May 6, 2024
05883e2
minor cleanup of code for extra derivatives in helfsurface:
gmao-rreichle May 16, 2024
aa8254c
Merge branch 'develop' into feature/boresca_helfand
biljanaorescanin May 20, 2024
4782ab7
Merge branch 'develop' into feature/boresca_helfand
biljanaorescanin Jun 2, 2024
60f8afb
Merge branch 'develop' into feature/boresca_helfand
gmao-rreichle Jun 17, 2024
d24b7cc
Merge branch 'develop' into feature/boresca_helfand
gmao-rreichle Sep 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 188 additions & 10 deletions ...GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1097,6 +1097,31 @@ subroutine SetServices ( GC, RC )
RC=STATUS )
VERIFY_(STATUS)

! for helfsurface
call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'pert_canopy_temperature' ,&
UNITS = 'K' ,&
SHORT_NAME = 'TC_pert' ,&
FRIENDLYTO = trim(COMP_NAME) ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'pert_canopy_specific_humidity' ,&
UNITS = 'kg kg-1' ,&
SHORT_NAME = 'QC_pert' ,&
FRIENDLYTO = trim(COMP_NAME) ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'canopy_specific_humidity' ,&
UNITS = 'kg kg-1' ,&
Expand Down Expand Up @@ -1352,6 +1377,51 @@ subroutine SetServices ( GC, RC )
RC=STATUS )
VERIFY_(STATUS)

! for helfsurface
call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'surface_heat_exchange_coeff_pert_tc_helfsurface',&
UNITS = 'kg m-2 s-1' ,&
SHORT_NAME = 'VKH_pert_tc' ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'surface_heat_exchange_coeff_pert_qc_helfsurface',&
UNITS = 'kg m-2 s-1' ,&
SHORT_NAME = 'VKH_pert_qc' ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'tc difference, optional in helfsurface' ,&
UNITS = 'kg m-2 s-1' ,&
SHORT_NAME = 'DTC' ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'qc difference, optional in helfsurface' ,&
UNITS = 'kg m-2 s-1' ,&
SHORT_NAME = 'DQC' ,&
DIMS = MAPL_DimsTileTile ,&
NUM_SUBTILES = NUM_SUBTILES ,&
VLOCATION = MAPL_VLocationNone ,&
RESTART = MAPL_RestartSkip ,&
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddInternalSpec(GC ,&
LONG_NAME = 'surface_momentum_exchange_coefficient',&
UNITS = 'kg m-2 s-1' ,&
Expand Down Expand Up @@ -3005,6 +3075,16 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
real, dimension(:,:), pointer :: DCH
real, dimension(:,:), pointer :: DCQ

! for helfsurface
real, dimension(:,:), pointer :: TC_pert
real, dimension(:,:), pointer :: QC_pert

real, allocatable :: VKH_pert_tc(:)
real, allocatable :: VKH_pert_qc(:)

real, dimension(:,:), pointer :: DQC
real, dimension(:,:), pointer :: DTC

! -----------------------------------------------------
! EXPORT Pointers
! -----------------------------------------------------
Expand Down Expand Up @@ -3086,10 +3166,17 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
real :: SCALE4ZVG
real :: SCALE4Z0_u
real :: MIN_VEG_HEIGHT

! for helfsurface
integer :: incl_HelfandMO_extra_derivs
real(ESMF_KIND_R8), parameter :: small_TC = 0.0001
real(ESMF_KIND_R8), parameter :: small_QC = 0.000001

type(CATCH_WRAP) :: wrap
type (T_CATCH_STATE), pointer :: CATCH_INTERNAL_STATE

incl_HelfandMO_extra_derivs = 0

!=============================================================================
! Begin...
! ------------------------------------------------------------------------------
Expand Down Expand Up @@ -3190,6 +3277,12 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,DCQ , 'DCQ' , RC=STATUS)
VERIFY_(STATUS)

! for helfsurface
call MAPL_GetPointer(INTERNAL,DQC , 'DQC' , RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,DTC , 'DTC' , RC=STATUS)
VERIFY_(STATUS)

! Pointers to outputs
!--------------------
Expand Down Expand Up @@ -3324,6 +3417,16 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
allocate(IWATER(NT),STAT=STATUS)
VERIFY_(STATUS)

! for helfsurface
allocate(TC_pert(NT,NUM_SUBTILES),STAT=STATUS)
VERIFY_(STATUS)
allocate(QC_pert(NT,NUM_SUBTILES),STAT=STATUS)
VERIFY_(STATUS)
allocate(VKH_pert_tc(NT),STAT=STATUS)
VERIFY_(STATUS)
allocate(VKH_pert_qc(NT),STAT=STATUS)
VERIFY_(STATUS)

! Vegetation types used to index into tables
!--------------------------------------------

Expand Down Expand Up @@ -3438,6 +3541,29 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
CALL helfsurface( UWINDLMTILE,VWINDLMTILE,TA,TC(:,N),QA,QC(:,N),PSL,PSMB,Z0T(:,N),lai, &
Copy link
Contributor

@weiyuan-jiang weiyuan-jiang Apr 1, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this intentional to call this function three times? or it should be
if (incl_HelfandMo_extra_derives ==0) then
call once
else if(incl_HelfandMO_extra_derivs ==1) then
call your addition ( twice)
endif

IWATER,DZE,niter,nt,RHOH,VKH,VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS, &
t2m,q2m,u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m,CHOOSEZ0)

! for helfsurface
if(CATCH_INTERNAL_STATE%CHOOSEMOSFC==1 .and. incl_HelfandMO_extra_derivs ==1) then
Copy link
Contributor

@weiyuan-jiang weiyuan-jiang Apr 1, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CATCH_INTERNAL_STATE%CHOOSEMOSFC==1 is redundant here

TC_pert(:,N)=TC(:,N)+small_TC
CALL helfsurface( UWINDLMTILE,VWINDLMTILE,TA,TC_pert(:,N),QA,QC(:,N),PSL,PSMB,Z0T(:,N),lai, &
IWATER,DZE,niter,nt,RHOH,VKH,VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS, &
t2m,q2m,u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m,CHOOSEZ0)
VKH_pert_tc =VKH

QC_pert(:,N)=QC(:,N)+small_QC
CALL helfsurface( UWINDLMTILE,VWINDLMTILE,TA,TC(:,N),QA,QC_pert(:,N),PSL,PSMB,Z0T(:,N),lai, &
IWATER,DZE,niter,nt,RHOH,VKH,VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS, &
t2m,q2m,u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m,CHOOSEZ0)
VKH_pert_qc=VKH

CALL helfsurface( UWINDLMTILE,VWINDLMTILE,TA,TC(:,N),QA,QC(:,N),PSL,PSMB,Z0T(:,N),lai, &
IWATER,DZE,niter,nt,RHOH,VKH,VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS, &
t2m,q2m,u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m,CHOOSEZ0)

DTC(:,N) = (VKH_pert_tc - VKH ) / small_TC
DQC(:,N) = (VKH_pert_qc - VKH ) / small_QC

endif

CM(:,N) = VKM
CH(:,N) = VKH
Expand Down Expand Up @@ -3533,6 +3659,11 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
deallocate(IWATER)
deallocate(PSMB)
deallocate(PSL)
! for helfsurface
deallocate(VKH_pert_tc)
deallocate(VKH_pert_qc)
deallocate(TC_pert )
deallocate(QC_pert )

! All done
! ------------------------------------------------------------------------------
Expand Down Expand Up @@ -3572,17 +3703,19 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC )
! Local derived type aliases
! ------------------------------------------------------------------------------

type(MAPL_MetaComp),pointer :: MAPL
type(ESMF_Alarm) :: ALARM
type(MAPL_MetaComp),pointer :: MAPL
type(ESMF_Alarm) :: ALARM

integer :: IM,JM
integer :: incl_Louis_extra_derivs
real :: SCALE4ZVG
real :: SCALE4Z0_u
real :: MIN_VEG_HEIGHT
type(ESMF_VM) :: VM
type (T_CATCH_STATE), pointer :: CATCH_INTERNAL_STATE
type (CATCH_WRAP) :: wrap
integer :: IM,JM
integer :: incl_Louis_extra_derivs
integer :: incl_HelfandMO_extra_derivs

real :: SCALE4ZVG
real :: SCALE4Z0_u
real :: MIN_VEG_HEIGHT
type(ESMF_VM) :: VM
type (T_CATCH_STATE), pointer :: CATCH_INTERNAL_STATE
type (CATCH_WRAP) :: wrap

! ------------------------------------------------------------------------------
! Begin: Get the target components name and
Expand Down Expand Up @@ -3611,6 +3744,10 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC )
VERIFY_(STATUS)

call MAPL_GetResource ( MAPL, incl_Louis_extra_derivs, Label="INCL_LOUIS_EXTRA_DERIVS:", DEFAULT=1, RC=STATUS)
VERIFY_(STATUS)

! for helfsurface
call MAPL_GetResource ( MAPL, incl_HelfandMO_extra_derivs, Label="INCL_HELFANDMO_EXTRA_DERIVS:", DEFAULT=1, RC=STATUS)
VERIFY_(STATUS)

call ESMF_VMGetCurrent(VM, rc=STATUS)
Expand Down Expand Up @@ -3809,6 +3946,9 @@ subroutine Driver ( RC )
real, dimension(:,:), pointer :: RBC002
real, dimension(:,:), pointer :: ROC001
real, dimension(:,:), pointer :: ROC002
! for helfsurface
real, dimension(:,:), pointer :: DQC
real, dimension(:,:), pointer :: DTC

! -----------------------------------------------------
! EXPORT Pointers
Expand Down Expand Up @@ -4340,6 +4480,10 @@ subroutine Driver ( RC )
call MAPL_GetPointer(INTERNAL,FR ,'FR' ,RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,DCQ ,'DCQ' ,RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,DCH ,'DCH' ,RC=STATUS); VERIFY_(STATUS)
! for helfsurface
call MAPL_GetPointer(INTERNAL,DTC ,'DTC' ,RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,DQC ,'DQC' ,RC=STATUS); VERIFY_(STATUS)

if (CATCH_INTERNAL_STATE%N_CONST_LAND4SNWALB /= 0) then
call MAPL_GetPointer(INTERNAL,RDU001 ,'RDU001' , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL,RDU002 ,'RDU002' , RC=STATUS); VERIFY_(STATUS)
Expand Down Expand Up @@ -4979,6 +5123,29 @@ subroutine Driver ( RC )
DHSDQA(:,N)=max(0.0,-MAPL_CP*DCH(:,N)*MAPL_VIREPS*TC(:,N)*(TC(:,N)-TA))
enddo
endif

! for helfsurface (adding derivatives for helfand MO scheme helfsurface)
if(CATCH_INTERNAL_STATE%CHOOSEMOSFC==1 .and. incl_HelfandMO_extra_derivs ==1) then

do N=1,NUM_SUBTILES

! getting DSHSBT - derivative of sensible heat flux w.r.t. Tc (ground temperature)
! note CH and CQ are VKH (values are assigned in RUN1 after the helfsurf call is returned)
DSHSBT(:,N)= MAPL_CP*( CH(:,N) +max(0.0,(TC(:,N)-TA) * DTC(:,N)))

! getting DHSDQA - cross derivative of sensible heat flux w.r.t. Qc (ground humidity)
DHSDQA(:,N)=max(0.0,(TC(:,N)-TA) * MAPL_CP * DQC(:,N))

! getting DEVSBT - derivative of latent heat flux w.r.t. Qc (ground humidity)
DEVSBT(:,N)=CQ(:,N)+max(0.0,(QC(:,N)-QA) * DQC(:,N))

! getting DEDTC - cross derivative of latent heat flux w.r.t. Tc (ground temperature)
DEDTC(:,N) =max(0.0,(QC(:,N)-QA) * DTC(:,N))

enddo ! N-loop (NUM_SUBTILES)

endif ! if Helfand

! BLWX = EMIS*MAPL_STFBOL*TA*TA*TA
! ALWX = -3.0*BLWX*TA
! BLWX = 4.0*BLWX
Expand Down Expand Up @@ -5959,6 +6126,11 @@ subroutine RUN0(gc, import, export, clock, rc)
real, pointer :: fr(:,:)=>null()
real, pointer :: DCQ(:,:)=>null()
real, pointer :: DCH(:,:)=>null()

! for helfsurface
real, pointer :: DTC(:,:)=>null()
real, pointer :: DQC(:,:)=>null()

!! -prognostic-variables-
real, pointer :: tc(:,:)=>null()
real, pointer :: qc(:,:)=>null()
Expand Down Expand Up @@ -6101,6 +6273,12 @@ subroutine RUN0(gc, import, export, clock, rc)
call MAPL_GetPointer(INTERNAL, catdef, 'CATDEF', rc=status)
VERIFY_(status)

! for helfsurface
call MAPL_GetPointer(INTERNAL, DTC, 'DTC', rc=status)
VERIFY_(status)
call MAPL_GetPointer(INTERNAL, DQC, 'DQC', rc=status)
VERIFY_(status)

! Number of tiles and a dummy real array
ntiles = size(HTSNNN1)
allocate(dummy(ntiles), stat=status)
Expand Down