diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 index 4de080553..7778c06c1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 @@ -3916,7 +3916,7 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) Iam=trim(COMP_NAME)//"::RUN1" - ! Get component's offline mode from its pvt internal state + ! Get component's offline mode from its private internal state call ESMF_UserCompGetInternalState(gc, 'CatchcnInternal', wrap, status) VERIFY_(status) catchcn_internal => wrap%ptr @@ -4481,7 +4481,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_Alarm) :: ALARM integer :: IM,JM - integer :: incl_Louis_extra_derivs real :: SCALE4Z0 @@ -4506,8 +4505,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_Get(MAPL, RUNALARM=ALARM, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource ( MAPL, incl_Louis_extra_derivs, Label="INCL_LOUIS_EXTRA_DERIVS:", DEFAULT=1, RC=STATUS) - VERIFY_(STATUS) call MAPL_GetResource ( MAPL, SCALE4Z0, Label="SCALE4Z0:", DEFAULT=0.5, RC=STATUS) VERIFY_(STATUS) @@ -6155,7 +6152,7 @@ subroutine Driver ( RC ) ALWN(:,N) = -3.0*BLWN(:,N)*TC(:,N) BLWN(:,N) = 4.0*BLWN(:,N) end do - if(catchcn_internal%CHOOSEMOSFC==0 .and. incl_Louis_extra_derivs ==1) then + if(catchcn_internal%CHOOSEMOSFC==0 .and. catchcn_internal%MOSFC_EXTRA_DERIVS_OFFL_LAND==1) then do N=1,NUM_SUBTILES DEVSBT(:,N)=CQ(:,N)+max(0.0,-DCQ(:,N)*MAPL_VIREPS*TC(:,N)*(QC(:,N)-QA)) DEDTC(:,N) =max(0.0,-DCQ(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(QC(:,N)-QA)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 index fa2a73b47..aa34a47ee 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 @@ -4445,7 +4445,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) type(MAPL_MetaComp),pointer :: MAPL type(ESMF_Alarm) :: ALARM integer :: IM,JM - integer :: incl_Louis_extra_derivs real :: SCALE4Z0 @@ -4470,8 +4469,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_Get(MAPL, RUNALARM=ALARM, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource ( MAPL, incl_Louis_extra_derivs, Label="INCL_LOUIS_EXTRA_DERIVS:", DEFAULT=1, RC=STATUS) - VERIFY_(STATUS) call MAPL_GetResource ( MAPL, SCALE4Z0, Label="SCALE4Z0:", DEFAULT=0.5, RC=STATUS) VERIFY_(STATUS) @@ -6185,7 +6182,7 @@ subroutine Driver ( RC ) ALWN(:,N) = -3.0*BLWN(:,N)*TC(:,N) BLWN(:,N) = 4.0*BLWN(:,N) end do - if(catchcn_internal%CHOOSEMOSFC==0 .and. incl_Louis_extra_derivs ==1) then + if(catchcn_internal%CHOOSEMOSFC==0 .and. catchcn_internal%MOSFC_EXTRA_DERIVS_OFFL_LAND==1) then do N=1,NUM_SUBTILES DEVSBT(:,N)=CQ(:,N)+max(0.0,-DCQ(:,N)*MAPL_VIREPS*TC(:,N)*(QC(:,N)-QA)) DEDTC(:,N) =max(0.0,-DCQ(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(QC(:,N)-QA)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index b48dbc164..6f0fa2747 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -1,6 +1,7 @@ #include "MAPL_Generic.h" #define DEALLOC_(A) if(associated(A))then;A=0;if(MAPL_ShmInitialized)then; call MAPL_DeAllocNodeArray(A,rc=STATUS);else; deallocate(A,stat=STATUS);endif;_VERIFY(STATUS);NULLIFY(A);endif + !============================================================================= module GEOS_CatchGridCompMod @@ -61,62 +62,62 @@ module GEOS_CatchGridCompMod !#-- use catch_wrap_stateMod -implicit none -private - -! !PUBLIC MEMBER FUNCTIONS: - -public SetServices - -! -!EOP - -integer,parameter :: FSAT=1 ! Saturated subtile -integer,parameter :: FTRN=2 ! Transition subtile -integer,parameter :: FWLT=3 ! Wilting subtile -integer,parameter :: FSNW=4 ! Snowcover subtile - -integer,parameter :: NUM_SUBTILES=4 - -! Vegetation type as follows: -! 1: BROADLEAF EVERGREEN TREES -! 2: BROADLEAF DECIDUOUS TREES -! 3: NEEDLELEAF TREES -! 4: GROUND COVER -! 5: BROADLEAF SHRUBS -! 6: DWARF TREES (TUNDRA) -!=================================================== -!ALT: we currently use only 6 types (see above) -! in the legacy code we used to have 8 -! (or 10 with the sea and land ice) with -! these additional entries -! 7: BARE SOIL -! 8: DESERT - -integer,parameter :: NTYPS = MAPL_NUMVEGTYPES - -! Veg-dep. vector SAI factor for scaling of rough length (now exp(-.5) ) -real, parameter :: SAI4ZVG(NTYPS) = (/ 0.60653, 0.60653, 0.60653, 1.0, 1.0, 1.0 /) -real, parameter :: HPBL = 1000. -!real, parameter :: MIN_VEG_HEIGHT = 0.01 -real, parameter :: Z0_BY_ZVEG = 0.13 -real, parameter :: D0_BY_ZVEG = 0.66 - -! Emissivity values from Wilber et al (1999, NATA-TP-1999-209362) -! Fu-Liou bands have been combined to Chou bands (though these are broadband only) -! IGBP veg types have been mapped to Sib-Mosaic types -! Details in ~suarez/Emiss on cerebus - -real, parameter :: EMSVEG(NTYPS) = (/ 0.99560, 0.99000, 0.99560, 0.99320, & - 0.99280, 0.99180 /) -real, parameter :: EMSBARESOIL = 0.94120 -real, parameter :: EMSSNO = 0.99999 - -! moved SURFLAY from catchment.F90 to enable run-time changes for off-line system -! - reichle, 29 Oct 2010 - -! real, parameter :: SURFLAY = 20. ! moved to GetResource in RUN2 LLT:12Jul2013 -! SURFLAY moved to internal state + implicit none + private + + ! !PUBLIC MEMBER FUNCTIONS: + + public SetServices + + ! + !EOP + + integer,parameter :: FSAT=1 ! Saturated subtile + integer,parameter :: FTRN=2 ! Transition subtile + integer,parameter :: FWLT=3 ! Wilting subtile + integer,parameter :: FSNW=4 ! Snowcover subtile + + integer,parameter :: NUM_SUBTILES=4 + + ! Vegetation type as follows: + ! 1: BROADLEAF EVERGREEN TREES + ! 2: BROADLEAF DECIDUOUS TREES + ! 3: NEEDLELEAF TREES + ! 4: GROUND COVER + ! 5: BROADLEAF SHRUBS + ! 6: DWARF TREES (TUNDRA) + !=================================================== + !ALT: we currently use only 6 types (see above) + ! in the legacy code we used to have 8 + ! (or 10 with the sea and land ice) with + ! these additional entries + ! 7: BARE SOIL + ! 8: DESERT + + integer,parameter :: NTYPS = MAPL_NUMVEGTYPES + + ! Veg-dep. vector SAI factor for scaling of rough length (now exp(-.5) ) + real, parameter :: SAI4ZVG(NTYPS) = (/ 0.60653, 0.60653, 0.60653, 1.0, 1.0, 1.0 /) + real, parameter :: HPBL = 1000. + !real, parameter :: MIN_VEG_HEIGHT = 0.01 + real, parameter :: Z0_BY_ZVEG = 0.13 + real, parameter :: D0_BY_ZVEG = 0.66 + + ! Emissivity values from Wilber et al (1999, NATA-TP-1999-209362) + ! Fu-Liou bands have been combined to Chou bands (though these are broadband only) + ! IGBP veg types have been mapped to Sib-Mosaic types + ! Details in ~suarez/Emiss on cerebus + + real, parameter :: EMSVEG(NTYPS) = (/ 0.99560, 0.99000, 0.99560, 0.99320, & + 0.99280, 0.99180 /) + real, parameter :: EMSBARESOIL = 0.94120 + real, parameter :: EMSSNO = 0.99999 + + ! moved SURFLAY from catchment.F90 to enable run-time changes for off-line system + ! - reichle, 29 Oct 2010 + + ! real, parameter :: SURFLAY = 20. ! moved to GetResource in RUN2 LLT:12Jul2013 + ! SURFLAY moved to internal state contains @@ -1366,7 +1367,7 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'surface_moisture_exchange_coffiecient',& + LONG_NAME = 'surface_moisture_exchange_coefficient',& UNITS = 'kg m-2 s-1' ,& SHORT_NAME = 'CQ' ,& DIMS = MAPL_DimsTileTile ,& @@ -1397,29 +1398,62 @@ subroutine SetServices ( GC, RC ) RESTART = RESTART ,& RC=STATUS ) VERIFY_(STATUS) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND == 1) then + + ! for *analytical* extra derivatives in louissurface + + call MAPL_AddInternalSpec(GC, & + SHORT_NAME = 'delCH_delTVA', & + LONG_NAME = 'partial_derivative_of_CH_wrt_virtual_Tair', & + UNITS = '1', & + DIMS = MAPL_DimsTileTile, & + NUM_SUBTILES = NUM_SUBTILES ,& + VLOCATION = MAPL_VLocationNone, & + RESTART = MAPL_RestartSkip ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC, & + SHORT_NAME = 'delCQ_delTVA', & + LONG_NAME = 'partial_derivative_of_CQ_wrt_virtual_Tair', & + UNITS = '1', & + DIMS = MAPL_DimsTileTile, & + NUM_SUBTILES = NUM_SUBTILES ,& + VLOCATION = MAPL_VLocationNone, & + RESTART = MAPL_RestartSkip ,& + RC=STATUS ) + VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC, & - SHORT_NAME = 'DCH', & - LONG_NAME = 'ch difference, optional in louissurface', & - UNITS = '1', & - DIMS = MAPL_DimsTileTile, & - NUM_SUBTILES = NUM_SUBTILES ,& - VLOCATION = MAPL_VLocationNone, & - RESTART = MAPL_RestartSkip ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC, & - SHORT_NAME = 'DCQ', & - LONG_NAME = 'cq difference, optional in louissurface', & - UNITS = '1', & - DIMS = MAPL_DimsTileTile, & - NUM_SUBTILES = NUM_SUBTILES ,& - VLOCATION = MAPL_VLocationNone, & - RESTART = MAPL_RestartSkip ,& - RC=STATUS ) - VERIFY_(STATUS) + elseif (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND >= 2) then + + ! for *numerical* extra derivatives in helfsurface and louissurface + + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'partial_derivative_of_CH_wrt_canopy_temperature', & + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'delCH_delTC' ,& + DIMS = MAPL_DimsTileTile ,& + NUM_SUBTILES = NUM_SUBTILES ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartSkip ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'partial_derivative_of_CQ_wrt_canopy_specific_humidity', & + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'delCQ_delQC' ,& + DIMS = MAPL_DimsTileTile ,& + NUM_SUBTILES = NUM_SUBTILES ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartSkip ,& + RC=STATUS ) + VERIFY_(STATUS) + end if + + !---------- GOSWIM snow impurity related variables ---------- if (CATCH_INTERNAL_STATE%N_CONST_LAND4SNWALB /= 0) then @@ -2920,7 +2954,6 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_TimerOn(MAPL,"INITIALIZE") - ! retrieve internal state call ESMF_UserCompGetInternalState ( GC, 'CatchInternal',wrap,status ) @@ -3076,7 +3109,7 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_Clock), intent(inout) :: CLOCK !The clock integer,optional, intent(out ) :: RC !Error code: -! !DESCRIPTION: Does the cds computation and roughness length +! !DESCRIPTION: Compute roughness length and exchange coefficients ("cds"), incl. derivatives !EOP ! ErrLog Variables @@ -3122,8 +3155,16 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) real, dimension(:,:), pointer :: CQ real, dimension(:,:), pointer :: FR real, dimension(:,:), pointer :: WW - real, dimension(:,:), pointer :: DCH - real, dimension(:,:), pointer :: DCQ + + ! for analytical extra derivatives (louissurface) + + real, dimension(:,:), pointer :: delCH_delTVA + real, dimension(:,:), pointer :: delCQ_delTVA + + ! for numerical extra derivatives (louissurface, helfsurface) + + real, dimension(:,:), pointer :: delCH_delTC + real, dimension(:,:), pointer :: delCQ_delQC ! ----------------------------------------------------- ! EXPORT Pointers @@ -3206,9 +3247,25 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) real :: SCALE4ZVG real :: SCALE4Z0_u real :: MIN_VEG_HEIGHT + + ! ------------------------------------- + ! + ! for numerical extra derivatives (louissurface, helfsurface) + + real, parameter :: MOSFC_pert_fac = 0.001 ! size of multiplicative pert for numerical derivatives + + ! Louis needs 2d arrays; Helfand would work with 1d arrays but use 2d arrays to avoid "if" statements + + real, dimension(:,:), allocatable :: DeltaTC, CHpert + real, dimension(:,:), allocatable :: DeltaQC, CQpert + + real, dimension(:,:), allocatable :: DummyZ0T, DummyCM + + ! ------------------------------------- type(CATCH_WRAP) :: wrap type (T_CATCH_STATE), pointer :: CATCH_INTERNAL_STATE + !============================================================================= ! Begin... @@ -3253,7 +3310,7 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) INTERNAL_ESMF_STATE=INTERNAL , & RC=STATUS ) VERIFY_(STATUS) - + call MAPL_GetResource ( MAPL, CHOOSEZ0, Label="CHOOSEZ0:", DEFAULT=3, RC=STATUS) VERIFY_(STATUS) call ESMF_VMGetCurrent(VM, rc=STATUS) @@ -3306,11 +3363,24 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL,WW , 'WW' , RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,DCH , 'DCH' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,DCQ , 'DCQ' , RC=STATUS) - VERIFY_(STATUS) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND == 1) then + + call MAPL_GetPointer(INTERNAL,delCH_delTVA , 'delCH_delTVA' , RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,delCQ_delTVA , 'delCQ_delTVA' , RC=STATUS) + VERIFY_(STATUS) + + elseif (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND >= 2) then + call MAPL_GetPointer(INTERNAL,delCQ_delQC , 'delCQ_delQC' , RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,delCH_delTC , 'delCH_delTC' , RC=STATUS) + VERIFY_(STATUS) + + end if + + ! Pointers to outputs !-------------------- @@ -3444,6 +3514,21 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) allocate(IWATER(NT),STAT=STATUS) VERIFY_(STATUS) + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND>=2) then + + ! allocate variables for numerical extra derivatives (louissurface, helfsurface) + + allocate(DeltaTC( NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + allocate(DeltaQC( NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + + allocate(CHpert( NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + allocate(CQpert( NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + + allocate(DummyZ0T(NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + allocate(DummyCM( NT,NUM_SUBTILES),STAT=STATUS); VERIFY_(STATUS) + + end if + ! Vegetation types used to index into tables !-------------------------------------------- @@ -3535,55 +3620,193 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) if(associated(Z0 )) Z0 = Z0T(:,N) if(associated(D0 )) D0 = D0T -! Compute the three surface exchange coefficients -!------------------------------------------------- +! Compute surface exchange coefficients +!--------------------------------------- + + call MAPL_TimerOn(MAPL,"-SURF") ! timer for computation of MOSFC exchange coeffs and derivs (Louis or Helfand) + + if (CATCH_INTERNAL_STATE%CHOOSEMOSFC.eq.0) then - call MAPL_TimerOn(MAPL,"-SURF") - if(CATCH_INTERNAL_STATE%CHOOSEMOSFC.eq.0) then - WW(:,N) = 0. - CM(:,N) = 0. + ! Louis surface turbulence + + WW(:,N) = 0. + CM(:,N) = 0. + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND==1) then + + ! analytical extra derivatives (default for Louis) + + call louissurface(3,N,UU,WW,PS,TA,TC,QA,QC,PCU,LAI,Z0T,DZE,CM,CN,RIB,ZT,ZQ,CH,CQ,UUU,UCN,RE,delCH_delTVA,delCQ_delTVA) + + else + + ! none .or. numerical extra derivatives + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND>=2) then + + ! Prep calculation of numerical extra derivatives. Start with calling louissurface with perturbed inputs and + ! save only the perturbed exchange coeffs. The final call with nominal inputs produces the unperturbed + ! exchange coeffs and other outputs (CN, RIB, ZT, ZQ, etc). + ! Must use properly initialized dummmies for Z0T and CM because these are intent(inout). + + ! perturb TC: send in (TC+DeltaTC), get back CHpert + + DeltaTC = MOSFC_pert_fac*TC + + DummyZ0T = Z0T + DummyCM = CM + + call louissurface( 3,N,UU,WW,PS,TA,TC+DeltaTC,QA,QC ,PCU,LAI,DummyZ0T,DZE,DummyCM,CN,RIB,ZT,ZQ,CHpert,CQ ,UUU,UCN,RE) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND==2) then + + ! perturb QC: send in (QC+DeltaQC), get back CQpert + + DeltaQC = MOSFC_pert_fac*QC + + DummyZ0T = Z0T + DummyCM = CM + + call louissurface(3,N,UU,WW,PS,TA,TC ,QA,QC+DeltaQC,PCU,LAI,DummyZ0T,DZE,DummyCM,CN,RIB,ZT,ZQ,CH ,CQpert,UUU,UCN,RE) + + end if + + end if + + ! Call with nominal inputs [after calls with perturbed inputs to obtain correct outputs (CN, RIB, ZT, ZQ, etc.)] + + call louissurface(3,N,UU,WW,PS,TA,TC,QA,QC,PCU,LAI,Z0T,DZE,CM,CN,RIB,ZT,ZQ,CH,CQ,UUU,UCN,RE) + + end if ! MOSFC_EXTRA_DERIVS_OFFL_LAND + + elseif (CATCH_INTERNAL_STATE%CHOOSEMOSFC.eq.1)then - call louissurface(3,N,UU,WW,PS,TA,TC,QA,QC,PCU,LAI,Z0T,DZE,CM,CN,RIB,ZT,ZQ,CH,CQ,UUU,UCN,RE,DCH,DCQ) + ! Helfand surface turbulence - elseif (CATCH_INTERNAL_STATE%CHOOSEMOSFC.eq.1)then - - niter = 6 ! number of internal iterations in the helfand MO surface layer routine - IWATER = 3 - - PSMB = PS * 0.01 ! convert to MB -! Approximate pressure at top of surface layer: hydrostatic, eqn of state using avg temp and press - PSL = PSMB * (1. - (DZE*MAPL_GRAV)/(MAPL_RGAS*(TA+TC(:,N)) ) ) / & - (1. + (DZE*MAPL_GRAV)/(MAPL_RGAS*(TA+TC(:,N)) ) ) - - 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) - - CM(:,N) = VKM - CH(:,N) = VKH - CQ(:,N) = VKH - - CN = (MAPL_KARMAN/ALOG(DZE/Z0T(:,N) + 1.0)) * (MAPL_KARMAN/ALOG(DZE/Z0T(:,N) + 1.0)) - ZT = Z0T(:,N) - ZQ = Z0T(:,N) - RE = 0. - UUU = UU - UCN = 0. - + niter = 6 ! number of internal iterations in the Helfand MO surface layer routine + IWATER = 3 + + PSMB = PS * 0.01 ! convert to MB + + ! Approximate pressure at top of surface layer: hydrostatic, eqn of state using avg temp and press + PSL = PSMB * (1. - (DZE*MAPL_GRAV)/(MAPL_RGAS*(TA+TC(:,N)) ) ) / & + (1. + (DZE*MAPL_GRAV)/(MAPL_RGAS*(TA+TC(:,N)) ) ) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND==2) then + + ! Prep calculation of numerical extra derivatives. Start with calling louissurface with perturbed inputs and + ! save only the perturbed exchange coeffs. The final call with nominal inputs produces the unperturbed + ! exchange coeffs and other outputs (CN, RIB, ZT, ZQ, etc). + ! Must use properly initialized dummmies for Z0T and CM because these are intent(inout). + + ! perturb TC: send in (TC+DeltaTC), get back CHpert + + DeltaTC( :,N) = MOSFC_pert_fac*TC(:,N) + + DummyZ0T(:,N) = Z0T(:,N) + + CALL helfsurface(UWINDLMTILE,VWINDLMTILE,TA,TC(:,N)+DeltaTC(:,N),QA,QC(:,N) ,PSL,PSMB,DummyZ0T(:,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) + + CHpert( :,N) = VKH + + ! perturb QC: send in (QC+DeltaQC), get back CQpert + + DeltaQC( :,N) = MOSFC_pert_fac*QC(:,N) + + DummyZ0T(:,N) = Z0T(:,N) + + CALL helfsurface(UWINDLMTILE,VWINDLMTILE,TA,TC(:,N) ,QA,QC(:,N)+DeltaQC(:,N),PSL,PSMB,DummyZ0T(:,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) + + CQpert( :,N) = VKH + + end if ! MOSFC_EXTRA_DERIVS_OFFL_LAND==2 + + ! Call with nominal inputs [after calls with perturbed inputs to obtain correct outputs (Z0T, [*]2m, [*]10m, etc.)] + + 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) + + CM(:,N) = VKM + CH(:,N) = VKH + CQ(:,N) = VKH + + CN = (MAPL_KARMAN/ALOG(DZE/Z0T(:,N) + 1.0)) * (MAPL_KARMAN/ALOG(DZE/Z0T(:,N) + 1.0)) + ZT = Z0T(:,N) + ZQ = Z0T(:,N) + RE = 0. + UUU = UU + UCN = 0. + ! Aggregate to tiles for MO only diagnostics !-------------------------------------------- - if(associated(MOU50M))MOU50M = MOU50M + U50M(:)*FR(:,N) - if(associated(MOV50M))MOV50M = MOV50M + V50M(:)*FR(:,N) - if(associated(MOT10M))MOT10M = MOT10M + T10M(:)*FR(:,N) - if(associated(MOQ10M))MOQ10M = MOQ10M + Q10M(:)*FR(:,N) - if(associated(MOU10M))MOU10M = MOU10M + U10M(:)*FR(:,N) - if(associated(MOV10M))MOV10M = MOV10M + V10M(:)*FR(:,N) - if(associated(MOT2M))MOT2M = MOT2M + T2M(:)*FR(:,N) - if(associated(MOQ2M))MOQ2M = MOQ2M + Q2M(:)*FR(:,N) - if(associated(MOU2M))MOU2M = MOU2M + U2M(:)*FR(:,N) - if(associated(MOV2M))MOV2M = MOV2M + V2M(:)*FR(:,N) + if(associated(MOU50M))MOU50M = MOU50M + U50M(:)*FR(:,N) + if(associated(MOV50M))MOV50M = MOV50M + V50M(:)*FR(:,N) + if(associated(MOT10M))MOT10M = MOT10M + T10M(:)*FR(:,N) + if(associated(MOQ10M))MOQ10M = MOQ10M + Q10M(:)*FR(:,N) + if(associated(MOU10M))MOU10M = MOU10M + U10M(:)*FR(:,N) + if(associated(MOV10M))MOV10M = MOV10M + V10M(:)*FR(:,N) + if(associated(MOT2M ))MOT2M = MOT2M + T2M( :)*FR(:,N) + if(associated(MOQ2M ))MOQ2M = MOQ2M + Q2M( :)*FR(:,N) + if(associated(MOU2M ))MOU2M = MOU2M + U2M( :)*FR(:,N) + if(associated(MOV2M ))MOV2M = MOV2M + V2M( :)*FR(:,N) + + endif ! CHOOSEMOSFC + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND==2) then + + ! finalize numerical derivatives + + delCH_delTC(:,N) = (CHpert(:,N) - CH(:,N)) / DeltaTC(:,N) + delCQ_delQC(:,N) = (CQpert(:,N) - CQ(:,N)) / DeltaQC(:,N) + + elseif (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND==3) then + + ! finalize numerical derivatives (valid for Louis only!) + + ! For Louis, the exchange coeffs depend only on the *virtual* temperature (not true for Helfand). + ! + ! This lets us compute the derivatives of the exchange coefficients w.r.t. both TC and QC from + ! just one additional call to louissurface() with perturbed TC. + ! + ! In the following, "del" indicates a derivative and "Delta" indicates a difference term. + ! + ! We have: + ! + ! (1) CH = CQ + ! + ! (2) TVC = TC*(1 + eps*QC) (virtual temperature; eps = MAPL_VIREPS) + ! + ! (2a) ==> delTVC_delQC = eps*TC + ! + ! (2b) ==> DeltaTVC = (TVCpert - TVC) = TCpert*(1 + eps*QC) - TC*(1 + eps*QC) = DeltaTC*(1 + eps*QC) + ! + ! (3) delCH_delTC = (CHpert - CH)/DeltaTC + ! + ! (4) delCH_delTVC = (CHpert - CH)/DeltaTVC = (CHpert - CH)/(DeltaTC*(1 + eps*QC)) + ! + ! Using (1)-(4), we have: + ! + ! delCQ_delQC = delCH_delQC using (1) + ! + ! = delCH_delTVC * delTVC_delQC using chain rule + ! + ! = (CHpert - CH)/(DeltaTC*(1 + eps*QC)) * delTVC_delQC using (4) + ! + ! = (CHpert - CH)/DeltaTC * 1/(1+eps*QC) * eps*TC using (2a) + ! + ! = delCH_delTC * 1/(1+eps*QC) * eps*TC using (3) + + delCH_delTC(:,N) = (CHpert(:,N) - CH(:,N)) / DeltaTC(:,N) + delCQ_delQC(:,N) = delCH_delTC(:,N) * MAPL_VIREPS*TC(:,N)/(1.+MAPL_VIREPS*QC(:,N)) + endif + call MAPL_TimerOff(MAPL,"-SURF") ! Aggregate to tile @@ -3653,9 +3876,17 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) deallocate(IWATER) deallocate(PSMB) deallocate(PSL) + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND>=2) then + deallocate(DeltaTC ) + deallocate(DeltaQC ) + deallocate(CHpert ) + deallocate(CQpert ) + deallocate(dummyZ0T) + deallocate(dummyCM ) + end if -! All done -! ------------------------------------------------------------------------------ + ! All done + ! ------------------------------------------------------------------------------ call MAPL_TimerOff ( MAPL, "RUN1" ) call MAPL_TimerOff ( MAPL, "TOTAL" ) @@ -3692,17 +3923,17 @@ 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 + + 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 @@ -3730,9 +3961,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_Get(MAPL, RUNALARM=ALARM, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource ( MAPL, incl_Louis_extra_derivs, Label="INCL_LOUIS_EXTRA_DERIVS:", DEFAULT=1, RC=STATUS) - VERIFY_(STATUS) - call ESMF_VMGetCurrent(VM, rc=STATUS) select case (CATCH_INTERNAL_STATE%Z0_FORMULATION) @@ -3918,8 +4146,10 @@ subroutine Driver ( RC ) real, dimension(:,:), pointer :: cm real, dimension(:,:), pointer :: cq real, dimension(:,:), pointer :: fr - real, dimension(:,:), pointer :: dcq - real, dimension(:,:), pointer :: dch + real, dimension(:,:), pointer :: delCQ_delTVA + real, dimension(:,:), pointer :: delCH_delTVA + real, dimension(:,:), pointer :: delCH_delTC + real, dimension(:,:), pointer :: delCQ_delQC real, dimension(:,:), pointer :: RDU001 real, dimension(:,:), pointer :: RDU002 real, dimension(:,:), pointer :: RDU003 @@ -4106,7 +4336,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:,:) :: evsbt real,pointer,dimension(:,:) :: devsbt real,pointer,dimension(:,:) :: DEDTC - real,pointer,dimension(:,:) :: DHSDQA + real,pointer,dimension(:,:) :: DHSDQC real,pointer,dimension(:,:) :: CFT real,pointer,dimension(:,:) :: RA real,pointer,dimension(:,:) :: CFQ @@ -4472,8 +4702,19 @@ subroutine Driver ( RC ) call MAPL_GetPointer(INTERNAL,CM ,'CM' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL,CQ ,'CQ' ,RC=STATUS); VERIFY_(STATUS) 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) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND == 1) then + + call MAPL_GetPointer(INTERNAL,delCQ_delTVA ,'delCQ_delTVA' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,delCH_delTVA ,'delCH_delTVA' ,RC=STATUS); VERIFY_(STATUS) + + elseif (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND >= 2) then + + call MAPL_GetPointer(INTERNAL,delCH_delTC ,'delCH_delTC' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,delCQ_delQC ,'delCQ_delQC' ,RC=STATUS); VERIFY_(STATUS) + + end if + 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) @@ -4683,7 +4924,7 @@ subroutine Driver ( RC ) allocate(EVSBT (NTILES,NUM_SUBTILES)) allocate(DEVSBT (NTILES,NUM_SUBTILES)) allocate(DEDTC (NTILES,NUM_SUBTILES)) - allocate(DHSDQA (NTILES,NUM_SUBTILES)) + allocate(DHSDQC (NTILES,NUM_SUBTILES)) allocate(CFT (NTILES,NUM_SUBTILES)) allocate(CFQ (NTILES,NUM_SUBTILES)) allocate(TCO (NTILES,NUM_SUBTILES)) @@ -5106,44 +5347,140 @@ subroutine Driver ( RC ) RDC = max(VGRDA(VEG)*min(1., LAI/VGRDB(VEG)),0.001) RHO = PS/(MAPL_RGAS*(TA*(1+MAPL_VIREPS*QA))) - DEDTC=0.0 - DHSDQA=0.0 - if(CATCH_INTERNAL_STATE%CATCH_OFFLINE /=0) then + !-------------------------------------------------------------------------------------------------------- + ! + ! MOSFC variable names and description: + ! + !-------------------------------------------------------------------------------------------------------- + ! GEOS_CatchGridComp.F90 | catchment.F90 | dimension | Description + !-------------------------------------------------------------------------------------------------------- + ! TA | TM | NT | surface (lowest model level) air temperature + ! QA | QM | NT | surface (lowest model level) air spec humidity + !-------------------------------------------------------------------------------------------------------- + ! TC | TC | NT-by-NSBT | canopy (air) temperature + ! QC | QA (!) | NT-by-NSBT | canopy (air) specific humidity + ! CH | - | NT-by-NSBT | exchange coeff for heat + ! CQ | - | NT-by-NSBT | exchange coeff for humidity + ! EVSBT | ETURB | NT-by-NSBT | evaporation + ! DEVSBT | DEDQA | NT-by-NSBT | deriv of evap w.r.t. canopy spec humidity + ! DEDTC | DEDTC | NT-by-NSBT | deriv of evap w.r.t. canopy temperature + ! SHSBT | HSTURB | NT-by-NSBT | sensible heat flux (SH) + ! DHSDQC (formerly DHSDQA) | DHSDQA | NT-by-NSBT | deriv of SH w.r.t. canopy spec humidity + ! DSHSBT | DHSDTC | NT-by-NSBT | deriv of SH w.r.t. canopy temperature + !-------------------------------------------------------------------------------------------------------- + ! *SBT = sub-tile (?) + ! NT = number of tiles + ! NSBT = number of subtiles (per tile) + ! + ! For land, CH = CQ in Helfand and Louis. + ! + ! + ! MOSFC equations: + ! + ! EVSBT = CQ * (QC - QA) + ! SHSBT = Cp * CH * (TC - TA) [ Cp = MAPL_CP ] + ! + ! Derivatives obtained via product rule. See equations below. + ! + ! For analytical derivatives, additionally use the following identities: + ! + ! virtual TC: TVC = TC*(1 + eps)*QC [ eps = MAPL_VIREPS ] + ! virtual TA: TVA = TA*(1 + eps)*QA + ! + ! delTVC_delQC = TC*eps + ! delTVC_delTC = (1 + eps) + ! + ! delCQ_delQC = delCQ_delTVC * delTVC_delQC + ! delCH_delTC = delCH_delTVC * delTVC_delTC + ! + ! CQ=CQ(Ri) where Ri is proportional to deltaTVA=TVA-TVC --> produces a minus sign + ! + ! delCQ_delTVC = -1*delCQ_delTVA + ! delCH_delTVC = -1*delCH_delTVA + ! + !-------------------------------------------------------------------------------------------------- + ! reichle, 9/9/2024 + !-------------------------------------------------------------------------------------------------- + + ! initialize derivatives that may not be filled later + + DEDTC =0.0 + DHSDQC=0.0 + + if (CATCH_INTERNAL_STATE%CATCH_OFFLINE /=0) then + + ! Catchment in offline (land-only) mode + do N=1,NUM_SUBTILES - CFT (:,N) = 1.0 - CFQ (:,N) = 1.0 - SHSBT (:,N) = MAPL_CP*CH(:,N)*(TC(:,N)-TA) - EVSBT (:,N) = CQ(:,N)*(QC(:,N)-QA) - DSHSBT(:,N) = MAPL_CP*CH(:,N) - DEVSBT(:,N) = CQ(:,N) - BLWN(:,N) = EMIS*MAPL_STFBOL*TC(:,N)*TC(:,N)*TC(:,N) - ALWN(:,N) = -3.0*BLWN(:,N)*TC(:,N) - BLWN(:,N) = 4.0*BLWN(:,N) + + CFT (:,N) = 1.0 + CFQ (:,N) = 1.0 + + SHSBT(:,N) = MAPL_CP*CH(:,N)*(TC(:,N)-TA) + EVSBT(:,N) = CQ(:,N)*(QC(:,N)-QA) + + BLWN( :,N) = EMIS*MAPL_STFBOL*TC(:,N)*TC(:,N)*TC(:,N) + ALWN( :,N) = -3.0*BLWN(:,N)*TC(:,N) + BLWN( :,N) = 4.0*BLWN(:,N) + end do - if(CATCH_INTERNAL_STATE%CHOOSEMOSFC==0 .and. incl_Louis_extra_derivs ==1) then + + select case (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND) + + case (0) ! ignore derivatives of exchange coeffs w.r.t. canopy temp and specific humidity + do N=1,NUM_SUBTILES - DEVSBT(:,N)=CQ(:,N)+max(0.0,-DCQ(:,N)*MAPL_VIREPS*TC(:,N)*(QC(:,N)-QA)) - DEDTC(:,N) =max(0.0,-DCQ(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(QC(:,N)-QA)) - DSHSBT(:,N)=MAPL_CP*(CH(:,N)+max(0.0,-DCH(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(TC(:,N)-TA))) - DHSDQA(:,N)=max(0.0,-MAPL_CP*DCH(:,N)*MAPL_VIREPS*TC(:,N)*(TC(:,N)-TA)) - enddo - endif - ! BLWX = EMIS*MAPL_STFBOL*TA*TA*TA - ! ALWX = -3.0*BLWX*TA - ! BLWX = 4.0*BLWX + DEVSBT(:,N) = CQ(:,N) + DSHSBT(:,N) = MAPL_CP* CH(:,N) + end do + + case (1) ! Louis only: analytical derivatives of exchange coeffs w.r.t. canopy temp and specific humidity + + _ASSERT( CATCH_INTERNAL_STATE%CHOOSEMOSFC==0, 'must use Louis scheme for MOSFC analytical derivatives' ) + + do N=1,NUM_SUBTILES + DEVSBT(:,N) = CQ(:,N) + max( 0.0, -delCQ_delTVA(:,N)* MAPL_VIREPS*TC(:,N) *(QC(:,N)-QA) ) + DEDTC( :,N) = max( 0.0, -delCQ_delTVA(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(QC(:,N)-QA) ) + DSHSBT(:,N) = MAPL_CP*( CH(:,N) + max( 0.0, -delCH_delTVA(:,N)*(1.+MAPL_VIREPS*QC(:,N))*(TC(:,N)-TA) ) ) + DHSDQC(:,N) = max( 0.0, -MAPL_CP*delCH_delTVA(:,N)* MAPL_VIREPS*TC(:,N) *(TC(:,N)-TA) ) + end do + + case (2,3) ! numerical derivatives of exchange coeffs w.r.t. canopy temp and specific humidity + + do N=1,NUM_SUBTILES + DEVSBT(:,N) = CQ(:,N) + max( 0.0, delCQ_delQC( :,N)* (QC(:,N)-QA) ) + DEDTC( :,N) = max( 0.0, delCH_delTC( :,N)* (QC(:,N)-QA) ) + DSHSBT(:,N) = MAPL_CP*( CH(:,N) + max( 0.0, delCH_delTC( :,N)* (TC(:,N)-TA) ) ) + DHSDQC(:,N) = max( 0.0, MAPL_CP*delCQ_delQC( :,N)* (TC(:,N)-TA) ) + end do + + case default + + _ASSERT(.false., 'unknown MOSFC_EXTRA_DERIVS_OFFL_LAND') + + end select + else + + ! GCM: Catchment coupled to atmosphere + do N=1,NUM_SUBTILES + CFT (:,N) = (CH(:,N)/CTATM) CFQ (:,N) = (CQ(:,N)/CQATM) + SHSBT (:,N) = (SH + DSH *(TC(:,N)-THATM))*CFT(:,N) EVSBT (:,N) = (EVAP+ DEVAP*(QC(:,N)-QHATM))*CFQ(:,N) - DSHSBT(:,N) = DSH *CFT(:,N) - DEVSBT(:,N) = DEVAP*CFQ(:,N) - ALWN(:,N)=ALW - BLWN(:,N)=BLW + DSHSBT(:,N) = DSH * CFT(:,N) + DEVSBT(:,N) = DEVAP* CFQ(:,N) + + ALWN (:,N) = ALW + BLWN (:,N) = BLW + end do - end if + + end if ! Catchment offline ! Compute DQS; make sure QC is between QA and QSAT; compute RA. ! @@ -5152,7 +5489,19 @@ subroutine Driver ( RC ) ! Some 40 lines below, duplicate code was present within #ifdef LAND_UPD block ! (later changed to "if (LAND_FIX)") and was removed in Jan 2022. ! - reichle, 14 Jan 2022. - + ! + ! reichle, 9/9/2024: + ! + ! WHY IS QC RESET *AFTER* IT WAS USED TO CALCULATE THE EXCHANGE COEFFS (AND DERIVS) ABOVE??? + ! + ! For reference, the following comment was copied from from LDASsa m3-16, specifically, from + ! reichle-LDASsa_m3-16_6/src/Components/GEOSlana_GridComp/process_cat.F90 (Lines 330-333): + ! + ! ! compute surface exchange coefficients etc BEFORE possibly resetting + ! ! profile of Qair-QAx-Qsat(surf) -- for consistency with two-stage + ! ! run-method in GEOS_CatchGridComp.F90 + ! ! reichle+qliu, 9 Oct 2008 + do N=1,NUM_SUBTILES DQS(:,N) = GEOS_DQSAT ( TC(:,N), PS, QSAT=QSAT(:,N), PASCALS=.true., RAMP=0.0 ) QC (:,N) = min(max(QA(:),QSAT(:,N)),QC(:,N)) @@ -5162,15 +5511,15 @@ subroutine Driver ( RC ) QC(:,FSNW) = QSAT(:,FSNW) - ! -------------------------------------------------------------------------- + ! -------------------------------------------------------------------------- ! get total solid precip ! -------------------------------------------------------------------------- SLDTOT = SNO+ICE+FRZR - ! protect the forcing from unsavory values, as per practice in offline - ! driver - ! -------------------------------------------------------------------------- + ! protect the forcing from unsavory values, as per practice in offline + ! driver + ! -------------------------------------------------------------------------- _ASSERT(count(PLS<0.)==0,'needs informative message') _ASSERT(count(PCU<0.)==0,'needs informative message') @@ -5218,25 +5567,25 @@ subroutine Driver ( RC ) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEDTC (:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, SHSBT (:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, DHSDQA(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, DHSDQC(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DSHSBT(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, EVSBT (:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEDTC (:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, SHSBT (:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, DHSDQA(:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, DHSDQC(:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DSHSBT(:,FTRN), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, EVSBT (:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEDTC (:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, SHSBT (:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, DHSDQA(:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, DHSDQC(:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DSHSBT(:,FWLT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, EVSBT (:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEDTC (:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, SHSBT (:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, DHSDQA(:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, DHSDQC(:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DSHSBT(:,FSNW), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, TA, mask=mask, rc=status); VERIFY_(STATUS) @@ -5627,13 +5976,13 @@ subroutine Driver ( RC ) UUU ,& EVSBT(:,FSAT), DEVSBT(:,FSAT), DEDTC(:,FSAT) ,& - SHSBT(:,FSAT), DHSDQA(:,FSAT), DSHSBT(:,FSAT),& + SHSBT(:,FSAT), DHSDQC(:,FSAT), DSHSBT(:,FSAT),& EVSBT(:,FTRN), DEVSBT(:,FTRN), DEDTC(:,FTRN) ,& - SHSBT(:,FTRN), DHSDQA(:,FTRN), DSHSBT(:,FTRN),& + SHSBT(:,FTRN), DHSDQC(:,FTRN), DSHSBT(:,FTRN),& EVSBT(:,FWLT), DEVSBT(:,FWLT), DEDTC(:,FWLT) ,& - SHSBT(:,FWLT), DHSDQA(:,FWLT), DSHSBT(:,FWLT),& + SHSBT(:,FWLT), DHSDQC(:,FWLT), DSHSBT(:,FWLT),& EVSBT(:,FSNW), DEVSBT(:,FSNW), DEDTC(:,FSNW) ,& - SHSBT(:,FSNW), DHSDQA(:,FSNW), DSHSBT(:,FSNW),& + SHSBT(:,FSNW), DHSDQC(:,FSNW), DSHSBT(:,FSNW),& TA , QA ,& @@ -6115,7 +6464,7 @@ subroutine Driver ( RC ) deallocate(EVSBT ) deallocate(DEVSBT ) deallocate(DEDTC ) - deallocate(DHSDQA ) + deallocate(DHSDQC ) deallocate(CFT ) deallocate(CFQ ) deallocate(TCO ) @@ -6181,13 +6530,16 @@ subroutine RUN0(gc, import, export, clock, rc) real, pointer :: ps(:)=>null() !! INTERNAL pointers - !! -asnow-emis-ww-fr- + !! -asnow-emis-ww-fr-D[xx] real, pointer :: asnow(:)=>null() real, pointer :: emis(:)=>null() real, pointer :: ww(:,:)=>null() real, pointer :: fr(:,:)=>null() - real, pointer :: DCQ(:,:)=>null() - real, pointer :: DCH(:,:)=>null() + real, pointer :: delCQ_delTVA(:,:)=>null() + real, pointer :: delCH_delTVA(:,:)=>null() + real, pointer :: delCH_delTC(:,:)=>null() + real, pointer :: delCQ_delQC(:,:)=>null() + !! -prognostic-variables- real, pointer :: tc(:,:)=>null() real, pointer :: qc(:,:)=>null() @@ -6267,10 +6619,23 @@ subroutine RUN0(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(INTERNAL, ww, 'WW', rc=status) VERIFY_(status) - call MAPL_GetPointer(INTERNAL, DCQ, 'DCQ', rc=status) - VERIFY_(status) - call MAPL_GetPointer(INTERNAL, DCH, 'DCH', rc=status) - VERIFY_(status) + + if (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND == 1) then + + call MAPL_GetPointer(INTERNAL, delCQ_delTVA, 'delCQ_delTVA', rc=status) + VERIFY_(status) + call MAPL_GetPointer(INTERNAL, delCH_delTVA, 'delCH_delTVA', rc=status) + VERIFY_(status) + + elseif (CATCH_INTERNAL_STATE%MOSFC_EXTRA_DERIVS_OFFL_LAND >= 2) then + + call MAPL_GetPointer(INTERNAL, delCH_delTC, 'delCH_delTC', rc=status) + VERIFY_(status) + call MAPL_GetPointer(INTERNAL, delCQ_delQC, 'delCQ_delQC', rc=status) + VERIFY_(status) + + end if + call MAPL_GetPointer(INTERNAL, tc, 'TC', rc=status) VERIFY_(status) call MAPL_GetPointer(INTERNAL, qc, 'QC', rc=status) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 index cb9af81e0..4c50fedcd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 @@ -12,7 +12,9 @@ module catch_wrap_stateMod type T_CATCH_STATE type (ESMF_FieldBundle) :: Bundle + ! logical :: LDAS_CORRECTOR + ! ! CATCH_OFFLINE: ! 0: DEFAULT for GCM, (WW,CH,CM,CQ,FR) are required in Catchment restart file ! 1: DEFAULT for GEOSldas, (WW,CH,CM,CQ,FR) are optional @@ -21,10 +23,11 @@ module catch_wrap_stateMod ! see also GEOSldas repo: src/Applications/LDAS_App/GEOSldas_LDAS.rc integer :: CATCH_OFFLINE ! - ! resource parameters from GEOS_SurfaceGridComp.rc integer :: CATCH_SPINUP + ! + ! some (but not all) resource parameters from GEOS_SurfaceGridComp.rc: integer :: USE_ASCATZ0, Z0_FORMULATION, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB - integer :: CHOOSEMOSFC, SNOW_ALBEDO_INFO + integer :: CHOOSEMOSFC, MOSFC_EXTRA_DERIVS_OFFL_LAND, SNOW_ALBEDO_INFO real :: SURFLAY real :: FWETC, FWETL logical :: USE_FWET_FOR_RUNOFF @@ -36,6 +39,7 @@ module catch_wrap_stateMod end type CATCH_WRAP type, extends(T_CATCH_STATE) :: T_CATCHCN_STATE + ! resource parameters from GEOS_SurfaceGridComp.rc: integer :: ATM_CO2, PRESCRIBE_DVG real :: CO2 integer :: CO2_YEAR_IN @@ -50,21 +54,79 @@ module catch_wrap_stateMod subroutine surface_params_to_wrap_state(statePtr, scf, rc) - ! *********************************************** ! - ! see GEOS_SurfaceGridComp.rc for documentation ! - ! *********************************************** ! + ! obtain resource variables from "SURFRC" file; they are ultimately stored in CATCH_INTERNAL_STATE or CATCHCN_INTERNAL_STATE - class(T_CATCH_STATE), pointer, intent(inout) :: statePtr - type(ESMF_Config), intent(inout) :: SCF - integer, optional, intent(out) :: rc + class(T_CATCH_STATE), pointer, intent(inout) :: statePtr + type(ESMF_Config), intent(inout) :: SCF + integer, optional, intent( out) :: rc + real :: FWETC_default, FWETL_default - integer:: status + integer:: status, ii - call MAPL_GetResource( SCF, statePtr%SURFLAY, label='SURFLAY:', DEFAULT=50., __RC__ ) - call MAPL_GetResource( SCF, statePtr%USE_ASCATZ0, label='USE_ASCATZ0:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%CHOOSEMOSFC, label='CHOOSEMOSFC:', DEFAULT=1, __RC__ ) - call MAPL_GetResource( SCF, statePtr%USE_FWET_FOR_RUNOFF, label='USE_FWET_FOR_RUNOFF:', DEFAULT=.FALSE., __RC__ ) - call MAPL_GetResource( SCF, statePtr%Z0_FORMULATION, label='Z0_FORMULATION:', DEFAULT=4, __RC__ ) + ! ************************************************* ! + ! For documentation, see GEOS_SurfaceGridComp.rc. ! + ! ************************************************* ! + + call MAPL_GetResource( SCF, statePtr%SURFLAY, label='SURFLAY:', DEFAULT=50., __RC__ ) + call MAPL_GetResource( SCF, statePtr%USE_ASCATZ0, label='USE_ASCATZ0:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%CHOOSEMOSFC, label='CHOOSEMOSFC:', DEFAULT=1, __RC__ ) + + ! MOSFC_EXTRA_DERIVS_OFFL_LAND: Resource parameter for *offline* (LDAS) mode. + ! + ! Over *land*, use derivatives of exchange coeffs w.r.t. temp. & humidity. + ! + ! 0 : None, default for Helfand. + ! 1 : Analytical derivs, default for Louis, *not* available for Helfand. + ! 2 : Numerical derivs. + ! 3 : Numerical derivs, via virtual temp., *not* available for Helfand, same as 2 but faster than 2. + ! + ! Runtimes: Helfand takes ~10 times longer than Louis. In offline mode, Helfand consumes + ! about as much CPU as Catchment. Option 2 triples the runtime of the MOSFC scheme. + ! Option 3 doubles the runtime of the Louis scheme. + + if (statePtr%CATCH_OFFLINE==0) then + + statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND = 0 ! must be 0 for GCM + + else + + ! offline (LDAS) mode; default for MOSFC_EXTRA_DERIVS_OFFL_LAND depends on CHOOSEMOSFC (Louis or Helfand) + + if (statePtr%CHOOSEMOSFC==0) then + + ! Louis + call MAPL_GetResource( SCF, statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND, label='MOSFC_EXTRA_DERIVS_OFFL_LAND:', DEFAULT=1, __RC__ ) + ! make sure parameter value is allowed + ii = statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND ; _ASSERT(ii>=0 .and. ii<=3, 'unknown MOSFC_EXTRA_DERIVS_OFFL_LAND for Louis ') + + elseif (statePtr%CHOOSEMOSFC==1) then + + ! Helfand + call MAPL_GetResource( SCF, statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND, label='MOSFC_EXTRA_DERIVS_OFFL_LAND:', DEFAULT=0, __RC__ ) + ! make sure parameter value is allowed (analytical derivs not implemented for Helfand) + ii = statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND ; _ASSERT(ii==0 .or. ii==2, 'unknown MOSFC_EXTRA_DERIVS_OFFL_LAND for Helfand') + + else + + _ASSERT(.FALSE.,'unknown CHOOSEMOSFC') + + end if + + end if + + ! for CatchCN, must have MOSFC_EXTRA_DERIVS_OFFL_LAND<=1 (numerical derivatives not yet implemented for CatchCN) + + select type (statePtr) + type is (T_CATCHCN_STATE) ! CATCHCN + + _ASSERT( statePtr%MOSFC_EXTRA_DERIVS_OFFL_LAND<=1, 'selected choice for MOSFC_EXTRA_DERIVS_OFFL_LAND not yet implemented for CatchCN') + + end select + + ! ------------------------- + + call MAPL_GetResource( SCF, statePtr%USE_FWET_FOR_RUNOFF, label='USE_FWET_FOR_RUNOFF:', DEFAULT=.FALSE., __RC__ ) + call MAPL_GetResource( SCF, statePtr%Z0_FORMULATION, label='Z0_FORMULATION:', DEFAULT=4, __RC__ ) if (.NOT. statePtr%USE_FWET_FOR_RUNOFF) then FWETC_default = 0.02 @@ -73,22 +135,23 @@ subroutine surface_params_to_wrap_state(statePtr, scf, rc) FWETC_default = 0.005 ! NOT ready for science! FWETL_default = 0.025 ! NOT ready for science! endif - call MAPL_GetResource( SCF, statePtr%FWETC, label='FWETC:', DEFAULT=FWETC_default, __RC__ ) - call MAPL_GetResource( SCF, statePtr%FWETL, label='FWETL:', DEFAULT=FWETL_default, __RC__ ) - call MAPL_GetResource( SCF, statePtr%SNOW_ALBEDO_INFO, label='SNOW_ALBEDO_INFO:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%RUN_IRRIG, label='RUN_IRRIG:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%IRRIG_METHOD, label='IRRIG_METHOD:', DEFAULT=0, __RC__ ) + + call MAPL_GetResource( SCF, statePtr%FWETC, label='FWETC:', DEFAULT=FWETC_default, __RC__ ) + call MAPL_GetResource( SCF, statePtr%FWETL, label='FWETL:', DEFAULT=FWETL_default, __RC__ ) + call MAPL_GetResource( SCF, statePtr%SNOW_ALBEDO_INFO, label='SNOW_ALBEDO_INFO:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%RUN_IRRIG, label='RUN_IRRIG:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%IRRIG_METHOD, label='IRRIG_METHOD:', DEFAULT=0, __RC__ ) select type (statePtr) type is (T_CATCHCN_STATE) ! CATCHCN - call MAPL_GetResource( SCF, statePtr%DTCN, label='DTCN:', DEFAULT=5400., __RC__ ) - call MAPL_GetResource( SCF, statePtr%ATM_CO2, label='ATM_CO2:', DEFAULT=2, __RC__ ) - call MAPL_GetResource( SCF, statePtr%PRESCRIBE_DVG, label='PRESCRIBE_DVG:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%CO2, label='CO2:', DEFAULT=350.e-6, __RC__ ) - call MAPL_GetResource( SCF, statePtr%CO2_YEAR_IN, label='CO2_YEAR:', DEFAULT=-9999, __RC__ ) + call MAPL_GetResource( SCF, statePtr%DTCN, label='DTCN:', DEFAULT=5400., __RC__ ) + call MAPL_GetResource( SCF, statePtr%ATM_CO2, label='ATM_CO2:', DEFAULT=2, __RC__ ) + call MAPL_GetResource( SCF, statePtr%PRESCRIBE_DVG, label='PRESCRIBE_DVG:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%CO2, label='CO2:', DEFAULT=350.e-6, __RC__ ) + call MAPL_GetResource( SCF, statePtr%CO2_YEAR_IN, label='CO2_YEAR:', DEFAULT=-9999, __RC__ ) end select diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index ff9d8fca8..1c05fae75 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -29,13 +29,20 @@ # # #################################################################################### -# ---- Surface layer turbulence scheme +# ---- Monin-Obukhov (MO) surface layer turbulence scheme # # 0 : Louis (MERRA, Fortuna-DAS, SMAP NRv4/4.1/5/7.2/8.1/9.1/10.0) -# 1 : Helfand Monin-Obukhov (Fortuna-AR5, Ganymed, Heracles, Icarus-3_2, MERRA-2) +# 1 : Helfand (Fortuna-AR5, Ganymed, Heracles, Icarus-3_2, MERRA-2) # -# GEOSagcm=>CHOOSEMOSFC: 1 -# GEOSldas=>CHOOSEMOSFC: 0 +# Note: For *offline* simulations, optional use of extra MO derivatives is supported +# through rc parameter MOSFC_EXTRA_DERIVS_LAND (see catch_wrap_state.F90). +# +# GEOSagcm=>CHOOSEMOSFC: 1 +# GEOSagcm=>MOSFC_EXTRA_DERIVS_OFFL_LAND: 0 +# +# GEOSldas=>CHOOSEMOSFC: 0 +# GEOSldas=>MOSFC_EXTRA_DERIVS_OFFL_LAND: 1 + # ---- Thickness of surface layer for soil moisture [mm] #