diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOS_SurfaceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOS_SurfaceGridComp.F90 index 2f2a2fbd4..ec1fc9805 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOS_SurfaceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOS_SurfaceGridComp.F90 @@ -2697,6 +2697,24 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'depth_to_water_table_from_surface',& + UNITS = 'm' ,& + SHORT_NAME = 'WATERTABLED' ,& + DIMS = MAPL_DimsHorzOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'change_in_free_surface_water_reservoir_on_peat',& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'FSWCHANGE' ,& + DIMS = MAPL_DimsHorzOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + IF(LSM_CHOICE > 1) THEN call MAPL_AddExportSpec(GC ,& LONG_NAME = 'CN_exposed_leaf-area_index',& @@ -5052,24 +5070,26 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) real, pointer, dimension(:,:) :: T2MWET => NULL() ! GOSWIM (internal/export variables from catch/catchcn) - real, pointer, dimension(:,:,:) :: RDU001 => NULL() - real, pointer, dimension(:,:,:) :: RDU002 => NULL() - real, pointer, dimension(:,:,:) :: RDU003 => NULL() - real, pointer, dimension(:,:,:) :: RDU004 => NULL() - real, pointer, dimension(:,:,:) :: RDU005 => NULL() - real, pointer, dimension(:,:,:) :: RBC001 => NULL() - real, pointer, dimension(:,:,:) :: RBC002 => NULL() - real, pointer, dimension(:,:,:) :: ROC001 => NULL() - real, pointer, dimension(:,:,:) :: ROC002 => NULL() - real, pointer, dimension(:,:) :: RMELTDU001 => NULL() - real, pointer, dimension(:,:) :: RMELTDU002 => NULL() - real, pointer, dimension(:,:) :: RMELTDU003 => NULL() - real, pointer, dimension(:,:) :: RMELTDU004 => NULL() - real, pointer, dimension(:,:) :: RMELTDU005 => NULL() - real, pointer, dimension(:,:) :: RMELTBC001 => NULL() - real, pointer, dimension(:,:) :: RMELTBC002 => NULL() - real, pointer, dimension(:,:) :: RMELTOC001 => NULL() - real, pointer, dimension(:,:) :: RMELTOC002 => NULL() + real, pointer, dimension(:,:,:) :: RDU001 => NULL() + real, pointer, dimension(:,:,:) :: RDU002 => NULL() + real, pointer, dimension(:,:,:) :: RDU003 => NULL() + real, pointer, dimension(:,:,:) :: RDU004 => NULL() + real, pointer, dimension(:,:,:) :: RDU005 => NULL() + real, pointer, dimension(:,:,:) :: RBC001 => NULL() + real, pointer, dimension(:,:,:) :: RBC002 => NULL() + real, pointer, dimension(:,:,:) :: ROC001 => NULL() + real, pointer, dimension(:,:,:) :: ROC002 => NULL() + real, pointer, dimension(:,:) :: RMELTDU001 => NULL() + real, pointer, dimension(:,:) :: RMELTDU002 => NULL() + real, pointer, dimension(:,:) :: RMELTDU003 => NULL() + real, pointer, dimension(:,:) :: RMELTDU004 => NULL() + real, pointer, dimension(:,:) :: RMELTDU005 => NULL() + real, pointer, dimension(:,:) :: RMELTBC001 => NULL() + real, pointer, dimension(:,:) :: RMELTBC002 => NULL() + real, pointer, dimension(:,:) :: RMELTOC001 => NULL() + real, pointer, dimension(:,:) :: RMELTOC002 => NULL() + real, pointer, dimension(:,:) :: WATERTABLED => NULL() + real, pointer, dimension(:,:) :: FSWCHANGE => NULL() ! CN model real, pointer, dimension(:,:) :: CNLAI => NULL() @@ -5329,6 +5349,8 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) real, pointer, dimension(:) :: RMELTBC002TILE => NULL() real, pointer, dimension(:) :: RMELTOC001TILE => NULL() real, pointer, dimension(:) :: RMELTOC002TILE => NULL() + real, pointer, dimension(:) :: WATERTABLEDTILE => NULL() + real, pointer, dimension(:) :: FSWCHANGETILE => NULL() real, pointer, dimension(:) :: CNLAITILE => NULL() real, pointer, dimension(:) :: CNTLAITILE => NULL() @@ -6141,6 +6163,9 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_GetPointer(EXPORT , RMELTBC002 , 'RMELTBC002', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT , RMELTOC001 , 'RMELTOC001', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT , RMELTOC002 , 'RMELTOC002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT , WATERTABLED, 'WATERTABLED', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT , FSWCHANGE , 'FSWCHANGE', RC=STATUS); VERIFY_(STATUS) + IF(LSM_CHOICE > 1) THEN call MAPL_GetPointer(EXPORT , CNLAI , 'CNLAI' , RC=STATUS); VERIFY_(STATUS) @@ -6714,6 +6739,8 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) call MKTILE(RMELTBC002 ,RMELTBC002TILE ,NT,RC=STATUS); VERIFY_(STATUS) call MKTILE(RMELTOC001 ,RMELTOC001TILE ,NT,RC=STATUS); VERIFY_(STATUS) call MKTILE(RMELTOC002 ,RMELTOC002TILE ,NT,RC=STATUS); VERIFY_(STATUS) + call MKTILE(WATERTABLED,WATERTABLEDTILE,NT,RC=STATUS); VERIFY_(STATUS) + call MKTILE(FSWCHANGE ,FSWCHANGETILE ,NT,RC=STATUS); VERIFY_(STATUS) IF (LSM_CHOICE > 1) THEN call MKTILE(CNLAI ,CNLAITILE ,NT,RC=STATUS); VERIFY_(STATUS) @@ -7532,15 +7559,17 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) if(associated(ROC002)) call MAPL_LocStreamTransform( LOCSTREAM,ROC002(:,:,N) ,ROC002TILE(:,N), RC=STATUS); VERIFY_(STATUS) END DO - if(associated(RMELTDU001))call MAPL_LocStreamTransform( LOCSTREAM,RMELTDU001 ,RMELTDU001TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTDU002))call MAPL_LocStreamTransform( LOCSTREAM,RMELTDU002 ,RMELTDU002TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTDU003))call MAPL_LocStreamTransform( LOCSTREAM,RMELTDU003 ,RMELTDU003TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTDU004))call MAPL_LocStreamTransform( LOCSTREAM,RMELTDU004 ,RMELTDU004TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTDU005))call MAPL_LocStreamTransform( LOCSTREAM,RMELTDU005 ,RMELTDU005TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTBC001))call MAPL_LocStreamTransform( LOCSTREAM,RMELTBC001 ,RMELTBC001TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTBC002))call MAPL_LocStreamTransform( LOCSTREAM,RMELTBC002 ,RMELTBC002TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTOC001))call MAPL_LocStreamTransform( LOCSTREAM,RMELTOC001 ,RMELTOC001TILE, RC=STATUS); VERIFY_(STATUS) - if(associated(RMELTOC002))call MAPL_LocStreamTransform( LOCSTREAM,RMELTOC002 ,RMELTOC002TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTDU001 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTDU001 ,RMELTDU001TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTDU002 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTDU002 ,RMELTDU002TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTDU003 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTDU003 ,RMELTDU003TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTDU004 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTDU004 ,RMELTDU004TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTDU005 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTDU005 ,RMELTDU005TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTBC001 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTBC001 ,RMELTBC001TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTBC002 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTBC002 ,RMELTBC002TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTOC001 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTOC001 ,RMELTOC001TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(RMELTOC002 ))call MAPL_LocStreamTransform(LOCSTREAM,RMELTOC002 ,RMELTOC002TILE, RC=STATUS); VERIFY_(STATUS) + if(associated(WATERTABLED))call MAPL_LocStreamTransform(LOCSTREAM,WATERTABLED,WATERTABLEDTILE,RC=STATUS); VERIFY_(STATUS) + if(associated(FSWCHANGE ))call MAPL_LocStreamTransform(LOCSTREAM,FSWCHANGE ,FSWCHANGETILE, RC=STATUS); VERIFY_(STATUS) if(associated(CNLAI)) then call MAPL_LocStreamTransform( LOCSTREAM,CNLAI ,CNLAITILE , RC=STATUS) @@ -8074,6 +8103,8 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) if(associated(RMELTBC002TILE )) deallocate(RMELTBC002TILE ) if(associated(RMELTOC001TILE )) deallocate(RMELTOC001TILE ) if(associated(RMELTOC002TILE )) deallocate(RMELTOC002TILE ) + if(associated(WATERTABLEDTILE)) deallocate(WATERTABLEDTILE) + if(associated(FSWCHANGETILE )) deallocate(FSWCHANGETILE ) if(associated(CNLAITILE )) deallocate(CNLAITILE ) if(associated(CNTLAITILE )) deallocate(CNTLAITILE ) if(associated(CNSAITILE )) deallocate(CNSAITILE ) @@ -8409,6 +8440,10 @@ subroutine DOTYPE(type,RC) VERIFY_(STATUS) call MAPL_GetPointer(GEX(type), dum, 'RMELTOC002' , ALLOC=associated(RMELTOC002TILE ), notFoundOK=.true., RC=STATUS) VERIFY_(STATUS) + call MAPL_GetPointer(GEX(type), dum, 'WATERTABLED', ALLOC=associated(WATERTABLEDTILE ),notFoundOK=.true., RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(GEX(type), dum, 'FSWCHANGE' , ALLOC=associated(FSWCHANGETILE ) , notFoundOK=.true., RC=STATUS) + VERIFY_(STATUS) IF (LSM_CHOICE > 1) THEN call MAPL_GetPointer(GEX(type), dum, 'CNLAI' , ALLOC=associated(CNLAITILE ), notFoundOK=.true., RC=STATUS) @@ -8977,6 +9012,8 @@ subroutine DOTYPE(type,RC) if(associated(RMELTBC002TILE)) call FILLOUT_TILE(GEX(type), 'RMELTBC002' , RMELTBC002TILE , XFORM, RC=STATUS);VERIFY_(STATUS) if(associated(RMELTOC001TILE)) call FILLOUT_TILE(GEX(type), 'RMELTOC001' , RMELTOC001TILE , XFORM, RC=STATUS);VERIFY_(STATUS) if(associated(RMELTOC002TILE)) call FILLOUT_TILE(GEX(type), 'RMELTOC002' , RMELTOC002TILE , XFORM, RC=STATUS);VERIFY_(STATUS) + if(associated(WATERTABLEDTILE))call FILLOUT_TILE(GEX(type), 'WATERTABLED', WATERTABLEDTILE, XFORM, RC=STATUS);VERIFY_(STATUS) + if(associated(FSWCHANGETILE)) call FILLOUT_TILE(GEX(type), 'FSWCHANGE' , FSWCHANGETILE , XFORM, RC=STATUS);VERIFY_(STATUS) if(associated(CNLAITILE)) then call FILLOUT_TILE(GEX(type), 'CNLAI' , CNLAITILE , XFORM, RC=STATUS) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 index fd981de57..bfca66893 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 @@ -920,6 +920,8 @@ subroutine SetServices ( GC, RC ) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTBC002', CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC001', CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC002', CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WATERTABLED',CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'FSWCHANGE', CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) if (DO_GOSWIM /= 0) then call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RDU001', CHILD_ID = CATCH(1), RC=STATUS) ; VERIFY_(STATUS) @@ -1290,6 +1292,8 @@ subroutine SetServices ( GC, RC ) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTBC002', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC001', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC002', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WATERTABLED',CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'FSWCHANGE', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) if (DO_GOSWIM /= 0) then call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RDU001', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 index 1e385acb2..e790bcd5f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 @@ -984,6 +984,8 @@ subroutine SetServices ( GC, RC ) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTBC002', CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC001', CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RMELTOC002', CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WATERTABLED',CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'FSWCHANGE' , CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) if (DO_GOSWIM /= 0) then call MAPL_AddExportSpec ( GC, SHORT_NAME = 'RDU001', CHILD_ID = CATCHCN, RC=STATUS) ; VERIFY_(STATUS) 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 c4f373144..32173c29d 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 @@ -61,8 +61,9 @@ module GEOS_CatchCNCLM40GridCompMod use MAPL_ConstantsMod,only: Tzero => MAPL_TICE, pi => MAPL_PI use clm_time_manager, only: get_days_per_year, get_step_size use pftvarcon, only: noveg - USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist, irrigation_rate, & - gndtmp + USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist, & + catch_calc_zbar, catch_calc_watertabled, irrigation_rate, & + gndtmp implicit none private @@ -3732,7 +3733,25 @@ subroutine SetServices ( GC, RC ) SHORT_NAME = 'RMELTOC002' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'depth_to_water_table_from_surface',& + UNITS = 'm' ,& + SHORT_NAME = 'WATERTABLED' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'change_in_free_surface_water_reservoir_on_peat',& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'FSWCHANGE' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) VERIFY_(STATUS) !EOS @@ -4892,6 +4911,8 @@ subroutine Driver ( RC ) real, pointer, dimension(:) :: RMELTOC001 real, pointer, dimension(:) :: RMELTOC002 real, pointer, dimension(:) :: IRRIGRATE + real, pointer, dimension(:) :: WATERTABLED + real, pointer, dimension(:) :: FSWCHANGE ! -------------------------------------------------------------------------- ! Local pointers for tile variables @@ -4920,7 +4941,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: ghflxsno, ghflxtskin real,pointer,dimension(:) :: SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 real,pointer,dimension(:) :: LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW - real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1 + real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1, FSW_CHANGE real,pointer,dimension(:) :: WCHANGE, ECHANGE, HSNACC, EVACC, SHACC real,pointer,dimension(:) :: SNOVR, SNOVF, SNONR, SNONF real,pointer,dimension(:) :: VSUVR, VSUVF @@ -5528,15 +5549,17 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,CNLAI41, 'CNLAI41' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNLAI42, 'CNLAI42' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNLAI43, 'CNLAI43' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU001,'RMELTDU001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU002,'RMELTDU002', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU003,'RMELTDU003', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU004,'RMELTDU004', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU005,'RMELTDU005', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTBC001,'RMELTBC001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTBC002,'RMELTBC002', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTOC001,'RMELTOC001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTOC002,'RMELTOC002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU001,'RMELTDU001', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU002,'RMELTDU002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU003,'RMELTDU003', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU004,'RMELTDU004', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU005,'RMELTDU005', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTBC001,'RMELTBC001', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTBC002,'RMELTBC002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTOC001,'RMELTOC001', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTOC002,'RMELTOC002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WATERTABLED ,'WATERTABLED', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FSWCHANGE , 'FSWCHANGE', RC=STATUS); VERIFY_(STATUS) IF (RUN_IRRIG /= 0) call MAPL_GetPointer(EXPORT,IRRIGRATE ,'IRRIGRATE' , RC=STATUS); VERIFY_(STATUS) NTILES = size(PS) @@ -5871,6 +5894,7 @@ subroutine Driver ( RC ) allocate(fveg2 (NTILES)) allocate(FICE1 (NTILES)) allocate(SLDTOT (NTILES)) + allocate(FSW_CHANGE(NTILES)) allocate(SHSBT (NTILES,NUM_SUBTILES)) allocate(DSHSBT (NTILES,NUM_SUBTILES)) @@ -6418,9 +6442,9 @@ subroutine Driver ( RC ) ! gkw: obtain catchment area fractions and soil moisture ! ------------------------------------------------------ -call catch_calc_soil_moist( ntiles, veg1, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & - srfexc, rzexc, catdef, car1, car2, car4, sfmc, rzmc, prmc ) + call catch_calc_soil_moist( ntiles, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & + ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, bf1, bf2, & + srfexc, rzexc, catdef, car1, car2, car4, sfmc, rzmc, prmc ) ! obtain saturated canopy resistance following Farquhar, CLM4 implementation @@ -6546,9 +6570,11 @@ subroutine Driver ( RC ) ! soil temperatures ! ----------------- - zbar = -sqrt(1.e-20+catdef(n)/bf1(n))+bf2(n) + + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTMP) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) HT(:)=GHTCNT(:,N) - CALL GNDTMP(poros(n),zbar,ht,frice,tp,soilice) + CALL GNDTMP(poros(n),-1.*zbar,ht,frice,tp,soilice) ! note minus sign for zbar ! At the CatchCNGridComp level, tp1, tp2, .., tp6 are export variables in units of Kelvin, ! - rreichle & borescan, 6 Nov 2020 @@ -6567,7 +6593,6 @@ subroutine Driver ( RC ) ! baseflow ! -------- - zbar = sqrt(1.e-20+catdef(n)/bf1(n))-bf2(n) bflow(n) = (1.-frice)*1000.* & cond(n)*exp(-(bf3(n)-ashift)-gnu(n)*zbar)/gnu(n) IF(catdef(n) >= cdcr1(n)) bflow(n) = 0. @@ -7294,9 +7319,9 @@ subroutine Driver ( RC ) IF ((RUN_IRRIG /= 0).AND.(ntiles >0)) THEN - CALL CATCH_CALC_SOIL_MOIST ( & - NTILES,VEG1,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & + CALL CATCH_CALC_SOIL_MOIST ( & + NTILES,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, CAR1, CAR2, CAR4, sfmc, rzmc, prmc) call irrigation_rate (IRRIG_METHOD, & @@ -7557,7 +7582,7 @@ subroutine Driver ( RC ) TSURF ,& SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 ,& LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW ,& - TCSORIG1, TPSN1IN1, TPSN1OUT1 ,& + TCSORIG1, TPSN1IN1, TPSN1OUT1, FSW_CHANGE ,& TC1_0=TC1_0, TC2_0=TC2_0, TC4_0=TC4_0 ,& QA1_0=QA1_0, QA2_0=QA2_0, QA4_0=QA4_0 ,& RCONSTIT=RCONSTIT, RMELT=RMELT, TOTDEPOS=TOTDEPOS, LHACC=LHACC) @@ -7781,6 +7806,10 @@ subroutine Driver ( RC ) if(associated(RMELTBC002)) RMELTBC002 = RMELT(:,7) if(associated(RMELTOC001)) RMELTOC001 = RMELT(:,8) if(associated(RMELTOC002)) RMELTOC002 = RMELT(:,9) + if(associated(FSWCHANGE)) FSWCHANGE = FSW_CHANGE + if(associated(WATERTABLED)) then + WATERTABLED = catch_calc_watertabled( BF1, BF2, CDCR2, POROS, WPWET, CATDEF ) + endif if(associated(TPSN1OUT)) then where(WESNN(1,:)>0.) @@ -7965,7 +7994,8 @@ subroutine Driver ( RC ) deallocate(TOTDEPOS ) deallocate(RMELT ) deallocate(FICE1 ) - deallocate(SLDTOT ) + deallocate(SLDTOT ) + deallocate(FSW_CHANGE) deallocate( btran ) deallocate( wgt ) deallocate( bt1 ) @@ -8482,11 +8512,11 @@ subroutine RUN0(gc, import, export, clock, rc) rzexccp = rzexc call catch_calc_soil_moist( & ! intent(in) - ntiles, nint(veg1), dzsf, vgwmax, cdcr1, cdcr2, & + ntiles, dzsf, vgwmax, cdcr1, cdcr2, & psis, bee, poros, wpwet, & ars1, ars2, ars3, & ara1, ara2, ara3, ara4, & - arw1, arw2, arw3, arw4, & + arw1, arw2, arw3, arw4, bf1, bf2, & ! intent(inout) ! from process_cat srfexccp, rzexccp, catdefcp, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/CLM45/CN_DriverMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/CLM45/CN_DriverMod.F90 index 8cad3d43a..a9769a60f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/CLM45/CN_DriverMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/CLM45/CN_DriverMod.F90 @@ -29,7 +29,7 @@ module CN_DriverMod #else use CNDecompCascadeMod_CENTURY, only : init_decompcascade #endif - use catch_constants, only: DZTC=>CATCH_DZTC, DZGT=>CATCH_DZGT + use catch_constants, only: DZTSURF=>CATCH_DZTSURF, DZGT=>CATCH_DZGT use SurfParams, only: LAND_FIX ! use update_model_para4cn, only : LocalTileID, upd_tileid ! useful for debugging @@ -230,7 +230,7 @@ subroutine CN_Driver(istep,nch,nveg,nzone,daylength, & real, pointer :: psisat(:,:) !soil water potential at saturation for CN code (MPa) real, pointer :: psiwilt(:) !root-zone soil water potential at wilting point (MPa) real, pointer :: soilpsi(:,:) !soil water potential in each soil layer (MPa) - real, pointer :: h2osoi_liq(:,:) !column liquid water (kg/m2) (new) + real, pointer :: h2osoi_liq(:,:) !column liquid water (kg/m2) (new) real, pointer :: wf(:) !soil water as frac. of whc for top 0.05 m real, pointer :: wf2(:) !soil water as frac. of whc for top 0.17 m real, pointer :: qflx_drain(:) !sub-surface runoff (mm H2O /s) @@ -639,7 +639,7 @@ subroutine CN_Driver(istep,nch,nveg,nzone,daylength, & ! ---------------- t_soisno(n,1) = tp1(nc) ! soil layer temperature (K) t_grnd(n) = tgw(nc,nz) ! ground surface temperature (K) - tsoi17(n) = (DZTC*tgw(nc,nz)+(DZGT(1)-DZTC)*tp1(nc)+(0.17-DZGT(1))*tp2(nc))/0.17 ! soil temperature in top 17cm of soil (Kelvin) + tsoi17(n) = (DZTSURF*tgw(nc,nz)+(DZGT(1)-DZTSURF)*tp1(nc)+(0.17-DZGT(1))*tp2(nc))/0.17 ! soil temperature in top 17cm of soil (Kelvin) ! fzeng: tgw is for the top 5cm; tp1 is for the 2nd 5cm; tp2 is for the next 10cm ! see Koster et al., 2000, JGR ! The depths are hard coded here. Improve this? 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 77dcc8b6a..5caaa0425 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 @@ -61,8 +61,10 @@ module GEOS_CatchCNCLM45GridCompMod use MAPL_ConstantsMod,only: Tzero => MAPL_TICE, pi => MAPL_PI use clm_time_manager, only: get_days_per_year, get_step_size use pftvarcon, only: noveg - USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist, irrigation_rate, & - gndtmp + USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist, & + catch_calc_zbar, catch_calc_watertabled, irrigation_rate, & + gndtmp + use update_model_para4cn, only : upd_curr_date_time implicit none @@ -3669,9 +3671,28 @@ subroutine SetServices ( GC, RC ) SHORT_NAME = 'RMELTOC002' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) + RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'depth_to_water_table_from_surface',& + UNITS = 'm' ,& + SHORT_NAME = 'WATERTABLED' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'change_in_free_surface_water_reservoir_on_peat',& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'FSWCHANGE' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + !EOS call MAPL_TimerAdd(GC, name="RUN1" ,RC=STATUS) @@ -4838,6 +4859,8 @@ subroutine Driver ( RC ) real, pointer, dimension(:) :: RMELTOC001 real, pointer, dimension(:) :: RMELTOC002 real, pointer, dimension(:) :: IRRIGRATE + real, pointer, dimension(:) :: WATERTABLED + real, pointer, dimension(:) :: FSWCHANGE ! -------------------------------------------------------------------------- ! Local pointers for tile variables @@ -4866,7 +4889,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: ghflxsno, ghflxtskin real,pointer,dimension(:) :: SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 real,pointer,dimension(:) :: LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW - real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1 + real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1, FSW_CHANGE real,pointer,dimension(:) :: WCHANGE, ECHANGE, HSNACC, EVACC, SHACC real,pointer,dimension(:) :: SNOVR, SNOVF, SNONR, SNONF real,pointer,dimension(:) :: VSUVR, VSUVF @@ -5390,116 +5413,116 @@ subroutine Driver ( RC ) ! EXPORT POINTERS ! ----------------------------------------------------- - call MAPL_GetPointer(EXPORT,EVAPOUT,'EVAPOUT',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SUBLIM,'SUBLIM',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SHOUT, 'SHOUT' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RUNOFF, 'RUNOFF' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVPINT, 'EVPINT' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVPSOI, 'EVPSOI' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVPVEG, 'EVPVEG' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVPICE, 'EVPICE' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WAT10CM,'WAT10CM',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WATSOI, 'WATSOI' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ICESOI, 'ICESOI' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVPSNO, 'EVPSNO' ,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,BFLOW, 'BASEFLOW',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RUNSURF,'RUNSURF',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SMELT, 'SMELT' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,HLWUP, 'HLWUP' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SWNDSRF,'SWNDSRF',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LWNDSRF,'LWNDSRF',ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,HLATN, 'HLATN' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,QINFIL, 'QINFIL' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,AR1, 'AR1' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,AR2, 'AR2' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RZEQ, 'RZEQ' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,GHFLX, 'GHFLX' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPSURF, 'TPSURF' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPSN1, 'TPSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPUST, 'TPUNST' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPSAT, 'TPSAT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPWLT, 'TPWLT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ASNOW, 'ASNOW' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SHSNOW, 'SHSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,AVETSNOW,'AVETSNOW', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,FRSAT, 'FRSAT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,FRUST, 'FRUST' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,FRWLT, 'FRWLT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP1, 'TP1' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP2, 'TP2' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP3, 'TP3' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP4, 'TP4' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP5, 'TP5' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TP6, 'TP6' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EMIS, 'EMIS' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ALBVR, 'ALBVR' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ALBVF, 'ALBVF' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ALBNR, 'ALBNR' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ALBNF, 'ALBNF' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DELTS, 'DELTS' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DELQS, 'DELQS' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TST , 'TST' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,QST , 'QST' ,ALLOC=.true.,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LST , 'LST' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WET1 , 'WET1' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WET2 , 'WET2' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WET3 , 'WET3' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WCSF , 'WCSF' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WCRZ , 'WCRZ' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,WCPR , 'WCPR' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,ACCUM, 'ACCUM' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SNOMAS,'SNOWMASS', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SNOWDP, 'SNOWDP' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,EVLAND, 'EVLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,PRLAND, 'PRLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SNOLAND, 'SNOLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DRPARLAND, 'DRPARLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DFPARLAND, 'DFPARLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LHSNOW, 'LHSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SWNETSNOW1, 'SWNETSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LWUPSNOW, 'LWUPSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LWDNSNOW, 'LWDNSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TCSORIG, 'TCSORIG' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPSN1IN, 'TPSN1IN' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TPSN1OUT, 'TPSN1OUT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LHLAND, 'LHLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SHLAND, 'SHLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SWLAND, 'SWLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SWDOWNLAND, 'SWDOWNLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,LWLAND, 'LWLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,GHLAND, 'GHLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,GHSNOW, 'GHSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,GHTSKIN,'GHTSKIN', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SMLAND, 'SMLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TWLAND, 'TWLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TELAND, 'TELAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,TSLAND, 'TSLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DWLAND, 'DWLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,DHLAND, 'DHLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SPLAND, 'SPLAND' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SPWATR, 'SPWATR' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SPSNOW, 'SPSNOW' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNLAI, 'CNLAI' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNTLAI, 'CNTLAI' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNSAI, 'CNSAI' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNTOTC, 'CNTOTC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNVEGC, 'CNVEGC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNFROOTC,'CNFROOTC', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNNPP, 'CNNPP' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNGPP, 'CNGPP' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNSR, 'CNSR' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNNEE, 'CNNEE' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNXSMR, 'CNXSMR' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNADD, 'CNADD' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNLOSS, 'CNLOSS' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNBURN, 'CNBURN' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,PARABS, 'PARABS' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,PARINC, 'PARINC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SCSAT, 'SCSAT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SCUNS, 'SCUNS' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,BTRANT, 'BTRANT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,SIF, 'SIF' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNCO2, 'CNCO2' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVAPOUT , 'EVAPOUT',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SUBLIM , 'SUBLIM' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SHOUT , 'SHOUT' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RUNOFF , 'RUNOFF' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVPINT , 'EVPINT' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVPSOI , 'EVPSOI' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVPVEG , 'EVPVEG' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVPICE , 'EVPICE' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WAT10CM , 'WAT10CM',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WATSOI , 'WATSOI' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ICESOI , 'ICESOI' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVPSNO , 'EVPSNO' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,BFLOW , 'BASEFLOW',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RUNSURF , 'RUNSURF',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SMELT , 'SMELT' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,HLWUP , 'HLWUP' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SWNDSRF , 'SWNDSRF',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LWNDSRF , 'LWNDSRF',ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,HLATN , 'HLATN' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,QINFIL , 'QINFIL' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,AR1 , 'AR1' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,AR2 , 'AR2' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RZEQ , 'RZEQ' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,GHFLX , 'GHFLX' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPSURF , 'TPSURF' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPSN1 , 'TPSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPUST , 'TPUNST' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPSAT , 'TPSAT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPWLT , 'TPWLT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ASNOW , 'ASNOW' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SHSNOW , 'SHSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,AVETSNOW , 'AVETSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FRSAT , 'FRSAT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FRUST , 'FRUST' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FRWLT , 'FRWLT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP1 , 'TP1' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP2 , 'TP2' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP3 , 'TP3' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP4 , 'TP4' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP5 , 'TP5' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TP6 , 'TP6' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EMIS , 'EMIS' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ALBVR , 'ALBVR' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ALBVF , 'ALBVF' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ALBNR , 'ALBNR' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ALBNF , 'ALBNF' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DELTS , 'DELTS' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DELQS , 'DELQS' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TST , 'TST' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,QST , 'QST' ,ALLOC=.true., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LST , 'LST' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WET1 , 'WET1' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WET2 , 'WET2' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WET3 , 'WET3' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WCSF , 'WCSF' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WCRZ , 'WCRZ' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WCPR , 'WCPR' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,ACCUM , 'ACCUM' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SNOMAS , 'SNOWMASS' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SNOWDP , 'SNOWDP' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,EVLAND , 'EVLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,PRLAND , 'PRLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SNOLAND , 'SNOLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DRPARLAND , 'DRPARLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DFPARLAND , 'DFPARLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LHSNOW , 'LHSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SWNETSNOW1 , 'SWNETSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LWUPSNOW , 'LWUPSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LWDNSNOW , 'LWDNSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TCSORIG , 'TCSORIG' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPSN1IN , 'TPSN1IN' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TPSN1OUT , 'TPSN1OUT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LHLAND , 'LHLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SHLAND , 'SHLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SWLAND , 'SWLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SWDOWNLAND , 'SWDOWNLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,LWLAND , 'LWLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,GHLAND , 'GHLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,GHSNOW , 'GHSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,GHTSKIN , 'GHTSKIN' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SMLAND , 'SMLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TWLAND , 'TWLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TELAND , 'TELAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,TSLAND , 'TSLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DWLAND , 'DWLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,DHLAND , 'DHLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SPLAND , 'SPLAND' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SPWATR , 'SPWATR' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SPSNOW , 'SPSNOW' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNLAI , 'CNLAI' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNTLAI , 'CNTLAI' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNSAI , 'CNSAI' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNTOTC , 'CNTOTC' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNVEGC , 'CNVEGC' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNFROOTC , 'CNFROOTC' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNNPP , 'CNNPP' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNGPP , 'CNGPP' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNSR , 'CNSR' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNNEE , 'CNNEE' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNXSMR , 'CNXSMR' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNADD , 'CNADD' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNLOSS , 'CNLOSS' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNBURN , 'CNBURN' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,PARABS , 'PARABS' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,PARINC , 'PARINC' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SCSAT , 'SCSAT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SCUNS , 'SCUNS' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,BTRANT , 'BTRANT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,SIF , 'SIF' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNCO2 , 'CNCO2' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNFIRE_CNT , 'CNFIRE_CNT' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNSOM_CLOSS , 'CNSOM_CLOSS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNNDEPLOY , 'CNNDEPLOY' , RC=STATUS); VERIFY_(STATUS) @@ -5526,18 +5549,21 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,CNFUELC , 'CNFUELC' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNTOTLITC , 'CNTOTLITC' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,CNCWDC , 'CNCWDC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNROOT , 'CNROOT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,CNFSEL, 'CNFSEL' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU001,'RMELTDU001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU002,'RMELTDU002', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU003,'RMELTDU003', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU004,'RMELTDU004', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTDU005,'RMELTDU005', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTBC001,'RMELTBC001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTBC002,'RMELTBC002', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTOC001,'RMELTOC001', RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT,RMELTOC002,'RMELTOC002', RC=STATUS); VERIFY_(STATUS) - IF (RUN_IRRIG /= 0) call MAPL_GetPointer(EXPORT,IRRIGRATE ,'IRRIGRATE' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNROOT , 'CNROOT' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,CNFSEL , 'CNFSEL' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU001 ,'RMELTDU001' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU002 ,'RMELTDU002' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU003 ,'RMELTDU003' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU004 ,'RMELTDU004' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTDU005 ,'RMELTDU005' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTBC001 ,'RMELTBC001' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTBC002 ,'RMELTBC002' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTOC001 ,'RMELTOC001' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,RMELTOC002 ,'RMELTOC002' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WATERTABLED ,'WATERTABLED' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FSWCHANGE ,'FSWCHANGE' , RC=STATUS); VERIFY_(STATUS) + + IF (RUN_IRRIG /= 0) call MAPL_GetPointer(EXPORT,IRRIGRATE ,'IRRIGRATE' , RC=STATUS); VERIFY_(STATUS) NTILES = size(PS) @@ -5865,6 +5891,7 @@ subroutine Driver ( RC ) allocate(fveg2 (NTILES)) allocate(FICE1 (NTILES)) allocate(SLDTOT (NTILES)) + allocate(FSW_CHANGE(NTILES)) allocate(SHSBT (NTILES,NUM_SUBTILES)) allocate(DSHSBT (NTILES,NUM_SUBTILES)) @@ -6476,10 +6503,12 @@ subroutine Driver ( RC ) ! gkw: obtain catchment area fractions and soil moisture ! ------------------------------------------------------ -call catch_calc_soil_moist( ntiles, veg1, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & - srfexc, rzexc, catdef, car1, car2, car4, sfmc, rzmc, prmc, & - SWSRF1OUT=SWSRF1, SWSRF2OUT=SWSRF2, SWSRF4OUT=SWSRF4 ) + call catch_calc_soil_moist( ntiles, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & + ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, bf1, bf2, & + srfexc, rzexc, catdef, car1, car2, car4, sfmc, rzmc, prmc, & + SWSRF1OUT=SWSRF1, SWSRF2OUT=SWSRF2, SWSRF4OUT=SWSRF4 ) + + ! obtain saturated canopy resistance following Farquhar, CLM4 implementation @@ -6621,9 +6650,11 @@ subroutine Driver ( RC ) ! soil temperatures ! ----------------- - zbar = -sqrt(1.e-20+catdef(n)/bf1(n))+bf2(n) + + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTMP) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) HT(:)=GHTCNT(:,N) - CALL GNDTMP(poros(n),zbar,ht,frice,tp,soilice) + CALL GNDTMP(poros(n),-1.*zbar,ht,frice,tp,soilice) ! note minus sign for zbar ! At the CatchCNGridComp level, tp1, tp2, .., tp6 are export variables in units of Kelvin, ! - rreichle & borescan, 6 Nov 2020 @@ -6642,7 +6673,6 @@ subroutine Driver ( RC ) ! baseflow ! -------- - zbar = sqrt(1.e-20+catdef(n)/bf1(n))-bf2(n) bflow(n) = (1.-frice)*1000.* & cond(n)*exp(-(bf3(n)-ashift)-gnu(n)*zbar)/gnu(n) IF(catdef(n) >= cdcr1(n)) bflow(n) = 0. @@ -7540,9 +7570,9 @@ subroutine Driver ( RC ) IF ((RUN_IRRIG /= 0).AND.(ntiles >0)) THEN - CALL CATCH_CALC_SOIL_MOIST ( & - NTILES,VEG1,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & + CALL CATCH_CALC_SOIL_MOIST ( & + NTILES,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, CAR1, CAR2, CAR4, sfmc, rzmc, prmc) call irrigation_rate (IRRIG_METHOD, & @@ -7802,7 +7832,7 @@ subroutine Driver ( RC ) TSURF ,& SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 ,& LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW ,& - TCSORIG1, TPSN1IN1, TPSN1OUT1 ,& + TCSORIG1, TPSN1IN1, TPSN1OUT1, FSW_CHANGE ,& TC1_0=TC1_0, TC2_0=TC2_0, TC4_0=TC4_0 ,& QA1_0=QA1_0, QA2_0=QA2_0, QA4_0=QA4_0 ,& RCONSTIT=RCONSTIT, RMELT=RMELT, TOTDEPOS=TOTDEPOS, LHACC=LHACC) @@ -8031,6 +8061,10 @@ subroutine Driver ( RC ) if(associated(RMELTBC002)) RMELTBC002 = RMELT(:,7) if(associated(RMELTOC001)) RMELTOC001 = RMELT(:,8) if(associated(RMELTOC002)) RMELTOC002 = RMELT(:,9) + if(associated(FSWCHANGE)) FSWCHANGE = FSW_CHANGE + if(associated(WATERTABLED)) then + WATERTABLED = catch_calc_watertabled( BF1, BF2, CDCR2, POROS, WPWET, CATDEF ) + endif if(associated(TPSN1OUT)) then where(WESNN(1,:)>0.) @@ -8215,7 +8249,8 @@ subroutine Driver ( RC ) deallocate(TOTDEPOS ) deallocate(RMELT ) deallocate(FICE1 ) - deallocate(SLDTOT ) + deallocate(SLDTOT ) + deallocate(FSW_CHANGE) deallocate( btran ) deallocate( wgt ) deallocate( bt1 ) @@ -8776,11 +8811,11 @@ subroutine RUN0(gc, import, export, clock, rc) rzexccp = rzexc call catch_calc_soil_moist( & ! intent(in) - ntiles, nint(veg1), dzsf, vgwmax, cdcr1, cdcr2, & + ntiles, dzsf, vgwmax, cdcr1, cdcr2, & psis, bee, poros, wpwet, & ars1, ars2, ars3, & ara1, ara2, ara3, ara4, & - arw1, arw2, arw3, arw4, & + arw1, arw2, arw3, arw4, bf1, bf2, & ! intent(inout) ! from process_cat srfexccp, rzexccp, catdefcp, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchcn_iau.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchcn_iau.F90 index 89481e2bd..5587dc238 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchcn_iau.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchcn_iau.F90 @@ -23,8 +23,9 @@ module catchcn_iau ! *********************************************************************** subroutine apply_catchcn_iau( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1, bf2, & TG1_INC, TG2_INC, TG4_INC, & TC1_INC, TC2_INC, TC4_INC, QC1_INC, QC2_INC, QC4_INC, & CAPAC_INC, CATDEF_INC, RZEXC_INC, SRFEXC_INC, & @@ -39,12 +40,12 @@ subroutine apply_catchcn_iau( NTILES, & ! CATCHMENT MODEL PARAMETERS - integer, dimension( NTILES), intent(in) :: VEG real, dimension( NTILES), intent(in) :: DZSF, VGWMAX, CDCR1, CDCR2 real, dimension( NTILES), intent(in) :: PSIS, BEE, POROS, WPWET real, dimension( NTILES), intent(in) :: ARS1, ARS2, ARS3 real, dimension( NTILES), intent(in) :: ARA1, ARA2, ARA3, ARA4 real, dimension( NTILES), intent(in) :: ARW1, ARW2, ARW3, ARW4 + real, dimension( NTILES), intent(in) :: bf1, bf2 ! CATCHMENT-CN MODEL PROGNOSTIC INCREMENTS @@ -95,8 +96,9 @@ subroutine apply_catchcn_iau( NTILES, & ! make sure that updated prognostics are OK call check_catchcn_progn( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1, bf2, & TG1, TG2, TG4, TC1, TC2, TC4, QC1, QC2, QC4, & CAPAC, CATDEF, RZEXC, SRFEXC, & GHTCNT, WESNN, HTSNNN, SNDZN ) @@ -106,8 +108,9 @@ end subroutine apply_catchcn_iau ! *********************************************************************** subroutine check_catchcn_progn( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1, bf2, & TG1, TG2, TG4, TC1, TC2, TC4, QC1, QC2, QC4, & CAPAC, CATDEF, RZEXC, SRFEXC, & GHTCNT, WESNN, HTSNNN, SNDZN ) @@ -132,12 +135,12 @@ subroutine check_catchcn_progn( NTILES, & ! CATCHMENT MODEL PARAMETERS - integer, dimension( NTILES), intent(in) :: VEG real, dimension( NTILES), intent(in) :: DZSF, VGWMAX, CDCR1, CDCR2 real, dimension( NTILES), intent(in) :: PSIS, BEE, POROS, WPWET real, dimension( NTILES), intent(in) :: ARS1, ARS2, ARS3 real, dimension( NTILES), intent(in) :: ARA1, ARA2, ARA3, ARA4 real, dimension( NTILES), intent(in) :: ARW1, ARW2, ARW3, ARW4 + real, dimension( NTILES), intent(in) :: bf1, bf2 ! CATCHMENT-CN MODEL PROGNOSTICS @@ -219,9 +222,9 @@ subroutine check_catchcn_progn( NTILES, & ! lower bound on catdef, - reichle, 3 Apr 2012 call catch_calc_soil_moist( & - NTILES,veg,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + NTILES,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & + ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, & ar1, ar2, ar4 ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 index 2b218d726..fc7ac11ba 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 @@ -90,20 +90,21 @@ MODULE CATCHMENT_CN_MODEL N_sm => CATCH_N_ZONES, & SATCAPFR => CATCH_SATCAPFR, & PHIGT => CATCH_PHIGT, & - DZTC => CATCH_DZTC, & + DZTSURF => CATCH_DZTSURF, & DZGT => CATCH_DZGT, & - FSN => CATCH_FSN + FSN => CATCH_FSN, & + PEATCLSM_POROS_THRESHOLD, & + PEATCLSM_ZBARMAX_4_SYSOIL - USE SURFPARAMS, ONLY: CSOIL_2, RSWILT, & LAND_FIX, FLWALPHA - + USE lsm_routines, only : & INTERC, RZDRAIN, BASE, PARTITION, RZEQUIL,& - gndtp0, gndtmp, & - catch_calc_soil_moist, & + gndtp0, gndtmp, & + catch_calc_soil_moist, catch_calc_zbar, & catch_calc_wtotl, dampen_tc_oscillations, & - SRUNOFF + SRUNOFF USE SIBALB_COEFF, ONLY: coeffsib @@ -168,7 +169,7 @@ SUBROUTINE CATCHCN ( & EVACC, SHACC, TSURF, & SH_SNOW, AVET_SNOW, WAT_10CM, TOTWAT_SOIL, TOTICE_SOIL, & LH_SNOW, LWUP_SNOW, LWDOWN_SNOW, NETSW_SNOW, & - TCSORIG, TPSN1IN, TPSN1OUT, & + TCSORIG, TPSN1IN, TPSN1OUT,FSW_CHANGE, & TC1_0, TC2_0, TC4_0, QA1_0, QA2_0, QA4_0, EACC_0, & RCONSTIT, RMELT, TOTDEPOS, LHACC & ) @@ -230,11 +231,12 @@ SUBROUTINE CATCHCN ( & HSNACC, EVACC, SHACC REAL, INTENT(OUT), DIMENSION(:) :: GHFLUXSNO, GHTSKIN - REAL, INTENT(OUT), DIMENSION(:) :: SH_SNOW, AVET_SNOW, & + REAL, INTENT(OUT), DIMENSION(:) :: SH_SNOW, AVET_SNOW, & WAT_10CM, TOTWAT_SOIL, TOTICE_SOIL - REAL, INTENT(OUT), DIMENSION(:) :: LH_SNOW, LWUP_SNOW, & + REAL, INTENT(OUT), DIMENSION(:) :: LH_SNOW, LWUP_SNOW, & LWDOWN_SNOW, NETSW_SNOW - REAL, INTENT(OUT), DIMENSION(:) :: TCSORIG, TPSN1IN, TPSN1OUT + REAL, INTENT(OUT), DIMENSION(:) :: TCSORIG, TPSN1IN, TPSN1OUT, & + FSW_CHANGE REAL, INTENT(OUT), DIMENSION(:), OPTIONAL :: LHACC @@ -253,7 +255,7 @@ SUBROUTINE CATCHCN ( & RC, SATCAP, SNWFRC, POTFRC, ESNFRC, EVSNOW, SHFLUXS, HLWUPS, & HFTDS1, HFTDS2, HFTDS4, DHFT1, DHFT2, DHFT4, TPSNB, & QSATTC, DQSDTC, SWSRF1, SWSRF2, SWSRF4, AR4, & - FCAN, THRUL, THRUC, RZEQOL, frice, srfmx, & + FCAN, THRUL_VOL, THRUC_VOL, RZEQOL, frice, srfmx, & srfmn, RCST1, RCST2, EVAPFR, RDCX, EVAP1, EVAP2, & EVAP4, SHFLUX1, SHFLUX2, SHFLUX4, HLWUP1, HLWUP2, HLWUP4, & GHFLUX1, GHFLUX2, GHFLUX4, RZI, TC1SF, TC2SF, TC4SF, ar1old, & @@ -280,7 +282,7 @@ SUBROUTINE CATCHCN ( & REAL, DIMENSION(N_SM) :: T1, AREA, tkgnd, fhgnd - REAL :: TG1SN, TG2SN, TG4SN, DTG1SN,DTG2SN,DTG4SN, ZBAR, THETAF, & + REAL :: TG1SN, TG2SN, TG4SN, DTG1SN,DTG2SN,DTG4SN, ZBAR, THETAF, & XFICE, FH21, FH21W, FH21I, FH21D, DFH21W, DFH21I, DFH21D, & EVSN, SHFLS, HUPS, HCORR, SWNET0, HLWDWN0, TMPSNW, HLWTC, & DHLWTC, HSTURB, DHSDEA, DHSDTC, ESATTC, ETURB, DEDEA, DEDTC, & @@ -622,7 +624,7 @@ SUBROUTINE CATCHCN ( & !**** DETERMINE INITIAL VALUE OF RZEQ: CALL RZEQUIL ( & - NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET, & + NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET,POROS, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4, & RZEQOL & @@ -638,6 +640,7 @@ SUBROUTINE CATCHCN ( & CALL PARTITION ( & NCH,DTSTEP,DZSF,RZEXC, RZEQOL,VGWMAX,CDCR1,CDCR2, & PSIS,BEE,poros,WPWET, & + bf1, bf2, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4,BUG, & SRFEXC,CATDEF,RUNSRF, & @@ -685,16 +688,17 @@ SUBROUTINE CATCHCN ( & else phi=PHIGT end if - ZBAR=-SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTP0) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) THETAF=.5 DO LAYER=1,6 HT(LAYER)=GHTCNT(LAYER,N) ENDDO - CALL GNDTP0( & - T1,phi,ZBAR,THETAF, & - HT, & - fh21w,fH21i,fh21d,dfh21w,dfh21i,dfh21D,tp & + CALL GNDTP0( & + T1,phi,-1.*ZBAR,THETAF, & ! note minus sign for zbar + HT, & + fh21w,fH21i,fh21d,dfh21w,dfh21i,dfh21D,tp & ) HFTDS1(N)=-FH21W @@ -729,9 +733,7 @@ SUBROUTINE CATCHCN ( & ENDDO - !**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !**** 1. SATURATED FRACTION DO N=1,NCH @@ -750,9 +752,7 @@ SUBROUTINE CATCHCN ( & SHFLUX1, HLWUP1, GHFLUX1 & ) - !**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !**** 2. SUBSATURATED BUT UNSTRESSED FRACTION CALL RSURFP2 ( NCH, RSSAT, SWSRF2, QSAT2, QA2, PSUR, WPWET, RSURF ) @@ -764,20 +764,19 @@ SUBROUTINE CATCHCN ( & SWNETF, HLWDWN, ALW2, BLW2, & QM, CSOIL, CCANOP, PSUR, & HFTDS2, DHFT2, RD, RSURF, POTFRC, & - TG2SF, TC2SF, QA2, & + TG2SF, TC2SF, QA2, & EVAP2, EVROOT2, EVSURF2, EVINT2, & SHFLUX2, HLWUP2, GHFLUX2 & ) -!**** - !**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !**** 3. WILTING FRACTION + ! MB: For ET calculation, AR4 surface wetness is set to WPWET + WHERE (POROS > PEATCLSM_POROS_THRESHOLD) SWSRF4 = WPWET ! PEATCLSM + CALL RSURFP2 ( NCH, RSSAT, SWSRF4, QSAT4, QA4, PSUR, WPWET, RSURF ) - CALL FLUXES ( & NCH, FVEG, DTSTEP, QSAT4, DQS4, & ETURB4, DEDQA4X, DEDTC4X, HSTURB4, DHSDQA4X, DHSDTC4X, & @@ -785,12 +784,11 @@ SUBROUTINE CATCHCN ( & SWNETF, HLWDWN, ALW4, BLW4, & QM, CSOIL, CCANOP, PSUR, & HFTDS4, DHFT4, RD, RSURF, POTFRC, & - TG4SF, TC4SF, QA4, & + TG4SF, TC4SF, QA4, & EVAP4, EVROOT4, EVSURF4, EVINT4, & SHFLUX4, HLWUP4, GHFLUX4 & ) - !**** -------------------------------------------------------- !**** B. SNOW-COVERED FRACTION. @@ -801,7 +799,12 @@ SUBROUTINE CATCHCN ( & T1(1) = TG1(N)-TF T1(2) = TG2(N)-TF T1(3) = TG4(N)-TF - AREA(1)= AR1(N) + ! MB: avoid division by zero (AR1=0) in PEATCLSM equations + IF(POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + AREA(1)= amax1(AR1(N),2.E-20) + ELSE + AREA(1) = AR1(N) + END IF AREA(2)= AR2(N) AREA(3)= AR4(N) pr = trainc(n)+trainl(n)+tsnow(n)+tice(n)+tfrzr(n) @@ -818,7 +821,7 @@ SUBROUTINE CATCHCN ( & tkgnd(2)=1.8 tkgnd(3)=1.8 raddn=hlwdwn(n)+swnets(n) - zc1=-(DZTC*0.5) + zc1=-(DZTSURF*0.5) hups=0.0 !**** 1. RUN SNOW MODEL: @@ -958,9 +961,9 @@ SUBROUTINE CATCHCN ( & TG4SN=TPSNB(N) ENDIF - TG1(N)=TG1SF(N)*(1-AREASC)+TG1SN*AREASC - TG2(N)=TG2SF(N)*(1-AREASC)+TG2SN*AREASC - TG4(N)=TG4SF(N)*(1-AREASC)+TG4SN*AREASC + TG1(N)=TG1SF(N)*(1.-AREASC)+TG1SN*AREASC + TG2(N)=TG2SF(N)*(1.-AREASC)+TG2SN*AREASC + TG4(N)=TG4SF(N)*(1.-AREASC)+TG4SN*AREASC @@ -972,9 +975,9 @@ SUBROUTINE CATCHCN ( & DTC2SN=TPSN1(N)-TC2(N) DTC4SN=TPSN1(N)-TC4(N) - TC1(N)=TC1SF(N)*(1-AREASC)+TPSN1(N)*AREASC - TC2(N)=TC2SF(N)*(1-AREASC)+TPSN1(N)*AREASC - TC4(N)=TC4SF(N)*(1-AREASC)+TPSN1(N)*AREASC + TC1(N)=TC1SF(N)*(1.-AREASC)+TPSN1(N)*AREASC + TC2(N)=TC2SF(N)*(1.-AREASC)+TPSN1(N)*AREASC + TC4(N)=TC4SF(N)*(1.-AREASC)+TPSN1(N)*AREASC ! TC1(N)=TC1SF(N)*(1-AREASC)+TC1_ORIG(N)*AREASC ! TC2(N)=TC2SF(N)*(1-AREASC)+TC2_ORIG(N)*AREASC @@ -1050,17 +1053,19 @@ SUBROUTINE CATCHCN ( & else phi=PHIGT end if - ZBAR=-SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTMP) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) THETAF=.5 DO LAYER=1,6 HT(LAYER)=GHTCNT(LAYER,N) ENDDO FH21=-GHFLUX(N) - CALL GNDTMP( & - phi,zbar, & - ht, & - xfice,tp, soilice,DTS=dtstep,THETAF=thetaf,FH21=fh21) + CALL GNDTMP( & + phi, -1.*zbar, & ! note minus sign for zbar + ht, & + xfice, tp, soilice, & + DTS=dtstep, THETAF=thetaf, FH21=fh21) DO LAYER=1,6 GHTCNT(LAYER,N)=HT(LAYER) @@ -1088,7 +1093,8 @@ SUBROUTINE CATCHCN ( & NCH, DTSTEP, & EVAPFR, SATCAP, TG1, RA1, RC, & AR1,AR2,AR4,CDCR1, ESATFR, & - RZEQOL,SRFMN,WPWET,VGWMAX, & + RZEQOL,SRFMN,WPWET,VGWMAX, POROS, & + BF1, BF2, ARS1, ARS2, ARS3, & CAPAC, RZEXC, CATDEF, SRFEXC, & ESOI, EVEG, EINT, & ECORR & @@ -1114,7 +1120,8 @@ SUBROUTINE CATCHCN ( & CALL RZDRAIN ( & NCH,DTSTEP,VGWMAX,SATCAP,RZEQOL,AR1,WPWET, & - tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & + tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros, & + BF1, BF2, ARS1, ARS2, ARS3, BUG, & CAPAC,RZEXC,SRFEXC,CATDEF,RUNSRF & ) @@ -1127,9 +1134,9 @@ SUBROUTINE CATCHCN ( & !**** COMPUTE BASEFLOW FROM TOPMODEL EQUATIONS CALL BASE ( & - NCH, DTSTEP,BF1, BF2, BF3, CDCR1, FRICE, COND, GNU, & - CATDEF, & - BFLOW & + NCH, DTSTEP,BF1, BF2, BF3, CDCR1, FRICE, COND, GNU,AR1, POROS,& + ARS1, ARS2, ARS3, & + CATDEF, BFLOW & ) ! --------------------------------------------------------------------- @@ -1144,7 +1151,7 @@ SUBROUTINE CATCHCN ( & NCH, DTSTEP, FWETC, FWETL, TRAINLX, TRAINCX, SMELT, & SATCAP, BUG, & CAPAC, & - THRUL, THRUC & + THRUL_VOL, THRUC_VOL & ) IF (BUG) THEN @@ -1153,9 +1160,11 @@ SUBROUTINE CATCHCN ( & !**** DETERMINE SURFACE RUNOFF AND INFILTRATION RATES: - CALL SRUNOFF ( NCH,DTSTEP,UFW4RO, FWETC, FWETL, & - AR1,ar2,ar4,THRUL, THRUC,frice,tp1,srfmx,BUG, & - SRFEXC,RUNSRF, & + CALL SRUNOFF ( NCH, DTSTEP, UFW4RO, FWETC, FWETL, & + AR1, AR2, AR4, THRUL_VOL, THRUC_VOL, & + FRICE, TP1, SRFMX, BUG, & + VGWMAX, RZEQOL, POROS, & + SRFEXC, RZEXC, RUNSRF, & QINFIL & ) @@ -1168,7 +1177,7 @@ SUBROUTINE CATCHCN ( & !**** RECOMPUTE RZEXC: CALL RZEQUIL ( & - NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET, & + NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET,POROS, & ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & RZEQ & ) @@ -1194,8 +1203,8 @@ SUBROUTINE CATCHCN ( & ! note revised interface - reichle, 3 Apr 2012 CALL CATCH_CALC_SOIL_MOIST ( & - nch,ityp1,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & + nch,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, & AR1, AR2, AR4, & sfmc, rzmc, prmc, & @@ -1276,6 +1285,12 @@ SUBROUTINE CATCHCN ( & WCHANGE(N) = (WTOT(N)-WTOT_ORIG(N))/DTSTEP ECHANGE(N) = (ENTOT(N)-ENTOT_ORIG(N))/DTSTEP + !FSW_CHANGE IS THE CHANGE IN THE FREE-STANDING WATER, RELEVANT FOR PEATLAND ONLY + FSW_CHANGE(N) = 0. + IF(POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + pr = trainc(n)+trainl(n)+tsnow(n)+tice(n)+tfrzr(n) + FSW_CHANGE(N) = PR - EVAP(N) - RUNOFF(N) - WCHANGE(N) + ENDIF ! Perform check on sum of AR1 and AR2, to avoid calculation of negative ! wilting fraction due to roundoff, outside of catchment: @@ -1902,7 +1917,7 @@ SUBROUTINE FLUXES ( & DTC=DTCDRY+WETFRC*(DTCWET-DTCDRY) DTG=DTGDRY+WETFRC*(DTGWET-DTGDRY) DEA=DEADRY+WETFRC*(DEAWET-DEADRY) - ENBAL1=(1-FVEG(N))*SWNET(N) + FVEG(N)*(HLWTC + DHLWTC*DTC) + & + ENBAL1=(1.-FVEG(N))*SWNET(N) + FVEG(N)*(HLWTC + DHLWTC*DTC) + & (1.-FVEG(N))*HLWDWN(N) - (HLWTG + DHLWTC*DTG) - & DHDTG*(TG(N)-TC(N)) - & RHOTERM*ALHE*(ESATTC(N)+DESDTC(N)*(TG(N)-TCOLD) & @@ -2149,7 +2164,8 @@ SUBROUTINE WUPDAT ( & NCH, DTSTEP, & EVAP, SATCAP, TC, RA, RC, & AR1,AR2,AR4,CDCR1, ESATFR, & - RZEQ,SRFMN,WPWET,VGWMAX, & + RZEQ,SRFMN,WPWET,VGWMAX, POROS, & + BF1, BF2, ARS1, ARS2, ARS3, & CAPAC, RZEXC, CATDEF, SRFEXC, & EVROOT, EVSURF, EVINT, & ECORR & @@ -2163,7 +2179,8 @@ SUBROUTINE WUPDAT ( & INTEGER, INTENT(IN) :: NCH REAL, INTENT(IN) :: DTSTEP REAL, INTENT(IN), DIMENSION(NCH) :: EVAP, SATCAP, TC, RA, RC, AR1, & - AR2, AR4, CDCR1, ESATFR, RZEQ, SRFMN, WPWET, VGWMAX + AR2, AR4, CDCR1, ESATFR, RZEQ, SRFMN, WPWET, VGWMAX, & + POROS, BF1, BF2, ARS1, ARS2, ARS3 REAL, INTENT(INOUT), DIMENSION(NCH) :: CAPAC, RZEXC, CATDEF, & SRFEXC, EVROOT, EVSURF, EVINT @@ -2172,11 +2189,11 @@ SUBROUTINE WUPDAT ( & INTEGER N REAL EGRO,CNDSAT,CNDUNS,CNDV,CNDS,WILT,EGROMX,RZEMAX + REAL :: ZBAR1,SYSOIL,ET_CATDEF,AR1eq !**** !**** ----------------------------------------------------------------- DO 100 N = 1, NCH - !**** !**** PARTITION EVAP BETWEEN INTERCEPTION AND GROUND RESERVOIRS. !**** @@ -2217,17 +2234,30 @@ SUBROUTINE WUPDAT ( & !**** REMOVE MOISTURE FROM RESERVOIRS: !**** - IF (CATDEF(N) .LT. CDCR1(N)) THEN + IF (CATDEF(N) .LT. CDCR1(N)) THEN CAPAC(N) = AMAX1(0., CAPAC(N) - EVINT(N)*DTSTEP) RZEXC(N) = RZEXC(N) - EVROOT(N)*(1.-ESATFR(N))*DTSTEP SRFEXC(N) = SRFEXC(N) - EVSURF(N)*(1.-ESATFR(N))*DTSTEP - CATDEF(N) = CATDEF(N) + (EVSURF(N) + EVROOT(N))*ESATFR(N)*DTSTEP -! 05.12.98: FIRST ATTEMPT TO INCLUDE BEDROCK - ELSE + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN + CATDEF(N) = CATDEF(N) + (EVSURF(N) + EVROOT(N))*ESATFR(N)*DTSTEP + ! 05.12.98: FIRST ATTEMPT TO INCLUDE BEDROCK + ELSE + ! PEAT + ! MB: accounting for water ponding on AR1 + ! same approach as for RZFLW (see subroutine RZDRAIN for + ! comments) + ZBAR1 = catch_calc_zbar( BF1(N), BF2(N), CATDEF(N) ) + SYSOIL = (2.*bf1(N)*amin1(amax1(zbar1,0.),PEATCLSM_ZBARMAX_4_SYSOIL) + 2.*bf1(N)*bf2(N))/1000. + SYSOIL = amin1(SYSOIL,poros(N)) + ET_CATDEF = SYSOIL*(EVSURF(N) + EVROOT(N))*ESATFR(N)/(1.*AR1(N)+SYSOIL*(1.-AR1(N))) + AR1eq = (1.+ars1(N)*(catdef(N)))/(1.+ars2(N)*(catdef(N))+ars3(N)*(catdef(N))**2) + CATDEF(N) = CATDEF(N) + (1.-AR1eq)*ET_CATDEF + ENDIF + ELSE CAPAC(N) = AMAX1(0., CAPAC(N) - EVINT(N)*DTSTEP) RZEXC(N) = RZEXC(N) - EVROOT(N)*DTSTEP SRFEXC(N) = SRFEXC(N) - EVSURF(N)*DTSTEP - ENDIF + ENDIF !**** 100 CONTINUE @@ -2300,7 +2330,6 @@ END SUBROUTINE RSURFP2 !**** ----------------------------------------------------------------- !**** ///////////////////////////////////////////////////////////////// !**** ----------------------------------------------------------------- - ! ******************************************************************* subroutine catchcn_calc_tsurf( NTILES, tc1, tc2, tc4, wesnn, htsnn, & ar1, ar2, ar4, tsurf ) @@ -2389,11 +2418,11 @@ end subroutine catchcn_calc_tsurf_excl_snow ! ******************************************************************* - subroutine catchcn_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & - psis, bee, poros, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & - srfexc, rzexc, catdef, tc1, tc2, tc4, tg1, tg2, tg4, & - wesnn, htsnn, ghtcnt, & + subroutine catchcn_calc_etotl( NTILES, dzsf, vgwmax, cdcr1, cdcr2, & + psis, bee, poros, wpwet,bf1, bf2, & + ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & + srfexc, rzexc, catdef, tc1, tc2, tc4, tg1, tg2, tg4, & + wesnn, htsnn, ghtcnt, & etotl ) ! compute total energy stored in land tiles @@ -2407,10 +2436,9 @@ subroutine catchcn_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & integer, intent(in) :: NTILES - integer, dimension( NTILES), intent(in) :: vegcls real, dimension( NTILES), intent(in) :: dzsf real, dimension( NTILES), intent(in) :: vgwmax - real, dimension( NTILES), intent(in) :: cdcr1, cdcr2 + real, dimension( NTILES), intent(in) :: cdcr1, cdcr2, bf1, bf2 real, dimension( NTILES), intent(in) :: psis, bee, poros, wpwet real, dimension( NTILES), intent(in) :: ars1, ars2, ars3 real, dimension( NTILES), intent(in) :: ara1, ara2, ara3, ara4 @@ -2443,9 +2471,9 @@ subroutine catchcn_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & rzexc_tmp = rzexc ! rzexc is "inout" in catch_calc_soil_moist() catdef_tmp = catdef ! catdef is "inout" in catch_calc_soil_moist() - call catch_calc_soil_moist( & - NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & + call catch_calc_soil_moist( & + NTILES, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & + ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4,bf1, bf2,& srfexc_tmp, rzexc_tmp, catdef_tmp, ar1, ar2, ar4 ) ! compute snow-free tsurf 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 53837a54d..ebb728495 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 @@ -45,7 +45,7 @@ module GEOS_CatchGridCompMod SNWALB_NIRMAX => CATCH_SNWALB_NIRMAX, & SLOPE => CATCH_SNWALB_SLOPE - USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist + USE lsm_routines, ONLY : sibalb, catch_calc_soil_moist, catch_calc_watertabled !#for_ldas_coupling use catch_incr @@ -2678,6 +2678,24 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'depth_to_water_table_from_surface',& + UNITS = 'm' ,& + SHORT_NAME = 'WATERTABLED' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'change_in_free_surface_water_reservoir_on_peat',& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'FSWCHANGE' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + !EOS call MAPL_TimerAdd(GC, name="INITIALIZE",RC=STATUS) @@ -3909,6 +3927,8 @@ subroutine Driver ( RC ) real, pointer, dimension(:) :: RMELTBC002 real, pointer, dimension(:) :: RMELTOC001 real, pointer, dimension(:) :: RMELTOC002 + real, pointer, dimension(:) :: WATERTABLED + real, pointer, dimension(:) :: FSWCHANGE ! -------------------------------------------------------------------------- ! Local pointers for tile variables @@ -3937,7 +3957,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: ghflxsno, ghflxtskin real,pointer,dimension(:) :: SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 real,pointer,dimension(:) :: LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW - real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1 + real,pointer,dimension(:) :: TCSORIG1, TPSN1IN1, TPSN1OUT1, FSW_CHANGE real,pointer,dimension(:) :: WCHANGE, ECHANGE, HSNACC, EVACC, SHACC real,pointer,dimension(:) :: SNOVR, SNOVF, SNONR, SNONF real,pointer,dimension(:) :: VSUVR, VSUVF @@ -4448,6 +4468,8 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,RMELTBC002,'RMELTBC002', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,RMELTOC001,'RMELTOC001', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,RMELTOC002,'RMELTOC002', RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,WATERTABLED,'WATERTABLED',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT,FSWCHANGE, 'FSWCHANGE', RC=STATUS); VERIFY_(STATUS) NTILES = size(PS) @@ -4516,6 +4538,7 @@ subroutine Driver ( RC ) allocate(SUMEV (NTILES)) allocate(FICE1 (NTILES)) allocate(SLDTOT (NTILES)) ! total solid precip + allocate(FSW_CHANGE(NTILES)) allocate(SHSBT (NTILES,NUM_SUBTILES)) allocate(DSHSBT (NTILES,NUM_SUBTILES)) @@ -5342,8 +5365,8 @@ subroutine Driver ( RC ) SNDZN_INCR (3,:) = SNDZN3_INCR call apply_catch_incr(NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & - ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, bf1, bf2, & TCFSAT_INCR, TCFTRN_INCR, TCFWLT_INCR, QCFSAT_INCR, QCFTRN_INCR, QCFWLT_INCR, & CAPAC_INCR, CATDEF_INCR, RZEXC_INCR, SRFEXC_INCR, & GHTCNT_INCR, WESNN_INCR, HTSNNN_INCR, SNDZN_INCR, & @@ -5436,7 +5459,7 @@ subroutine Driver ( RC ) ENTOT,WTOT, WCHANGE, ECHANGE, HSNACC, EVACC, SHACC ,& SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1 ,& LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW ,& - TCSORIG1, TPSN1IN1, TPSN1OUT1 ,& + TCSORIG1, TPSN1IN1, TPSN1OUT1,FSW_CHANGE ,& lonbeg,lonend,latbeg,latend ,& TC1_0=TC1_0, TC2_0=TC2_0, TC4_0=TC4_0 ,& QA1_0=QA1_0, QA2_0=QA2_0, QA4_0=QA4_0 ,& @@ -5621,6 +5644,10 @@ subroutine Driver ( RC ) if(associated(RMELTBC002)) RMELTBC002 = RMELT(:,7) if(associated(RMELTOC001)) RMELTOC001 = RMELT(:,8) if(associated(RMELTOC002)) RMELTOC002 = RMELT(:,9) + if(associated(FSWCHANGE )) FSWCHANGE = FSW_CHANGE + if(associated(WATERTABLED )) then + WATERTABLED = catch_calc_watertabled( BF1, BF2, CDCR2, POROS, WPWET, CATDEF ) + endif if(associated(TPSN1)) then where(WESNN(1,:)>0.) @@ -5784,6 +5811,7 @@ subroutine Driver ( RC ) deallocate(RMELT ) deallocate(FICE1 ) deallocate(SLDTOT ) + deallocate(FSW_CHANGE) RETURN_(ESMF_SUCCESS) @@ -5861,6 +5889,8 @@ subroutine RUN0(gc, import, export, clock, rc) real, pointer :: arw2(:)=>null() real, pointer :: arw3(:)=>null() real, pointer :: arw4(:)=>null() + real, pointer :: bf1(:)=>null() + real, pointer :: bf2(:)=>null() !! Miscellaneous integer :: ntiles @@ -5954,6 +5984,10 @@ subroutine RUN0(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(INTERNAL, arw4, 'ARW4', rc=status) VERIFY_(status) + call MAPL_GetPointer(INTERNAL, bf1, 'BF1', rc=status) + VERIFY_(status) + call MAPL_GetPointer(INTERNAL, bf2, 'BF2', rc=status) + VERIFY_(status) call MAPL_GetPointer(INTERNAL, srfexc, 'SRFEXC', rc=status) VERIFY_(status) call MAPL_GetPointer(INTERNAL, rzexc, 'RZEXC', rc=status) @@ -6007,11 +6041,11 @@ subroutine RUN0(gc, import, export, clock, rc) rzexccp = rzexc call catch_calc_soil_moist( & ! intent(in) - ntiles, nint(ity), dzsf, vgwmax, cdcr1, cdcr2, & + ntiles, dzsf, vgwmax, cdcr1, cdcr2, & psis, bee, poros, wpwet, & ars1, ars2, ars3, & ara1, ara2, ara3, ara4, & - arw1, arw2, arw3, arw4, & + arw1, arw2, arw3, arw4,bf1, bf2, & ! intent(inout) ! from process_cat srfexccp, rzexccp, catdefcp, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catch_incr.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catch_incr.F90 index 3649b719a..e53013cf7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catch_incr.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catch_incr.F90 @@ -21,9 +21,10 @@ module catch_incr ! *********************************************************************** - subroutine apply_catch_incr( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + subroutine apply_catch_incr( NTILES, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1, bf2, & TC1_INC, TC2_INC, TC4_INC, QC1_INC, QC2_INC, QC4_INC, & CAPAC_INC, CATDEF_INC, RZEXC_INC, SRFEXC_INC, & GHTCNT_INC, WESNN_INC, HTSNNN_INC, SNDZN_INC, & @@ -37,12 +38,12 @@ subroutine apply_catch_incr( NTILES, & ! CATCHMENT MODEL PARAMETERS - integer, dimension( NTILES), intent(in) :: VEG real, dimension( NTILES), intent(in) :: DZSF, VGWMAX, CDCR1, CDCR2 real, dimension( NTILES), intent(in) :: PSIS, BEE, POROS, WPWET real, dimension( NTILES), intent(in) :: ARS1, ARS2, ARS3 real, dimension( NTILES), intent(in) :: ARA1, ARA2, ARA3, ARA4 real, dimension( NTILES), intent(in) :: ARW1, ARW2, ARW3, ARW4 + real, dimension( NTILES), intent(in) :: bf1, bf2 ! CATCHMENT MODEL PROGNOSTIC INCREMENTS @@ -87,8 +88,9 @@ subroutine apply_catch_incr( NTILES, & ! make sure that updated prognostics are OK call check_catch_progn( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & - ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1,bf2, & TC1, TC2, TC4, QC1, QC2, QC4, & CAPAC, CATDEF, RZEXC, SRFEXC, & GHTCNT, WESNN, HTSNNN, SNDZN ) @@ -98,8 +100,9 @@ end subroutine apply_catch_incr ! *********************************************************************** subroutine check_catch_progn( NTILES, & - VEG, DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & + DZSF, VGWMAX, CDCR1, CDCR2, PSIS, BEE, POROS, WPWET, & ARS1, ARS2, ARS3, ARA1, ARA2, ARA3, ARA4, ARW1, ARW2, ARW3, ARW4, & + bf1,bf2, & TC1, TC2, TC4, QC1, QC2, QC4, & CAPAC, CATDEF, RZEXC, SRFEXC, & GHTCNT, WESNN, HTSNNN, SNDZN ) @@ -124,12 +127,12 @@ subroutine check_catch_progn( NTILES, & ! CATCHMENT MODEL PARAMETERS - integer, dimension( NTILES), intent(in) :: VEG real, dimension( NTILES), intent(in) :: DZSF, VGWMAX, CDCR1, CDCR2 real, dimension( NTILES), intent(in) :: PSIS, BEE, POROS, WPWET real, dimension( NTILES), intent(in) :: ARS1, ARS2, ARS3 real, dimension( NTILES), intent(in) :: ARA1, ARA2, ARA3, ARA4 real, dimension( NTILES), intent(in) :: ARW1, ARW2, ARW3, ARW4 + real, dimension( NTILES), intent(in) :: bf1,bf2 ! CATCHMENT MODEL PROGNOSTICS @@ -210,9 +213,9 @@ subroutine check_catch_progn( NTILES, & ! lower bound on catdef, - reichle, 3 Apr 2012 call catch_calc_soil_moist( & - NTILES,veg,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + NTILES,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & + ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, & ar1, ar2, ar4 ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 index d5c2da131..3c743903f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 @@ -92,7 +92,9 @@ MODULE CATCHMENT_MODEL N_sm => CATCH_N_ZONES, & SATCAPFR => CATCH_SATCAPFR, & PHIGT => CATCH_PHIGT, & - DZTC => CATCH_DZTC + DZTSURF => CATCH_DZTSURF, & + PEATCLSM_POROS_THRESHOLD, & + PEATCLSM_ZBARMAX_4_SYSOIL USE SURFPARAMS, ONLY: & LAND_FIX, ASTRFR, STEXP, RSWILT, & @@ -101,10 +103,10 @@ MODULE CATCHMENT_MODEL USE lsm_routines, only : & INTERC, RZDRAIN, BASE, PARTITION, RZEQUIL,& gndtp0, gndtmp, & - catch_calc_soil_moist, & + catch_calc_soil_moist, catch_calc_zbar, & catch_calc_wtotl, dampen_tc_oscillations, & SRUNOFF - + USE SIBALB_COEFF, ONLY: coeffsib USE STIEGLITZSNOW, ONLY: & @@ -157,7 +159,8 @@ SUBROUTINE CATCHMENT ( & EVACC, SHACC, & SH_SNOW, AVET_SNOW, WAT_10CM, TOTWAT_SOIL, TOTICE_SOIL, & LH_SNOW, LWUP_SNOW, LWDOWN_SNOW, NETSW_SNOW, & - TCSORIG, TPSN1IN, TPSN1OUT,lonbeg,lonend,latbeg,latend, & + TCSORIG, TPSN1IN, TPSN1OUT,FSW_CHANGE , & + lonbeg,lonend,latbeg,latend, & TC1_0, TC2_0, TC4_0, QA1_0, QA2_0, QA4_0, EACC_0, & RCONSTIT, RMELT, TOTDEPOS, LHACC) @@ -218,11 +221,12 @@ SUBROUTINE CATCHMENT ( & HSNACC, EVACC, SHACC REAL, INTENT(OUT), DIMENSION(NCH) :: GHFLUXSNO, GHTSKIN - REAL, INTENT(OUT), DIMENSION(NCH) :: SH_SNOW, AVET_SNOW, & + REAL, INTENT(OUT), DIMENSION(NCH) :: SH_SNOW, AVET_SNOW, & WAT_10CM, TOTWAT_SOIL, TOTICE_SOIL - REAL, INTENT(OUT), DIMENSION(NCH) :: LH_SNOW, LWUP_SNOW, & + REAL, INTENT(OUT), DIMENSION(NCH) :: LH_SNOW, LWUP_SNOW, & LWDOWN_SNOW, NETSW_SNOW - REAL, INTENT(OUT), DIMENSION(NCH) :: TCSORIG, TPSN1IN, TPSN1OUT + REAL, INTENT(OUT), DIMENSION(NCH) :: TCSORIG, TPSN1IN, TPSN1OUT, & + FSW_CHANGE REAL, INTENT(OUT), DIMENSION(NCH), OPTIONAL :: LHACC @@ -242,7 +246,8 @@ SUBROUTINE CATCHMENT ( & RC, SATCAP, SNWFRC, POTFRC, ESNFRC, EVSNOW, SHFLUXS, HLWUPS, & HFTDS1, HFTDS2, HFTDS4, DHFT1, DHFT2, DHFT4, TPSNB, & QSATTC, DQSDTC, SWSRF1, SWSRF2, SWSRF4, AR4, RX11, RX21, RX12, & - RX14, RX24, RX22, EIRFRC, FCAN, THRUL, THRUC,RZEQOL, frice, srfmx, & + RX14, RX24, RX22, EIRFRC, FCAN, THRUL_VOL, THRUC_VOL, & + RZEQOL, frice, srfmx, & srfmn, RCST, EVAPFR, RCUN, PAR, PDIR, RDCX, EVAP1, EVAP2, & EVAP4, SHFLUX1, SHFLUX2, SHFLUX4, HLWUP1, HLWUP2, HLWUP4, & GHFLUX1, GHFLUX2, GHFLUX4, RZI, TC1SF, TC2SF, TC4SF, ar1old, & @@ -610,7 +615,7 @@ SUBROUTINE CATCHMENT ( & !**** DETERMINE INITIAL VALUE OF RZEQ: CALL RZEQUIL ( & - NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET, & + NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET,POROS, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4, & RZEQOL & @@ -630,6 +635,7 @@ SUBROUTINE CATCHMENT ( & CALL PARTITION ( & NCH,DTSTEP,DZSF,RZEXC, RZEQOL,VGWMAX,CDCR1,CDCR2, & PSIS,BEE,poros,WPWET, & + bf1, bf2, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4,BUG, & SRFEXC,CATDEF,RUNSRF, & @@ -683,9 +689,12 @@ SUBROUTINE CATCHMENT ( & end if if (LAND_FIX) then - ZBAR =-SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) ! zbar bug fix, - reichle, 16 Nov 2015 - else - ZBAR=-SQRT(1.e-20+catdef(n)/bf1(n))-bf2(n) + ! zbar bug fix, - reichle, 16 Nov 2015 + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTP0) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) + else + ! zbar minus sign applied in call to GNDTP0 + ZBAR = SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) ! old bug is wrong sign for bf2 here end if THETAF=.5 @@ -693,10 +702,10 @@ SUBROUTINE CATCHMENT ( & HT(LAYER)=GHTCNT(LAYER,N) ENDDO - CALL GNDTP0( & - T1,phi,ZBAR,THETAF, & - HT, & - fh21w,fH21i,fh21d,dfh21w,dfh21i,dfh21D,tp & + CALL GNDTP0( & + T1,phi,-1.*ZBAR,THETAF, & ! note minus sign for zbar + HT, & + fh21w,fH21i,fh21d,dfh21w,dfh21i,dfh21D,tp & ) HFTDS1(N)=-FH21W @@ -770,7 +779,7 @@ SUBROUTINE CATCHMENT ( & !**** 3. WILTING FRACTION !CC print*,'energy4' CALL ENERGY4 ( & - NCH, DTSTEP, ITYP, UM, RCST, & + NCH, DTSTEP, ITYP,POROS, UM, RCST, & ETURB4, DEDQA4X, DEDTC4X, HSTURB4, DHSDQA4X, DHSDTC4X, & QM, RA4, SWNETF, HLWDWN, PSUR, & RDCX, HFTDS4, DHFT4, QSAT4, DQS4, ALW4, BLW4, & @@ -828,7 +837,12 @@ SUBROUTINE CATCHMENT ( & T1(1) = TC1(N)-TF T1(2) = TC2(N)-TF T1(3) = TC4(N)-TF - AREA(1)= AR1(N) + ! MB: avoid division by zero (AR1=0) in PEATCLSM equations + IF(POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + AREA(1)= amax1(AR1(N),2.E-20) + ELSE + AREA(1) = AR1(N) + END IF AREA(2)= AR2(N) AREA(3)= AR4(N) pr = trainc(n)+trainl(n)+tsnow(n)+tice(n)+tfrzr(n) @@ -845,7 +859,7 @@ SUBROUTINE CATCHMENT ( & tkgnd(2)=1.8 tkgnd(3)=1.8 raddn=hlwdwn(n)+swnets(n) - zc1=-(DZTC*0.5) + zc1=-(DZTSURF*0.5) hups=0.0 !**** 1. RUN SNOW MODEL: @@ -990,9 +1004,9 @@ SUBROUTINE CATCHMENT ( & - TC1(N)=TC1SF(N)*(1-AREASC)+TC1SN*AREASC - TC2(N)=TC2SF(N)*(1-AREASC)+TC2SN*AREASC - TC4(N)=TC4SF(N)*(1-AREASC)+TC4SN*AREASC + TC1(N)=TC1SF(N)*(1.-AREASC)+TC1SN*AREASC + TC2(N)=TC2SF(N)*(1.-AREASC)+TC2SN*AREASC + TC4(N)=TC4SF(N)*(1.-AREASC)+TC4SN*AREASC EVSNOW(N)=EVSN esno(n)=evsnow(n)*asnow(n)*DTSTEP ! to have esno in mm/20min (03-17-99) @@ -1028,7 +1042,7 @@ SUBROUTINE CATCHMENT ( & GHFLUX(N)=(1.-ASNOW(N))* & (GHFLUX1(N)*AR1(N)+GHFLUX2(N)*AR2(N)+GHFLUX4(N)*AR4(N)) & +ASNOW(N)*GHFLUXS(N) - GHTSKIN(N)=(1.-ASNOW(N))* & + GHTSKIN(N)=(1.-ASNOW(N))* & (GHFLUX1(N)*AR1(N)+GHFLUX2(N)*AR2(N)+GHFLUX4(N)*AR4(N)) & -ASNOW(N)*ghfluxsno(N) ENDDO @@ -1049,9 +1063,12 @@ SUBROUTINE CATCHMENT ( & phi=PHIGT end if if (LAND_FIX) then - ZBAR =-SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) ! zbar bug fix, - reichle, 16 Nov 2015 + ! zbar bug fix, - reichle, 16 Nov 2015 + ! zbar function - reichle, 29 Jan 2022 (minus sign applied in call to GNDTMP) + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) else - ZBAR=-SQRT(1.e-20+catdef(n)/bf1(n))-bf2(n) + ! zbar minus sign applied in call to GNDTMP + ZBAR = SQRT(1.e-20+catdef(n)/bf1(n))+bf2(n) ! old bug is wrong sign for bf2 here end if THETAF=.5 DO LAYER=1,6 @@ -1059,10 +1076,11 @@ SUBROUTINE CATCHMENT ( & ENDDO FH21=-GHFLUX(N) - CALL GNDTMP( & - phi,zbar, & - ht, & - xfice,tp, soilice,DTS=dtstep,THETAF=thetaf,FH21=fh21) + CALL GNDTMP( & + phi, -1.*zbar, & ! note minus sign for zbar + ht, & + xfice, tp, soilice, & + DTS=dtstep, THETAF=thetaf, FH21=fh21) DO LAYER=1,6 GHTCNT(LAYER,N)=HT(LAYER) @@ -1107,7 +1125,8 @@ SUBROUTINE CATCHMENT ( & CALL WUPDAT ( & NCH, DTSTEP, EVAPFR, SATCAP, TC1, RA1, RC, & RX11,RX21,RX12,RX22,RX14,RX24, & - AR1,AR2,AR4,CDCR1,EIRFRC,RZEQOL,srfmn,WPWET,VGWMAX, & + AR1,AR2,AR4,CDCR1,EIRFRC,RZEQOL,srfmn,WPWET,VGWMAX,POROS, & + BF1, BF2, ARS1, ARS2, ARS3, & CAPAC, RZEXC, CATDEF, SRFEXC, & EINT, ESOI, EVEG & ) @@ -1122,7 +1141,8 @@ SUBROUTINE CATCHMENT ( & CALL RZDRAIN ( & NCH,DTSTEP,VGWMAX,SATCAP,RZEQOL,AR1,WPWET, & - tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & + tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros, & + BF1, BF2, ARS1, ARS2, ARS3, BUG, & CAPAC,RZEXC,SRFEXC,CATDEF,RUNSRF & ) @@ -1135,9 +1155,9 @@ SUBROUTINE CATCHMENT ( & !**** COMPUTE BASEFLOW FROM TOPMODEL EQUATIONS CALL BASE ( & - NCH, DTSTEP,BF1, BF2, BF3, CDCR1, FRICE, COND, GNU, & - CATDEF, & - BFLOW & + NCH, DTSTEP,BF1, BF2, BF3, CDCR1, FRICE, COND, GNU,AR1, POROS,& + ARS1, ARS2, ARS3, & + CATDEF, BFLOW & ) ! --------------------------------------------------------------------- @@ -1152,7 +1172,7 @@ SUBROUTINE CATCHMENT ( & NCH, DTSTEP, FWETC, FWETL, TRAINLX, TRAINCX, SMELT, & SATCAP, BUG, & CAPAC, & - THRUL, THRUC & + THRUL_VOL, THRUC_VOL & ) IF (BUG) THEN @@ -1161,9 +1181,11 @@ SUBROUTINE CATCHMENT ( & !**** DETERMINE SURFACE RUNOFF AND INFILTRATION RATES: - CALL SRUNOFF ( NCH,DTSTEP,UFW4RO, FWETC, FWETL, & - AR1,ar2,ar4,THRUL, THRUC,frice,tp1,srfmx,BUG, & - SRFEXC,RUNSRF, & + CALL SRUNOFF ( NCH, DTSTEP, UFW4RO, FWETC, FWETL, & + AR1, AR2, AR4, THRUL_VOL, THRUC_VOL, & + FRICE, TP1, SRFMX, BUG, & + VGWMAX, RZEQOL, POROS, & + SRFEXC, RZEXC, RUNSRF, & QINFIL & ) @@ -1176,7 +1198,7 @@ SUBROUTINE CATCHMENT ( & !**** RECOMPUTE RZEXC: CALL RZEQUIL ( & - NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET, & + NCH, CATDEF, VGWMAX,CDCR1,CDCR2,WPWET,POROS, & ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & RZEQ & ) @@ -1202,8 +1224,8 @@ SUBROUTINE CATCHMENT ( & ! note revised interface - reichle, 3 Apr 2012 CALL CATCH_CALC_SOIL_MOIST ( & - nch,ityp,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4, & + nch,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4,arw1,arw2,arw3,arw4,bf1,bf2, & srfexc,rzexc,catdef, & AR1, AR2, AR4, & sfmc, rzmc, prmc, & @@ -1284,6 +1306,12 @@ SUBROUTINE CATCHMENT ( & WCHANGE(N) = (WTOT(N)-WTOT_ORIG(N))/DTSTEP ECHANGE(N) = (ENTOT(N)-ENTOT_ORIG(N))/DTSTEP + !FSW_CHANGE IS THE CHANGE IN THE FREE-STANDING WATER, RELEVANT FOR PEATLAND ONLY + FSW_CHANGE(N) = 0. + IF(POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + pr = trainc(n)+trainl(n)+tsnow(n)+tice(n)+tfrzr(n) + FSW_CHANGE(N) = PR - EVAP(N) - RUNOFF(N) - WCHANGE(N) + ENDIF ! Perform check on sum of AR1 and AR2, to avoid calculation of negative ! wilting fraction due to roundoff, outside of catchment: @@ -1812,6 +1840,7 @@ SUBROUTINE energy1 ( & DEDEA(CHNO) = DEDQA(CHNO) * EPSILON / PSUR(CHNO) DHSDEA(CHNO) = DHSDQA(CHNO) * EPSILON / PSUR(CHNO) + 100 CONTINUE !**** @@ -2142,7 +2171,7 @@ END SUBROUTINE energy2 !**** ----------------------------------------------------------------- !**** SUBROUTINE energy4 ( & - NCH, DTSTEP, ITYP, UM, RCIN, & + NCH, DTSTEP, ITYP, POROS,UM, RCIN, & ETURB, DEDQA, DEDTC, HSTURB, DHSDQA, DHSDTC, & QM, RA, SWNET, HLWDWN, PSUR, & RDC, HFTDS, DHFTDS, & @@ -2161,7 +2190,7 @@ SUBROUTINE energy4 ( & REAL, INTENT(IN), DIMENSION(NCH) :: UM, RCIN, ETURB, HSTURB, QM, RA, & SWNET, HLWDWN, PSUR, RDC, HFTDS, DHFTDS, QSATTC, DQSDTC, & ALWRAD, BLWRAD, EMAXRT, CSOIL, SWSRF, POTFRC, WPWET, DEDQA, & - DEDTC, DHSDQA, DHSDTC + DEDTC, DHSDQA, DHSDTC, POROS LOGICAL, INTENT(IN) :: BUG REAL, INTENT(INOUT), DIMENSION(NCH) :: TC, QA @@ -2172,7 +2201,7 @@ SUBROUTINE energy4 ( & INTEGER ChNo, N REAL, DIMENSION(NCH) :: DEDEA, DHSDEA, EM, ESATTC, DESDTC, EA, RC, & - DRCDTC, DRCDEA + DRCDTC, DRCDEA, SWSRF4 REAL DELTC, DELEA !**** @@ -2200,6 +2229,15 @@ SUBROUTINE energy4 ( & DEDEA(CHNO) = DEDQA(CHNO) * EPSILON / PSUR(CHNO) DHSDEA(CHNO) = DHSDQA(CHNO) * EPSILON / PSUR(CHNO) + IF (POROS(CHNO) < PEATCLSM_POROS_THRESHOLD) THEN + ! mineral soil + SWSRF4(CHNO) = SWSRF(CHNO) + ELSE + ! PEAT + ! MB: For ET calculation, AR4 surface wetness is set to WPWET + SWSRF4(CHNO) = WPWET(CHNO) + ENDIF + 100 CONTINUE !**** @@ -2217,7 +2255,7 @@ SUBROUTINE energy4 ( & ENDDO CALL RSURFP2 ( & - NCH, UM, RDC, SWSRF, ESATTC, EA, WPWET, & + NCH, UM, RDC, SWSRF4, ESATTC, EA, WPWET, & RC, & RX1, RX2 & ) @@ -2566,9 +2604,10 @@ END SUBROUTINE TMPFAC !**** [ BEGIN WUPDAT ] !**** SUBROUTINE WUPDAT ( & - NCH, DTSTEP, EVAP, SATCAP, TC, RA, RC, & + NCH, DTSTEP, EVAP, SATCAP, TC, RA, RC, & RX11,RX21,RX12,RX22,RX14,RX24, AR1,AR2,AR4,CDCR1, & - EIRFRC,RZEQ,srfmn,WPWET,VGWMAX, & + EIRFRC,RZEQ,srfmn,WPWET,VGWMAX, POROS, & + BF1, BF2, ARS1, ARS2, ARS3, & CAPAC, RZEXC, CATDEF, SRFEXC, & EINT, ESOI, EVEG & ) @@ -2583,7 +2622,7 @@ SUBROUTINE WUPDAT ( & REAL, INTENT(IN) :: DTSTEP REAL, INTENT(IN), DIMENSION(NCH) :: EVAP, SATCAP, TC, RA, RC, RX11, & RX21, RX12, RX22, RX14, RX24, AR1, AR2, AR4, CDCR1, EIRFRC, & - RZEQ, srfmn, WPWET, VGWMAX + RZEQ, srfmn, WPWET, VGWMAX, POROS, BF1, BF2, ARS1, ARS2, ARS3 REAL, INTENT(INOUT), DIMENSION(NCH) :: CAPAC, CATDEF, RZEXC, SRFEXC @@ -2592,6 +2631,7 @@ SUBROUTINE WUPDAT ( & INTEGER CHNO REAL EGRO, CNDSAT, CNDUNS, ESATFR, cndv, cnds, WILT, egromx,rzemax + REAL :: ZBAR1,SYSOIL,ET_CATDEF,AR1eq !**** !**** ----------------------------------------------------------------- @@ -2672,11 +2712,25 @@ SUBROUTINE WUPDAT ( & !**** REMOVE MOISTURE FROM RESERVOIRS: !**** - IF (CATDEF(CHNO) .LT. CDCR1(CHNO)) THEN + IF (CATDEF(CHNO) .LT. CDCR1(CHNO)) THEN CAPAC(CHNO) = AMAX1(0., CAPAC(CHNO) - EINT(CHNO)) RZEXC(CHNO) = RZEXC(CHNO) - EVEG(CHNO)*(1.-ESATFR) SRFEXC(CHNO) = SRFEXC(CHNO) - ESOI(CHNO)*(1.-ESATFR) - CATDEF(CHNO) = CATDEF(CHNO) + (ESOI(CHNO) + EVEG(CHNO))*ESATFR + + IF (POROS(CHNO) < PEATCLSM_POROS_THRESHOLD) THEN + CATDEF(CHNO) = CATDEF(CHNO) + (ESOI(CHNO) + EVEG(CHNO))*ESATFR + ELSE + ! PEAT + ! MB: accounting for water ponding on AR1 + ! same approach as for RZFLW (see subroutine RZDRAIN for + ! comments) + ZBAR1 = catch_calc_zbar( BF1(CHNO), BF2(CHNO), CATDEF(CHNO) ) + SYSOIL = (2.*bf1(CHNO)*amin1(amax1(zbar1,0.),PEATCLSM_ZBARMAX_4_SYSOIL) + 2.*bf1(CHNO)*bf2(CHNO))/1000. + SYSOIL = amin1(SYSOIL,poros(CHNO)) + ET_CATDEF = SYSOIL*(ESOI(CHNO) + EVEG(CHNO))*ESATFR/(1.*AR1(CHNO)+SYSOIL*(1.-AR1(CHNO))) + AR1eq = (1.+ars1(chno)*(catdef(chno)))/(1.+ars2(chno)*(catdef(chno))+ars3(chno)*(catdef(chno))**2) + CATDEF(CHNO) = CATDEF(CHNO) + (1.-AR1eq)*ET_CATDEF + ENDIF ! 05.12.98: first attempt to include bedrock ELSE CAPAC(CHNO) = AMAX1(0., CAPAC(CHNO) - EINT(CHNO)) @@ -2947,7 +3001,7 @@ end subroutine catch_calc_tsurf_excl_snow ! ******************************************************************* subroutine catch_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & - psis, bee, poros, wpwet, & + psis, bee, poros, wpwet, bf1, bf2, & ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & srfexc, rzexc, catdef, tc1, tc2, tc4, wesnn, htsnn, ghtcnt, & etotl ) @@ -2966,7 +3020,7 @@ subroutine catch_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & integer, dimension( NTILES), intent(in) :: vegcls real, dimension( NTILES), intent(in) :: dzsf real, dimension( NTILES), intent(in) :: vgwmax - real, dimension( NTILES), intent(in) :: cdcr1, cdcr2 + real, dimension( NTILES), intent(in) :: cdcr1, cdcr2, bf1, bf2 real, dimension( NTILES), intent(in) :: psis, bee, poros, wpwet real, dimension( NTILES), intent(in) :: ars1, ars2, ars3 real, dimension( NTILES), intent(in) :: ara1, ara2, ara3, ara4 @@ -2998,9 +3052,9 @@ subroutine catch_calc_etotl( NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, & rzexc_tmp = rzexc ! rzexc is "inout" in catch_calc_soil_moist() catdef_tmp = catdef ! catdef is "inout" in catch_calc_soil_moist() - call catch_calc_soil_moist( & - NTILES, vegcls, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, & + call catch_calc_soil_moist( & + NTILES, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & + ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, arw2, arw3, arw4, bf1, bf2, & srfexc_tmp, rzexc_tmp, catdef_tmp, ar1, ar2, ar4 ) ! compute snow-free tsurf diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_constants.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_constants.f90 index e5c0c7b2f..92d07809f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_constants.f90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_constants.f90 @@ -70,9 +70,11 @@ module catch_constants ! --------------------------------------------------------------------------- ! - ! layer depth associated with snow-free land temperatures + ! layer depth associated with snow-free land surface soil temperatures ! - ! Note: DZTC = .05 is a hardwired setting of the depth of the bottom of + ! ================== + ! Note by Randy Koster & Rolf Reichle when CSOIL=200. was still used (~2018): + ! DZTC = .05 is a hardwired setting of the depth of the bottom of ! the surface soil layer. It should be made a parameter that is tied to ! the heat capacity CSOIL, which had been set to either CSOIL_1 or ! CSOIL_2 based on vegetation type. For now we leave @@ -87,8 +89,13 @@ module catch_constants ! are other impacts in wet climates regarding the effect of ! the depth of the water table on the thermal conductivity; these impacts ! are presumably very small. - - REAL, PARAMETER, PUBLIC :: CATCH_DZTC = 0.05 ! m layer depth for tc1, tc2, tc4 + ! ================== + ! + ! DZTSURF (formerly DZTC) is the layer depth associated w/ surface soil temperatures: + ! Catchment: tc1, tc2, tc4 (CSOIL also includes veg canopy) + ! CatchmentCN: tg1, tg2, tg4 (tc[X] are separate canopy temperatures) + + REAL, PARAMETER, PUBLIC :: CATCH_DZTSURF = 0.05 ! m layer depth for tc[X] or tg[X] ! --------------------------------------------------------------------------- ! @@ -127,6 +134,25 @@ module catch_constants REAL, PARAMETER, PUBLIC :: CATCH_SATCAPFR = 0.2 ! SATCAP = SATCAPFR * LAI + + ! peatCLSM implementation smahanam 3-16-2021 + ! + ! Use of peat-specific hydrology (PEATCLSM) is triggered by a porosity threshold. + ! Porosity of peat tiles depends on bcs version. + ! + ! bcs version | source of peat info | porosity + ! ----------------------------------------------------------------- + ! NLv3, NLv4 | HWSD | poros=0.80 + ! NLv5 | PEATMAP | poros=0.93 + ! + ! - reichle, 26 Jan 2022 + + REAL, PARAMETER, PUBLIC :: PEATCLSM_POROS_THRESHOLD = 0.90 ! [m3/m3] + + ! max zbar for specific yield calc in PEATCLSM + + REAL, PARAMETER, PUBLIC :: PEATCLSM_ZBARMAX_4_SYSOIL = 0.45 ! [m] + contains ! **************************************************************************************** @@ -142,39 +168,42 @@ subroutine echo_catch_constants(logunit) write (logunit,*) write (logunit,*) 'echo_catch_constants():' write (logunit,*) - write (logunit,*) 'CATCH_N_PFAFS = ', CATCH_N_PFAFS - write (logunit,*) 'CATCH_N_PFAFSROUTE = ', CATCH_N_PFAFSROUTE - write (logunit,*) 'CATCH_N_SNOW = ', CATCH_N_SNOW - write (logunit,*) 'CATCH_N_GT = ', CATCH_N_GT - write (logunit,*) 'CATCH_N_ZONES = ', CATCH_N_ZONES - write (logunit,*) 'CATCH_SNWALB_RHOFS = ', CATCH_SNWALB_RHOFS - write (logunit,*) 'CATCH_SNWALB_VISMAX = ', CATCH_SNWALB_VISMAX - write (logunit,*) 'CATCH_SNWALB_NIRMAX = ', CATCH_SNWALB_NIRMAX - write (logunit,*) 'CATCH_SNWALB_SLOPE = ', CATCH_SNWALB_SLOPE - write (logunit,*) 'CATCH_MAXSNDEPTH = ', CATCH_MAXSNDEPTH - write (logunit,*) 'CATCH_DZ1MAX = ', CATCH_DZ1MAX - write (logunit,*) 'CATCH_DZTC = ', CATCH_DZTC - write (logunit,*) 'CATCH_DZGT = ', CATCH_DZGT - write (logunit,*) 'CATCH_PHIGT = ', CATCH_PHIGT - write (logunit,*) 'CATCH_ALHMGT = ', CATCH_ALHMGT - write (logunit,*) 'CATCH_FSN = ', CATCH_FSN - write (logunit,*) 'CATCH_SHR = ', CATCH_SHR - write (logunit,*) 'CATCH_SCONST = ', CATCH_SCONST - write (logunit,*) 'CATCH_CSOIL_1 = ', CATCH_CSOIL_1 - write (logunit,*) 'CATCH_C_CANOP = ', CATCH_C_CANOP - write (logunit,*) 'CATCH_SATCAPFR = ', CATCH_SATCAPFR + write (logunit,*) 'CATCH_N_PFAFS = ', CATCH_N_PFAFS + write (logunit,*) 'CATCH_N_PFAFSROUTE = ', CATCH_N_PFAFSROUTE + write (logunit,*) 'CATCH_N_SNOW = ', CATCH_N_SNOW + write (logunit,*) 'CATCH_N_GT = ', CATCH_N_GT + write (logunit,*) 'CATCH_N_ZONES = ', CATCH_N_ZONES + write (logunit,*) 'CATCH_SNWALB_RHOFS = ', CATCH_SNWALB_RHOFS + write (logunit,*) 'CATCH_SNWALB_VISMAX = ', CATCH_SNWALB_VISMAX + write (logunit,*) 'CATCH_SNWALB_NIRMAX = ', CATCH_SNWALB_NIRMAX + write (logunit,*) 'CATCH_SNWALB_SLOPE = ', CATCH_SNWALB_SLOPE + write (logunit,*) 'CATCH_MAXSNDEPTH = ', CATCH_MAXSNDEPTH + write (logunit,*) 'CATCH_DZ1MAX = ', CATCH_DZ1MAX + write (logunit,*) 'CATCH_DZTSURF = ', CATCH_DZTSURF + write (logunit,*) 'CATCH_DZGT = ', CATCH_DZGT + write (logunit,*) 'CATCH_PHIGT = ', CATCH_PHIGT + write (logunit,*) 'CATCH_ALHMGT = ', CATCH_ALHMGT + write (logunit,*) 'CATCH_FSN = ', CATCH_FSN + write (logunit,*) 'CATCH_SHR = ', CATCH_SHR + write (logunit,*) 'CATCH_SCONST = ', CATCH_SCONST + write (logunit,*) 'CATCH_CSOIL_1 = ', CATCH_CSOIL_1 + write (logunit,*) 'CATCH_C_CANOP = ', CATCH_C_CANOP + write (logunit,*) 'CATCH_SATCAPFR = ', CATCH_SATCAPFR + write (logunit,*) + write (logunit,*) 'PEATCLSM_POROS_THRESHOLD = ', PEATCLSM_POROS_THRESHOLD + write (logunit,*) 'PEATCLSM_ZBARMAX_4_SYSOIL = ', PEATCLSM_ZBARMAX_4_SYSOIL write (logunit,*) write (logunit,*) 'Constants from SURFPARAMS:' write (logunit,*) - write (logunit,*) 'LAND_FIX = ', LAND_FIX - write (logunit,*) 'CSOIL_2 = ', CSOIL_2 - write (logunit,*) 'WEMIN = ', WEMIN - write (logunit,*) 'AICEV = ', AICEV - write (logunit,*) 'AICEN = ', AICEN - write (logunit,*) 'FLWALPHA = ', FLWALPHA - write (logunit,*) 'ASTRFR = ', ASTRFR - write (logunit,*) 'STEXP = ', STEXP - write (logunit,*) 'RSWILT = ', RSWILT + write (logunit,*) 'LAND_FIX = ', LAND_FIX + write (logunit,*) 'CSOIL_2 = ', CSOIL_2 + write (logunit,*) 'WEMIN = ', WEMIN + write (logunit,*) 'AICEV = ', AICEV + write (logunit,*) 'AICEN = ', AICEN + write (logunit,*) 'FLWALPHA = ', FLWALPHA + write (logunit,*) 'ASTRFR = ', ASTRFR + write (logunit,*) 'STEXP = ', STEXP + write (logunit,*) 'RSWILT = ', RSWILT write (logunit,*) write (logunit,*) 'end echo_catch_constants()' write (logunit,*) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index de7fe8ff2..4f4835c2e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -31,12 +31,14 @@ MODULE lsm_routines N_SNOW => CATCH_N_SNOW, & N_GT => CATCH_N_GT, & RHOFS => CATCH_SNWALB_RHOFS, & - DZTC => CATCH_DZTC, & + DZTSURF => CATCH_DZTSURF, & DZGT => CATCH_DZGT, & PHIGT => CATCH_PHIGT, & FSN => CATCH_FSN, & SHR => CATCH_SHR, & - N_SM => CATCH_N_ZONES + N_SM => CATCH_N_ZONES, & + PEATCLSM_POROS_THRESHOLD, & + PEATCLSM_ZBARMAX_4_SYSOIL USE SURFPARAMS, ONLY: & LAND_FIX, CSOIL_2, WEMIN, AICEV, AICEN, & @@ -53,10 +55,16 @@ MODULE lsm_routines PRIVATE PUBLIC :: INTERC, SRUNOFF, RZDRAIN, BASE, PARTITION, RZEQUIL, gndtp0 - PUBLIC :: SIBALB, catch_calc_soil_moist, catch_calc_subtile2tile + PUBLIC :: SIBALB, catch_calc_soil_moist, catch_calc_zbar, catch_calc_watertabled + PUBLIC :: catch_calc_subtile2tile PUBLIC :: gndtmp, catch_calc_tp, catch_calc_wtotl, catch_calc_ght, catch_calc_FT PUBLIC :: dampen_tc_oscillations, irrigation_rate + INTERFACE catch_calc_zbar + MODULE PROCEDURE catch_calc_zbar_scalar + MODULE PROCEDURE catch_calc_zbar_vector + END INTERFACE catch_calc_zbar + ! --------------------------------------------------------------------------- ! ! ***** Do not define *public* Catchment model constants here. ***** @@ -94,7 +102,7 @@ SUBROUTINE INTERC ( & TRAINL, TRAINC, SMELT, & ! [kg m-2 s-1] SATCAP,BUG, & CAPAC, & - THRUL, THRUC & ! [kg m-2] !!!! + THRUL_VOL, THRUC_VOL & ! [kg m-2] !!!! ) !**** !**** THIS ROUTINE USES THE PRECIPITATION FORCING TO DETERMINE @@ -116,10 +124,12 @@ SUBROUTINE INTERC ( & REAL, INTENT(INOUT), DIMENSION(NCH) :: CAPAC - REAL, INTENT(OUT), DIMENSION(NCH) :: THRUC, THRUL ! [kg m-2] ("volume" units) !!! - - INTEGER CHNO - REAL WETINT, WATADD, CAVAIL, THRU1, THRU2, XTCORR,SMPERS + REAL, INTENT(OUT), DIMENSION(NCH) :: THRUL_VOL, THRUC_VOL ! [kg m-2] ("volume" units) !!! + + ! -------------------------- + + INTEGER :: CHNO + REAL :: WETINT, WATADD, CAVAIL, THRU1, THRU2, XTCORR, SMPERS, THRUL, THRUC !**** !**** ------------------------------------------------------------------ @@ -161,7 +171,7 @@ SUBROUTINE INTERC ( & THRU2=XTCORR*WATADD - THRUL(CHNO)=THRU1+THRU2 + THRUL=THRU1+THRU2 CAPAC(CHNO)=CAPAC(CHNO)+WATADD-THRU1-THRU2 @@ -198,13 +208,14 @@ SUBROUTINE INTERC ( & THRU2=XTCORR*WATADD - THRUC(CHNO)=THRU1+THRU2 + THRUC=THRU1+THRU2 + CAPAC(CHNO)=CAPAC(CHNO)+WATADD-THRU1-THRU2 !**** - IF (THRUL(CHNO)+THRUC(CHNO) .LT. -1.e-8) WRITE(*,*) 'THRU= ', & - THRUL(CHNO), THRUC(CHNO), TRAINC(CHNO), TRAINL(CHNO), SMELT(CHNO) - THRUL(CHNO)=AMAX1(0., THRUL(CHNO)) - THRUC(CHNO)=AMAX1(0., THRUC(CHNO)) + IF (THRUL+THRUC .LT. -1.e-8) WRITE(*,*) 'THRU= ', & + THRUL, THRUC, TRAINC(CHNO), TRAINL(CHNO), SMELT(CHNO) + THRUL_VOL(CHNO)=AMAX1(0., THRUL) + THRUC_VOL(CHNO)=AMAX1(0., THRUC) 100 CONTINUE !**** @@ -219,13 +230,14 @@ END SUBROUTINE INTERC !**** /////////////////////////////////////////////////// !**** =================================================== - SUBROUTINE SRUNOFF ( & - NCH,DTSTEP,UFW4RO, FWETC, FWETL, AR1,ar2,ar4, & - THRUL_VOL, THRUC_VOL, & ! [kg m-2] ("volume" units) !!! - frice,tp1,srfmx, BUG, & - SRFEXC, & - RUNSRF, & ! [kg m-2 s-1] (flux units) - QINFIL & ! [kg m-2 s-1] (flux units) + SUBROUTINE SRUNOFF ( & + NCH, DTSTEP, UFW4RO, FWETC, FWETL, AR1, AR2, AR4, & + THRUL_VOL, THRUC_VOL, & ! [kg m-2] ("volume" units) !!! + FRICE, TP1, SRFMX, BUG, & + VGWMAX, RZEQ, POROS, & + SRFEXC, RZEXC, & + RUNSRF, & ! [kg m-2 s-1] (flux units) + QINFIL & ! [kg m-2 s-1] (flux units) ) !**** NOTE: Input throughfall is in volume units, as are calcs throughout this subroutine [kg m-2] @@ -238,19 +250,26 @@ SUBROUTINE SRUNOFF ( & INTEGER, INTENT(IN) :: NCH REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL LOGICAL, INTENT(IN) :: UFW4RO - REAL, INTENT(IN), DIMENSION(NCH) :: AR1, ar2, ar4, frice, tp1, srfmx + REAL, INTENT(IN), DIMENSION(NCH) :: AR1, AR2, AR4, FRICE, TP1, SRFMX + REAL, INTENT(IN), DIMENSION(NCH) :: VGWMAX, RZEQ, POROS REAL, INTENT(IN), DIMENSION(NCH) :: THRUL_VOL, THRUC_VOL ! [kg m-2] LOGICAL, INTENT(IN) :: BUG - - REAL, INTENT(INOUT), DIMENSION(NCH) :: SRFEXC + REAL, INTENT(INOUT), DIMENSION(NCH) :: SRFEXC, RZEXC REAL, INTENT(INOUT), DIMENSION(NCH) :: RUNSRF ! [kg m-2 s-1] - REAL, INTENT(OUT), DIMENSION(NCH) :: QINFIL ! [kg m-2 s-1] - INTEGER N - REAL deficit,srun0,frun,qin, qinfil_l, qinfil_c, qcapac, excess_infil, srunc, srunl, ptotal + ! --------------------------- + + INTEGER :: N + REAL :: deficit,srun0,frun,qin, qinfil_l, qinfil_c, qcapac, excess_infil + REAL :: srunc, srunl, ptotal, excess, totcapac, watadd REAL, DIMENSION(NCH) :: THRUL, THRUC + ! constants for PEATCLSM piecewise linear relationship between surface runoff and AR1 + + REAL, PARAMETER :: SRUN_AR1_MIN = 0.5 + REAL, PARAMETER :: SRUN_AR1_SLOPE = 10. + !**** - - - - - - - - - - - - - - - - - - - - - - - - - ! calculations throughout srunoff() are in "volume" units [kg m-2] @@ -263,69 +282,177 @@ SUBROUTINE SRUNOFF ( & DO N=1,NCH if(.not.UFW4RO) then - + PTOTAL=THRUL(N) + THRUC(N) - frun=AR1(N) - srun0=PTOTAL*frun - - !**** Comment out this line in order to allow moisture - !**** to infiltrate soil: - ! if(tp1(n) .lt. 0.) srun0=ptotal - - if(ptotal-srun0 .gt. srfmx(n)-srfexc(n)) & - srun0=ptotal-(srfmx(n)-srfexc(n)) - - if (srun0 .gt. ptotal) srun0=ptotal - - RUNSRF(N)=RUNSRF(N)+srun0 - QIN=PTOTAL-srun0 - + + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN + ! Non-peatland + frun=AR1(N) + srun0=PTOTAL*frun + + !**** Comment out this line in order to allow moisture + !**** to infiltrate soil: + ! if(tp1(n) .lt. 0.) srun0=ptotal + + if(ptotal-srun0 .gt. srfmx(n)-srfexc(n)) & + srun0=ptotal-(srfmx(n)-srfexc(n)) + + if (srun0 .gt. ptotal) srun0=ptotal + + RUNSRF(N)=RUNSRF(N)+srun0 + QIN=PTOTAL-srun0 + SRFEXC(N)=SRFEXC(N)+QIN + ELSE + ! Peatland + ! MB: no Hortonian surface runoff + !Note (rdk, from discussion w/MB; email 01/04/2021): + ! forcing all the rain to fall onto + ! the soil (1-AR1 fraction) rather than + ! onto both the soil and the free water surface is a simple shortcut; + ! the key function of this code is to retain all rainwater in the system + ! and *only* to produce surface runoff when the + ! ground is already ridiculously wet. + ! This prevents problems (e.g., numerical instabilities) found in + ! discharge calculations elsewhere in the code. + + srun0 = 0. + ! handling numerical instability due to exceptional snow melt events at some pixels + ! avoid AR1 to increase much higher than > 0.5 by enabling runoff + !Added ramping to avoid potential oscillations (rdk, 09/18/20) + IF (AR1(N)>SRUN_AR1_MIN) srun0=PTOTAL*amin1(1.,(ar1(n)-SRUN_AR1_MIN)*SRUN_AR1_SLOPE) + + ! MB: even no surface runoff when srfmx is exceeded (activating macro-pore flow) + ! Rewrote code to determine excess over capacity all at once (rdk, 09/18/20) + + totcapac=(srfmx(n)-srfexc(n))+(vgwmax(n)-(rzeq(n)+rzexc(n))) + watadd=ptotal-srun0 + if (watadd .gt. totcapac) then + excess=watadd-totcapac + srun0=srun0+excess + srfexc(n)=srfmx(n) + rzexc(n)=vgwmax(n)-rzeq(n) + elseif(watadd .gt. srfmx(n)-srfexc(n)) then + excess=watadd-(srfmx(n)-srfexc(n)) + srfexc(n)=srfmx(n) + rzexc(n)=rzexc(n)+excess + else + srfexc(n)=srfexc(n)+watadd + endif + ! MB: check if VGWMAX is exceeded + !IF(RZEQ(N) + RZEXC(N) .GT. (VGWMAX(N))) THEN + ! srun0 = srun0 + RZEQ(N)+RZEXC(N)-VGWMAX(N) + ! RZEXC(N)=VGWMAX(N)-RZEQ(N) + ! ENDIF + !(Commented out following lines to retain water balance -- rdk, 9/18/20) + !if (srun0 .gt. ptotal) then + ! srun0=ptotal + ! endif + RUNSRF(N)=RUNSRF(N)+srun0 + QIN=PTOTAL-srun0 + !SRFEXC(N)=amin1(SRFEXC(N)+QIN,srfmx(n)) + ENDIF + endif if(UFW4RO) then !**** Compute runoff from large-scale and convective storms separately: - - deficit=srfmx(n)-srfexc(n) - srunl=AR1(n)*THRUL(n) - qinfil_l=(1.-ar1(n))*THRUL(n) - qcapac=deficit*FWETL - - if(qinfil_l .gt. qcapac) then - excess_infil=qinfil_l-qcapac - srunl=srunl+excess_infil - qinfil_l=qinfil_l-excess_infil - endif - - srunc=AR1(n)*THRUC(n) - qinfil_c=(1.-ar1(n))*THRUC(n) - qcapac=deficit*FWETC - - if(qinfil_c .gt. qcapac) then - excess_infil=qinfil_c-qcapac - srunc=srunc+excess_infil - qinfil_c=qinfil_c-excess_infil + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN + !non-peatland + deficit=srfmx(n)-srfexc(n) + srunl=AR1(n)*THRUL(n) + qinfil_l=(1.-ar1(n))*THRUL(n) + qcapac=deficit*FWETL + + if(qinfil_l .gt. qcapac) then + excess_infil=qinfil_l-qcapac + srunl=srunl+excess_infil + qinfil_l=qinfil_l-excess_infil + endif + + srunc=AR1(n)*THRUC(n) + qinfil_c=(1.-ar1(n))*THRUC(n) + qcapac=deficit*FWETC + + if(qinfil_c .gt. qcapac) then + excess_infil=qinfil_c-qcapac + srunc=srunc+excess_infil + qinfil_c=qinfil_c-excess_infil + endif + + !**** Comment out this line in order to allow moisture + !**** to infiltrate soil: + ! if(tp1(n) .lt. 0.) srun0=ptotal + + if (srunl .gt. THRUL(n)) srunl=THRUL(n) + + if (srunc .gt. THRUC(n)) srunc=THRUC(n) + + RUNSRF(N)=RUNSRF(N)+srunl+srunc + QIN=THRUL(n)+THRUC(n)-(srunl+srunc) + SRFEXC(N)=SRFEXC(N)+QIN + + else + ! peatland + ! MB: no Hortonian surface runoff + !Note (rdk, from discussion w/MB; email 01/04/2021): + ! forcing all the rain to fall onto + ! the soil (1-AR1 fraction) rather than + ! onto both the soil and the free water surface is a simple shortcut; + ! the key function of this code is to retain all rainwater in the system + ! and *only* to produce surface runoff when the + ! ground is already ridiculously wet. + ! This prevents problems (e.g., numerical instabilities) found in + ! discharge calculations elsewhere in the code. + + srunl = 0. + srunc = 0. + ! handling numerical instability due to exceptional snow melt events at some pixels + ! avoid AR1 to increase much higher than > 0.5 by enabling runoff + IF (AR1(N)>SRUN_AR1_MIN) THEN + !Added ramping to avoid potential oscillations (rdk, 09/18/20) + srunl = THRUL(n)*amin1(1.,(ar1(n)-SRUN_AR1_MIN)*SRUN_AR1_SLOPE) + srunc = THRUC(n)*amin1(1.,(ar1(n)-SRUN_AR1_MIN)*SRUN_AR1_SLOPE) + ENDIF + PTOTAL = THRUL(N) + THRUC(N) + SRUN0 = srunl + srunc + ! MB: even no surface runoff when srfmx is exceeded (activating macro-pore flow) + ! Rewrote code to determine excess over capacity all at once (rdk, 09/18/20) + totcapac=(srfmx(n)-srfexc(n))+(vgwmax(n)-(rzeq(n)+rzexc(n))) + watadd=ptotal-srun0 + if (watadd .gt. totcapac) then + excess=watadd-totcapac + srun0=srun0+excess + srfexc(n)=srfmx(n) + rzexc(n)=vgwmax(n)-rzeq(n) + elseif(watadd .gt. srfmx(n)-srfexc(n)) then + excess=watadd-(srfmx(n)-srfexc(n)) + srfexc(n)=srfmx(n) + rzexc(n)=rzexc(n)+excess + else + srfexc(n)=srfexc(n)+watadd + endif + !if (ptotal-srun0 .gt. srfmx(n)-srfexc(n)) then + ! excess=(ptotal-srun0)-(srfmx(n)-srfexc(n)) + ! rzexc(n)=rzexc(n) + excess + ! ptotal=ptotal-excess + ! endif + ! MB: check if VGWMAX is exceeded + !IF(RZEQ(N) + RZEXC(N) .GT. (VGWMAX(N))) THEN + ! srun0 = srun0 + RZEQ(N)+RZEXC(N)-VGWMAX(N) + ! RZEXC(N)=VGWMAX(N)-RZEQ(N) + ! ENDIF + RUNSRF(N)=RUNSRF(N)+srun0 + QIN=PTOTAL-srun0 + ! SRFEXC(N)=amin1(SRFEXC(N)+QIN,srfmx(n)) endif - - !**** Comment out this line in order to allow moisture - !**** to infiltrate soil: - ! if(tp1(n) .lt. 0.) srun0=ptotal - - if (srunl .gt. THRUL(n)) srunl=THRUL(n) - - if (srunc .gt. THRUC(n)) srunc=THRUC(n) - - RUNSRF(N)=RUNSRF(N)+srunl+srunc - QIN=THRUL(n)+THRUC(n)-(srunl+srunc) endif - SRFEXC(N)=SRFEXC(N)+QIN - ! convert outputs to flux units [kg m-2 s-1] RUNSRF(N)=RUNSRF(N)/DTSTEP QINFIL(N)=QIN/DTSTEP - + END DO RETURN @@ -336,12 +463,15 @@ END SUBROUTINE SRUNOFF !**** ///////////////////////////////////////////////////////////////// !**** ----------------------------------------------------------------- !**** - SUBROUTINE RZDRAIN ( & - NCH,DTSTEP,VGWMAX,SATCAP,RZEQ,AR1,WPWET, & - tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & - CAPAC,RZEXC,SRFEXC,CATDEF, & - RUNSRF & ! [kg m-2 s-1] - ) + + + SUBROUTINE RZDRAIN ( & + NCH, DTSTEP, VGWMAX, SATCAP, RZEQ, AR1, WPWET, & + TSA1, TSA2, TSB1, TSB2, ATAU, BTAU, CDCR2, POROS, & + BF1, BF2, ARS1, ARS2, ARS3, BUG, & + CAPAC, RZEXC, SRFEXC, CATDEF, & + RUNSRF & ! [kg m-2 s-1] + ) !----------------------------------------------------------------- ! defines drainage timescales: @@ -352,19 +482,23 @@ SUBROUTINE RZDRAIN ( & !----------------------------------------------------------------- IMPLICIT NONE - INTEGER, INTENT(IN) :: NCH REAL, INTENT(IN) :: DTSTEP REAL, INTENT(IN), DIMENSION(NCH) :: VGWMAX, SATCAP, RZEQ, AR1, wpwet - REAL, INTENT(IN), DIMENSION(NCH) :: tsa1, tsa2, tsb1, tsb2, atau, btau, CDCR2, poros + REAL, INTENT(IN), DIMENSION(NCH) :: TSA1, TSA2, TSB1, TSB2, ATAU, BTAU, CDCR2, POROS + REAL, INTENT(IN), DIMENSION(NCH) :: BF1, BF2, ARS1, ARS2, ARS3 LOGICAL, INTENT(IN) :: BUG - REAL, INTENT(INOUT), DIMENSION(NCH) :: RZEXC, SRFEXC, CATDEF, CAPAC - REAL, INTENT(INOUT), DIMENSION(NCH) :: RUNSRF ! [kg m-2 s-1] + REAL, INTENT(INOUT), DIMENSION(NCH) :: RZEXC, SRFEXC, CATDEF, CAPAC + REAL, INTENT(INOUT), DIMENSION(NCH) :: RUNSRF ! [kg m-2 s-1] + + ! -------------------- INTEGER N - REAL srflw,rzflw,FLOW,EXCESS,TSC0,tsc2,rzave,rz0,wanom,rztot, & - rzx,btaux,ax,bx,rzdif, rzavemin + REAL srflw,rzflw,FLOW,EXCESS,TSC0,tsc2,rzave,rz0,wanom,rztot, & + rzx,btaux,ax,bx,rzdif, rzavemin,ZBAR1,SYSOIL,RZFLW_CATDEF, & + EXCESS_CATDEF, CATDEF_PEAT_THRESHOLD, RZFLW_AR1, AR1eq + !**** - - - - - - - - - - - - - - - - - - - - - - - - - @@ -440,36 +574,99 @@ SUBROUTINE RZDRAIN ( & IF (CATDEF(N)-RZFLW .GT. CDCR2(N)) then RZFLW=CATDEF(N)-CDCR2(N) - end if + end if - CATDEF(N)=CATDEF(N)-RZFLW - RZEXC(N)=RZEXC(N)-RZFLW + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) then + ! mineral soil + CATDEF(N)=CATDEF(N)-RZFLW + RZEXC(N)=RZEXC(N)-RZFLW + else + !MB2021: use AR1eq, equilibrium assumption between water level in soil hummocks and surface water level in hollows + AR1eq = (1.+ars1(n)*(catdef(n)))/(1.+ars2(n)*(catdef(n))+ars3(n)*(catdef(n))**2) + ! PEAT + ! MB: accounting for water ponding on AR1 + ! RZFLOW is partitioned into two flux components: (1) going in/out ponding water volume and (1) going in/out unsaturated soil storage + ! Specific yield of ponding water surface fraction is 1.0 + ! calculate SYSOIL (see Dettmann and Bechtold, VZJ, 2016, for detailed theory) + ! SYSOIL in CLSM can be derived from first derivative of + ! f_catdef(zbar) = ((zbar + bf2)^2 +1.0E-20)*bf1 + ! division by 1000 to convert from m to mm gives (Note: catdef in PEATCLSM remains + ! the soil profile deficit, i.e. does not include the ponding water storage). + ! SYSOIL = (2*bf1*zbar + 2*bf1*bf2)/1000 + ! Note: zbar defined here positive below ground. + ! For the SYSOIL estimation zbar must be constrained to 0.0 to 0.45 m, + ! to avoid extrapolation errors due to the non-optimal + ! (linear) approximation with the bf1-bf2-CLSM function, + ! theoretical SYSOIL curve levels off approximately at 0 m and 0.45 m. + ZBAR1 = catch_calc_zbar( BF1(N), BF2(N), CATDEF(N) ) + SYSOIL = (2.*bf1(n)*amin1(amax1(zbar1,0.),PEATCLSM_ZBARMAX_4_SYSOIL) + 2.*bf1(n)*bf2(n))/1000. + SYSOIL = amin1(SYSOIL,poros(n)) + ! Calculate fraction of RZFLW removed/added to catdef + RZFLW_CATDEF = (1.-AR1eq)*SYSOIL*RZFLW/(1.*AR1eq+SYSOIL*(1.-AR1eq)) + CATDEF(N)=CATDEF(N)-RZFLW_CATDEF + ! MB: remove all RZFLW from RZEXC because the other part + ! flows into the surface water storage (microtopgraphy) + RZEXC(N)=RZEXC(N)-RZFLW + + ENDIF + !**** REMOVE ANY EXCESS FROM MOISTURE RESERVOIRS: IF(CAPAC(N) .GT. SATCAP(N)) THEN RZEXC(N)=RZEXC(N)+CAPAC(N)-SATCAP(N) CAPAC(N)=SATCAP(N) - ENDIF + ENDIF IF(RZEQ(N) + RZEXC(N) .GT. VGWMAX(N)) THEN EXCESS=RZEQ(N)+RZEXC(N)-VGWMAX(N) RZEXC(N)=VGWMAX(N)-RZEQ(N) - CATDEF(N)=CATDEF(N)-EXCESS + + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN + CATDEF(N)=CATDEF(N)-EXCESS + ELSE + ! PEAT + ! MB: like for RZFLW --> EXCESS_CATDEF is the fraction in/out of catdef + EXCESS_CATDEF=(1.-AR1eq)*SYSOIL*EXCESS/(1.*AR1eq+SYSOIL*(1.-AR1eq)) + CATDEF(N)=CATDEF(N)-EXCESS_CATDEF ENDIF + ENDIF + + IF (POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + ! MB: CATDEF Threshold at zbar=0 + ! water table not allowed to rise higher (numerically instable) + ! zbar<0 only occurred due to extreme infiltration rates + ! (noticed this only snow melt events, very few locations and times) + ! (--> NOTE: PEATCLSM has no Hortonian runoff for zbar > 0) + CATDEF_PEAT_THRESHOLD = ((BF2(N))**2.0)*BF1(N) + IF(CATDEF(N) .LT. CATDEF_PEAT_THRESHOLD) THEN + ! RUNSRF(N)=RUNSRF(N) + (CATDEF_PEAT_THRESHOLD - CATDEF(N)) + ! runoff from AR1 for zbar>0 + ! RZFLW_AR1 = RZFLW - RZFLW_CATDEF + (CATDEF_PEAT_THRESHOLD - CATDEF(N)) + ! AR1=0.5 at zbar=0 + ! SYsurface=0.5 at zbar=0 + ! RUNSRF(N) = RUNSRF(N) + amax1(0.0, RZFLW_AR1 - 0.5*1000.*ZBAR1) + ! + ! revised (rdk, 1/04/2021): take excess water from both + ! soil and free standing water, the latter assumed to cover area AR1=0.5 + RUNSRF(N) = RUNSRF(N) + (CATDEF_PEAT_THRESHOLD-CATDEF(N) + 0.5*1000.*(-ZBAR1))/DTSTEP + CATDEF(N)=CATDEF_PEAT_THRESHOLD + ENDIF + ENDIF - IF(CATDEF(N) .LT. 0.) THEN + IF(CATDEF(N) .LT. 0.) THEN ! bug fix: RUNSRF in flux units [kg m-2 s-1] for consistency with partition() ! reichle, 6 Feb 2022 - RUNSRF(N)=RUNSRF(N)-CATDEF(N)/DTSTEP + RUNSRF(N)=RUNSRF(N)-CATDEF(N)/DTSTEP CATDEF(N)=0. - ENDIF + ENDIF 100 ENDDO RETURN END SUBROUTINE RZDRAIN + !**** ----------------------------------------------------------------- !**** ///////////////////////////////////////////////////////////////// !**** ----------------------------------------------------------------- @@ -477,6 +674,7 @@ END SUBROUTINE RZDRAIN SUBROUTINE BASE ( & NCH,DTSTEP,BF1,BF2,BF3,CDCR1,FRICE,COND,GNU, & + AR1, POROS, ARS1, ARS2, ARS3, & CATDEF, & BFLOW & ) @@ -487,7 +685,7 @@ SUBROUTINE BASE ( & INTEGER, INTENT(IN) :: NCH REAL, INTENT(IN) :: DTSTEP REAL, INTENT(IN), DIMENSION(NCH) :: BF1, BF2, BF3, CDCR1, FRICE, COND, & - GNU + GNU, AR1, POROS, ars1, ars2, ars3 REAL, INTENT(INOUT), DIMENSION(NCH) :: CATDEF @@ -495,29 +693,67 @@ SUBROUTINE BASE ( & INTEGER N - REAL ZBAR, ashift + REAL ZBAR, ashift, CFRICE,Ksz_zero,m_Ivanov,v_slope,Ta,dztmp,SYSOIL,BFLOW_CATDEF,ICERAMP,AR1eq data ashift/0./ - DO N=1,NCH - ! note intentionally opposite sign w.r.t. zbar defined above, - reichle, 16 Nov 2015 - ZBAR=SQRT(1.e-20+catdef(n)/bf1(n))-bf2(n) - BFLOW(N)=(1.-FRICE(N))*1000.* & - cond(n)*exp(-(bf3(n)-ashift)-gnu(n)*zbar)/gnu(n) -! *1000 is to convert from m/s to mm/s - IF (CATDEF(N) .GE. CDCR1(N)) BFLOW(N)=0. -!#ifdef LAND_UPD - IF (LAND_FIX) THEN - bflow(n)=amin1(1000.*cond(n),bflow(n)) - ELSE -!#else - bflow(n)=amin1(cond(n),bflow(n)) - END IF -!#endif - CATDEF(N)=CATDEF(N)+BFLOW(N)*dtstep - ENDDO + ZBAR = catch_calc_zbar( bf1(n), bf2(n), catdef(n) ) + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN + BFLOW(N)=(1.-FRICE(N))*1000.* & + cond(n)*exp(-(bf3(n)-ashift)-gnu(n)*zbar)/gnu(n) + ! *1000 is to convert from m/s to mm/s + IF (CATDEF(N) .GE. CDCR1(N)) BFLOW(N)=0. + !#ifdef LAND_UPD + IF (LAND_FIX) THEN + bflow(n)=amin1(1000.*cond(n),bflow(n)) + ELSE + !#else + bflow(n)=amin1(cond(n),bflow(n)) + END IF + !#endif + CATDEF(N)=CATDEF(N)+BFLOW(N)*dtstep + ELSE + ! PEAT + ! MB: + IF (FRICE(N) .GE. 0.9999) THEN + CFRICE = 1. + ELSE + CFRICE = FRICE(N) + ENDIF + ! BFLOW in mm/s + ! based on Ivanov 1981 + ! Ksz0 in m/s + ! m_Ivanov [-] value depends on unit of Ksz0 and z + ! v_slope in m^(-1) + Ksz_zero=10. + m_Ivanov=3.0 + v_slope = 1.5e-05 + ! Ta in m2/s, BFLOW in mm/s + Ta = (Ksz_zero*(1.+100.*amax1(0.,ZBAR))**(1.-m_Ivanov))/(100.*(m_Ivanov-1.)) + BFLOW(N)=v_slope*Ta*1000. + ! handling numerical instability due to extreme snow melt events on partly frozen ground + ! --> allow BFLOW/DISCHARGE for zbar .LE. 0.05 + ICERAMP= AMAX1(0., AMIN1(1., ZBAR/0.05)) + ICERAMP= 1.-ICERAMP*CFRICE + BFLOW(N)=ICERAMP*BFLOW(N) + + ! MB: Remove water from CATDEF and surface water storage + IF (BFLOW(N) .NE. 0.0) THEN + ! PEAT + ! MB: accounting for water ponding on AR1 + ! same approach as for RZFLW (see subroutine RZDRAIN for + ! comments) + SYSOIL = (2.*bf1(N)*amin1(amax1(zbar,0.),PEATCLSM_ZBARMAX_4_SYSOIL) + 2.*bf1(N)*bf2(N))/1000. + SYSOIL = amin1(SYSOIL,poros(n)) + !MB2021: use AR1eq, equilibrium assumption between water level in soil hummocks and surface water level in hollows + AR1eq = (1.+ars1(n)*(catdef(n)))/(1.+ars2(n)*(catdef(n))+ars3(n)*(catdef(n))**2) + BFLOW_CATDEF = (1.-AR1eq)*SYSOIL*BFLOW(N)/(1.*AR1eq+SYSOIL*(1.-AR1eq)) + CATDEF(N)=CATDEF(N)+BFLOW_CATDEF*dtstep + ENDIF + ENDIF + ENDDO RETURN END SUBROUTINE BASE @@ -532,6 +768,7 @@ END SUBROUTINE BASE SUBROUTINE PARTITION ( & NCH,DTSTEP,DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & PSIS,BEE,poros,WPWET, & + BF1, BF2, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4,BUG, & srfexc,catdef, & @@ -547,7 +784,7 @@ SUBROUTINE PARTITION ( & REAL, INTENT(IN) :: DTSTEP REAL, INTENT(IN), DIMENSION(NCH) :: DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & + PSIS,BEE,poros,WPWET,BF1,BF2, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4 @@ -568,6 +805,7 @@ SUBROUTINE PARTITION ( & ARG1, EXPARG1, ARG2, EXPARG2, ARG3, EXPARG3 !, surflay LOGICAL :: LSTRESS + REAL :: ZBAR, ARREST DATA LSTRESS/.FALSE./ !,surflay/20./ @@ -705,10 +943,22 @@ SUBROUTINE PARTITION ( & ENDIF + IF (POROS(N) >= PEATCLSM_POROS_THRESHOLD) THEN + ! peat + ! MB: AR4 (wilting fraction) for peatland depending on water table depth + !ZBAR defined here positive below ground and in meter + ZBAR = catch_calc_zbar( BF1(N), BF2(N), CATDEF(N) ) + AR4(N)=amax1(0.,amin1(1.0,(ZBAR-0.30)/(1.0))) + ARREST = 1.0 - AR1(N) + AR4(N)=amin1(ARREST,AR4(N)) + AR2(N)=1.0-AR4(n)-AR1(N) + ENDIF + RZI(N)=RZEQYI SWSRF1(N)=1. !mjs: changed .001 temporarily because of large bee. + IF (POROS(N) < PEATCLSM_POROS_THRESHOLD) THEN SWSRF2(N)=AMIN1(1., AMAX1(0.01, RZEQYI)) SWSRF4(N)=AMIN1(1., AMAX1(0.01, WILT)) @@ -720,6 +970,18 @@ SUBROUTINE PARTITION ( & SWSRF2(N)=((SWSRF2(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) SWSRF4(N)=((SWSRF4(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) + ELSE + + ! PEAT + ! MB: for peatlands integrate across surface soil moisture distribution + ! coefficients fitted for equilibrium conditions + ! SWSRF2 and SWSRF4 as wetness (not moisture) + ! MB: bug April 2018, AMIN1 function due to problems when spin up from + ! scratch (i.e. dry soil at time=0) + SWSRF2(N)=0.79437 - 0.99996*AMIN1(1.5,ZBAR) + 0.68801*(AMIN1(1.5,ZBAR))**2 + & + 0.04186*(AMIN1(1.5,ZBAR))**3 - 0.15042*(AMIN1(1.5,ZBAR))**4 + SWSRF4(N)=SWSRF2(N) + ENDIF ! srfmx is the maximum amount of water that can be added to the surface layer ! The choice of defining SWSRF4 like SWSRF2 needs to be better examined. @@ -782,7 +1044,7 @@ END SUBROUTINE PARTITION !**** ----------------------------------------------------------------- SUBROUTINE RZEQUIL ( & - NCH,CATDEF,VGWMAX,CDCR1,CDCR2,WPWET, & + NCH,CATDEF,VGWMAX,CDCR1,CDCR2,WPWET,POROS, & ars1,ars2,ars3,ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4, & RZEQ & @@ -793,7 +1055,7 @@ SUBROUTINE RZEQUIL ( & INTEGER, INTENT(IN) :: NCH REAL, INTENT(IN), DIMENSION(NCH) :: CATDEF, VGWMAX, CDCR1, CDCR2, & WPWET, ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, & - arw2, arw3, arw4 + arw2, arw3, arw4, POROS REAL, INTENT(OUT), DIMENSION(NCH) :: RZEQ @@ -895,7 +1157,7 @@ subroutine gndtp0(t1,phi,zbar,thetaf,ht,fh21w,fh21i,fh21d, & ! calculate the boundaries, based on the layer thicknesses(DZGT) - zb(1)=-DZTC + zb(1)=-DZTSURF zb(2)=zb(1)-DZGT(1) shc(1)=shr0*(1.-phi)*DZGT(1) zc(1)=0.5*(zb(1)+zb(2)) @@ -958,7 +1220,7 @@ subroutine gndtp0(t1,phi,zbar,thetaf,ht,fh21w,fh21i,fh21d, & ! *********************************** ! lets get the thermal conductivity for the layers - a1=1-phi + a1=1.-phi tk1=1.01692+a1*(0.89865+1.06211*a1) xw=phi*(1.-fice(1)) a2=phi-xw @@ -980,7 +1242,7 @@ subroutine gndtp0(t1,phi,zbar,thetaf,ht,fh21w,fh21i,fh21d, & xklh(1)=(tksat-tkdry)*xwi + tkdry xklhw=tksat - denom=-(DZTC*0.5)-zc(1) + denom=-(DZTSURF*0.5)-zc(1) fh21w=-xklhw *(t1(1)-TF-tp(1))/denom fh21i=-xklh(1)*(t1(2)-TF-tp(1))/denom fh21d=-xklh(1)*(t1(3)-TF-tp(1))/denom @@ -1474,9 +1736,9 @@ END SUBROUTINE SIBALB ! ================================================================================ subroutine catch_calc_soil_moist( & - NTILES,vegcls,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & + NTILES,dzsf,vgwmax,cdcr1,cdcr2,psis,bee,poros,wpwet, & ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & + ara3,ara4,arw1,arw2,arw3,arw4,bf1, bf2, & srfexc,rzexc,catdef, & ar1, ar2, ar4, & sfmc, rzmc, prmc, & @@ -1528,11 +1790,10 @@ subroutine catch_calc_soil_moist( & implicit none integer, intent(in) :: NTILES - integer, dimension(NTILES), intent(in) :: vegcls real, dimension(NTILES), intent(in) :: dzsf,vgwmax,cdcr1,cdcr2 real, dimension(NTILES), intent(in) :: wpwet,poros,psis - real, dimension(NTILES), intent(in) :: bee,ars1 + real, dimension(NTILES), intent(in) :: bee,ars1,bf1, bf2 real, dimension(NTILES), intent(in) :: ars2,ars3,ara1,ara2,ara3 real, dimension(NTILES), intent(in) :: ara4,arw1,arw2,arw3,arw4 @@ -1547,7 +1808,7 @@ subroutine catch_calc_soil_moist( & real, dimension(NTILES), intent(out), optional :: sfmcun real, dimension(NTILES), intent(out), optional :: rzmcun real, dimension(NTILES), intent(out), optional :: prmcun - real, dimension(NTILES), intent(out), optional :: swsrf1out, swsrf2out, swsrf4out + real, dimension(NTILES), intent(out), optional :: swsrf1out, swsrf2out, swsrf4out ! ---------------------------- ! @@ -1594,7 +1855,7 @@ subroutine catch_calc_soil_moist( & call rzequil( & NTILES, catdef, vgwmax, & - cdcr1, cdcr2, wpwet, & + cdcr1, cdcr2, wpwet, poros, & ars1, ars2, ars3, ara1, ara2, ara3, ara4, & arw1, arw2, arw3, arw4, & rzeq) @@ -1632,7 +1893,7 @@ subroutine catch_calc_soil_moist( & call partition( & NTILES,dtstep_dummy,dzsf,rzexc, & rzeq,vgwmax,cdcr1,cdcr2, & - psis,bee,poros,wpwet, & + psis,bee,poros,wpwet,bf1, bf2, & ars1,ars2,ars3, & ara1,ara2,ara3,ara4, & arw1,arw2,arw3,arw4,.false., & @@ -1704,6 +1965,62 @@ subroutine catch_calc_soil_moist( & end subroutine catch_calc_soil_moist + ! ******************************************************************* + + ! Calculate zbar for Catchment[CN] model. + ! + ! Convention: zbar positive below ground (downward). + ! + ! This convention applies to water calculations, incl. subroutines RZDRAIN(), + ! WUPDAT(), BASE(), PEATCLSM + ! + ! WARNING: + ! Opposite convention applies when zbar is used in ground heat + ! diffusion model, incl. subroutines GNDTP0(), GNDTMP(), GNDTMP_CN(). + ! + ! - reichle, 29 Jan 2022 + + function catch_calc_zbar_scalar( bf1, bf2, catdef ) result(zbar) + + implicit none + + real, intent(in) :: bf1, bf2, catdef + real :: zbar + + zbar = SQRT(1.e-20 + catdef/bf1) - bf2 + + end function catch_calc_zbar_scalar + + ! -------------------------- + + function catch_calc_zbar_vector( bf1, bf2, catdef ) result(zbar) + + ! vector version of catch_calc_zbar + + implicit none + + real, dimension(:), intent(in) :: bf1, bf2, catdef + real, dimension(size(bf1)) :: zbar + + zbar = SQRT(1.e-20 + catdef/bf1) - bf2 + + end function catch_calc_zbar_vector + + ! ******************************************************************* + + function catch_calc_watertabled( bf1, bf2, cdcr2, poros, wpwet, catdef ) result(wtd) + + ! calculate water table depth [m] + + implicit none + + real, dimension(:), intent(in) :: bf1, bf2, cdcr2, poros, wpwet, catdef + real, dimension(size(bf1)) :: wtd + + wtd = MIN( catch_calc_zbar(BF1,BF2,CATDEF), CDCR2/(1.-WPWET)/POROS/1000. ) + + end function catch_calc_watertabled + ! ******************************************************************* subroutine catch_calc_subtile2tile( NTILES, ar1, ar2, asnow, subtile_data, tile_data ) @@ -1781,18 +2098,18 @@ subroutine gndtmp(phi,zbar,ht,xfice,tp,FICE,dts,thetaf,fh21) REAL, INTENT(INOUT), DIMENSION(*) :: HT ! dimension(N_GT) - REAL, INTENT(OUT) :: XFICE + REAL, INTENT(OUT) :: XFICE REAL, INTENT(OUT), DIMENSION(*) :: TP, FICE ! dimension(N_GT) ! ---------------------------------- INTEGER :: L, LSTART, K - REAL, DIMENSION(N_GT) :: ZC, SHC, XKLH - REAL, DIMENSION(N_GT+1) :: FH, ZB + REAL, DIMENSION(N_GT) :: ZC, SHC, XKLH + REAL, DIMENSION(N_GT+1) :: FH, ZB REAL :: SHW0, SHI0, SHR0, WS, XW, A1, TK1, A2, TK2, TK3, TKSAT REAL :: XWI, XD1, XD2, XKTH, TKDRY LOGICAL :: full_version - + !data dz/0.0988,0.1952,0.3859,0.7626,1.5071,10.0/ !DATA PHI/0.45/, FSN/3.34e+8/, SHR/2.4E6/ @@ -1826,13 +2143,13 @@ subroutine gndtmp(phi,zbar,ht,xfice,tp,FICE,dts,thetaf,fh21) !---------------------------------- ! calculate the boundaries, based on the layer thicknesses(DZGT) - - zb(1)=-DZTC ! Bottom of surface layer, which is handled outside - ! this routine. + + zb(1)=-DZTSURF ! Bottom of surface layer, which is handled outside + ! this routine. do l=1,N_GT zb(l+1)=zb(l)-DZGT(l) - shc(l)=shr0*(1-phi)*DZGT(l) + shc(l)=shr0*(1.-phi)*DZGT(l) enddo ! evaluates the temperatures in the soil layers based on the heat values. @@ -1880,128 +2197,128 @@ subroutine gndtmp(phi,zbar,ht,xfice,tp,FICE,dts,thetaf,fh21) zc(l)=0.5*(zb(l)+zb(l+1)) enddo - ! evaluates: layer thermal conductivities - ! ***************************************** - ! from farouki(cold regions sci and tech, 5, 1981, - ! 67-75) the values for tk1,tk2,tk3 are as follows: - ! tk2=2.2**(phi(l,ibv)-xw), tk3=.57**xw, and - ! tk1=3**(1-phi(l,ibv) for %sand<.5, and - ! for computation purposes i have fit these eqs. - ! to 2nd order polynomials. - ! *********************************** - ! input: - ! sklh - soil heat conductivities of layers - ! zb - soil layer boundaries, m - ! zc - soil layer centers, m - ! DZGT - layer thickness, m - ! w - soil water content, m - ! phi - soil porosity, dimensionless - ! q - % sand, silt, clay, peat - ! fice - fraction of ice in layers - ! output: - ! xklh - thermal conductivity, w m-2 k-1 - ! *********************************** - - ! lets get the thermal conductivity for the layers - - do k=1,N_GT - - a1=1-phi ! ROCK FRACTION - tk1=1.01692+a1*(0.89865+1.06211*a1) - xw=phi*(1.-fice(k)) ! FOR SATURATED SOIL, XW HERE IS - ! THE LIQUID WATER FRACTION - a2=phi-xw ! FOR SATURATED SOIL, THE ICE FRACTION - - tk2=1.00543+a2*(0.723371+.464342*a2) - tk3=0.998899+xw*(-0.548043+0.120291*xw) - tksat=tk1*tk2*tk3 - - xwi=1.0 - if (zbar .le. zb(k+1))then - xwi=thetaf - elseif (zbar .ge. zb(k+1) .and. zbar .le. zb(k))then - xd1=zb(k)-zbar - xd2=zbar-zb(k+1) - xwi=((xd1*thetaf)+xd2)/(xd1+xd2) - endif +! evaluates: layer thermal conductivities +! ***************************************** +! from farouki(cold regions sci and tech, 5, 1981, +! 67-75) the values for tk1,tk2,tk3 are as follows: +! tk2=2.2**(phi(l,ibv)-xw), tk3=.57**xw, and +! tk1=3**(1-phi(l,ibv) for %sand<.5, and +! for computation purposes i have fit these eqs. +! to 2nd order polynomials. +! *********************************** +! input: +! sklh - soil heat conductivities of layers +! zb - soil layer boundaries, m +! zc - soil layer centers, m +! DZGT - layer thickness, m +! w - soil water content, m +! phi - soil porosity, dimensionless +! q - % sand, silt, clay, peat +! fice - fraction of ice in layers +! output: +! xklh - thermal conductivity, w m-2 k-1 +! *********************************** + +! lets get the thermal conductivity for the layers + + do k=1,N_GT + + a1=1.-phi ! ROCK FRACTION + tk1=1.01692+a1*(0.89865+1.06211*a1) + xw=phi*(1.-fice(k)) ! FOR SATURATED SOIL, XW HERE IS + ! THE LIQUID WATER FRACTION + a2=phi-xw ! FOR SATURATED SOIL, THE ICE FRACTION + + tk2=1.00543+a2*(0.723371+.464342*a2) + tk3=0.998899+xw*(-0.548043+0.120291*xw) + tksat=tk1*tk2*tk3 + + xwi=1.0 + if (zbar .le. zb(k+1))then + xwi=thetaf + elseif (zbar .ge. zb(k+1) .and. zbar .le. zb(k))then + xd1=zb(k)-zbar + xd2=zbar-zb(k+1) + xwi=((xd1*thetaf)+xd2)/(xd1+xd2) + endif + + xwi=min(xwi,1.) + tkdry=0.226 ! = .039*0.45^(-2.2), from Farouki, p. 71 + xklh(k)=(tksat-tkdry)*xwi + tkdry + + enddo + +! evaluates heat flux between layers due to heat diffussion +! *********************************** +! input: +! zb - soil layer boundaries, m +! zc - soil layer centers, m +! DZGT - layer thickness, m +! fice - fraction of ice in layers +! tp - temperature of layers, c +! shw - specific heat of water +! shi - specific heat of ice +! fsn - heat of fusion J/m +! output: +! fh - heat flux between layers +! *********************************** - xwi=min(xwi,1.) - tkdry=0.226 ! = .039*0.45^(-2.2), from Farouki, p. 71 - xklh(k)=(tksat-tkdry)*xwi + tkdry - - enddo - - ! evaluates heat flux between layers due to heat diffussion - ! *********************************** - ! input: - ! zb - soil layer boundaries, m - ! zc - soil layer centers, m - ! DZGT - layer thickness, m - ! fice - fraction of ice in layers - ! tp - temperature of layers, c - ! shw - specific heat of water - ! shi - specific heat of ice - ! fsn - heat of fusion J/m - ! output: - ! fh - heat flux between layers - ! *********************************** - - - ! total heat flux is via diffusion along the temperature gradient - fh(N_GT+1)=0. - fh(1)=fh21 - do k=2,N_GT - ! THIS xkth is NEW (ie., Agnes corrected) - it should be fixed in all - ! codes I'm using - xkth=((zb(k)-zc(k-1))*xklh(k-1)+(zc(k)-zb(k))*xklh(k)) & - /(zc(k)-zc(k-1)) - fh(k)=-xkth*(tp(k-1)-tp(k))/(zc(k-1)-zc(k)) - enddo - - ! update the heat contents in the model layers; ht(l) - ! IF THERE'S SNOW THIS WILL HAVE TO BE MODIFIED L=1,N - - do k=1,N_GT - ht(k)=ht(k)+(fh(k+1)-fh(k))*dts - enddo - - ! evaluates the temperatures in the soil layers based on the heat - ! values. - ! *********************************** - ! input: - ! xw - water in soil layers, m - ! ht - heat in soil layers - ! fsn - heat of fusion of water J/m - ! shc - specific heat capacity of soil - ! shi - specific heat capacity of ice - ! shw - specific heat capcity of water - ! snowd - snow depth, equivalent water m - ! output: - ! tp - temperature of layers, c - ! fice - fraction of ice of layers - ! pre - extra precipitation, i.e. snowmelt, m s-1 - ! snowd - snow depth after melting, equivalent water m. - ! *********************************** - ! determine fraction of ice and temp of soil layers based on layer - ! heat and water content - do 1000 k=1,N_GT - ws=phi*DZGT(k) ! saturated water content -! xl=l -! xw=(1/(7-xl))*ws - xw=0.5*ws ! For calculations here, assume soil + +! total heat flux is via diffusion along the temperature gradient + fh(N_GT+1)=0. + fh(1)=fh21 + do k=2,N_GT +! THIS xkth is NEW (ie., Agnes corrected) - it should be fixed in all +! codes I'm using + xkth=((zb(k)-zc(k-1))*xklh(k-1)+(zc(k)-zb(k))*xklh(k)) & + /(zc(k)-zc(k-1)) + fh(k)=-xkth*(tp(k-1)-tp(k))/(zc(k-1)-zc(k)) + enddo + +! update the heat contents in the model layers; ht(l) +! IF THERE'S SNOW THIS WILL HAVE TO BE MODIFIED L=1,N + + do k=1,N_GT + ht(k)=ht(k)+(fh(k+1)-fh(k))*dts + enddo + +! evaluates the temperatures in the soil layers based on the heat +! values. +! *********************************** +! input: +! xw - water in soil layers, m +! ht - heat in soil layers +! fsn - heat of fusion of water J/m +! shc - specific heat capacity of soil +! shi - specific heat capacity of ice +! shw - specific heat capcity of water +! snowd - snow depth, equivalent water m +! output: +! tp - temperature of layers, c +! fice - fraction of ice of layers +! pre - extra precipitation, i.e. snowmelt, m s-1 +! snowd - snow depth after melting, equivalent water m. +! *********************************** +! determine fraction of ice and temp of soil layers based on layer +! heat and water content + do 1000 k=1,N_GT + ws=phi*DZGT(k) ! saturated water content +! xl=l +! xw=(1/(7-xl))*ws + xw=0.5*ws ! For calculations here, assume soil ! is half full of water. - fice(k)=AMAX1( 0., AMIN1( 1., -ht(k)/(fsn*xw) ) ) - - IF(FICE(K) .EQ. 1.) THEN - tp(k)=(ht(k)+xw*fsn)/(shc(k)+xw*shi0) - ELSEIF(FICE(K) .EQ. 0.) THEN - tp(k)=ht(k)/(shc(k)+xw*shw0) - ELSE - TP(K)=0. - ENDIF + fice(k)=AMAX1( 0., AMIN1( 1., -ht(k)/(fsn*xw) ) ) + + IF(FICE(K) .EQ. 1.) THEN + tp(k)=(ht(k)+xw*fsn)/(shc(k)+xw*shi0) + ELSEIF(FICE(K) .EQ. 0.) THEN + tp(k)=ht(k)/(shc(k)+xw*shw0) + ELSE + TP(K)=0. + ENDIF + + 1000 continue - 1000 continue - end if ! full version of gndtmp() ! determine the value of xfice @@ -2016,13 +2333,32 @@ subroutine gndtmp(phi,zbar,ht,xfice,tp,FICE,dts,thetaf,fh21) do l=lstart,N_GT xfice=xfice+fice(l) - enddo - xfice=xfice/((N_GT+1)-lstart) + enddo + + IF (phi < PEATCLSM_POROS_THRESHOLD) THEN + xfice=xfice/((N_GT+1)-lstart) + ELSE + !PEAT + !MB: only first layer for total runoff reduction + xfice=AMIN1(1.0,fice(1)) + ENDIF Return end subroutine gndtmp + ! ******************************************************************* + + + + + + + + + + + ! ******************************************************************* subroutine catch_calc_tp( NTILES, poros, ghtcnt, tp, fice ) @@ -2081,7 +2417,7 @@ subroutine catch_calc_tp( NTILES, poros, ghtcnt, tp, fice ) do k=1,N_gt - shc(k) = shr0*(1-phi)*DZGT(k) + shc(k) = shr0*(1.-phi)*DZGT(k) ws = phi*DZGT(k) ! PORE SPACE IN LAYER @@ -2205,7 +2541,7 @@ subroutine catch_calc_ght( dzgt, poros, tp, fice, ghtcnt ) phi=PHIGT end if - shc = shr0*(1-phi)*DZGT + shc = shr0*(1.-phi)*DZGT ws = phi*DZGT ! PORE SPACE IN LAYER diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/SurfParams.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/SurfParams.F90 index 3e2398649..1ed10273f 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/SurfParams.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/SurfParams.F90 @@ -2,23 +2,23 @@ !========================================================== module SurfParams - + ! Justin, 12 Apr 2018 - Created - Replaces LAND_UPD ifdefs , added SurfParams_init, ! called from GEOS_LandGridCompMod, GEOS_LandiceGridCompMod ! jkolassa, 11 Jun 2020 - added LSM_CHOICE as input; introduced separate land parameter collections ! use MAPL_ExceptionHandling implicit none - + private - + public :: SurfParams_init - + ! --------------------------------------------------------------------------- ! Switch for catchment physics parameters ! --------------------------------------------------------------------------- - - + + ! Variables set by SurfParams_init: REAL, PUBLIC, SAVE :: CSOIL_2 ! J/K - heat capacity associated w/ tsurf REAL, PUBLIC, SAVE :: WEMIN ! kg/m^2 minimum SWE in areal fraction @@ -26,113 +26,113 @@ module SurfParams REAL, PUBLIC, SAVE :: FLWALPHA ! SRFLW multiplier with SRFLW < 0 REAL, PUBLIC, SAVE :: ASTRFR, STEXP ! stress parameters in energy2 REAL, PUBLIC, SAVE :: RSWILT ! parameters in rsurfp - + LOGICAL, PUBLIC, SAVE :: LAND_FIX ! Used for fixes and init changes that ! are still default in Icarus GCM - contains - +contains + ! Call to get "constants" that really are variables changeable during land/landice initialization - + subroutine SurfParams_init(LAND_PARAMS,LSM_CHOICE, rc) - + implicit none - + CHARACTER(*), INTENT(IN) :: LAND_PARAMS INTEGER, INTENT(IN) :: LSM_CHOICE INTEGER, OPTIONAL, INTENT(OUT) :: rc - + LOGICAL, SAVE :: init_called = .FALSE. ! Flag if SurfParams_init has been called - + ! --------------------------------------------------------------------------- - + if (init_called) then ! already called !write (*,*) "SurfParams_init being called again" return end if - + if (LSM_CHOICE==1) then - - select case (LAND_PARAMS) - case ("Icarus") ! "Old" LDASsa physics, current default for Icarus GCM - LAND_FIX = .FALSE. - CSOIL_2 = 200. - WEMIN = 26. - AICEV = 0.149 - AICEN = 19.851 - FLWALPHA = 1. ! i.e., FLWALPHA unchanged - ASTRFR = 0.333 - STEXP = 1. - RSWILT = 500. - - case ("V24_C05") ! V24_C05 changes, default for LDAS m4-17-0 - LAND_FIX = .TRUE. - CSOIL_2 = 70000. ! Post H5_0 - WEMIN = 13. - AICEV = 0.107 - AICEN = 19.893 - FLWALPHA = 0.01 - ASTRFR = 1. - STEXP = 2. - RSWILT = 2000. - - case ("NRv7.2") ! SMAP NRv7.2 changes, default for after LDAS m4-17-6 - LAND_FIX = .TRUE. - CSOIL_2 = 70000. ! Post H5_0 - WEMIN = 13. - AICEV = 0.149 - AICEN = 19.851 - FLWALPHA = 0.04 - ASTRFR = 0.333 ! reverted - STEXP = 1. ! reverted - RSWILT = 500. ! reverted - - case DEFAULT - _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') - end select - + + select case (LAND_PARAMS) + case ("Icarus") ! "Old" LDASsa physics, current default for Icarus GCM + LAND_FIX = .FALSE. + CSOIL_2 = 200. + WEMIN = 26. + AICEV = 0.149 + AICEN = 19.851 + FLWALPHA = 1. ! i.e., FLWALPHA unchanged + ASTRFR = 0.333 + STEXP = 1. + RSWILT = 500. + + case ("V24_C05") ! V24_C05 changes, default for LDAS m4-17-0 + LAND_FIX = .TRUE. + CSOIL_2 = 70000. ! Post H5_0 + WEMIN = 13. + AICEV = 0.107 + AICEN = 19.893 + FLWALPHA = 0.01 + ASTRFR = 1. + STEXP = 2. + RSWILT = 2000. + + case ("NRv7.2") ! SMAP NRv7.2 changes, default for after LDAS m4-17-6 + LAND_FIX = .TRUE. + CSOIL_2 = 70000. ! Post H5_0 + WEMIN = 13. + AICEV = 0.149 + AICEN = 19.851 + FLWALPHA = 0.04 + ASTRFR = 0.333 ! reverted + STEXP = 1. ! reverted + RSWILT = 500. ! reverted + + case DEFAULT + _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') + end select + else if (LSM_CHOICE==2) then - - select case (LAND_PARAMS) - case ("CN_CLM40") ! parameters to reproduce Fanwei Zeng's Catchment-CN.4.0 runs (e0004s_transientCO2_05) done with build /gpfsm/dnb31/fzeng/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 - LAND_FIX = .TRUE. - CSOIL_2 = 70000. ! Post H5_0 - WEMIN = 13. - AICEV = 0.149 - AICEN = 19.851 - FLWALPHA = 1. - ASTRFR = 0.333 ! reverted - STEXP = 1. ! reverted - RSWILT = 1500. - - case DEFAULT - _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') - end select - - else if (LSM_CHOICE==3) then - select case (LAND_PARAMS) - - case ("CN_CLM45") ! parameters to reproduce Eunjee Lee's Catchment-CN4.5 fire carbon emission simulations - LAND_FIX = .TRUE. - CSOIL_2 = 70000. ! Post H5_0 - WEMIN = 13. - AICEV = 0.107 - AICEN = 19.893 - FLWALPHA = 0.005 - ASTRFR = 0.333 ! reverted - STEXP = 1. ! reverted - RSWILT = 2000. - - case DEFAULT - _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') - end select + + select case (LAND_PARAMS) + case ("CN_CLM40") ! parameters to reproduce Fanwei Zeng's Catchment-CN.4.0 runs (e0004s_transientCO2_05) done with build /gpfsm/dnb31/fzeng/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 + LAND_FIX = .TRUE. + CSOIL_2 = 70000. ! Post H5_0 + WEMIN = 13. + AICEV = 0.149 + AICEN = 19.851 + FLWALPHA = 1. + ASTRFR = 0.333 ! reverted + STEXP = 1. ! reverted + RSWILT = 1500. + + case DEFAULT + _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') + end select + + else if (LSM_CHOICE==3) then + select case (LAND_PARAMS) + + case ("CN_CLM45") ! parameters to reproduce Eunjee Lee's Catchment-CN4.5 fire carbon emission simulations + LAND_FIX = .TRUE. + CSOIL_2 = 70000. ! Post H5_0 + WEMIN = 13. + AICEV = 0.107 + AICEN = 19.893 + FLWALPHA = 0.005 + ASTRFR = 0.333 ! reverted + STEXP = 1. ! reverted + RSWILT = 2000. + + case DEFAULT + _ASSERT(.FALSE.,'LAND_PARAMS not valid or incompatible with LSM_CHOICE') + end select else - _ASSERT(.FALSE.,'land model choice not valid') + _ASSERT(.FALSE.,'land model choice not valid') end if ! LSM_CHOICE - + init_called = .TRUE. _RETURN(_SUCCESS) - + end subroutine SurfParams_init - + endmodule SurfParams diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 2675f127d..3d25d1671 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -185,7 +185,7 @@ PROGRAM mkCatchParam GridnameT='til/'//trim(Gridname) endif - if(process_peat) PEATSOURCE = 'PEATMAP' + if(use_PEATMAP) PEATSOURCE = 'PEATMAP' if(jpl_height) VEGZSOURCE = 'JPL' if(n_threads == 1) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 74e8c80c0..5c776d1f0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3758,7 +3758,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) logical :: regrid,write_file INTEGER, allocatable, dimension (:) :: soil_class_top,soil_class_com REAL :: sf,factor,wp_wetness,fac_count - logical :: file_exists + logical :: CatchParamsNC_file_exists REAL, ALLOCATABLE, DIMENSION (:,:) :: parms4file ! PEAT-clsm modification ! Below parameters are from Table 2 of: @@ -4036,7 +4036,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) end do end do - if(process_peat) then + if(use_PEATMAP) then print *, 'PMAP_THRESH : ', pmap_thresh allocate(pmapr (1:i_highd,1:j_highd)) status = NF_OPEN ('data/CATCH/PEATMAP_mask.nc4', NF_NOWRITE, ncid) @@ -4274,7 +4274,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) allocate(btau_2cm(1:n_SoilClasses)) allocate(a_wpsurf(1:n_SoilClasses)) allocate(a_porosurf(1:n_SoilClasses)) - if(process_peat) then + if(use_PEATMAP) then fname = trim(c_data)//'SoilClasses-SoilHyd-TauParam.peatmap' else fname = trim(c_data)//'SoilClasses-SoilHyd-TauParam.dat' @@ -4375,7 +4375,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) !$OMP data_vec4,data_vec5,data_vec6,cF_lim, & !$OMP table_map,soil_class_top,soil_class_com, & !$OMP soc_vec,poc_vec,ncells_top,ncells_top_pro,& -!$OMP ncells_sub_pro,process_peat) & +!$OMP ncells_sub_pro,use_PEATMAP) & !$OMP PRIVATE(n,i,j,k,icount,t_count,i1,i2,ss_clay, & !$OMP ss_sand,ss_clay_all,ss_sand_all, & !$OMP ss_oc_all,cFamily,factor,o_cl,o_clp,ktop, & @@ -4439,7 +4439,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) if (sum(cFamily) == 0.) o_cl = 1 if (sum(cFamily) > 0.) o_cl = maxloc(cFamily, dim = 1) - if (process_peat) then + if (use_PEATMAP) then ! if 50% or more of the tile surface is covered with peat, we assume the tile is peat if (cFamily(4)/real(i) > 0.5) then o_cl = 4 @@ -4604,9 +4604,9 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) ! call process_peatmap (nx, ny, gfiler, pmap) - inquire(file='clsm/catch_params.nc4', exist=file_exists) + inquire(file='clsm/catch_params.nc4', exist=CatchParamsNC_file_exists) - if(file_exists) then + if(CatchParamsNC_file_exists) then status = NF_OPEN ('clsm/catch_params.nc4', NF_WRITE, ncid) ; VERIFY_(STATUS) allocate (parms4file (1:maxcat, 1:10)) endif @@ -4664,7 +4664,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) fac_surf = soil_class_top(n) fac = soil_class_com(n) - if(process_peat) then + if(use_PEATMAP) then ! the maximum peat soil depth is set to the value Michel used to derive parameters (1334.) if (fac_surf == 253) soildepth(n) = 5000. ! max(soildepth(n),5000.) ! reseet subsurface tro peat if surface soil type is peat @@ -4699,8 +4699,8 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) endif end do write (11,'(a)')' ' - write (11,'(a)')'FMT=i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4' - write (11,'(a)')'TileIndex PfafID SoilClassTop SoilClassProfile BEE PSIS POROS Ks_at_SURF WPWET SoilDepth %Grav %OCTop %OCProf %Sand_top %Clay_top %Sand_prof %Clay_prof WPWET_SURF POROS_SURF' + write (11,'(a)')'FMT=i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4,f8.4' + write (11,'(a)')'TileIndex PfafID SoilClassTop SoilClassProfile BEE PSIS POROS Ks_at_SURF WPWET SoilDepth %Grav %OCTop %OCProf %Sand_top %Clay_top %Sand_prof %Clay_prof WPWET_SURF POROS_SURF PMAP' close (10, status = 'keep') close (11, status = 'keep') close (12, status = 'keep') @@ -4713,7 +4713,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler) atau_2cm,btau_2cm) deallocate (soildepth, grav_vec,soc_vec,poc_vec,& ncells_top,ncells_top_pro,ncells_sub_pro,soil_class_top,soil_class_com) - if(file_exists) then + if(CatchParamsNC_file_exists) then status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'BEE' ) ,(/1/),(/maxcat/), parms4file (:, 1)) ; VERIFY_(STATUS) status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'COND' ) ,(/1/),(/maxcat/), parms4file (:, 2)) ; VERIFY_(STATUS) status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'POROS') ,(/1/),(/maxcat/), parms4file (:, 3)) ; VERIFY_(STATUS) @@ -4852,82 +4852,84 @@ END FUNCTION center_pix_int0 ! -------------------------------------------------------------------------------------- - SUBROUTINE process_peatmap (nc, nr, gfiler, pmap) - - implicit none - integer , parameter :: N_lon_pm = 43200, N_lat_pm = 21600 - integer, intent (in) :: nc, nr - real, pointer, dimension (:), intent (inout) :: pmap - character(*), intent (in) :: gfiler - integer :: i,j, status, varid, ncid - integer :: NTILES - REAL, ALLOCATABLE, dimension (:) :: count_pix - REAL, ALLOCATABLE, dimension (:,:) :: data_grid, pm_grid - INTEGER, ALLOCATABLE, dimension (:,:) :: tile_id - character*100 :: fout - - ! Reading number of tiles - ! ----------------------- - - open (20, file = 'clsm/catchment.def', form = 'formatted', status = 'old', action = 'read') - - read (20, *) NTILES - - close (20, status = 'keep') - - ! READ PEATMAP source data files and regrid - ! ----------------------------------------- - - status = NF_OPEN ('data/CATCH/PEATMAP_mask.nc4', NF_NOWRITE, ncid) - - allocate (pm_grid (1 : NC , 1 : NR)) - allocate (data_grid (1 : N_lon_pm, 1 : N_lat_pm)) - - status = NF_INQ_VARID (ncid,'PEATMAP',VarID) ; VERIFY_(STATUS) - status = NF_GET_VARA_REAL (ncid,VarID, (/1,1/),(/N_lon_pm, N_lat_pm/), data_grid) ; VERIFY_(STATUS) - - call RegridRasterReal(data_grid, pm_grid) - - status = NF_CLOSE(ncid) - - ! Grid to tile - ! ------------ - - ! Reading tile-id raster file - - allocate(tile_id(1:nc,1:nr)) - - open (10,file=trim(gfiler)//'.rst',status='old',action='read', & - form='unformatted',convert='little_endian') - - do j=1,nr - read(10)tile_id(:,j) - end do - - close (10,status='keep') - - allocate (pmap (1:NTILES)) - allocate (count_pix (1:NTILES)) - - pmap = 0. - count_pix = 0. - - do j = 1,nr - do i = 1, nc - if((tile_id(i,j).gt.0).and.(tile_id(i,j).le.NTILES)) then - if(pm_grid(i,j) > 0.) pmap (tile_id(i,j)) = pmap (tile_id(i,j)) + pm_grid(i,j) - count_pix (tile_id(i,j)) = count_pix (tile_id(i,j)) + 1. - endif - end do - end do - - where (count_pix > 0.) pmap = pmap/count_pix - - deallocate (count_pix) - deallocate (pm_grid) - deallocate (tile_id) - - END SUBROUTINE process_peatmap +! this subroutine seems obsolete, commented out for now - reichle, 9 Feb 2022 + +! SUBROUTINE process_peatmap (nc, nr, gfiler, pmap) +! +! implicit none +! integer , parameter :: N_lon_pm = 43200, N_lat_pm = 21600 +! integer, intent (in) :: nc, nr +! real, pointer, dimension (:), intent (inout) :: pmap +! character(*), intent (in) :: gfiler +! integer :: i,j, status, varid, ncid +! integer :: NTILES +! REAL, ALLOCATABLE, dimension (:) :: count_pix +! REAL, ALLOCATABLE, dimension (:,:) :: data_grid, pm_grid +! INTEGER, ALLOCATABLE, dimension (:,:) :: tile_id +! character*100 :: fout +! +! ! Reading number of tiles +! ! ----------------------- +! +! open (20, file = 'clsm/catchment.def', form = 'formatted', status = 'old', action = 'read') +! +! read (20, *) NTILES +! +! close (20, status = 'keep') +! +! ! READ PEATMAP source data files and regrid +! ! ----------------------------------------- +! +! status = NF_OPEN ('data/CATCH/PEATMAP_mask.nc4', NF_NOWRITE, ncid) +! +! allocate (pm_grid (1 : NC , 1 : NR)) +! allocate (data_grid (1 : N_lon_pm, 1 : N_lat_pm)) +! +! status = NF_INQ_VARID (ncid,'PEATMAP',VarID) ; VERIFY_(STATUS) +! status = NF_GET_VARA_REAL (ncid,VarID, (/1,1/),(/N_lon_pm, N_lat_pm/), data_grid) ; VERIFY_(STATUS) +! +! call RegridRasterReal(data_grid, pm_grid) +! +! status = NF_CLOSE(ncid) +! +! ! Grid to tile +! ! ------------ +! +! ! Reading tile-id raster file +! +! allocate(tile_id(1:nc,1:nr)) +! +! open (10,file=trim(gfiler)//'.rst',status='old',action='read', & +! form='unformatted',convert='little_endian') +! +! do j=1,nr +! read(10)tile_id(:,j) +! end do +! +! close (10,status='keep') +! +! allocate (pmap (1:NTILES)) +! allocate (count_pix (1:NTILES)) +! +! pmap = 0. +! count_pix = 0. +! +! do j = 1,nr +! do i = 1, nc +! if((tile_id(i,j).gt.0).and.(tile_id(i,j).le.NTILES)) then +! if(pm_grid(i,j) > 0.) pmap (tile_id(i,j)) = pmap (tile_id(i,j)) + pm_grid(i,j) +! count_pix (tile_id(i,j)) = count_pix (tile_id(i,j)) + 1. +! endif +! end do +! end do +! +! where (count_pix > 0.) pmap = pmap/count_pix +! +! deallocate (count_pix) +! deallocate (pm_grid) +! deallocate (tile_id) +! +! END SUBROUTINE process_peatmap ! ==================================================================== @@ -6398,7 +6400,7 @@ SUBROUTINE map_country_codes (NC, NR, gfiler) read (20, *) maxcat - ! READ PEATMAP source data files and regrid + ! READ country code source data files and regrid ! ----------------------------------------- status = NF_OPEN ('data/CATCH/GADM_Country_and_USStates_codes_1km.nc4', NF_NOWRITE, ncid) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index d517bfe14..300ba87e2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -43,8 +43,8 @@ module rmTinyCatchParaMod public :: Get_MidTime, Time_Interp_Fac, compute_stats, c_data public :: ascat_r0, jpl_canoph, NC_VarID, init_bcs_config INTEGER, PARAMETER, public:: SRTM_maxcat = 291284 - logical, public, save :: process_peat = .true. - logical, public, save :: jpl_height = .true. + logical, public, save :: use_PEATMAP = .true. + logical, public, save :: jpl_height = .true. character*8, public, save :: LAIBCS = 'MODGEO' character*4, public, save :: SOILBCS = 'HWSD' character*6, public, save :: MODALB = 'MODIS2' @@ -85,43 +85,43 @@ SUBROUTINE init_bcs_config (LBSV) SOILBCS = 'NGDC' MODALB = 'MODIS1' GNU = 2.17 - process_peat = .false. - jpl_height = .false. + use_PEATMAP = .false. + jpl_height = .false. case ("GM4", "ICA") LAIBCS = 'GSWP2' SOILBCS = 'NGDC' MODALB = 'MODIS2' - process_peat = .false. - jpl_height = .false. + use_PEATMAP = .false. + jpl_height = .false. case ("NL3") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - process_peat = .false. - jpl_height = .false. + use_PEATMAP = .false. + jpl_height = .false. case ("NL4") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - process_peat = .false. - jpl_height = .true. + use_PEATMAP = .false. + jpl_height = .true. case ("NL5") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - process_peat = .true. - jpl_height = .true. + use_PEATMAP = .true. + jpl_height = .true. case ("DEV") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - process_peat = .true. - jpl_height = .true. + use_PEATMAP = .true. + jpl_height = .true. end select @@ -3584,7 +3584,7 @@ SUBROUTINE create_model_para_woesten (Maskfile) !c------------------------------------------------------------------------- - if(process_peat) then + if(use_PEATMAP) then fname = trim(c_data)//'SoilClasses-SoilHyd-TauParam.peatmap' else fname = trim(c_data)//'SoilClasses-SoilHyd-TauParam.dat' @@ -3608,7 +3608,7 @@ SUBROUTINE create_model_para_woesten (Maskfile) write (fout,'(i2.2,i2.2,i4.4)')nint(a_sand(n)),nint(a_clay(n)),nint(100*a_oc(n)) if(n == n_SoilClasses) then - if(process_peat) then + if(use_PEATMAP) then open (120,file=trim(losfile)//trim(fout)//'.peat', & form='formatted',status='old') else @@ -3763,7 +3763,7 @@ SUBROUTINE create_model_para_woesten (Maskfile) write(*,*)'Warnning 1: pfafstetter mismatched' stop endif - if((process_peat).and.(soil_class_top(n) == 253)) then + if((use_PEATMAP).and.(soil_class_top(n) == 253)) then meanlu = 9.3 stdev = 0.12 minlu = 8.5 @@ -3834,7 +3834,7 @@ SUBROUTINE create_model_para_woesten (Maskfile) !$OMP taberr1,taberr2,normerr1,normerr2, & !$OMP taberr3,taberr4,normerr3,normerr4, & !$OMP gwatdep,gwan,grzexcn,gfrc,soil_class_com, & -!$OMP n_threads, low_ind, upp_ind, process_peat )& +!$OMP n_threads, low_ind, upp_ind, use_PEATMAP ) & !$OMP PRIVATE(k,li,ui,n,i,watdep,wan,rzexcn,frc,ST,AC, & !$OMP COESKEW,profdep) @@ -3884,7 +3884,7 @@ SUBROUTINE create_model_para_woesten (Maskfile) if(soil_class_com(n) == 253) then ! Michel Bechtold paper - PEATCLSM_fitting_CLSM_params.R produced these data values. - if(process_peat) then + if(use_PEATMAP) then ars1(n) = -7.9514018e-03 ars2(n) = 6.2297356e-02 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_Catch.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_Catch.F90 index 31ee0586b..f79225031 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_Catch.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_Catch.F90 @@ -2,9 +2,19 @@ #include "MAPL_Generic.h" program Scale_Catch + use MAPL - use lsm_routines, ONLY: catch_calc_tp, catch_calc_ght - USE CATCH_CONSTANTS, ONLY: N_GT=>CATCH_N_GT, DZGT=>CATCH_DZGT + + use LSM_ROUTINES, ONLY: & + catch_calc_soil_moist, & + catch_calc_tp, & + catch_calc_ght + + USE CATCH_CONSTANTS, ONLY: & + N_GT => CATCH_N_GT, & + DZGT => CATCH_DZGT, & + PEATCLSM_POROS_THRESHOLD + implicit none character(256) :: fname1, fname2, fname3 @@ -82,45 +92,9 @@ program Scale_Catch type(catch_rst) catch(3) - interface - subroutine calc_soil_moist( & - ncat,vegcls,dzsf,vgwmax,cdcr1,cdcr2,wpwet,poros, & - psis,bee,ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & - srfexc,rzexc,catdef, & - sfmc, rzmc, prmc, & - werror, sfmcun, rzmcun, prmcun ) - - implicit none - - integer, parameter :: KSNGL=4 - integer, intent(in) :: ncat - integer, dimension(ncat), intent(in) :: vegcls - - real(KIND=KSNGL), dimension(ncat), intent(in) :: dzsf,vgwmax,cdcr1,cdcr2 - real(KIND=KSNGL), dimension(ncat), intent(in) :: wpwet,poros,psis - real(KIND=KSNGL), dimension(ncat), intent(in) :: bee,ars1 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ars2,ars3,ara1,ara2,ara3 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ara4,arw1,arw2,arw3,arw4 - - real(KIND=KSNGL), dimension(ncat), intent(inout) :: srfexc, rzexc, catdef - - real(KIND=KSNGL), dimension(ncat), intent(out) :: sfmc, rzmc, prmc - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: werror - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: sfmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: rzmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: prmcun - end subroutine calc_soil_moist - end interface - - real, allocatable, dimension(:) :: sfmc, rzmc, prmc, werror, sfmcun, rzmcun, prmcun, dzsf - integer, allocatable, dimension(:) :: vegcls - real, allocatable, dimension(:) :: vegdyn + real, allocatable, dimension(:) :: dzsf, ar1, ar2, ar4 real, allocatable, dimension(:,:) :: TP_IN, GHT_IN, FICE, GHT_OUT, TP_OUT - real, allocatable, dimension(:) :: swe_in,depth_in,areasc_in,areasc_out, depth_out - + real, allocatable, dimension(:) :: swe_in, depth_in, areasc_in, areasc_out, depth_out type(Netcdf4_fileformatter) :: formatter(3) type(Filemetadata) :: cfg(3) @@ -213,9 +187,6 @@ end subroutine calc_soil_moist call readcatch ( 20,catch(new) ) end if - allocate( vegcls(ntiles) ) - vegcls(:) = catch(new)%ity(:) - ! Create Scaled Catch ! ------------------- sca = 3 @@ -266,27 +237,26 @@ end subroutine calc_soil_moist catch(sca)%catdef = catch(old)%catdef * (catch(new)%cdcr1 / catch(old)%cdcr1) end where -! Sanity Check +! Sanity Check (catch_calc_soil_moist() forces consistency betw. srfexc, rzexc, catdef) ! ------------ print *, 'Performing Sanity Check ...' allocate ( dzsf(ntiles) ) - allocate ( sfmc(ntiles) ) - allocate ( rzmc(ntiles) ) - allocate ( prmc(ntiles) ) - allocate ( werror(ntiles) ) - allocate ( sfmcun(ntiles) ) - allocate ( rzmcun(ntiles) ) - allocate ( prmcun(ntiles) ) + allocate ( ar1( ntiles) ) + allocate ( ar2( ntiles) ) + allocate ( ar4( ntiles) ) dzsf = SURFLAY - call calc_soil_moist( ntiles,vegcls,dzsf, & - catch(sca)%vgwmax,catch(sca)%cdcr1,catch(sca)%cdcr2,catch(sca)%wpwet,catch(sca)%poros, & - catch(sca)%psis,catch(sca)%bee,catch(sca)%ars1,catch(sca)%ars2,catch(sca)%ars3,catch(sca)%ara1,catch(sca)%ara2, & - catch(sca)%ara3,catch(sca)%ara4,catch(sca)%arw1,catch(sca)%arw2,catch(sca)%arw3,catch(sca)%arw4, & - catch(sca)%srfexc,catch(sca)%rzexc,catch(sca)%catdef, & - sfmc, rzmc, prmc, werror, sfmcun, rzmcun, prmcun ) - + call catch_calc_soil_moist( ntiles, dzsf, & + catch(sca)%vgwmax, catch(sca)%cdcr1, catch(sca)%cdcr2, & + catch(sca)%psis, catch(sca)%bee, catch(sca)%poros, catch(sca)%wpwet, & + catch(sca)%ars1, catch(sca)%ars2, catch(sca)%ars3, & + catch(sca)%ara1, catch(sca)%ara2, catch(sca)%ara3, catch(sca)%ara4, & + catch(sca)%arw1, catch(sca)%arw2, catch(sca)%arw3, catch(sca)%arw4, & + catch(sca)%bf1, catch(sca)%bf2, & + catch(sca)%srfexc, catch(sca)%rzexc, catch(sca)%catdef, & + ar1, ar2, ar4 ) + n = count( catch(sca)%catdef .ne. catch(new)%catdef ) write(6,300) n,100*n/ntiles n = count( catch(sca)%srfexc .ne. catch(new)%srfexc ) @@ -393,6 +363,15 @@ end subroutine calc_soil_moist endif + ! PEATCLSM - ensure low CATDEF on peat tiles where "old" restart is not also peat + ! ------------------------------------------------------------------------------- + + where ( (catch(old)%poros < PEATCLSM_POROS_THRESHOLD) .and. (catch(sca)%poros >= PEATCLSM_POROS_THRESHOLD) ) + catch(sca)%catdef = 25. + catch(sca)%rzexc = 0. + catch(sca)%srfexc = 0. + end where + ! Write Scaled Catch ! ------------------ if (filetype ==0) then @@ -413,9 +392,12 @@ end subroutine calc_soil_moist stop contains - subroutine allocatch (ntiles,catch) - integer ntiles - type(catch_rst) catch + + subroutine allocatch (ntiles,catch) + + integer ntiles + + type(catch_rst) catch allocate( catch% bf1(ntiles) ) allocate( catch% bf2(ntiles) ) @@ -479,11 +461,11 @@ subroutine allocatch (ntiles,catch) end subroutine allocatch subroutine readcatch_nc4 (catch,formatter, rc) - type(catch_rst) catch - type(Netcdf4_fileformatter) :: formatter - integer, optional, intent(out) :: rc - integer :: status - character(256) :: Iam = "readcatch_nc4" + type(catch_rst) catch + type(Netcdf4_fileformatter) :: formatter + integer, optional, intent(out) :: rc + integer :: status + character(256) :: Iam = "readcatch_nc4" call MAPL_VarRead(formatter,"BF1",catch%bf1, __RC__) call MAPL_VarRead(formatter,"BF2",catch%bf2, __RC__) @@ -744,525 +726,3 @@ subroutine writecatch (unit,catch) end subroutine writecatch end program - - subroutine calc_soil_moist( & - ncat,vegcls,dzsf,vgwmax,cdcr1,cdcr2,wpwet,poros, & - psis,bee,ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & - srfexc,rzexc,catdef, & - sfmc, rzmc, prmc, & - werror, sfmcun, rzmcun, prmcun ) - - ! Calculate diagnostic soil moisture content from prognostic - ! excess/deficit variables. - ! - ! On input, also check validity of prognostic excess/deficit variables - ! and modify if necessary. Perturbed or updated excess/deficit variables - ! in data assimilation integrations may be unphysical. - ! Optional output "werror" contains excess or missing water related - ! to inconsistency. - ! - ! Optional outputs "smfcun", "rzmcun", "prmcun" are surface, - ! root zone, and profile moisture content for unsaturated areas only, - ! ie. excluding the saturated area of the catchment. - ! - ! NOTE: When calling with optional output arguments, use keywords - ! unless arguments are in proper order! - ! - ! Example: - ! (don't want "werror" as output, but want "*mcun" output) - ! - ! call calc_soil_moist( & - ! ncat, ... & - ! sfmc, rzmc, prmc, & - ! sfmcun=sfmc_unsat, & - ! rzmcun=rzmc_unsat, & - ! prmcun=prmc_unsat ) - ! - ! replaces moisture_sep_22_2003.f (and older moisture.f) - ! - ! koster+reichle, Feb 5, 2004 - ! - ! revised - koster+reichle, Mar 19, 2004 - ! - ! added optional *un output - koster+reichle, Apr 6, 2004 - ! - ! ---------------------------------------------------------------- - - - implicit none - - integer, parameter :: KSNGL=4 - integer, intent(in) :: ncat - integer, dimension(ncat), intent(in) :: vegcls - - real(KIND=KSNGL), dimension(ncat), intent(in) :: dzsf,vgwmax,cdcr1,cdcr2 - real(KIND=KSNGL), dimension(ncat), intent(in) :: wpwet,poros,psis - real(KIND=KSNGL), dimension(ncat), intent(in) :: bee,ars1 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ars2,ars3,ara1,ara2,ara3 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ara4,arw1,arw2,arw3,arw4 - - real(KIND=KSNGL), dimension(ncat), intent(inout) :: srfexc, rzexc, catdef - - real(KIND=KSNGL), dimension(ncat), intent(out) :: sfmc, rzmc, prmc - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: werror - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: sfmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: rzmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: prmcun - - ! ---------------------------- - ! - ! local variables - - integer :: n - - real(KIND=KSNGL), parameter :: dtstep_dummy = -9999. - - real(KIND=KSNGL), dimension(ncat) :: rzeq, runsrf_dummy, catdef_dummy - real(KIND=KSNGL), dimension(ncat) :: ar1, ar2, ar4, prmc_orig - real(KIND=KSNGL), dimension(ncat) :: srfmn, srfmx, swsrf1, swsrf2, swsrf4, rzi - - - ! -------------------------------------------------------------------- - ! - ! compute soil water storage upon input [mm] - - do n=1,ncat - prmc_orig(n) = & - (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) - enddo - - ! ----------------------------------- - ! - ! check limits of catchment deficit - ! - ! increased minimum catchment deficit from 0.01 to 1. to make the - ! check work with perturbed parameters and initial condition - ! reichle, 16 May 01 - ! - ! IT REALLY SHOULD WORK WITH catdef > 0 (rather than >1.) ???? - ! reichle, 5 Feb 2004 - - do n=1,ncat - catdef(n)=max(1.,min(cdcr2(n),catdef(n))) - end do - - ! ------------------------------------------------------------------ - ! - ! check limits of root zone excess - ! - ! calculate root zone equilibrium moisture for given catchment deficit - - call rzequil( & - ncat, vegcls, catdef, vgwmax, & - cdcr1, cdcr2, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, & - arw1, arw2, arw3, arw4, & - rzeq) - - ! assume srfexc=0 and constrain rzexc appropriately - ! (iteration would be needed to contrain srfexc and rzexc simultaneously) - - do n=1,ncat - rzexc(n)=max(wpwet(n)*vgwmax(n)-rzeq(n),min(vgwmax(n)-rzeq(n),rzexc(n))) - end do - - ! this translates into: - ! - ! wilting level < rzmc < porosity - ! - ! or more precisely: wpwet*vgwmax < rzeq+rzexc < vgwmax - ! - ! NOTE: root zone moisture is not allowed to drop below wilting level - - ! ----------------------------------------------------------------- - ! - ! Call partition() for computation of surface moisture content. - ! - ! Call to partition() also checks limits of surface excess. - ! - ! Call partition with dtstep_dummy: - ! In partition, dtstep is only used for a correction that - ! puts water into runsrf (for which runsrf_dummy is used here). - ! Also use catdef_dummy because partition() updates catdef - ! whenever srfexc exceeds physical bounds, but this is not desired here. - - runsrf_dummy = 0. - catdef_dummy = catdef - - call partition( & - ncat,dtstep_dummy,vegcls,dzsf,rzexc, & - rzeq,vgwmax,cdcr1,cdcr2, & - psis,bee,poros,wpwet, & - ars1,ars2,ars3, & - ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4,.false., & - srfexc,catdef_dummy,runsrf_dummy, & - ar1, ar2, ar4,srfmx,srfmn, & - swsrf1,swsrf2,swsrf4,rzi & - ) - - ! compute surface, root zone, and profile soil moisture - - do n=1,ncat - - sfmc(n) = poros(n) * & - (swsrf1(n)*ar1(n) + swsrf2(n)*ar2(n) + swsrf4(n)*ar4(n)) - - rzmc(n) = (rzeq(n)+rzexc(n)+srfexc(n))*poros(n)/vgwmax(n) - - ! compute revised soil water storage [mm] - - prmc(n) = & - (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) - - ! compute error in soil water storage [mm] (if argument is present) - - if (present(werror)) werror(n)=(prmc(n)-prmc_orig(n)) - - ! convert to volumetric soil moisture - ! note: dzpr = (cdcr2/(1-wpwet)) / poros - - prmc(n) = prmc(n)*poros(n) / (cdcr2(n)/(1.-wpwet(n))) - - - ! check for negative soil moisture - - if ( (sfmc(n)<.0) .or. (rzmc(n)<.0) .or. (prmc(n)<.0) ) then - - write (*,*) 'FOUND NEGATIVE SOIL MOISTURE CONTENT.... stopping' - write (*,*) n, sfmc(n), rzmc(n), prmc(n) - stop - end if - - ! compute moisture content in unsaturated areas [m3/m3] (if arg present) - - if (ar1(n)<1.) then - - if (present(prmcun)) prmcun(n)=(prmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - if (present(rzmcun)) rzmcun(n)=(rzmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - if (present(sfmcun)) sfmcun(n)=(sfmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - - else - - if (present(prmcun)) prmcun(n)=poros(n) - if (present(rzmcun)) rzmcun(n)=poros(n) - if (present(sfmcun)) sfmcun(n)=poros(n) - - end if - - enddo - - return - - end subroutine calc_soil_moist - - SUBROUTINE PARTITION ( & - NCH,DTSTEP,ITYP,DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2,& - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4,BUG, & - srfexc,catdef,runsrf, & - AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1,SWSRF2,SWSRF4,RZI & - ) - - IMPLICIT NONE - -! ------------------------------------------------------------------- - INTEGER, INTENT(IN) :: NCH - INTEGER, INTENT(IN), DIMENSION(NCH) :: ITYP - - REAL, INTENT(IN) :: DTSTEP - REAL, INTENT(IN), DIMENSION(NCH) :: DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4 - - LOGICAL, INTENT(IN) :: BUG -! ------------------------------------------------------------------- - REAL, INTENT(INOUT), DIMENSION(NCH) :: srfexc,catdef,runsrf -! ------------------------------------------------------------------- - REAL, INTENT(OUT), DIMENSION(NCH) :: AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1, SWSRF2, SWSRF4, RZI -! ------------------------------------------------------------------- - INTEGER :: N - - REAL :: cor, A150, W150, WMIN, AX, WMNEW, WRZ, TERM1, TERM2, TERM3, & - AREA0, AREA1, AREA2, AREA3, AREA4, ASCALE, WILT, D1, D2, CDI, & - DELTA1, DELTA2, DELTA4, MULTAR, CATDEFX, RZEQX, RZEQW, FACTOR, & - X0, RZEQY, CATDEFW, AR1W, ASUM, RZEQYI, RZEQWI, RZEQXI, AR20, & - ARG1, EXPARG1, ARG2, EXPARG2, ARG3, EXPARG3 !, surflay - - LOGICAL :: LSTRESS - - - DATA LSTRESS/.FALSE./ !,surflay/20./ - -!**** -!**** -------------------------------------------------- - -!rr next line for debugging, sep 23, 2003, reichle -!rr -!rr write (*,*) 'entering partition()' - - DO N=1,NCH - - WILT=WPWET(N) - WRZ=RZEXC(N)/VGWMAX(N) - CATDEFX=AMIN1( CATDEF(N) , CDCR1(N) ) - -! CDI DEFINES IF THE SHAPE PARAMETER IS ADJUSTED IN ONE OR TWO SEGMENTS - if (ara1(n) .ne. ara3(n)) then - cdi=(ara4(n)-ara2(n))/(ara1(n)-ara3(n)) - else - cdi=0. - endif - - AR1(N)= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFX) & - /(1.+ars2(n)*CATDEFX+ars3(n)*CATDEFX*CATDEFX))) - - if (CATDEFX .ge. cdi) then - ax=ara3(n)*CATDEFX+ara4(n) - else - ax=ara1(n)*CATDEFX+ara2(n) - endif - - WMIN=AMIN1(1.,AMAX1(0.,arw4(n)+(1.-arw4(n))* & - (1.+arw1(n)*CATDEFX) & - /(1.+arw2(n)*CATDEFX+arw3(n)*CATDEFX*CATDEFX))) - -!**** CRITICAL VALUE 1: AVERAGE MOISTURE IN ROOT ZONE AT WMIN -!**** ASSOCIATED WITH CATDEF. - ARG1=AMAX1(-40., AMIN1(40., -AX*(1.-WMIN))) - EXPARG1=EXP(ARG1) - RZEQX=(WMIN-1.-(2./AX))*EXPARG1 + WMIN + (2./AX) - RZEQXI=AX*EXPARG1 * & - ( -1. -(2./AX) - (2./(AX*AX)) + WMIN + (WMIN/AX) ) & - + WMIN + 2./AX - AR20=1.+(-AX-1.+AX*WMIN)*EXPARG1 - RZEQXI=RZEQXI/(AR20+1.E-20) - -!**** CRITICAL VALUE 2: AVERAGE MOISTURE IN ROOT ZONE WHEN WMIN -!**** IS EXACTLY AT WILTING POINT. - ARG2=AMAX1(-40., AMIN1(40., -AX*(1.-WILT))) - EXPARG2=EXP(ARG2) - RZEQW=(WILT-1.-(2./AX))*EXPARG2 + WILT + (2./AX) - RZEQWI=AX*EXPARG2 * & - ( -1. -(2./AX) - (2./(AX*AX)) + WILT + (WILT/AX) ) & - + WILT + 2./AX - AR20=1.+(-AX-1.+AX*WILT)*EXPARG2 - RZEQWI=RZEQWI/(AR20+1.E-20) - -!**** SITUATION 1: CATDEF LE CDCR1 - IF(CATDEF(N) .LE. CDCR1(N)) THEN - RZEQY=RZEQX+WRZ - RZEQYI=RZEQXI+WRZ - WMNEW=WMIN+WRZ - ARG3=AMAX1(-40., AMIN1(40., -AX*(1.-WMNEW))) - EXPARG3=EXP(ARG3) - AREA1=(1.+AX-AX*WMIN)*EXPARG1 - AREA2=(1.+AX-AX*WMNEW)*EXPARG3 - IF(WMNEW .GE. WILT) THEN - AR1(N)=AR1(N)+AREA2-AREA1 - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - IF(WMNEW .LT. WILT) THEN - AREA3=(1.+AX-AX*WILT)*EXPARG2 - AR1(N)=AR1(N)+AREA3-AREA1 - AR2(N)=1.-AR1(N) - FACTOR=(RZEQX+WRZ-WILT)/(RZEQW-WILT) - AR1(N)=AR1(N)*FACTOR - AR2(N)=AR2(N)*FACTOR - AR4(N)=1.-FACTOR - ENDIF - ENDIF - -!**** SITUATION 2: CATDEF GT CDCR1 - IF(CATDEF(N) .GT. CDCR1(N)) THEN - FACTOR=(CDCR2(N)-CATDEF(N))/(CDCR2(N)-CDCR1(N)) - RZEQY=WILT+(RZEQX-WILT)*FACTOR+WRZ - RZEQYI=WILT+(RZEQXI-WILT)*FACTOR+WRZ - - IF(RZEQY .LT. WILT) THEN - IF(RZEQY .LT. WILT-.001) THEN -!rr WRITE(*,*) 'RZEXC WAY TOO LOW! N=',N,' RZEQY=',RZEQY -!rr WRITE(*,*) 'SRFEXC=',SRFEXC(N),'RZEXC=',RZEXC(N), -!rr & 'CATDEF=',CATDEF(N) -! ELSE -! WRITE(*,*) 'RZEXC TOO LOW N=',N - ENDIF - RZEQY=WILT - RZEQYI=WILT - ENDIF - - IF(RZEQY .GE. RZEQX) THEN ! RZEXC BRINGS MOISTURE ABOVE CDCR1 POINT - WMNEW=WMIN+(RZEQY-RZEQX) - ARG3=AMAX1(-40., AMIN1(40., -AX*(1.-WMNEW))) - EXPARG3=EXP(ARG3) - AREA1=(1.+AX-AX*WMIN)*EXPARG1 - AREA2=(1.+AX-AX*WMNEW)*EXPARG3 - AR1(N)=AR1(N)+(AREA2-AREA1) - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - - IF(RZEQY .LT. RZEQX .AND. RZEQY .GE. RZEQW) THEN - CATDEFW=CDCR2(N)+((RZEQW-WILT)/(RZEQX-WILT))*(CDCR1(N)-CDCR2(N)) - AR1W= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFW) & - /(1.+ars2(n)*CATDEFW+ars3(n)*CATDEFW*CATDEFW))) - FACTOR=(RZEQY-RZEQW)/(RZEQX-RZEQW) - AR1(N)=AR1W+FACTOR*(AR1(N)-AR1W) - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - - IF(RZEQY .LT. RZEQW) THEN - CATDEFW=CDCR2(N)+((RZEQW-WILT)/(RZEQX-WILT))*(CDCR1(N)-CDCR2(N)) - AR1W= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFW) & - /(1.+ars2(n)*CATDEFW+ars3(n)*CATDEFW*CATDEFW))) - AR1(N)=AR1W - AR2(N)=1.-AR1(N) - FACTOR=(RZEQY-WILT)/(RZEQW-WILT) - AR1(N)=AR1(N)*FACTOR - AR2(N)=AR2(N)*FACTOR - AR4(N)=1.-FACTOR - ENDIF - - ENDIF - - RZI(N)=RZEQYI - - SWSRF1(N)=1. -!mjs: changed .001 temporarily because of large bee. - SWSRF2(N)=AMIN1(1., AMAX1(0.01, RZEQYI)) - SWSRF4(N)=AMIN1(1., AMAX1(0.01, WILT)) - -!**** EXTRAPOLATION OF THE SURFACE WETNESSES - -! 1st step: surface wetness in the unstressed fraction without considering -! the surface excess; we just assume an equilibrium profile from -! the middle of the root zone to the surface. - - SWSRF2(N)=((SWSRF2(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) - SWSRF4(N)=((SWSRF4(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) - -! srfmx is the maximum amount of water that can be added to the surface layer -! The choice of defining SWSRF4 like SWSRF2 needs to be better examined. - srfmx(n)=ar2(n)*(1.-swsrf2(n))*(dzsf(n)*poros(n)) - srfmx(n)=srfmx(n)+ar4(n)*(1.-swsrf4(n))*(dzsf(n)*poros(n)) -!**** For calculation of srfmn, assume surface moisture associated with -!**** AR1 is constantly replenished by water table. - srfmn(n)=-(ar2(n)*swsrf2(n)+ar4(n)*swsrf4(n))*(dzsf(n)*poros(n)) - - if(srfexc(n).gt.srfmx(n)) then - cor=srfexc(n)-srfmx(n) ! The correction is here - srfexc(n)=srfmx(n) - catdef(n)=catdef(n)-cor - if(catdef(n).lt.0.) then - runsrf(n)=runsrf(n)-catdef(n)/dtstep - catdef(n)=0. - endif - else if(srfexc(n).lt.srfmn(n)) then - cor=srfexc(n)-srfmn(n) - catdef(n)=catdef(n)-cor - srfexc(n)=srfmn(n) - else - cor=0. - endif - - SWSRF2(N)=SWSRF2(N)+SRFEXC(N)/(dzsf(n)*poros(n)*(1.-ar1(n))+1.e-20) - SWSRF2(N)=AMIN1(1., AMAX1(1.E-5, SWSRF2(N))) - swsrf4(n)=swsrf4(n)+srfexc(n)/(dzsf(n)*poros(n)*(1.-ar1(n))+1.e-20) - SWSRF4(N)=AMIN1(1., AMAX1(1.E-5, SWSRF4(N))) - - IF (AR1(N) .ge. 1.-1.E-5) then - AR1(N)=1. - AR2(N)=0. - AR4(N)=0. - SWSRF2(N)=1. - SWSRF4(N)=wilt - ENDIF - - IF (AR1(N) .LT. 0.) then -!rr IF(AR1(N) .LT. -1.E-3) WRITE(*,*) 'AR1 TOO LOW: AR1=',AR1(N) - AR1(N)=0. - ENDIF - ar1(n)=amax1(0., amin1(1., ar1(n))) - ar2(n)=amax1(0., amin1(1., ar2(n))) - ar4(n)=amax1(0., amin1(1., ar4(n))) - asum=ar1(n)+ar2(n)+ar4(n) - if(asum .lt. .9999 .or. asum .gt. 1.0001) then - write(*,*) 'Areas do not add to 1: sum=',asum,'N=',n - endif - - - ENDDO - - - RETURN - END SUBROUTINE PARTITION - - SUBROUTINE RZEQUIL ( & - NCH,ITYP,CATDEF,VGWMAX,CDCR1,CDCR2,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4, & - RZEQ & - ) - - IMPLICIT NONE - - INTEGER, INTENT(IN) :: NCH - INTEGER, INTENT(IN), DIMENSION(NCH) :: ITYP - REAL, INTENT(IN), DIMENSION(NCH) :: CATDEF, VGWMAX, CDCR1, CDCR2, & - WPWET, ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, & - arw2, arw3, arw4 - - REAL, INTENT(OUT), DIMENSION(NCH) :: RZEQ - - INTEGER N - REAL AX,WMIN,ASCALE,cdi,wilt,catdefx,factor,ARG1,EXPARG1 - -! ---------------------------------------------------------------------- - - DO N=1,NCH - - WILT=WPWET(N) - CATDEFX=AMIN1( CATDEF(N) , CDCR1(N) ) - -! CDI DEFINES IF THE SHAPE PARAMETER IS ADJUSTED IN ONE OR TWO SEGMENTS - if (ara1(n) .ne. ara3(n)) then - cdi=(ara4(n)-ara2(n))/(ara1(n)-ara3(n)) - else - cdi=0. - endif - - if (CATDEFX .ge. cdi) then - ax=ara3(n)*CATDEFX+ara4(n) - else - ax=ara1(n)*CATDEFX+ara2(n) - endif - - WMIN=AMIN1(1.,AMAX1(0.,arw4(n)+(1.-arw4(n))*(1.+arw1(n)*CATDEFX) & - /(1.+arw2(n)*CATDEFX+arw3(n)*CATDEFX*CATDEFX))) - - ARG1=AMAX1(-40., AMIN1(40., -AX*(1.-WMIN))) - EXPARG1=EXP(ARG1) - RZEQ(N)=(WMIN-1.-(2./AX))*EXPARG1 + WMIN + (2./AX) - - IF(CATDEF(N) .GT. CDCR1(N)) THEN - FACTOR=(CDCR2(N)-CATDEF(N))/(CDCR2(N)-CDCR1(N)) - RZEQ(N)=WILT+(RZEQ(N)-WILT)*FACTOR - ENDIF - -! scaling: - RZEQ(N)=AMIN1(1.,AMAX1(0.,RZEQ(N))) - RZEQ(N)=RZEQ(N)*VGWMAX(N) - - ENDDO - - RETURN - END SUBROUTINE RZEQUIL diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_CatchCN.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_CatchCN.F90 index 27e204979..cd2bce354 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_CatchCN.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/Scale_CatchCN.F90 @@ -1,9 +1,20 @@ #define I_AM_MAIN #include "MAPL_Generic.h" + program Scale_CatchCN + use MAPL - use lsm_routines, ONLY: catch_calc_tp, catch_calc_ght - USE CATCH_CONSTANTS, ONLY: N_GT=>CATCH_N_GT, DZGT=>CATCH_DZGT + + use LSM_ROUTINES, ONLY: & + catch_calc_soil_moist, & + catch_calc_tp, & + catch_calc_ght + + USE CATCH_CONSTANTS, ONLY: & + N_GT => CATCH_N_GT, & + DZGT => CATCH_DZGT, & + PEATCLSM_POROS_THRESHOLD + implicit none character(256) :: fname1, fname2, fname3 @@ -16,9 +27,10 @@ program Scale_CatchCN integer :: iargc real :: SURFLAY ! (Ganymed-3 and earlier) SURFLAY=20.0 for Old Soil Params ! (Ganymed-4 and later ) SURFLAY=50.0 for New Soil Params + real :: WEMIN_IN, WEMIN_OUT character*256 :: arg(6) - integer, parameter :: nveg = 4 + integer, parameter :: nveg = 4 integer, parameter :: nzone = 3 integer :: VAR_COL, VAR_PFT integer, parameter :: VAR_COL_CLM40 = 40 ! number of CN column restart variables @@ -27,7 +39,6 @@ program Scale_CatchCN integer, parameter :: VAR_COL_CLM45 = 35 ! number of CN column restart variables integer, parameter :: VAR_PFT_CLM45 = 75 ! number of CN PFT variables per column - real :: WEMIN_IN, WEMIN_OUT logical :: clm45 = .false. integer :: un_dim3 @@ -109,47 +120,16 @@ program Scale_CatchCN type(catch_rst) catch(3) - interface - subroutine calc_soil_moist( & - ncat,dzsf,vgwmax,cdcr1,cdcr2,wpwet,poros, & - psis,bee,ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & - srfexc,rzexc,catdef, & - sfmc, rzmc, prmc, & - werror, sfmcun, rzmcun, prmcun ) - - implicit none - - integer, parameter :: KSNGL=4 - integer, intent(in) :: ncat - - real(KIND=KSNGL), dimension(ncat), intent(in) :: dzsf,vgwmax,cdcr1,cdcr2 - real(KIND=KSNGL), dimension(ncat), intent(in) :: wpwet,poros,psis - real(KIND=KSNGL), dimension(ncat), intent(in) :: bee,ars1 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ars2,ars3,ara1,ara2,ara3 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ara4,arw1,arw2,arw3,arw4 - - real(KIND=KSNGL), dimension(ncat), intent(inout) :: srfexc, rzexc, catdef - - real(KIND=KSNGL), dimension(ncat), intent(out) :: sfmc, rzmc, prmc - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: werror - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: sfmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: rzmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: prmcun - end subroutine calc_soil_moist - end interface - - real, allocatable, dimension(:) :: sfmc, rzmc, prmc, werror, sfmcun, rzmcun, prmcun, dzsf - real, allocatable, dimension(:) :: vegdyn + real, allocatable, dimension(:) :: dzsf, ar1, ar2, ar4 real, allocatable, dimension(:,:) :: TP_IN, GHT_IN, FICE, GHT_OUT, TP_OUT - real, allocatable, dimension(:) :: swe_in,depth_in,areasc_in,areasc_out, depth_out + real, allocatable, dimension(:) :: swe_in, depth_in, areasc_in, areasc_out, depth_out type(Netcdf4_fileformatter) :: formatter(3) type(Filemetadata) :: cfg(3) integer :: i, rc, filetype integer :: status + character(256) :: Iam = "Scale_CatchCN" + ! Usage ! ----- if (iargc() /= 6) then @@ -171,13 +151,13 @@ end subroutine calc_soil_moist ! ------------------------------- read(arg(3),'(a)') fname3 - call MAPL_NCIOGetFileType(fname1, filetype,rc=rc) + call MAPL_NCIOGetFileType(fname1, filetype, __RC__) if (filetype == 0) then - call formatter(1)%open(trim(fname1),pFIO_READ,rc=rc) - call formatter(2)%open(trim(fname2),pFIO_READ,rc=rc) - cfg(1)=formatter(1)%read(rc=rc) - cfg(2)=formatter(2)%read(rc=rc) + call formatter(1)%open(trim(fname1),pFIO_READ, __RC__) + call formatter(2)%open(trim(fname2),pFIO_READ, __RC__) + cfg(1)=formatter(1)%read(__RC__) + cfg(2)=formatter(2)%read(__RC__) ! else ! open(unit=10, file=trim(fname1), form='unformatted') ! open(unit=20, file=trim(fname2), form='unformatted') @@ -203,8 +183,8 @@ end subroutine calc_soil_moist if (filetype ==0) then - ntiles = cfg(1)%get_dimension('tile',rc=rc) - un_dim3 = cfg(1)%get_dimension('unknown_dim3',rc=rc) + ntiles = cfg(1)%get_dimension('tile', __RC__) + un_dim3 = cfg(1)%get_dimension('unknown_dim3', __RC__) if(un_dim3 == 105) then clm45 = .true. VAR_COL = VAR_COL_CLM45 @@ -249,6 +229,7 @@ end subroutine calc_soil_moist ! Create Scaled Catch ! ------------------- sca = 3 + catch(sca) = catch(new) ! 1) soil moisture prognostics @@ -295,27 +276,26 @@ end subroutine calc_soil_moist catch(sca)%catdef = catch(old)%catdef * (catch(new)%cdcr1 / catch(old)%cdcr1) end where -! Sanity Check +! Sanity Check (catch_calc_soil_moist() forces consistency betw. srfexc, rzexc, catdef) ! ------------ print *, 'Performing Sanity Check ...' allocate ( dzsf(ntiles) ) - allocate ( sfmc(ntiles) ) - allocate ( rzmc(ntiles) ) - allocate ( prmc(ntiles) ) - allocate ( werror(ntiles) ) - allocate ( sfmcun(ntiles) ) - allocate ( rzmcun(ntiles) ) - allocate ( prmcun(ntiles) ) + allocate ( ar1( ntiles) ) + allocate ( ar2( ntiles) ) + allocate ( ar4( ntiles) ) dzsf = SURFLAY - call calc_soil_moist( ntiles,dzsf, & - catch(sca)%vgwmax,catch(sca)%cdcr1,catch(sca)%cdcr2,catch(sca)%wpwet,catch(sca)%poros, & - catch(sca)%psis,catch(sca)%bee,catch(sca)%ars1,catch(sca)%ars2,catch(sca)%ars3,catch(sca)%ara1,catch(sca)%ara2, & - catch(sca)%ara3,catch(sca)%ara4,catch(sca)%arw1,catch(sca)%arw2,catch(sca)%arw3,catch(sca)%arw4, & - catch(sca)%srfexc,catch(sca)%rzexc,catch(sca)%catdef, & - sfmc, rzmc, prmc, werror, sfmcun, rzmcun, prmcun ) - + call catch_calc_soil_moist( ntiles, dzsf, & + catch(sca)%vgwmax, catch(sca)%cdcr1, catch(sca)%cdcr2, & + catch(sca)%psis, catch(sca)%bee, catch(sca)%poros, catch(sca)%wpwet, & + catch(sca)%ars1, catch(sca)%ars2, catch(sca)%ars3, & + catch(sca)%ara1, catch(sca)%ara2, catch(sca)%ara3, catch(sca)%ara4, & + catch(sca)%arw1, catch(sca)%arw2, catch(sca)%arw3, catch(sca)%arw4, & + catch(sca)%bf1, catch(sca)%bf2, & + catch(sca)%srfexc, catch(sca)%rzexc, catch(sca)%catdef, & + ar1, ar2, ar4 ) + n = count( catch(sca)%catdef .ne. catch(new)%catdef ) write(6,300) n,100*n/ntiles n = count( catch(sca)%srfexc .ne. catch(new)%srfexc ) @@ -422,12 +402,21 @@ end subroutine calc_soil_moist endif + ! PEATCLSM - ensure low CATDEF on peat tiles where "old" restart is not also peat + ! ------------------------------------------------------------------------------- + + where ( (catch(old)%poros < PEATCLSM_POROS_THRESHOLD) .and. (catch(sca)%poros >= PEATCLSM_POROS_THRESHOLD) ) + catch(sca)%catdef = 25. + catch(sca)%rzexc = 0. + catch(sca)%srfexc = 0. + end where + ! Write Scaled Catch ! ------------------ if (filetype ==0) then - cfg(3) = cfg(2) - call formatter(3)%create(fname3,rc=rc) - call formatter(3)%write(cfg(3),rc=rc) + cfg(3)=cfg(2) + call formatter(3)%create(fname3, __RC__) + call formatter(3)%write(cfg(3), __RC__) call writecatchcn_nc4 ( catch(sca), formatter(3) ,cfg(3) ) ! else ! call writecatchcn ( 30,catch(sca) ) @@ -442,11 +431,12 @@ end subroutine calc_soil_moist stop contains - subroutine allocatch (ntiles,catch) - - integer ntiles - type(catch_rst) catch + subroutine allocatch (ntiles,catch) + + integer ntiles + + type(catch_rst) catch allocate( catch% bf1(ntiles) ) allocate( catch% bf2(ntiles) ) @@ -528,13 +518,13 @@ end subroutine allocatch subroutine readcatchcn_nc4 (catch,formatter,cfg, rc) type(catch_rst) catch type(Filemetadata) :: cfg - type(Netcdf4_Fileformatter) :: formatter + type(Netcdf4_fileformatter) :: formatter integer, optional, intent(out) :: rc integer :: j, dim1,dim2 type(Variable), pointer :: myVariable character(len=:), pointer :: dname - character(256) :: Iam = "readcatchcn_nc4" integer :: status + character(256) :: Iam = "readcatchcn_nc4" call MAPL_VarRead(formatter,"BF1",catch%bf1, __RC__) call MAPL_VarRead(formatter,"BF2",catch%bf2, __RC__) @@ -970,522 +960,3 @@ end subroutine writecatchcn end program - subroutine calc_soil_moist( & - ncat,dzsf,vgwmax,cdcr1,cdcr2,wpwet,poros, & - psis,bee,ars1,ars2,ars3,ara1,ara2, & - ara3,ara4,arw1,arw2,arw3,arw4, & - srfexc,rzexc,catdef, & - sfmc, rzmc, prmc, & - werror, sfmcun, rzmcun, prmcun ) - - ! Calculate diagnostic soil moisture content from prognostic - ! excess/deficit variables. - ! - ! On input, also check validity of prognostic excess/deficit variables - ! and modify if necessary. Perturbed or updated excess/deficit variables - ! in data assimilation integrations may be unphysical. - ! Optional output "werror" contains excess or missing water related - ! to inconsistency. - ! - ! Optional outputs "smfcun", "rzmcun", "prmcun" are surface, - ! root zone, and profile moisture content for unsaturated areas only, - ! ie. excluding the saturated area of the catchment. - ! - ! NOTE: When calling with optional output arguments, use keywords - ! unless arguments are in proper order! - ! - ! Example: - ! (don't want "werror" as output, but want "*mcun" output) - ! - ! call calc_soil_moist( & - ! ncat, ... & - ! sfmc, rzmc, prmc, & - ! sfmcun=sfmc_unsat, & - ! rzmcun=rzmc_unsat, & - ! prmcun=prmc_unsat ) - ! - ! replaces moisture_sep_22_2003.f (and older moisture.f) - ! - ! koster+reichle, Feb 5, 2004 - ! - ! revised - koster+reichle, Mar 19, 2004 - ! - ! added optional *un output - koster+reichle, Apr 6, 2004 - ! - ! ---------------------------------------------------------------- - - - implicit none - - integer, parameter :: KSNGL=4 - integer, intent(in) :: ncat - - real(KIND=KSNGL), dimension(ncat), intent(in) :: dzsf,vgwmax,cdcr1,cdcr2 - real(KIND=KSNGL), dimension(ncat), intent(in) :: wpwet,poros,psis - real(KIND=KSNGL), dimension(ncat), intent(in) :: bee,ars1 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ars2,ars3,ara1,ara2,ara3 - real(KIND=KSNGL), dimension(ncat), intent(in) :: ara4,arw1,arw2,arw3,arw4 - - real(KIND=KSNGL), dimension(ncat), intent(inout) :: srfexc, rzexc, catdef - - real(KIND=KSNGL), dimension(ncat), intent(out) :: sfmc, rzmc, prmc - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: werror - - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: sfmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: rzmcun - real(KIND=KSNGL), dimension(ncat), intent(out), optional :: prmcun - - ! ---------------------------- - ! - ! local variables - - integer :: n - - real(KIND=KSNGL), parameter :: dtstep_dummy = -9999. - - real(KIND=KSNGL), dimension(ncat) :: rzeq, runsrf_dummy, catdef_dummy - real(KIND=KSNGL), dimension(ncat) :: ar1, ar2, ar4, prmc_orig - real(KIND=KSNGL), dimension(ncat) :: srfmn, srfmx, swsrf1, swsrf2, swsrf4, rzi - - - ! -------------------------------------------------------------------- - ! - ! compute soil water storage upon input [mm] - - do n=1,ncat - prmc_orig(n) = & - (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) - enddo - - ! ----------------------------------- - ! - ! check limits of catchment deficit - ! - ! increased minimum catchment deficit from 0.01 to 1. to make the - ! check work with perturbed parameters and initial condition - ! reichle, 16 May 01 - ! - ! IT REALLY SHOULD WORK WITH catdef > 0 (rather than >1.) ???? - ! reichle, 5 Feb 2004 - - do n=1,ncat - catdef(n)=max(1.,min(cdcr2(n),catdef(n))) - end do - - ! ------------------------------------------------------------------ - ! - ! check limits of root zone excess - ! - ! calculate root zone equilibrium moisture for given catchment deficit - - call rzequil( & - ncat, catdef, vgwmax, & - cdcr1, cdcr2, wpwet, & - ars1, ars2, ars3, ara1, ara2, ara3, ara4, & - arw1, arw2, arw3, arw4, & - rzeq) - - ! assume srfexc=0 and constrain rzexc appropriately - ! (iteration would be needed to contrain srfexc and rzexc simultaneously) - - do n=1,ncat - rzexc(n)=max(wpwet(n)*vgwmax(n)-rzeq(n),min(vgwmax(n)-rzeq(n),rzexc(n))) - end do - - ! this translates into: - ! - ! wilting level < rzmc < porosity - ! - ! or more precisely: wpwet*vgwmax < rzeq+rzexc < vgwmax - ! - ! NOTE: root zone moisture is not allowed to drop below wilting level - - ! ----------------------------------------------------------------- - ! - ! Call partition() for computation of surface moisture content. - ! - ! Call to partition() also checks limits of surface excess. - ! - ! Call partition with dtstep_dummy: - ! In partition, dtstep is only used for a correction that - ! puts water into runsrf (for which runsrf_dummy is used here). - ! Also use catdef_dummy because partition() updates catdef - ! whenever srfexc exceeds physical bounds, but this is not desired here. - - runsrf_dummy = 0. - catdef_dummy = catdef - - call partition( & - ncat,dtstep_dummy,dzsf,rzexc, & - rzeq,vgwmax,cdcr1,cdcr2, & - psis,bee,poros,wpwet, & - ars1,ars2,ars3, & - ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4,.false., & - srfexc,catdef_dummy,runsrf_dummy, & - ar1, ar2, ar4,srfmx,srfmn, & - swsrf1,swsrf2,swsrf4,rzi & - ) - - ! compute surface, root zone, and profile soil moisture - - do n=1,ncat - - sfmc(n) = poros(n) * & - (swsrf1(n)*ar1(n) + swsrf2(n)*ar2(n) + swsrf4(n)*ar4(n)) - - rzmc(n) = (rzeq(n)+rzexc(n)+srfexc(n))*poros(n)/vgwmax(n) - - ! compute revised soil water storage [mm] - - prmc(n) = & - (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) - - ! compute error in soil water storage [mm] (if argument is present) - - if (present(werror)) werror(n)=(prmc(n)-prmc_orig(n)) - - ! convert to volumetric soil moisture - ! note: dzpr = (cdcr2/(1-wpwet)) / poros - - prmc(n) = prmc(n)*poros(n) / (cdcr2(n)/(1.-wpwet(n))) - - - ! check for negative soil moisture - - if ( (sfmc(n)<.0) .or. (rzmc(n)<.0) .or. (prmc(n)<.0) ) then - - write (*,*) 'FOUND NEGATIVE SOIL MOISTURE CONTENT.... stopping' - write (*,*) n, sfmc(n), rzmc(n), prmc(n) - stop - end if - - ! compute moisture content in unsaturated areas [m3/m3] (if arg present) - - if (ar1(n)<1.) then - - if (present(prmcun)) prmcun(n)=(prmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - if (present(rzmcun)) rzmcun(n)=(rzmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - if (present(sfmcun)) sfmcun(n)=(sfmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) - - else - - if (present(prmcun)) prmcun(n)=poros(n) - if (present(rzmcun)) rzmcun(n)=poros(n) - if (present(sfmcun)) sfmcun(n)=poros(n) - - end if - - enddo - - return - - end subroutine calc_soil_moist - - SUBROUTINE PARTITION ( & - NCH,DTSTEP,DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4,BUG, & - srfexc,catdef,runsrf, & - AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1,SWSRF2,SWSRF4,RZI & - ) - - IMPLICIT NONE - -! ------------------------------------------------------------------- - INTEGER, INTENT(IN) :: NCH - - REAL, INTENT(IN) :: DTSTEP - REAL, INTENT(IN), DIMENSION(NCH) :: DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4 - - LOGICAL, INTENT(IN) :: BUG -! ------------------------------------------------------------------- - REAL, INTENT(INOUT), DIMENSION(NCH) :: srfexc,catdef,runsrf -! ------------------------------------------------------------------- - REAL, INTENT(OUT), DIMENSION(NCH) :: AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1, SWSRF2, SWSRF4, RZI -! ------------------------------------------------------------------- - INTEGER :: N - - REAL :: cor, A150, W150, WMIN, AX, WMNEW, WRZ, TERM1, TERM2, TERM3, & - AREA0, AREA1, AREA2, AREA3, AREA4, ASCALE, WILT, D1, D2, CDI, & - DELTA1, DELTA2, DELTA4, MULTAR, CATDEFX, RZEQX, RZEQW, FACTOR, & - X0, RZEQY, CATDEFW, AR1W, ASUM, RZEQYI, RZEQWI, RZEQXI, AR20, & - ARG1, EXPARG1, ARG2, EXPARG2, ARG3, EXPARG3 !, surflay - - LOGICAL :: LSTRESS - - - DATA LSTRESS/.FALSE./ !,surflay/20./ - -!**** -!**** -------------------------------------------------- - -!rr next line for debugging, sep 23, 2003, reichle -!rr -!rr write (*,*) 'entering partition()' - - DO N=1,NCH - - WILT=WPWET(N) - WRZ=RZEXC(N)/VGWMAX(N) - CATDEFX=AMIN1( CATDEF(N) , CDCR1(N) ) - -! CDI DEFINES IF THE SHAPE PARAMETER IS ADJUSTED IN ONE OR TWO SEGMENTS - if (ara1(n) .ne. ara3(n)) then - cdi=(ara4(n)-ara2(n))/(ara1(n)-ara3(n)) - else - cdi=0. - endif - - AR1(N)= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFX) & - /(1.+ars2(n)*CATDEFX+ars3(n)*CATDEFX*CATDEFX))) - - if (CATDEFX .ge. cdi) then - ax=ara3(n)*CATDEFX+ara4(n) - else - ax=ara1(n)*CATDEFX+ara2(n) - endif - - WMIN=AMIN1(1.,AMAX1(0.,arw4(n)+(1.-arw4(n))* & - (1.+arw1(n)*CATDEFX) & - /(1.+arw2(n)*CATDEFX+arw3(n)*CATDEFX*CATDEFX))) - -!**** CRITICAL VALUE 1: AVERAGE MOISTURE IN ROOT ZONE AT WMIN -!**** ASSOCIATED WITH CATDEF. - ARG1=AMAX1(-40., AMIN1(40., -AX*(1.-WMIN))) - EXPARG1=EXP(ARG1) - RZEQX=(WMIN-1.-(2./AX))*EXPARG1 + WMIN + (2./AX) - RZEQXI=AX*EXPARG1 * & - ( -1. -(2./AX) - (2./(AX*AX)) + WMIN + (WMIN/AX) ) & - + WMIN + 2./AX - AR20=1.+(-AX-1.+AX*WMIN)*EXPARG1 - RZEQXI=RZEQXI/(AR20+1.E-20) - -!**** CRITICAL VALUE 2: AVERAGE MOISTURE IN ROOT ZONE WHEN WMIN -!**** IS EXACTLY AT WILTING POINT. - ARG2=AMAX1(-40., AMIN1(40., -AX*(1.-WILT))) - EXPARG2=EXP(ARG2) - RZEQW=(WILT-1.-(2./AX))*EXPARG2 + WILT + (2./AX) - RZEQWI=AX*EXPARG2 * & - ( -1. -(2./AX) - (2./(AX*AX)) + WILT + (WILT/AX) ) & - + WILT + 2./AX - AR20=1.+(-AX-1.+AX*WILT)*EXPARG2 - RZEQWI=RZEQWI/(AR20+1.E-20) - -!**** SITUATION 1: CATDEF LE CDCR1 - IF(CATDEF(N) .LE. CDCR1(N)) THEN - RZEQY=RZEQX+WRZ - RZEQYI=RZEQXI+WRZ - WMNEW=WMIN+WRZ - ARG3=AMAX1(-40., AMIN1(40., -AX*(1.-WMNEW))) - EXPARG3=EXP(ARG3) - AREA1=(1.+AX-AX*WMIN)*EXPARG1 - AREA2=(1.+AX-AX*WMNEW)*EXPARG3 - IF(WMNEW .GE. WILT) THEN - AR1(N)=AR1(N)+AREA2-AREA1 - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - IF(WMNEW .LT. WILT) THEN - AREA3=(1.+AX-AX*WILT)*EXPARG2 - AR1(N)=AR1(N)+AREA3-AREA1 - AR2(N)=1.-AR1(N) - FACTOR=(RZEQX+WRZ-WILT)/(RZEQW-WILT) - AR1(N)=AR1(N)*FACTOR - AR2(N)=AR2(N)*FACTOR - AR4(N)=1.-FACTOR - ENDIF - ENDIF - -!**** SITUATION 2: CATDEF GT CDCR1 - IF(CATDEF(N) .GT. CDCR1(N)) THEN - FACTOR=(CDCR2(N)-CATDEF(N))/(CDCR2(N)-CDCR1(N)) - RZEQY=WILT+(RZEQX-WILT)*FACTOR+WRZ - RZEQYI=WILT+(RZEQXI-WILT)*FACTOR+WRZ - - IF(RZEQY .LT. WILT) THEN - IF(RZEQY .LT. WILT-.001) THEN -!rr WRITE(*,*) 'RZEXC WAY TOO LOW! N=',N,' RZEQY=',RZEQY -!rr WRITE(*,*) 'SRFEXC=',SRFEXC(N),'RZEXC=',RZEXC(N), -!rr & 'CATDEF=',CATDEF(N) -! ELSE -! WRITE(*,*) 'RZEXC TOO LOW N=',N - ENDIF - RZEQY=WILT - RZEQYI=WILT - ENDIF - - IF(RZEQY .GE. RZEQX) THEN ! RZEXC BRINGS MOISTURE ABOVE CDCR1 POINT - WMNEW=WMIN+(RZEQY-RZEQX) - ARG3=AMAX1(-40., AMIN1(40., -AX*(1.-WMNEW))) - EXPARG3=EXP(ARG3) - AREA1=(1.+AX-AX*WMIN)*EXPARG1 - AREA2=(1.+AX-AX*WMNEW)*EXPARG3 - AR1(N)=AR1(N)+(AREA2-AREA1) - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - - IF(RZEQY .LT. RZEQX .AND. RZEQY .GE. RZEQW) THEN - CATDEFW=CDCR2(N)+((RZEQW-WILT)/(RZEQX-WILT))*(CDCR1(N)-CDCR2(N)) - AR1W= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFW) & - /(1.+ars2(n)*CATDEFW+ars3(n)*CATDEFW*CATDEFW))) - FACTOR=(RZEQY-RZEQW)/(RZEQX-RZEQW) - AR1(N)=AR1W+FACTOR*(AR1(N)-AR1W) - AR2(N)=1.-AR1(N) - AR4(N)=0. - ENDIF - - IF(RZEQY .LT. RZEQW) THEN - CATDEFW=CDCR2(N)+((RZEQW-WILT)/(RZEQX-WILT))*(CDCR1(N)-CDCR2(N)) - AR1W= AMIN1(1.,AMAX1(0.,(1.+ars1(n)*CATDEFW) & - /(1.+ars2(n)*CATDEFW+ars3(n)*CATDEFW*CATDEFW))) - AR1(N)=AR1W - AR2(N)=1.-AR1(N) - FACTOR=(RZEQY-WILT)/(RZEQW-WILT) - AR1(N)=AR1(N)*FACTOR - AR2(N)=AR2(N)*FACTOR - AR4(N)=1.-FACTOR - ENDIF - - ENDIF - - RZI(N)=RZEQYI - - SWSRF1(N)=1. -!mjs: changed .001 temporarily because of large bee. - SWSRF2(N)=AMIN1(1., AMAX1(0.01, RZEQYI)) - SWSRF4(N)=AMIN1(1., AMAX1(0.01, WILT)) - -!**** EXTRAPOLATION OF THE SURFACE WETNESSES - -! 1st step: surface wetness in the unstressed fraction without considering -! the surface excess; we just assume an equilibrium profile from -! the middle of the root zone to the surface. - - SWSRF2(N)=((SWSRF2(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) - SWSRF4(N)=((SWSRF4(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) - -! srfmx is the maximum amount of water that can be added to the surface layer -! The choice of defining SWSRF4 like SWSRF2 needs to be better examined. - srfmx(n)=ar2(n)*(1.-swsrf2(n))*(dzsf(n)*poros(n)) - srfmx(n)=srfmx(n)+ar4(n)*(1.-swsrf4(n))*(dzsf(n)*poros(n)) -!**** For calculation of srfmn, assume surface moisture associated with -!**** AR1 is constantly replenished by water table. - srfmn(n)=-(ar2(n)*swsrf2(n)+ar4(n)*swsrf4(n))*(dzsf(n)*poros(n)) - - if(srfexc(n).gt.srfmx(n)) then - cor=srfexc(n)-srfmx(n) ! The correction is here - srfexc(n)=srfmx(n) - catdef(n)=catdef(n)-cor - if(catdef(n).lt.0.) then - runsrf(n)=runsrf(n)-catdef(n)/dtstep - catdef(n)=0. - endif - else if(srfexc(n).lt.srfmn(n)) then - cor=srfexc(n)-srfmn(n) - catdef(n)=catdef(n)-cor - srfexc(n)=srfmn(n) - else - cor=0. - endif - - SWSRF2(N)=SWSRF2(N)+SRFEXC(N)/(dzsf(n)*poros(n)*(1.-ar1(n))+1.e-20) - SWSRF2(N)=AMIN1(1., AMAX1(1.E-5, SWSRF2(N))) - swsrf4(n)=swsrf4(n)+srfexc(n)/(dzsf(n)*poros(n)*(1.-ar1(n))+1.e-20) - SWSRF4(N)=AMIN1(1., AMAX1(1.E-5, SWSRF4(N))) - - IF (AR1(N) .ge. 1.-1.E-5) then - AR1(N)=1. - AR2(N)=0. - AR4(N)=0. - SWSRF2(N)=1. - SWSRF4(N)=wilt - ENDIF - - IF (AR1(N) .LT. 0.) then -!rr IF(AR1(N) .LT. -1.E-3) WRITE(*,*) 'AR1 TOO LOW: AR1=',AR1(N) - AR1(N)=0. - ENDIF - ar1(n)=amax1(0., amin1(1., ar1(n))) - ar2(n)=amax1(0., amin1(1., ar2(n))) - ar4(n)=amax1(0., amin1(1., ar4(n))) - asum=ar1(n)+ar2(n)+ar4(n) - if(asum .lt. .9999 .or. asum .gt. 1.0001) then - write(*,*) 'Areas do not add to 1: sum=',asum,'N=',n - endif - - - ENDDO - - - RETURN - END SUBROUTINE PARTITION - - SUBROUTINE RZEQUIL ( & - NCH,CATDEF,VGWMAX,CDCR1,CDCR2,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4, & - RZEQ & - ) - - IMPLICIT NONE - - INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN), DIMENSION(NCH) :: CATDEF, VGWMAX, CDCR1, CDCR2, & - WPWET, ars1, ars2, ars3, ara1, ara2, ara3, ara4, arw1, & - arw2, arw3, arw4 - - REAL, INTENT(OUT), DIMENSION(NCH) :: RZEQ - - INTEGER N - REAL AX,WMIN,ASCALE,cdi,wilt,catdefx,factor,ARG1,EXPARG1 - -! ---------------------------------------------------------------------- - - DO N=1,NCH - - WILT=WPWET(N) - CATDEFX=AMIN1( CATDEF(N) , CDCR1(N) ) - -! CDI DEFINES IF THE SHAPE PARAMETER IS ADJUSTED IN ONE OR TWO SEGMENTS - if (ara1(n) .ne. ara3(n)) then - cdi=(ara4(n)-ara2(n))/(ara1(n)-ara3(n)) - else - cdi=0. - endif - - if (CATDEFX .ge. cdi) then - ax=ara3(n)*CATDEFX+ara4(n) - else - ax=ara1(n)*CATDEFX+ara2(n) - endif - - WMIN=AMIN1(1.,AMAX1(0.,arw4(n)+(1.-arw4(n))*(1.+arw1(n)*CATDEFX) & - /(1.+arw2(n)*CATDEFX+arw3(n)*CATDEFX*CATDEFX))) - - ARG1=AMAX1(-40., AMIN1(40., -AX*(1.-WMIN))) - EXPARG1=EXP(ARG1) - RZEQ(N)=(WMIN-1.-(2./AX))*EXPARG1 + WMIN + (2./AX) - - IF(CATDEF(N) .GT. CDCR1(N)) THEN - FACTOR=(CDCR2(N)-CATDEF(N))/(CDCR2(N)-CDCR1(N)) - RZEQ(N)=WILT+(RZEQ(N)-WILT)*FACTOR - ENDIF - -! scaling: - RZEQ(N)=AMIN1(1.,AMAX1(0.,RZEQ(N))) - RZEQ(N)=RZEQ(N)*VGWMAX(N) - - ENDDO - - RETURN - END SUBROUTINE RZEQUIL -