diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt index acaec2953..73ed82643 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt @@ -7,6 +7,7 @@ set (alldirs GEOSlana_GridComp GEOSroute_GridComp GEOSigni_GridComp + GEOSirrigation_GridComp ) esma_add_library (${this} 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 69d73008e..8a137bef9 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 @@ -14,12 +14,14 @@ module GEOS_LandGridCompMod ! to determine relevant time-dependent land-surface characteristics. ! All parameters calculated in VegdynGridComp are required by CatchGridComp. ! Furthermore, several exports of the Vegdyn routines are also exports -! from the Land composite, for use in other modules, such as the case -! for lai and grn needed in radiation. Vegdyn will be updated first. -! Then the catchment call will be issued. The composite exports -! consist of the union of the catchment exports with a subset of the -! vegdyn exports. All imports and exports are on the prescribed tile -! grid in the (IM, JM)=(NTILES, 1) convention. +! from the Land composite, for use in other modules. For example, +! lai and grn are needed in radiation. Vegdyn will be updated first. +! Then the Catch[CN] Phase=1 call will be issued. IrrigationGridComp +! computes the irrigation rate IMPORT required by Catch[CN], +! followed by Catch[CN] Phase=2 (incl. application of irrigation). +! The composite exports consist of the union of the catchment exports with a +! subset of the vegdyn and Irrigation exports. All imports and exports are +! on the prescribed tile grid in the (IM, JM)=(NTILES, 1) convention. ! ! !USES: @@ -27,11 +29,12 @@ module GEOS_LandGridCompMod use ESMF use MAPL - use GEOS_VegdynGridCompMod, only : VegdynSetServices => SetServices - use GEOS_CatchGridCompMod, only : CatchSetServices => SetServices - use GEOS_CatchCNGridCompMod, only : CatchCNSetServices => SetServices - use GEOS_IgniGridCompMod, only : IgniSetServices => SetServices -! use GEOS_RouteGridCompMod, only : RouteSetServices => SetServices + use GEOS_VegdynGridCompMod, only : VegdynSetServices => SetServices + use GEOS_CatchGridCompMod, only : CatchSetServices => SetServices + use GEOS_CatchCNGridCompMod, only : CatchCNSetServices => SetServices + use GEOS_IgniGridCompMod, only : IgniSetServices => SetServices + use GEOS_IrrigationGridCompMod, only : IrrigationSetServices => SetServices + ! use GEOS_RouteGridCompMod, only : RouteSetServices => SetServices implicit none private @@ -45,8 +48,8 @@ module GEOS_LandGridCompMod integer :: VEGDYN - integer, allocatable :: CATCH(:), ROUTE (:), CATCHCN (:) - integer :: LSM_CHOICE, RUN_ROUTE, DO_GOSWIM + integer, allocatable :: CATCH(:), ROUTE (:), CATCHCN (:), IRRIGATION(:) + integer :: LSM_CHOICE, RUN_ROUTE, DO_GOSWIM, RUN_IRRIG integer :: IGNI logical :: DO_FIRE_DANGER @@ -65,10 +68,11 @@ subroutine SetServices ( GC, RC ) type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component integer, optional :: RC ! return code -! !DESCRIPTION: The SetServices for the Physics GC needs to register its +! !DESCRIPTION: The SetServices for the Land GC needs to register its ! Initialize and Run. It uses the MAPL\_Generic construct for defining ! state specs and couplings among its children. In addition, it creates the -! children GCs (VegDyn, Catch, CatchCN, Route) and runs their respective SetServices. +! children GCs (VegDyn, Catch, CatchCN, Irrigation, Route) and runs their +! respective SetServices. !EOP @@ -102,7 +106,7 @@ subroutine SetServices ( GC, RC ) call ESMF_GridCompGet(GC ,& NAME=COMP_NAME ,& CONFIG=CF ,& - RC=STATUS ) + RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // 'SetServices' @@ -114,7 +118,7 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_GetResource ( MAPL, NUM_LDAS_ENSEMBLE, Label="NUM_LDAS_ENSEMBLE:", DEFAULT=1, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource ( MAPL, ens_id_width, Label="ENS_ID_WIDTH:", DEFAULT=0, RC=STATUS) + call MAPL_GetResource ( MAPL, ens_id_width, Label="ENS_ID_WIDTH:", DEFAULT=0, RC=STATUS) VERIFY_(STATUS) tmp = '' @@ -129,9 +133,9 @@ subroutine SetServices ( GC, RC ) call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS ) VERIFY_(STATUS) - call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run1, RC=STATUS ) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run1, RC=STATUS ) VERIFY_(STATUS) - call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run2, RC=STATUS ) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run2, RC=STATUS ) VERIFY_(STATUS) call ESMF_ConfigGetAttribute ( CF, NUM_CATCH, Label="NUM_CATCH_ENSEMBLES:", default=1, RC=STATUS) @@ -149,15 +153,20 @@ subroutine SetServices ( GC, RC ) ! and Runoff Routing Model (0: OFF, 1: ON) ! ------------------------------------------------------- - call MAPL_GetResource ( MAPL, LSM_CHOICE, Label="LSM_CHOICE:", DEFAULT=1, RC=STATUS) + call MAPL_GetResource (MAPL, LSM_CHOICE, Label="LSM_CHOICE:", DEFAULT=1, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource (MAPL, SURFRC, label = 'SURFRC:', default = 'GEOS_SurfaceGridComp.rc', RC=STATUS) ; VERIFY_(STATUS) - SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SURFRC, Label = 'SURFRC:', DEFAULT='GEOS_SurfaceGridComp.rc', RC=STATUS) + VERIFY_(STATUS) + + SCF = ESMF_ConfigCreate( rc=status) ; VERIFY_(STATUS) call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) - call MAPL_GetResource (SCF, RUN_ROUTE, label='RUN_ROUTE:', DEFAULT=0, __RC__ ) - call MAPL_GetResource (SCF, DO_GOSWIM, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) - call MAPL_GetResource (SCF, DO_FIRE_DANGER, label='FIRE_DANGER:', DEFAULT=.false., __RC__ ) - call ESMF_ConfigDestroy (SCF, __RC__) + + call MAPL_GetResource (SCF, RUN_ROUTE, label='RUN_ROUTE:', DEFAULT=0, __RC__ ) + call MAPL_GetResource (SCF, RUN_IRRIG, label='RUN_IRRIG:', DEFAULT=0, __RC__ ) + call MAPL_GetResource (SCF, DO_GOSWIM, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) + call MAPL_GetResource (SCF, DO_FIRE_DANGER, label='FIRE_DANGER:', DEFAULT=.false., __RC__ ) + + call ESMF_ConfigDestroy(SCF, __RC__) SELECT CASE (LSM_CHOICE) @@ -193,7 +202,23 @@ subroutine SetServices ( GC, RC ) end do end if - END SELECT + END SELECT ! LSM_CHOICE + + if(RUN_IRRIG /= 0) then + allocate (IRRIGATION(NUM_CATCH), stat=status) + VERIFY_(STATUS) + if (NUM_CATCH == 1) then + IRRIGATION(1) = MAPL_AddChild(GC, NAME='IRRIGATION'//trim(tmp), SS=IrrigationSetServices, RC=STATUS) + VERIFY_(STATUS) + else + do I = 1, NUM_CATCH + WRITE(TMP,'(I3.3)') I + GCName = 'ens' // trim(TMP) // ':IRRIGATION' + IRRIGATION(I) = MAPL_AddChild(GC, NAME=GCName, SS=IrrigationSetServices, RC=STATUS) + VERIFY_(STATUS) + end do + end if + end if ! IF(RUN_ROUTE == 1) THEN ! if (NUM_CATCH == 1) then @@ -670,6 +695,11 @@ subroutine SetServices ( GC, RC ) CHILD_ID = CATCH(1), & RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, & + SHORT_NAME = 'SPIRRG', & + CHILD_ID = CATCH(1), & + RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'WESNN1', & CHILD_ID = CATCH(1), & @@ -1149,6 +1179,8 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'SPSNOW' , CHILD_ID = CATCHCN(1), RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'SPIRRG' , CHILD_ID = CATCHCN(1), RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WESNN1' , CHILD_ID = CATCHCN(1), RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WESNN2' , CHILD_ID = CATCHCN(1), RC=STATUS ) @@ -1349,7 +1381,18 @@ subroutine SetServices ( GC, RC ) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'ROC002', CHILD_ID = CATCHCN(1), RC=STATUS) ; VERIFY_(STATUS) endif - END SELECT + END SELECT ! LSM_CHOICE (Catch, CatchCN) + + ! ------------------------------------------------------------------------------------------------------------------- + + if (RUN_IRRIG /= 0) then + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'IRRG_RATE_SPR', CHILD_ID = IRRIGATION(1),RC=STATUS ) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'IRRG_RATE_DRP', CHILD_ID = IRRIGATION(1),RC=STATUS ) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'IRRG_RATE_FRW', CHILD_ID = IRRIGATION(1),RC=STATUS ) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'IRRG_RATE_PDY', CHILD_ID = IRRIGATION(1),RC=STATUS ) ; VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'IRRG_RATE_TOT', CHILD_ID = IRRIGATION(1),RC=STATUS ) ; VERIFY_(STATUS) + end if + ! These are from RUN1 of vegdyn and the first catchment instance call MAPL_AddExportSpec ( GC, & @@ -1377,6 +1420,7 @@ subroutine SetServices ( GC, RC ) CHILD_ID = VEGDYN,& RC=STATUS ) VERIFY_(STATUS) + ! IF(RUN_ROUTE == 1) THEN ! call MAPL_AddExportSpec ( GC, & ! SHORT_NAME = 'QOUTFLOW', & @@ -1385,7 +1429,7 @@ subroutine SetServices ( GC, RC ) ! VERIFY_(STATUS) ! ENDIF - + if (DO_FIRE_DANGER) then call MAPL_AddExportSpec ( GC, SHORT_NAME = 'FFMC', CHILD_ID = IGNI, __RC__ ) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'GFMC', CHILD_ID = IGNI, __RC__ ) @@ -1453,6 +1497,26 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) end if + if (RUN_IRRIG /= 0) then + call MAPL_AddConnectivity ( & + GC ,& + SHORT_NAME = (/'POROS ', 'WPWET ' ,& + 'VGWMAX ', 'WCRZ '/) ,& + SRC_ID = CATCH(I) ,& + DST_ID = IRRIGATION(I) ,& + RC = STATUS ) + VERIFY_(STATUS) + + call MAPL_AddConnectivity ( & + GC ,& + SHORT_NAME = (/'IRRG_RATE_SPR', 'IRRG_RATE_DRP' ,& + 'IRRG_RATE_FRW', 'IRRG_RATE_PDY'/) ,& + SRC_ID = IRRIGATION(I) ,& + DST_ID = CATCH(I) ,& + RC = STATUS ) + VERIFY_(STATUS) + end if + ! IF(RUN_ROUTE == 1) THEN ! call MAPL_AddConnectivity ( & ! GC ,& @@ -1485,6 +1549,26 @@ subroutine SetServices ( GC, RC ) RC = STATUS ) VERIFY_(STATUS) end if + + if (RUN_IRRIG /= 0) then + call MAPL_AddConnectivity ( & + GC ,& + SHORT_NAME = (/'POROS ', 'WPWET ' ,& + 'VGWMAX ', 'WCRZ '/) ,& + SRC_ID = CATCH(I) ,& + DST_ID = IRRIGATION(I) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddConnectivity ( & + GC ,& + SHORT_NAME = (/'IRRG_RATE_SPR', 'IRRG_RATE_DRP' ,& + 'IRRG_RATE_FRW', 'IRRG_RATE_PDY'/) ,& + SRC_ID = IRRIGATION(I) ,& + DST_ID = CATCH(I) ,& + RC=STATUS ) + VERIFY_(STATUS) + end if ! IF(RUN_ROUTE == 1) THEN ! call MAPL_AddConnectivity ( & @@ -1497,6 +1581,17 @@ subroutine SetServices ( GC, RC ) ! VERIFY_(STATUS) ! ENDIF END SELECT + + if (RUN_IRRIG /= 0) then + call MAPL_AddConnectivity ( & + GC ,& + SHORT_NAME = (/'LAI '/) ,& + SRC_ID = VEGDYN ,& + DST_ID = IRRIGATION(I) ,& + RC=STATUS ) + VERIFY_(STATUS) + end if + END DO @@ -1511,8 +1606,9 @@ subroutine SetServices ( GC, RC ) call MAPL_GenericSetServices(GC, RC=STATUS ) VERIFY_(STATUS) - if (allocated(CATCH)) deallocate(CATCH) - if (allocated(CATCHCN)) deallocate(CATCHCN) + if (allocated(CATCH)) deallocate(CATCH) + if (allocated(CATCHCN)) deallocate(CATCHCN) + if (allocated(IRRIGATION)) deallocate(IRRIGATION) RETURN_(ESMF_SUCCESS) @@ -1550,14 +1646,14 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Local derived type aliases - type (MAPL_MetaComp ), pointer :: MAPL - type (MAPL_MetaComp ), pointer :: CHILD_MAPL - type (MAPL_LocStream ) :: LOCSTREAM - type (ESMF_DELayout ) :: LAYOUT - type (ESMF_Config ) :: CF - type (ESMF_GridComp ), pointer :: GCS(:) + type (MAPL_MetaComp ), pointer :: MAPL + type (MAPL_MetaComp ), pointer :: CHILD_MAPL + type (MAPL_LocStream ) :: LOCSTREAM + type (ESMF_DELayout ) :: LAYOUT + type (ESMF_Config ) :: CF + type (ESMF_GridComp ), pointer :: GCS(:) - integer :: I + integer :: I !============================================================================= @@ -1577,7 +1673,7 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) call MAPL_TimerOn(MAPL,"INITIALIZE", RC=STATUS ); VERIFY_(STATUS) - call MAPL_TimerOn(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) + call MAPL_TimerOn(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) ! Get the land tilegrid and the child components !----------------------------------------------- @@ -1623,8 +1719,8 @@ subroutine Run1(GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_Clock), intent(inout) :: CLOCK ! The clock integer, optional, intent( out) :: RC ! Error code -! !DESCRIPTION: This first run method calls the children's -! first run methods. VEGDYN has only one, and it is called here. +! !DESCRIPTION: This first run method calls the children's first run methods. +! VEGDYN and Irrigation have only one run method, which is called here. !EOP ! ErrLog Variables @@ -1660,13 +1756,13 @@ subroutine Run1(GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) call MAPL_TimerOn(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) - call MAPL_TimerOn(MAPL,"RUN1", RC=STATUS ); VERIFY_(STATUS) + call MAPL_TimerOn(MAPL,"RUN1", RC=STATUS ); VERIFY_(STATUS) call MAPL_Get (MAPL, GCS=GCS, GIM=GIM, GEX=GEX, GCnames=GCnames,rc=STATUS) VERIFY_(STATUS) -! Call the children's RUN methods -!-------------------------------- +! Call the children's RUN methods (PHASE=1) +!------------------------------------------ DO I = 1, size(GCS) call MAPL_TimerOn(MAPL,trim(GCnames(i)), RC=STATUS ); VERIFY_(STATUS) @@ -1676,7 +1772,7 @@ subroutine Run1(GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_TimerOff(MAPL,trim(GCnames(i)), RC=STATUS ); VERIFY_(STATUS) END DO - call MAPL_TimerOff(MAPL,"RUN1", RC=STATUS ); VERIFY_(STATUS) + call MAPL_TimerOff(MAPL,"RUN1", RC=STATUS ); VERIFY_(STATUS) call MAPL_TimerOff(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) RETURN_(ESMF_SUCCESS) @@ -1734,13 +1830,13 @@ subroutine Run2(GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) call MAPL_TimerOn(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) - call MAPL_TimerOn(MAPL,"RUN2", RC=STATUS ); VERIFY_(STATUS) + call MAPL_TimerOn(MAPL,"RUN2", RC=STATUS ); VERIFY_(STATUS) call MAPL_Get (MAPL, GCS=GCS, GIM=GIM, GEX=GEX, GCnames=GCnames,rc=STATUS) VERIFY_(STATUS) -! Call the children's RUN methods -!-------------------------------- +! Call the children's RUN methods (PHASE=2) +!------------------------------------------ DO I=1,size(GCS) if (I == VEGDYN) cycle call MAPL_TimerOn(MAPL,trim(GCnames(i)), RC=STATUS ); VERIFY_(STATUS) @@ -1750,7 +1846,7 @@ subroutine Run2(GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_TimerOff(MAPL,trim(GCnames(i)), RC=STATUS ); VERIFY_(STATUS) END DO - call MAPL_TimerOff(MAPL,"RUN2", RC=STATUS ); VERIFY_(STATUS) + call MAPL_TimerOff(MAPL,"RUN2", RC=STATUS ); VERIFY_(STATUS) call MAPL_TimerOff(MAPL,"TOTAL", RC=STATUS ); VERIFY_(STATUS) RETURN_(ESMF_SUCCESS) 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 05b20561d..c888d20ea 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 @@ -671,6 +671,42 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_SPR' ,& + LONG_NAME = 'irrigation_flux_sprinkler' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_DRP' ,& + LONG_NAME = 'irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_FRW' ,& + LONG_NAME = 'irrigation_flux_furrow' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_PDY' ,& + LONG_NAME = 'irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + ! EXPORT STATE: call MAPL_AddExportSpec ( GC, SHORT_NAME = 'LST', CHILD_ID = CATCHCN, RC=STATUS ) VERIFY_(STATUS) @@ -844,6 +880,8 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'SPSNOW' , CHILD_ID = CATCHCN, RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec ( GC, SHORT_NAME = 'SPIRRG' , CHILD_ID = CATCHCN, RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WESNN1' , CHILD_ID = CATCHCN, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, SHORT_NAME = 'WESNN2' , CHILD_ID = CATCHCN, RC=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 cbebc22a3..07ae5ad2d 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 @@ -66,7 +66,7 @@ module GEOS_CatchCNCLM40GridCompMod 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, & - catch_calc_zbar, catch_calc_peatclsm_waterlevel, irrigation_rate, & + catch_calc_zbar, catch_calc_peatclsm_waterlevel, & gndtmp use catch_wrap_stateMod @@ -198,7 +198,7 @@ subroutine SetServices ( GC, RC ) ! Local Variables type(MAPL_MetaComp), pointer :: MAPL=>null() - integer :: OFFLINE_MODE, RUN_IRRIG, ATM_CO2, PRESCRIBE_DVG, N_CONST_LAND4SNWALB + integer :: OFFLINE_MODE, ATM_CO2, PRESCRIBE_DVG, N_CONST_LAND4SNWALB integer :: RESTART, SNOW_ALBEDO_INFO ! Begin... @@ -219,7 +219,6 @@ subroutine SetServices ( GC, RC ) call MAPL_GetResource ( MAPL, OFFLINE_MODE, Label="CATCHMENT_OFFLINE:", DEFAULT=0, _RC) call MAPL_GetResource ( MAPL, ATM_CO2, Label="ATM_CO2:", _RC) call MAPL_GetResource ( MAPL, N_CONST_LAND4SNWALB, Label="N_CONST_LAND4SNWALB:", _RC) - call MAPL_GetResource ( MAPL, RUN_IRRIG, Label="RUN_IRRIG:", _RC) call MAPL_GetResource ( MAPL, PRESCRIBE_DVG, Label="PRESCRIBE_DVG:", _RC) call MAPL_GetResource ( MAPL, SNOW_ALBEDO_INFO, Label="SNOW_ALBEDO_INFO:", _RC) @@ -770,6 +769,42 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_SPR' ,& + LONG_NAME = 'irrigation_flux_sprinkler' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_DRP' ,& + LONG_NAME = 'irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_FRW' ,& + LONG_NAME = 'irrigation_flux_furrow' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_PDY' ,& + LONG_NAME = 'irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + ! !INTERNAL STATE: ! if is_offline, some variables ( in the last) are not required @@ -2033,117 +2068,11 @@ subroutine SetServices ( GC, RC ) endif -! IRRIGATION MODEL INTERNAL - - IF (RUN_IRRIG /= 0) THEN - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'fraction_of_irrigated_cropland',& - UNITS = '1' ,& - SHORT_NAME = 'IRRIGFRAC' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'fraction_of_paddy_cropland',& - UNITS = '1' ,& - SHORT_NAME = 'PADDYFRAC' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'Maximum_LAI' ,& - UNITS = '1' ,& - SHORT_NAME = 'LAIMAX' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'Minimum_LAI' ,& - UNITS = '1' ,& - SHORT_NAME = 'LAIMIN' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_primary_type' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMPT' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_secondary_type' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMST' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_primary_fraction' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMPF' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_secondary_fraction' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMSF' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - ENDIF - - !EOS ! EXPORT STATE: - IF (RUN_IRRIG /= 0) THEN - call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'irrigation_rate' ,& - UNITS = 'kg m-2 s-1' ,& - SHORT_NAME = 'IRRIGRATE' ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - VERIFY_(STATUS) - ENDIF - + call MAPL_AddExportSpec(GC, & LONG_NAME = 'evaporation' ,& UNITS = 'kg m-2 s-1' ,& @@ -3199,6 +3128,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + SHORT_NAME = 'SPIRRG', & + LONG_NAME = 'Spurious_irrigation_flux', & + UNITS = 'kg m-2 s-1', & + DIMS = MAPL_DimsTileOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& LONG_NAME = 'vegetation_type' ,& UNITS = '1' ,& @@ -4585,7 +4523,10 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: GRN real, dimension(:), pointer :: ASCATZ0 real, dimension(:), pointer :: NDVI - + real, dimension(:), pointer :: IRRG_RATE_SPR + real, dimension(:), pointer :: IRRG_RATE_DRP + real, dimension(:), pointer :: IRRG_RATE_FRW + real, dimension(:), pointer :: IRRG_RATE_PDY real, dimension(:,:), pointer :: DUDP real, dimension(:,:), pointer :: DUSV real, dimension(:,:), pointer :: DUWT @@ -4702,15 +4643,7 @@ subroutine Driver ( RC ) real, dimension(:,:), pointer :: RBC002 real, dimension(:,:), pointer :: ROC001 real, dimension(:,:), pointer :: ROC002 - real, dimension(:), pointer :: IRRIGFRAC - real, dimension(:), pointer :: PADDYFRAC - real, dimension(:), pointer :: LAIMAX - real, dimension(:), pointer :: LAIMIN - real, dimension(:), pointer :: CLMPT - real, dimension(:), pointer :: CLMST - real, dimension(:), pointer :: CLMPF - real, dimension(:), pointer :: CLMSF - + ! ----------------------------------------------------- ! EXPORT Pointers ! ----------------------------------------------------- @@ -4806,6 +4739,7 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: SPLAND real, dimension(:), pointer :: SPWATR real, dimension(:), pointer :: SPSNOW + real, dimension(:), pointer :: SPIRRG real, dimension(:), pointer :: CNLAI real, dimension(:), pointer :: CNTLAI @@ -4868,7 +4802,6 @@ subroutine Driver ( RC ) real, pointer, dimension(:) :: RMELTBC002 real, pointer, dimension(:) :: RMELTOC001 real, pointer, dimension(:) :: RMELTOC002 - real, pointer, dimension(:) :: IRRIGRATE real, pointer, dimension(:) :: PEATCLSM_WATERLEVEL real, pointer, dimension(:) :: PEATCLSM_FSWCHANGE @@ -4877,7 +4810,7 @@ subroutine Driver ( RC ) ! -------------------------------------------------------------------------- INTEGER,pointer,dimension(:) :: CAT_ID - real,pointer,dimension(:) :: dzsf + real,pointer,dimension(:) :: dzsf_in_mm real,pointer,dimension(:) :: swnetfree real,pointer,dimension(:) :: swnetsnow real,pointer,dimension(:) :: qa1 @@ -4908,7 +4841,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: fveg1, fveg2 real,pointer,dimension(:) :: FICE1TMP real,pointer,dimension(:) :: SLDTOT - + real,pointer,dimension(:) :: PLS_SPR ! real*8,pointer,dimension(:) :: fsum real,pointer,dimension(:,:) :: ghtcnt @@ -4985,6 +4918,8 @@ subroutine Driver ( RC ) integer :: NTILES integer :: I, J, K, N + real, dimension(:), allocatable :: AR4, werror ! for catch_calc_soil_moist() after irrigation application + ! dummy variables for call to get snow temp real :: FICE @@ -5014,7 +4949,6 @@ subroutine Driver ( RC ) ! unadulterated TC's and QC's real, pointer :: TC1_0(:), TC2_0(:), TC4_0(:) real, pointer :: QA1_0(:), QA2_0(:), QA4_0(:) - real, pointer :: PLSIN(:) ! CATCHMENT_SPINUP integer :: CurrMonth, CurrDay, CurrHour, CurrMin, CurrSec @@ -5249,6 +5183,10 @@ subroutine Driver ( RC ) call MAPL_GetPointer(IMPORT,SSSV ,'SSSV' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,SSWT ,'SSWT' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,SSSD ,'SSSD' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_SPR,'IRRG_RATE_SPR',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_DRP,'IRRG_RATE_DRP',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_FRW,'IRRG_RATE_FRW',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_PDY,'IRRG_RATE_PDY',RC=STATUS); VERIFY_(STATUS) ! ----------------------------------------------------- ! INTERNAL Pointers @@ -5348,17 +5286,6 @@ subroutine Driver ( RC ) call MAPL_GetPointer(INTERNAL,ROC002 ,'ROC002' , RC=STATUS); VERIFY_(STATUS) endif - IF (catchcn_internal%RUN_IRRIG /= 0) THEN - call MAPL_GetPointer(INTERNAL,IRRIGFRAC ,'IRRIGFRAC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,PADDYFRAC ,'PADDYFRAC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,LAIMAX ,'LAIMAX' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,LAIMIN ,'LAIMIN' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMPT ,'CLMPT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMST ,'CLMST' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMPF ,'CLMPF' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMSF ,'CLMSF' , RC=STATUS); VERIFY_(STATUS) - ENDIF - ! ----------------------------------------------------- ! EXPORT POINTERS ! ----------------------------------------------------- @@ -5455,6 +5382,8 @@ subroutine Driver ( RC ) 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) + if(CATCHCN_INTERNAL%RUN_IRRIG /= 0) & + call MAPL_GetPointer(EXPORT,SPIRRG, 'SPIRRG' ,ALLOC=.true.,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) @@ -5512,8 +5441,7 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,RMELTOC002,'RMELTOC002', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,PEATCLSM_WATERLEVEL,'PEATCLSM_WATERLEVEL',RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,PEATCLSM_FSWCHANGE ,'PEATCLSM_FSWCHANGE', RC=STATUS); VERIFY_(STATUS) - IF (catchcn_internal%RUN_IRRIG /= 0) call MAPL_GetPointer(EXPORT,IRRIGRATE ,'IRRIGRATE' , RC=STATUS); VERIFY_(STATUS) - + NTILES = size(PS) allocate( ityp(ntiles,nveg,nzone) ) @@ -5670,7 +5598,7 @@ subroutine Driver ( RC ) allocate(FICESOUT(N_SNOW,NTILES)) allocate(TILEZERO (NTILES)) - allocate(DZSF (NTILES)) + allocate(DZSF_IN_MM (NTILES)) allocate(SWNETFREE(NTILES)) allocate(SWNETSNOW(NTILES)) allocate(VEG1 (NTILES)) @@ -5699,7 +5627,7 @@ subroutine Driver ( RC ) allocate(SFMC (NTILES)) allocate(RZMC (NTILES)) allocate(PRMC (NTILES)) - allocate(ENTOT (NTILES)) + allocate(ENTOT (NTILES)) allocate(ghflxsno (NTILES)) allocate(ghflxtskin(NTILES)) allocate(WTOT (NTILES)) @@ -5761,7 +5689,7 @@ subroutine Driver ( RC ) allocate(QA1_0 (NTILES)) allocate(QA2_0 (NTILES)) allocate(QA4_0 (NTILES)) - allocate(PLSIN (NTILES)) + allocate(PLS_SPR (NTILES)) call ESMF_VMGetCurrent ( VM, RC=STATUS ) @@ -5883,7 +5811,7 @@ subroutine Driver ( RC ) ! surface layer depth for soil moisture ! -------------------------------------------------------------------------- - DZSF( :) = catchcn_internal%SURFLAY + DZSF_IN_MM( :) = catchcn_internal%SURFLAY ! same as DZSF but in units of [mm] ! -------------------------------------------------------------------------- ! build arrays from internal state @@ -6348,7 +6276,7 @@ subroutine Driver ( RC ) ! gkw: obtain catchment area fractions and soil moisture ! ------------------------------------------------------ - call catch_calc_soil_moist( ntiles, dzsf, vgwmax, cdcr1, cdcr2, psis, bee, poros, wpwet, & + call catch_calc_soil_moist( ntiles, dzsf_in_mm, 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 ) @@ -7043,28 +6971,49 @@ subroutine Driver ( RC ) ! gkw: end of main CN block - PLSIN = PLS - + PLS_SPR = PLS ! PLS_SPR = large-scale precip plus sprinkler irrigation (if present, see below) + ! -------------------------------------------------------------------------- - ! Call Irrigation Model + ! Add irrigation model imports ! -------------------------------------------------------------------------- - IF ((catchcn_internal%RUN_IRRIG /= 0).AND.(ntiles >0)) THEN - - 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 (catchcn_internal%IRRIG_METHOD, & - NTILES, AGCM_HH, AGCM_MI, sofmin, lons, IRRIGFRAC, PADDYFRAC, & - CLMPT,CLMST, CLMPF, CLMSF, LAIMAX, LAIMIN, LAI0, & - POROS, WPWET, VGWMAX, RZMC, IRRIGRATE) + IF ((CATCHCN_INTERNAL%RUN_IRRIG /= 0)) THEN - PLSIN = PLS + IRRIGRATE + where (IRRG_RATE_SPR > 0) + PLS_SPR = PLS_SPR + IRRG_RATE_SPR + end where + where (IRRG_RATE_DRP > 0) + RZEXC = RZEXC + IRRG_RATE_DRP * DT + end where + where (IRRG_RATE_FRW > 0) + RZEXC = RZEXC + IRRG_RATE_FRW * DT + end where + where (IRRG_RATE_PDY > 0) + SRFEXC = SRFEXC + IRRG_RATE_PDY * DT + end where - ENDIF - + ! after application of irrigation water, make sure soil moisture prognostics + ! (srfexc, rzexc, catdef) remain valid; + ! werror accounts for excess irrigation that cannot be absorbed by the soil; + ! sfmc, rzmc, prmc here are dummies that are required to get werror + + allocate(ar4( NTILES)) + allocate(werror(NTILES)) + + call catch_calc_soil_moist( & + NTILES, dzsf_in_mm, 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, werror ) + + SPIRRG = -werror / DT ! add excess irrigation into spurious term + + deallocate(ar4) + deallocate(werror) + + ENDIF + + #ifdef DBG_CNLSM_INPUTS call MAPL_Get(MAPL, LocStream=LOCSTREAM, RC=STATUS) VERIFY_(STATUS) @@ -7082,12 +7031,12 @@ subroutine Driver ( RC ) ! Inputs - call MAPL_VarWrite(unit, tilegrid, PCU, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, PLS, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, SNO, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, ICE, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, FRZR, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, UUU, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, PCU, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, PLS_SPR, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, SNO, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, ICE, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, FRZR, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, UUU, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, EVSBT (:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) @@ -7170,7 +7119,7 @@ subroutine Driver ( RC ) call MAPL_VarWrite(unit, tilegrid, VEG2, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, FVEG1, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, FVEG2, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, DZSF, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, DZSF_IN_MM, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, BF1, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, BF2, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, BF3, mask=mask, rc=status); VERIFY_(STATUS) @@ -7250,8 +7199,8 @@ subroutine Driver ( RC ) if (ntiles > 0) then call CATCHCN ( NTILES, LONS, LATS, DT,catchcn_internal%USE_FWET_FOR_RUNOFF, & ! LONS, LATS are in [radians] !!! - catchcn_internal%FWETC, catchcn_internal%FWETL, cat_id, VEG1,VEG2,FVEG1,FVEG2,DZSF ,& - PCU , PLSIN , SNO, ICE, FRZR ,& + catchcn_internal%FWETC, catchcn_internal%FWETL, cat_id, VEG1,VEG2,FVEG1,FVEG2,DZSF_IN_MM ,& + PCU , PLS_SPR, SNO, ICE, FRZR ,& UUU ,& EVSBT(:,FSAT), DEVSBT(:,FSAT), DEDTC(:,FSAT) ,& @@ -7650,7 +7599,7 @@ subroutine Driver ( RC ) deallocate(SNDZN ) deallocate(FICESOUT ) deallocate(TILEZERO ) - deallocate(DZSF ) + deallocate(DZSF_IN_MM ) deallocate(SWNETFREE) deallocate(SWNETSNOW) deallocate(VEG1 ) @@ -7804,7 +7753,7 @@ subroutine Driver ( RC ) deallocate( ht ) deallocate( tp ) deallocate( soilice ) - deallocate (PLSIN) + deallocate (PLS_SPR ) call MAPL_TimerOff ( MAPL, "-CATCHCNCLM40" ) RETURN_(ESMF_SUCCESS) @@ -8005,7 +7954,7 @@ subroutine RUN0(gc, import, export, clock, rc) !! Miscellaneous integer :: ntiles, nv, nz real, allocatable :: dummy(:) - real, allocatable :: dzsf(:), ar1(:), ar2(:), wesnn(:,:) + real, allocatable :: dzsf_in_mm(:), ar1(:), ar2(:), wesnn(:,:) real, allocatable :: catdefcp(:), srfexccp(:), rzexccp(:) real, allocatable :: VEG1(:), VEG2(:) integer, allocatable :: ityp(:,:,:) @@ -8213,14 +8162,14 @@ subroutine RUN0(gc, import, export, clock, rc) emis = emis*(1.-asnow) + EMSSNO*asnow ! Compute FR - ! Step 1: set dzsf + ! Step 1: set dzsf_in_mm ! Step 2: compute ar1, ar2 via call to catch_calc_soil_moist() ! Step 3: compute fr ! -step-1- - allocate(dzsf(ntiles), stat=status) + allocate(dzsf_in_mm(ntiles), stat=status) VERIFY_(status) - dzsf = catchcn_internal%SURFLAY + dzsf_in_mm = catchcn_internal%SURFLAY ! -step-2- allocate(ar1(ntiles), stat=status) @@ -8240,7 +8189,7 @@ subroutine RUN0(gc, import, export, clock, rc) rzexccp = rzexc call catch_calc_soil_moist( & ! intent(in) - ntiles, dzsf, vgwmax, cdcr1, cdcr2, & + ntiles, dzsf_in_mm, vgwmax, cdcr1, cdcr2, & psis, bee, poros, wpwet, & ars1, ars2, ars3, & ara1, ara2, ara3, ara4, & @@ -8273,7 +8222,7 @@ subroutine RUN0(gc, import, export, clock, rc) if (allocated(srfexccp)) deallocate(srfexccp) if (allocated(rzexccp)) deallocate(rzexccp) if (allocated(dummy)) deallocate(dummy) - if (allocated(dzsf)) deallocate(dzsf) + if (allocated(dzsf_in_mm)) deallocate(dzsf_in_mm) if (allocated(ar1)) deallocate(ar1) if (allocated(ar2)) deallocate(ar2) if (allocated(wesnn)) deallocate(wesnn) 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 1d87167a4..bfd249ca3 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 @@ -65,7 +65,7 @@ module GEOS_CatchCNCLM45GridCompMod 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, & - catch_calc_zbar, catch_calc_peatclsm_waterlevel, irrigation_rate, & + catch_calc_zbar, catch_calc_peatclsm_waterlevel, & gndtmp use update_model_para4cn, only : upd_curr_date_time @@ -772,6 +772,42 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_SPR' ,& + LONG_NAME = 'irrigation_flux_sprinkler' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_DRP' ,& + LONG_NAME = 'irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_FRW' ,& + LONG_NAME = 'irrigation_flux_furrow' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_PDY' ,& + LONG_NAME = 'irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) ! !INTERNAL STATE: @@ -1976,117 +2012,11 @@ subroutine SetServices ( GC, RC ) endif -! IRRIGATION MODEL INTERNAL - - IF (RUN_IRRIG /= 0) THEN - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'fraction_of_irrigated_cropland',& - UNITS = '1' ,& - SHORT_NAME = 'IRRIGFRAC' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'fraction_of_paddy_cropland',& - UNITS = '1' ,& - SHORT_NAME = 'PADDYFRAC' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'Maximum_LAI' ,& - UNITS = '1' ,& - SHORT_NAME = 'LAIMAX' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'Minimum_LAI' ,& - UNITS = '1' ,& - SHORT_NAME = 'LAIMIN' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_primary_type' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMPT' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_secondary_type' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMST' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_primary_fraction' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMPF' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'CLM_secondary_fraction' ,& - UNITS = '1' ,& - SHORT_NAME = 'CLMSF' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - - ENDIF - !EOS ! EXPORT STATE: - IF (RUN_IRRIG /= 0) THEN - call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'irrigation_rate' ,& - UNITS = 'kg m-2 s-1' ,& - SHORT_NAME = 'IRRIGRATE' ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - VERIFY_(STATUS) - ENDIF - call MAPL_AddExportSpec(GC, & LONG_NAME = 'evaporation' ,& UNITS = 'kg m-2 s-1' ,& @@ -4550,6 +4480,11 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: ASCATZ0 real, dimension(:), pointer :: NDVI + real, dimension(:), pointer :: IRRG_RATE_SPR + real, dimension(:), pointer :: IRRG_RATE_DRP + real, dimension(:), pointer :: IRRG_RATE_FRW + real, dimension(:), pointer :: IRRG_RATE_PDY + real, dimension(:,:), pointer :: DUDP real, dimension(:,:), pointer :: DUSV real, dimension(:,:), pointer :: DUWT @@ -4679,14 +4614,6 @@ subroutine Driver ( RC ) real, dimension(:,:), pointer :: RBC002 real, dimension(:,:), pointer :: ROC001 real, dimension(:,:), pointer :: ROC002 - real, dimension(:), pointer :: IRRIGFRAC - real, dimension(:), pointer :: PADDYFRAC - real, dimension(:), pointer :: LAIMAX - real, dimension(:), pointer :: LAIMIN - real, dimension(:), pointer :: CLMPT - real, dimension(:), pointer :: CLMST - real, dimension(:), pointer :: CLMPF - real, dimension(:), pointer :: CLMSF real, dimension(:), pointer :: T2M10D real, dimension(:), pointer :: TPREC10D real, dimension(:), pointer :: TPREC60D @@ -4851,7 +4778,6 @@ subroutine Driver ( RC ) real, pointer, dimension(:) :: RMELTBC002 real, pointer, dimension(:) :: RMELTOC001 real, pointer, dimension(:) :: RMELTOC002 - real, pointer, dimension(:) :: IRRIGRATE real, pointer, dimension(:) :: PEATCLSM_WATERLEVEL real, pointer, dimension(:) :: PEATCLSM_FSWCHANGE @@ -4891,7 +4817,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: fveg1, fveg2 real,pointer,dimension(:) :: FICE1TMP real,pointer,dimension(:) :: SLDTOT - + real,pointer,dimension(:) :: PLS_IN ! real*8,pointer,dimension(:) :: fsum real,pointer,dimension(:,:) :: ghtcnt @@ -4996,7 +4922,6 @@ subroutine Driver ( RC ) ! unadulterated TC's and QC's real, pointer :: TC1_0(:), TC2_0(:), TC4_0(:) real, pointer :: QA1_0(:), QA2_0(:), QA4_0(:) - real, pointer :: PLSIN(:) ! CATCHMENT_SPINUP integer :: CurrMonth, CurrDay, CurrHour, CurrMin, CurrSec @@ -5268,6 +5193,10 @@ subroutine Driver ( RC ) call MAPL_GetPointer(IMPORT,SSSV ,'SSSV' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,SSWT ,'SSWT' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,SSSD ,'SSSD' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_SPR,'IRRG_RATE_SPR',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_DRP,'IRRG_RATE_DRP',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_FRW,'IRRG_RATE_FRW',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_PDY,'IRRG_RATE_PDY',RC=STATUS); VERIFY_(STATUS) ! ----------------------------------------------------- ! INTERNAL Pointers @@ -5384,17 +5313,7 @@ subroutine Driver ( RC ) call MAPL_GetPointer(INTERNAL,ROC002 ,'ROC002' , RC=STATUS); VERIFY_(STATUS) endif - IF (catchcn_internal%RUN_IRRIG /= 0) THEN - call MAPL_GetPointer(INTERNAL,IRRIGFRAC ,'IRRIGFRAC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,PADDYFRAC ,'PADDYFRAC' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,LAIMAX ,'LAIMAX' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,LAIMIN ,'LAIMIN' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMPT ,'CLMPT' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMST ,'CLMST' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMPF ,'CLMPF' , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,CLMSF ,'CLMSF' , RC=STATUS); VERIFY_(STATUS) - ENDIF - + ! ----------------------------------------------------- ! EXPORT POINTERS ! ----------------------------------------------------- @@ -5550,9 +5469,7 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,RMELTOC002 ,'RMELTOC002' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,PEATCLSM_WATERLEVEL,'PEATCLSM_WATERLEVEL' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,PEATCLSM_FSWCHANGE ,'PEATCLSM_FSWCHANGE' , RC=STATUS); VERIFY_(STATUS) - - IF (catchcn_internal%RUN_IRRIG /= 0) call MAPL_GetPointer(EXPORT,IRRIGRATE ,'IRRIGRATE' , RC=STATUS); VERIFY_(STATUS) - + NTILES = size(PS) allocate( ityp(ntiles,nveg,nzone) ) @@ -5790,7 +5707,7 @@ subroutine Driver ( RC ) allocate(QA1_0 (NTILES)) allocate(QA2_0 (NTILES)) allocate(QA4_0 (NTILES)) - allocate(PLSIN (NTILES)) + allocate(PLS_IN (NTILES)) call ESMF_VMGetCurrent ( VM, RC=STATUS ) @@ -7324,27 +7241,27 @@ subroutine Driver ( RC ) ! gkw: end of main CN block - PLSIN = PLS + PLS_IN = PLS ! -------------------------------------------------------------------------- - ! Call Irrigation Model + ! Add irrigation model imports ! -------------------------------------------------------------------------- - IF ((catchcn_internal%RUN_IRRIG /= 0).AND.(ntiles >0)) THEN - - 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 (catchcn_internal%IRRIG_METHOD, & - NTILES, AGCM_HH, AGCM_MI, AGCM_S, lons, IRRIGFRAC, PADDYFRAC, & - CLMPT,CLMST, CLMPF, CLMSF, LAIMAX, LAIMIN, LAI0, & - POROS, WPWET, VGWMAX, RZMC, IRRIGRATE) - - PLSIN = PLS + IRRIGRATE - - ENDIF + IF (catchcn_internal%RUN_IRRIG /= 0) THEN + where (IRRG_RATE_SPR > 0) + PLS_IN = PLS_IN + IRRG_RATE_SPR + end where + where (IRRG_RATE_DRP > 0) + RZEXC = RZEXC + IRRG_RATE_DRP * DT + end where + where (IRRG_RATE_PDY > 0) + SRFEXC = SRFEXC + IRRG_RATE_PDY * DT + end where + where (IRRG_RATE_FRW > 0) + RZEXC = RZEXC + IRRG_RATE_FRW * DT + end where + + ENDIF #ifdef DBG_CNLSM_INPUTS call MAPL_Get(MAPL, LocStream=LOCSTREAM, RC=STATUS) @@ -7364,7 +7281,7 @@ subroutine Driver ( RC ) ! Inputs call MAPL_VarWrite(unit, tilegrid, PCU, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, PLS, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, PLS_IN, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, SNO, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, ICE, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, FRZR, mask=mask, rc=status); VERIFY_(STATUS) @@ -7531,7 +7448,7 @@ subroutine Driver ( RC ) call CATCHCN ( NTILES, LONS, LATS, DT,catchcn_internal%USE_FWET_FOR_RUNOFF, & ! LONS, LATS are in [radians] !!! catchcn_internal%FWETC, catchcn_internal%FWETL, cat_id, VEG1,VEG2,FVEG1,FVEG2,DZSF ,& - PCU , PLSIN , SNO, ICE, FRZR ,& + PCU , PLS_IN , SNO, ICE, FRZR ,& UUU ,& EVSBT(:,FSAT), DEVSBT(:,FSAT), DEDTC(:,FSAT) ,& @@ -8135,7 +8052,7 @@ subroutine Driver ( RC ) deallocate( ht ) deallocate( tp ) deallocate( soilice ) - deallocate (PLSIN) + deallocate (PLS_IN) call MAPL_TimerOff ( MAPL, "-CATCHCNCLM45" ) RETURN_(ESMF_SUCCESS) 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 17da8116a..92c620a59 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 @@ -746,6 +746,42 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddImportSpec(GC ,& + LONG_NAME = 'irrigation_flux_sprinkler' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'IRRG_RATE_SPR' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + LONG_NAME = 'irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'IRRG_RATE_DRP' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + LONG_NAME = 'irrigation_flux_furrow' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'IRRG_RATE_FRW' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + LONG_NAME = 'irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'IRRG_RATE_PDY' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + ! !INTERNAL STATE: ! if is_offline, some variables ( in the last) are not required @@ -2633,6 +2669,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + SHORT_NAME = 'SPIRRG', & + LONG_NAME = 'Spurious_irrigation_flux', & + UNITS = 'kg m-2 s-1', & + DIMS = MAPL_DimsTileOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& LONG_NAME = 'vegetation_type' ,& UNITS = '1' ,& @@ -4064,6 +4109,11 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: ASCATZ0 real, dimension(:), pointer :: NDVI + real, dimension(:), pointer :: IRRG_RATE_SPR + real, dimension(:), pointer :: IRRG_RATE_DRP + real, dimension(:), pointer :: IRRG_RATE_FRW + real, dimension(:), pointer :: IRRG_RATE_PDY + real, dimension(:,:), pointer :: DUDP real, dimension(:,:), pointer :: DUSV real, dimension(:,:), pointer :: DUWT @@ -4256,6 +4306,7 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: SPLH real, dimension(:), pointer :: SPWATR real, dimension(:), pointer :: SPSNOW + real, dimension(:), pointer :: SPIRRG real, dimension(:), pointer :: WAT10CM real, dimension(:), pointer :: WATSOI @@ -4323,6 +4374,7 @@ subroutine Driver ( RC ) real,pointer,dimension(:) :: ALWX, BLWX real,pointer,dimension(:) :: FICE1TMP real,pointer,dimension(:) :: SLDTOT + real,pointer,dimension(:) :: PLS_SPR ! real*8,pointer,dimension(:) :: fsum @@ -4398,6 +4450,8 @@ subroutine Driver ( RC ) integer :: NTILES integer :: I, N + real, dimension(:), allocatable :: AR4, werror ! for catch_calc_soil_moist() after irrigation application + ! dummy variables for call to get snow temp real :: FICE @@ -4643,6 +4697,11 @@ subroutine Driver ( RC ) call MAPL_GetPointer(IMPORT,SSWT ,'SSWT' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,SSSD ,'SSSD' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_SPR,'IRRG_RATE_SPR',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_DRP,'IRRG_RATE_DRP',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_FRW,'IRRG_RATE_FRW',RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT,IRRG_RATE_PDY,'IRRG_RATE_PDY',RC=STATUS); VERIFY_(STATUS) + ! ----------------------------------------------------- ! INTERNAL Pointers ! ----------------------------------------------------- @@ -4824,6 +4883,8 @@ subroutine Driver ( RC ) call MAPL_GetPointer(EXPORT,SPLH, 'SPLH' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,SPWATR, 'SPWATR' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT,SPSNOW, 'SPSNOW' , RC=STATUS); VERIFY_(STATUS) + if(CATCH_INTERNAL_STATE%RUN_IRRIG /= 0) & + call MAPL_GetPointer(EXPORT,SPIRRG, 'SPIRRG' ,ALLOC=.true.,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) @@ -4944,7 +5005,7 @@ subroutine Driver ( RC ) allocate(RCONSTIT (NTILES,N_SNOW,N_constit)) allocate(TOTDEPOS (NTILES,N_constit)) allocate(RMELT (NTILES,N_constit)) - + allocate(PLS_SPR (NTILES)) debugzth = .false. ! -------------------------------------------------------------------------- @@ -5533,8 +5594,8 @@ subroutine Driver ( RC ) ! driver ! -------------------------------------------------------------------------- - _ASSERT(count(PLS<0.)==0,'needs informative message') - _ASSERT(count(PCU<0.)==0,'needs informative message') + _ASSERT(count(PLS <0.)==0,'needs informative message') + _ASSERT(count(PCU <0.)==0,'needs informative message') _ASSERT(count(SLDTOT<0.)==0,'needs informative message') LAI0 = max(0.0001 , LAI) @@ -5550,6 +5611,48 @@ subroutine Driver ( RC ) TILEZERO = 0.0 + PLS_SPR = PLS ! PLS_SPR = large-scale precip plus sprinkler irrigation (if present, see below) + + ! -------------------------------------------------------------------------- + ! Add irrigation model imports + ! -------------------------------------------------------------------------- + + if(CATCH_INTERNAL_STATE%RUN_IRRIG /= 0) then + + where (IRRG_RATE_SPR > 0) + PLS_SPR = PLS_SPR + IRRG_RATE_SPR + end where + where (IRRG_RATE_DRP > 0) + RZEXC = RZEXC + IRRG_RATE_DRP * DT + end where + where (IRRG_RATE_FRW > 0) + RZEXC = RZEXC + IRRG_RATE_FRW * DT + end where + where (IRRG_RATE_PDY > 0) + SRFEXC = SRFEXC + IRRG_RATE_PDY * DT + end where + + ! after application of irrigation water, make sure soil moisture prognostics + ! (srfexc, rzexc, catdef) remain valid; + ! werror accounts for excess irrigation that cannot be absorbed by the soil; + ! sfmc, rzmc, prmc here are dummies that are required to get werror + + allocate(ar4( NTILES)) + allocate(werror(NTILES)) + + call catch_calc_soil_moist( & + NTILES, dzsf_in_mm, 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, werror ) + + SPIRRG = -werror / DT ! add excess irrigation into spurious term + + deallocate(ar4) + deallocate(werror) + + endif + call MAPL_TimerOn ( MAPL, "-CATCH" ) #ifdef DBG_CATCH_INPUTS @@ -5568,12 +5671,12 @@ subroutine Driver ( RC ) unit = unit_i ! Inputs - call MAPL_VarWrite(unit, tilegrid, PCU, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, PLS, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, SNO, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, ICE, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, FRZR, mask=mask, rc=status); VERIFY_(STATUS) - call MAPL_VarWrite(unit, tilegrid, UUU, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, PCU, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, PLS_SPR, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, SNO, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, ICE, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, FRZR, mask=mask, rc=status); VERIFY_(STATUS) + call MAPL_VarWrite(unit, tilegrid, UUU, mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, EVSBT (:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) call MAPL_VarWrite(unit, tilegrid, DEVSBT(:,FSAT), mask=mask, rc=status); VERIFY_(STATUS) @@ -5984,7 +6087,7 @@ subroutine Driver ( RC ) DT,CATCH_INTERNAL_STATE%USE_FWET_FOR_RUNOFF ,& CATCH_INTERNAL_STATE%FWETC, CATCH_INTERNAL_STATE%FWETL,& cat_id, VEG, DZSF_in_mm ,& ! cat_id is set to no-data above !!! - PCU , PLS , SNO, ICE, FRZR ,& + PCU , PLS_SPR , SNO, ICE, FRZR ,& UUU ,& EVSBT(:,FSAT), DEVSBT(:,FSAT), DEDTC(:,FSAT) ,& @@ -6500,6 +6603,7 @@ subroutine Driver ( RC ) deallocate(RMELT ) deallocate(FICE1TMP ) deallocate(SLDTOT ) + deallocate(PLS_SPR ) deallocate(FSW_CHANGE) RETURN_(ESMF_SUCCESS) @@ -6534,7 +6638,7 @@ subroutine RUN0(gc, import, export, clock, rc) !! ESMF/MAPL variables type(MAPL_MetaComp), pointer :: MAPL - type(ESMF_State) :: INTERNAL + type(ESMF_State) :: INTERNAL !! IMPORT pointers real, pointer :: ity(:)=>null() diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/CMakeLists.txt new file mode 100644 index 000000000..5bc517e1e --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/CMakeLists.txt @@ -0,0 +1,13 @@ +esma_set_this () + +set (srcs + GEOS_IrrigationGridComp.F90 + irrigation_model.F90 + ) + +include_directories (${INC_ESMF}) + +esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL) +target_include_directories (${this} PUBLIC ${INC_ESMF}) + + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/GEOS_IrrigationGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/GEOS_IrrigationGridComp.F90 new file mode 100644 index 000000000..ef953a7ab --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/GEOS_IrrigationGridComp.F90 @@ -0,0 +1,899 @@ +! $Id$ + +#include "MAPL_Generic.h" + + +!============================================================================= +module GEOS_IrrigationGridCompMod + +!BOP + +! !MODULE: GEOS_Irrigation -- child to the "Land" gridded component. + +!DESCRIPTION: +! {\tt GEOS\_Irrigation} is a gridded component that calculates +! irrigation rates that are (optionally) used in the Catchment[CN] model.\\ +! +! Exports from this routine are the instantaneous values of the +! irrigation rates from 4 different irrigation methods on tilespace: +! 1) drip, 2) sprinkler, 3) furrow, and 4) flood. +! Catchment[CN] (optionally) uses these irrigation rates as inputs to +! the water budget calculations. All imports and exports are stored +! in tile space. Soil parameters and root zone soil moisture are +! imports from the Catchment[CN] gridded component.\\ +! +! Temporally and spatially varying irrigation model parameters are +! from a binary data file.\\ +! +! All internals are static parameters.\\ +! +! IMPORTS: POROS, WPWET, VGWMAX, WCRZ, LAI\\ +! +! EXPORTS: IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, IRRG_RATE_TOT\\ +! +! INTERNALS: IRRG_IRRIGFRAC, IRRG_PADDYFRAC, IRRG_CROPIRRIGFRAC, IRRG_DOY_PLANT, IRRG_DOY_HARVEST, +! IRRG_TYPE, IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW, IRRG_LAIMIN, IRRG_LAIMAX\\ +! +! OPTIONAL INTERNALS: SRATE, DRATE, FRATE\\ +! +! !USES: + + use ESMF + use MAPL + use IRRIGATION_MODULE + + implicit none + private + +! !PUBLIC MEMBER FUNCTIONS: + + public SetServices + + integer :: IRRG_METHOD, IRRG_TRIGGER + integer :: RUN_IRRIG + + type IRRIG_WRAP + type (irrig_params), pointer :: PTR => null() + end type IRRIG_WRAP + +contains + +!BOP + +! !IROUTINE: SetServices -- Sets ESMF services for this component + +! !INTERFACE: + + subroutine SetServices ( GC, RC ) + +! !ARGUMENTS: + + type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component + integer, optional :: RC ! return code + +! !DESCRIPTION: This version uses the MAPL\_GenericSetServices. This function sets +! the Initialize and Finalize services, as well as allocating +! our instance of a generic state and putting it in the +! gridded component (GC). Here we only need to set the run method and +! add the state variable specifications (also generic) to our instance +! of the generic state. This is the way our true state variables get into +! the ESMF\_State INTERNAL, which is in the MAPL\_MetaComp. + +!EOP + +!============================================================================= +! +! ErrLog Variables + + + character(len=ESMF_MAXSTR) :: IAm + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + +! Local derived type aliases + + type(MAPL_MetaComp),pointer :: MAPL=>null() + type(ESMF_Config) :: SCF + character(len=ESMF_MAXSTR) :: SURFRC + + type (irrigation_model), pointer :: IM => null() + type (IRRIG_wrap) :: wrap + +!============================================================================= + +! Begin... + +!------------------------------------------------------------ +! Get my name and set-up traceback handle +!------------------------------------------------------------ + + call ESMF_GridCompGet(GC ,& + NAME=COMP_NAME ,& + RC=STATUS ) + + VERIFY_(STATUS) + + Iam = trim(COMP_NAME) // 'SetServices' + +! ----------------------------------------------------------- +! Get the configuration +! ----------------------------------------------------------- + call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) + VERIFY_(STATUS) + +! ----------------------------------------------------------- +! Get runtime switches +! ----------------------------------------------------------- + + call MAPL_GetResource(MAPL, SURFRC, label='SURFRC:', default='GEOS_SurfaceGridComp.rc', RC=STATUS); VERIFY_(STATUS) + + SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) + + call ESMF_ConfigLoadFile( SCF, SURFRC, rc=status) ; VERIFY_(STATUS) + + call ESMF_ConfigGetAttribute(SCF, label='RUN_IRRIG:' , value=RUN_IRRIG , DEFAULT=0, __RC__ ) + call ESMF_ConfigGetAttribute(SCF, label='IRRG_TRIGGER:', value=IRRG_TRIGGER, DEFAULT=0, __RC__ ) + call ESMF_ConfigGetAttribute(SCF, label='IRRG_METHOD:' , value=IRRG_METHOD , DEFAULT=0, __RC__ ) + + call ESMF_ConfigDestroy (SCF, __RC__) + + ! Leave GEOSirrigation_GridComp if RUN_IRRIG == 0 + if(RUN_IRRIG == 0) then + RETURN_(ESMF_SUCCESS) + endif + + _ASSERT( 0 <= IRRG_TRIGGER .and. IRRG_TRIGGER <= 1, Iam // ' bad value for IRRG_TRIGGER' ) + _ASSERT( 0 <= IRRG_METHOD .and. IRRG_METHOD <= 3, Iam // ' bad value for IRRG_METHOD' ) + +! ----------------------------------------------------------- +! Set the the Initialize and Run entry point +! ----------------------------------------------------------- + + call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS) + VERIFY_(STATUS) + call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_RUN, Run, RC=STATUS) + VERIFY_(STATUS) + +! BOS + +! ----------------------------------------------------------- +! Internal State +! ----------------------------------------------------------- + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_IRRIGFRAC' ,& + LONG_NAME = 'fraction_of_irrigated_cropland' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_PADDYFRAC' ,& + LONG_NAME = 'fraction_of_paddy_cropland' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_CROPIRRIGFRAC' ,& + LONG_NAME = 'Crop_irrigated_fraction' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + UNGRIDDED_DIMS = (/IRRG_NCROPS/) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_DOY_PLANT' ,& + LONG_NAME = 'DOY_start_planting' ,& + UNITS = 'day' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + UNGRIDDED_DIMS = (/IRRG_NSEASONS, IRRG_NCROPS/) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_DOY_HARVEST' ,& + LONG_NAME = 'DOY_end_harvesting' ,& + UNITS = 'day' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + UNGRIDDED_DIMS = (/IRRG_NSEASONS, IRRG_NCROPS/) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_TYPE' ,& + LONG_NAME = 'Preferred_Irrig_method=(0)CONCURRENT_(1)SPRINKLER_(2)DRIP_(3)FLOOD_(negative)AVOID',& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + UNGRIDDED_DIMS = (/IRRG_NCROPS/) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_IRRIGFRAC_SPR' ,& + LONG_NAME = 'fraction_of_sprinkler_irrigation' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_IRRIGFRAC_DRP' ,& + LONG_NAME = 'fraction_of_drip_irrigation' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_IRRIGFRAC_FRW' ,& + LONG_NAME = 'fraction_of_flood_irrigation' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_LAIMIN' ,& + LONG_NAME = 'Minimum_LAI_irrigated_crops' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'IRRG_LAIMAX' ,& + LONG_NAME = 'Maximum_LAI_irrigated_crops' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + + ! NOTE: UNGRIDDED_DIMS for SRATE, DRATE, and FRATE internal specs depends on IRRG_TRIGGER + + if (IRRG_TRIGGER == 0) then + + ! UNGRIDDED_DIMS = 2: irrigated crops (sprinkler/drip/furrow) and paddy (in order) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'SRATE' ,& + LONG_NAME ='crop_specific_irrigation_flux_sprinkler',& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/2/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'DRATE' ,& + LONG_NAME = 'crop_specific_irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/2/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'FRATE' ,& + LONG_NAME = 'crop_specific_irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/2/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + elseif (IRRG_TRIGGER == 1) then + + ! UNGRIDDED_DIMS = 26 crops of crop calendar + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'SRATE' ,& + LONG_NAME ='crop_specific_irrigation_flux_sprinkler',& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/IRRG_NCROPS/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'DRATE' ,& + LONG_NAME = 'crop_specific_irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/IRRG_NCROPS/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC ,& + SHORT_NAME = 'FRATE' ,& + LONG_NAME = 'crop_specific_irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + FRIENDLYTO = trim(COMP_NAME) ,& + RESTART = MAPL_RestartOptional ,& + UNGRIDDED_DIMS = (/IRRG_NCROPS/) ,& + RC=STATUS ) + VERIFY_(STATUS) + + endif + +! ----------------------------------------------------------- +! Export state +! ----------------------------------------------------------- + + call MAPL_AddExportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_SPR' ,& + LONG_NAME = 'irrigation_flux_sprinkler' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_DRP' ,& + LONG_NAME = 'irrigation_flux_drip' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_FRW' ,& + LONG_NAME = 'irrigation_flux_furrow' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_PDY' ,& + LONG_NAME = 'irrigation_flux_paddy' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + SHORT_NAME = 'IRRG_RATE_TOT' ,& + LONG_NAME = 'irrigation_flux_total' ,& + UNITS = 'kg m-2 s-1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + +! ----------------------------------------------------------- +! Import states +! ----------------------------------------------------------- + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'POROS' ,& + LONG_NAME = 'soil_porosity' ,& + UNITS = 'm3 m-3' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'WPWET' ,& + LONG_NAME = 'soil_wilting_point_in_degree_of_saturation_units' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'VGWMAX' ,& + LONG_NAME = 'max_rootzone_water_content' ,& + UNITS = 'kg m-2' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'WCRZ' ,& + LONG_NAME = 'soil_moisture_rootzone' ,& + UNITS = 'm3 m-3' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddImportSpec(GC ,& + SHORT_NAME = 'LAI' ,& + LONG_NAME = 'leaf_area_index' ,& + UNITS = '1' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + +! Allocate this instance of the internal state and put it in wrapper. +! ------------------------------------------------------------------- + + allocate(IM, stat=status ) + VERIFY_(STATUS) + call IM%init_model (SURFRC) + wrap%ptr => IM%irrig_params + +! Save pointer to the wrapped internal state in the GC +! ---------------------------------------------------- + + call ESMF_UserCompSetInternalState ( GC, 'irrigation_state',wrap,status ) + VERIFY_(STATUS) + +! Clocks +!------- + + call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) + VERIFY_(STATUS) + call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) + VERIFY_(STATUS) + +!------------------------------------------------------------ +! Set generic init and final methods +!------------------------------------------------------------ + + call MAPL_GenericSetServices(GC, RC=STATUS) + VERIFY_(STATUS) + + RETURN_(ESMF_SUCCESS) + + end subroutine SetServices + + ! ----------------------------------------------------------- + ! INITIALIZE -- Initialize method for the irrigation component + ! ----------------------------------------------------------- + + subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC ) + + ! ARGUMENTS: + + type(ESMF_GridComp), intent(inout) :: GC + type(ESMF_State), intent(inout) :: IMPORT + type(ESMF_State), intent(inout) :: EXPORT + type(ESMF_Clock), intent(inout) :: CLOCK + integer, optional, intent( out) :: RC + + ! ErrLog Variables + + character(len=ESMF_MAXSTR) :: IAm="Initialize" + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + + ! Locals + + type (MAPL_MetaComp), pointer :: MAPL=>null() + type (ESMF_State ) :: INTERNAL + + ! INTERNAL pointers + + real, dimension(:), pointer :: IRRG_IRRIGFRAC + real, dimension(:), pointer :: IRRG_PADDYFRAC + real, dimension(:,:), pointer :: IRRG_CROPIRRIGFRAC + real, dimension(:,:), pointer :: IRRG_TYPE + real, dimension(:,:), pointer :: SRATE + real, dimension(:,:), pointer :: DRATE + real, dimension(:,:), pointer :: FRATE + +! EXPORT pointers + + real, dimension(:), pointer :: IRRG_RATE_SPR + real, dimension(:), pointer :: IRRG_RATE_DRP + real, dimension(:), pointer :: IRRG_RATE_FRW + real, dimension(:), pointer :: IRRG_RATE_PDY + real, dimension(:), pointer :: IRRG_RATE_TOT + + type(irrigation_model), pointer :: IM + type (IRRIG_wrap) :: wrap + +! Get the target components name and set-up traceback handle. +! ----------------------------------------------------------- + + call ESMF_GridCompGet(GC, name=COMP_NAME, RC=STATUS ) + VERIFY_(STATUS) + + Iam = trim(COMP_NAME) // "Initialize" + + call MAPL_GenericInitialize ( GC, import, export, clock, rc=status ) + VERIFY_(STATUS) + call ESMF_UserCompGetInternalState ( GC, 'irrigation_state',wrap,status ) + VERIFY_(STATUS) + allocate (IM) + IM = irrigation_model() + IM%irrig_params = wrap%ptr + + ! Get my internal MAPL_Generic state + ! ----------------------------------------------------------- + + call MAPL_GetObjectFromGC(GC, MAPL, STATUS) + VERIFY_(STATUS) + + call MAPL_TimerOn(MAPL,"INITIALIZE") + + call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) + VERIFY_(STATUS) + + ! get pointers to internal variables + ! ---------------------------------- + + call MAPL_GetPointer(INTERNAL, IRRG_IRRIGFRAC ,'IRRG_IRRIGFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_PADDYFRAC ,'IRRG_PADDYFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_CROPIRRIGFRAC ,'IRRG_CROPIRRIGFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_TYPE ,'IRRG_TYPE', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, SRATE ,'SRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, DRATE ,'DRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, FRATE ,'FRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + + ! get pointers to EXPORT variable + ! ------------------------------- + call MAPL_GetPointer(EXPORT, IRRG_RATE_SPR, 'IRRG_RATE_SPR', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_DRP, 'IRRG_RATE_DRP', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_FRW, 'IRRG_RATE_FRW', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_PDY, 'IRRG_RATE_PDY', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_TOT, 'IRRG_RATE_TOT', RC=STATUS) ; VERIFY_(STATUS) + + ! Update IRRG_IRRIGFRAC and IRRG_PADDYFRAC for applications that are run on regular tiles in + ! which IRRG_IRRIGFRAC and IRRG_PADDYFRAC in BCs are fractions. + ! The irrigation model would run on tiles with IRRG_IRRIGFRAC + IRRG_PADDYFRAC > IRRG_FRAC_THRES (default is 0.01). + + where (IRRG_IRRIGFRAC + IRRG_PADDYFRAC > IM%irrig_thres) + + ! uncomment the following block to assign the entire cell to the largest fraction: + + ! where (IRRG_PADDYFRAC >= IRRG_IRRIGFRAC) + ! IRRG_PADDYFRAC = 1. + ! IRRG_IRRIGFRAC = 0. + ! elsewhere + ! IRRG_PADDYFRAC = 0. + ! IRRG_IRRIGFRAC = 1. + ! endwhere + + elsewhere + IRRG_PADDYFRAC = 0. + IRRG_IRRIGFRAC = 0. + endwhere + + if (IRRG_TRIGGER == 0) then + + ! LAI based trigger: scale soil moisture to LAI seasonal cycle + ! ============================================================ + + call IM%update_irates (IRRG_RATE_SPR,IRRG_RATE_DRP,IRRG_RATE_FRW,IRRG_RATE_PDY, & + IRRG_IRRIGFRAC,IRRG_PADDYFRAC,SRATE,DRATE,FRATE) + + else + + ! crop calendar based irrigation + ! ============================== + + call IM%update_irates (IRRG_RATE_SPR,IRRG_RATE_DRP,IRRG_RATE_FRW,IRRG_RATE_PDY, & + IRRG_CROPIRRIGFRAC,SRATE,DRATE,FRATE) + + endif + + ! Scale computed IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, and IRRG_RATE_PDY to the total + ! irrigated tile fraction before exporting to Catchment[CN]. + + IRRG_RATE_SPR = IRRG_RATE_SPR * IRRG_IRRIGFRAC + IRRG_RATE_DRP = IRRG_RATE_DRP * IRRG_IRRIGFRAC + IRRG_RATE_FRW = IRRG_RATE_FRW * IRRG_IRRIGFRAC + IRRG_RATE_PDY = IRRG_RATE_PDY * IRRG_PADDYFRAC + + if(associated(IRRG_RATE_TOT)) IRRG_RATE_TOT = IRRG_RATE_SPR + IRRG_RATE_DRP + IRRG_RATE_FRW +IRRG_RATE_PDY + + call MAPL_TimerOff(MAPL,"INITIALIZE") + RETURN_(ESMF_SUCCESS) + + end subroutine INITIALIZE + +! ----------------------------------------------------------- +! RUN -- Run method for the irrigation component +! ----------------------------------------------------------- + + subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) + +! ----------------------------------------------------------- +! !ARGUMENTS: + + type(ESMF_GridComp), intent(inout) :: GC + type(ESMF_State), intent(inout) :: IMPORT + type(ESMF_State), intent(inout) :: EXPORT + type(ESMF_Clock), intent(inout) :: CLOCK + integer, optional, intent( out) :: RC + +! ErrLog Variables + + character(len=ESMF_MAXSTR) :: IAm + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + +! Locals + + type (MAPL_MetaComp), pointer :: MAPL=>null() + type (ESMF_State ) :: INTERNAL + +! INTERNAL pointers + + real, dimension(:), pointer :: IRRG_IRRIGFRAC + real, dimension(:), pointer :: IRRG_PADDYFRAC + real, dimension(:), pointer :: IRRG_IRRIGFRAC_SPR + real, dimension(:), pointer :: IRRG_IRRIGFRAC_DRP + real, dimension(:), pointer :: IRRG_IRRIGFRAC_FRW + real, dimension(:), pointer :: IRRG_LAIMIN + real, dimension(:), pointer :: IRRG_LAIMAX + real, dimension(:,:), pointer :: IRRG_CROPIRRIGFRAC + real, dimension(:,:), pointer :: IRRG_TYPE + real, dimension(:,:,:), pointer :: IRRG_DOY_PLANT + real, dimension(:,:,:), pointer :: IRRG_DOY_HARVEST + real, dimension(:,:), pointer :: SRATE + real, dimension(:,:), pointer :: DRATE + real, dimension(:,:), pointer :: FRATE + +! EXPORT pointers + + real, dimension(:), pointer :: IRRG_RATE_SPR + real, dimension(:), pointer :: IRRG_RATE_DRP + real, dimension(:), pointer :: IRRG_RATE_FRW + real, dimension(:), pointer :: IRRG_RATE_PDY + real, dimension(:), pointer :: IRRG_RATE_TOT + +! IMPORT pointers + + real, dimension(:), pointer :: POROS + real, dimension(:), pointer :: WPWET + real, dimension(:), pointer :: VGWMAX + real, dimension(:), pointer :: WCRZ + real, dimension(:), pointer :: LAI + +! Time attributes + + type(ESMF_Time) :: CURRENT_TIME + integer :: AGCM_YY, AGCM_MM, AGCM_DD, AGCM_MI, AGCM_S, AGCM_HH, dofyr + +! Others/ Locals + + type(irrigation_model), pointer :: IM + type (IRRIG_wrap) :: wrap + real, dimension(:), pointer :: lons + integer :: ntiles, n + real, dimension(:), allocatable :: local_hour, SMWP, SMSAT, SMREF, SMCNT, RZDEF + real :: DT, T1, T2 + + ! Get the target components name and set-up traceback handle. + ! ----------------------------------------------------------- + + call ESMF_GridCompGet(GC, name=COMP_NAME, RC=STATUS ) + VERIFY_(STATUS) + + Iam = trim(COMP_NAME) // "Run" + + call ESMF_UserCompGetInternalState ( GC, 'irrigation_state',wrap,status ) + VERIFY_(STATUS) + allocate (IM) + IM = irrigation_model() + IM%irrig_params = wrap%ptr + + ! Get my internal MAPL_Generic state + ! ----------------------------------------------------------- + + call MAPL_GetObjectFromGC(GC, MAPL, STATUS) + VERIFY_(STATUS) + + call MAPL_Get(MAPL, HEARTBEAT = DT, RC=STATUS) + VERIFY_(STATUS) + + call MAPL_TimerOn(MAPL,"RUN") + + call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) + VERIFY_(STATUS) + + ! get pointers to internal variables + ! ---------------------------------- + + call MAPL_GetPointer(INTERNAL, IRRG_IRRIGFRAC ,'IRRG_IRRIGFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_PADDYFRAC ,'IRRG_PADDYFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_CROPIRRIGFRAC ,'IRRG_CROPIRRIGFRAC', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_DOY_PLANT ,'IRRG_DOY_PLANT', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_DOY_HARVEST ,'IRRG_DOY_HARVEST', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_TYPE ,'IRRG_TYPE', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_IRRIGFRAC_SPR ,'IRRG_IRRIGFRAC_SPR', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_IRRIGFRAC_DRP ,'IRRG_IRRIGFRAC_DRP', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_IRRIGFRAC_FRW ,'IRRG_IRRIGFRAC_FRW', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_LAIMIN ,'IRRG_LAIMIN', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, IRRG_LAIMAX ,'IRRG_LAIMAX', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, SRATE ,'SRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, DRATE ,'DRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, FRATE ,'FRATE', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + + ! get pointers to EXPORT variable + ! ------------------------------- + call MAPL_GetPointer(EXPORT, IRRG_RATE_SPR ,'IRRG_RATE_SPR', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_DRP ,'IRRG_RATE_DRP', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_FRW ,'IRRG_RATE_FRW', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_PDY ,'IRRG_RATE_PDY', ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, IRRG_RATE_TOT ,'IRRG_RATE_TOT', RC=STATUS) ; VERIFY_(STATUS) + + + + ! get pointers to IMPORT variables + ! -------------------------------- + + call MAPL_GetPointer(IMPORT, POROS ,'POROS', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT, WPWET ,'WPWET', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT, VGWMAX ,'VGWMAX', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT, WCRZ ,'WCRZ', RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetPointer(IMPORT, LAI ,'LAI', RC=STATUS) ; VERIFY_(STATUS) + + ! Get time and parameters from local state + ! ---------------------------------------- + + call ESMF_ClockGet ( CLOCK, currTime=CURRENT_TIME, RC=STATUS ) + VERIFY_(STATUS) + + call ESMF_TimeGet ( CURRENT_TIME, YY = AGCM_YY, & + MM = AGCM_MM, & + DD = AGCM_DD, & + H = AGCM_HH, & + M = AGCM_MI, & + S = AGCM_S , & + dayOfYear = dofyr , & + rc=status ) + VERIFY_(STATUS) + + call MAPL_Get (MAPL, TILELONS = LONS, & ! longitude in [radians] + INTERNAL_ESMF_STATE = INTERNAL, RC=STATUS ) + VERIFY_(STATUS) + + + ! call irrigation model + ! --------------------- + + NTILES = SIZE (LONS) + + allocate (local_hour (1:NTILES)) + allocate (SMWP (1:NTILES)) + allocate (SMSAT (1:NTILES)) + allocate (SMREF (1:NTILES)) + allocate (SMCNT (1:NTILES)) + allocate (RZDEF (1:NTILES)) + + ! soil moisture state + SMWP = VGWMAX * WPWET ! RZ soil moisture content at wilting point [mm] + SMSAT = VGWMAX ! RZ soil moisture at saturation [mm] + SMCNT = (VGWMAX/POROS) * WCRZ ! actual RZ soil moisture content [mm] + + DO N = 1, NTILES + + ! local time [hour] + + local_hour(n) = AGCM_HH + AGCM_MI / 60. + AGCM_S / 3600. + 12.* (lons(n)/MAPL_PI) ! longitude in [radians] + IF (local_hour(n) >= 24.) local_hour(n) = local_hour(n) - 24. + IF (local_hour(n) < 0.) local_hour(n) = local_hour(n) + 24. + T1 = CEILING (local_hour(n)) - DT/3600. + T2 = FLOOR (local_hour(n) + 1) + DT/3600. + if((local_hour(n) >= T1).and.(local_hour(n) < T2))then + local_hour(n) = real(NINT(local_hour(n))) + end if + + ! The reference soil moisture content is set to lower tercile of the root zone soil + ! moisture range [mm] to be consistent with ASTRFR = 0.333 used in Catchment[CN]. + ! + ! Note on choice of SMREF: + ! Perhaps, soil field capacity (FIELDCAP) is the desired parameter here - the upper limit + ! of water content that the soil can hold for plants after excess water drained off downward + ! quickly. In future, could switch to FIELDCAP, which has already been derived + ! on tiles and is available in the irrigation_IMxJM_DL.dat file. + + SMREF (n) = VGWMAX (n) * (wpwet (n) + (1. - wpwet (n))/ 3.) + + ! rootzone moisture deficit to reach complete soil saturation for paddy [mm] + + RZDEF (n) = MAX(SMSAT(n) - SMCNT(n), 0.) + + END DO + + if (IRRG_TRIGGER == 0) then + + ! LAI based trigger: scale soil moisture to LAI seasonal cycle + ! ============================================================ + + call IM%run_model(IRRG_METHOD, local_hour, & + IRRG_IRRIGFRAC, IRRG_PADDYFRAC, IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW, & + SMWP,SMSAT,SMREF,SMCNT, LAI, IRRG_LAIMIN, IRRG_LAIMAX, RZDEF, & + IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, & + SRATE, DRATE, FRATE) + + else + + ! crop calendar based irrigation + ! ============================== + + call IM%run_model (dofyr,local_hour, & + IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW, & + IRRG_CROPIRRIGFRAC,IRRG_DOY_PLANT,IRRG_DOY_HARVEST,IRRG_TYPE , & + SMWP,SMSAT,SMREF,SMCNT, RZDEF, & + IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, & + SRATE, DRATE, FRATE) + + endif + + ! Scale computed IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, and IRRG_RATE_PDY to the total + ! irrigated tile fraction before exporting to Catchment[CN]. + + IRRG_RATE_SPR = IRRG_RATE_SPR * IRRG_IRRIGFRAC + IRRG_RATE_DRP = IRRG_RATE_DRP * IRRG_IRRIGFRAC + IRRG_RATE_FRW = IRRG_RATE_FRW * IRRG_IRRIGFRAC + IRRG_RATE_PDY = IRRG_RATE_PDY * IRRG_PADDYFRAC + + if(associated(IRRG_RATE_TOT)) IRRG_RATE_TOT = IRRG_RATE_SPR + IRRG_RATE_DRP + IRRG_RATE_FRW +IRRG_RATE_PDY + + deallocate (local_hour, SMWP, SMSAT, SMREF, SMCNT, RZDEF, IM) + + call MAPL_TimerOff(MAPL,"RUN") + RETURN_(ESMF_SUCCESS) + + end subroutine RUN + +end module GEOS_IrrigationGridCompMod + + +! ========================= EOF ================================================================================== diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/irrigation_model.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/irrigation_model.F90 new file mode 100644 index 000000000..f7fd51517 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/irrigation_model.F90 @@ -0,0 +1,766 @@ +#include "MAPL_Generic.h" + +MODULE IRRIGATION_MODULE + + USE MAPL + USE ESMF + + IMPLICIT NONE + + ! Irrigation Module + ! + ! First Version: March 21, 2021 (Sarith Mahanama) + ! Second Version: June 25, 2024 (Stefano Casirati) + ! + ! This module computes irrigation rates by 4 different methods: sprinkler, flood, furrow, and drip. + ! The computed irrigation rates (exports) are imports to Catchment[CN], where the irrigation water + ! is added to water balance as follows: + ! 1) sprinkler irrigation rate added to large scale precipitation (for irrigated tiles fractions); + ! 2) drip irrigation volume added to rootzone excess (for irrigated tiles fractions); + ! 3) furrow irrigation volume added to rootzone excess (for irrigated tiles fractions); + ! 4) flood irrigation volume added to surface excess (only for paddy tiles fractions). + + ! The model uses the rootzone soil moisture state at the local start time of irrigation to compute + ! irrigation rates for the day and maintains the same rate throughout the irrigation duration. + ! + ! Sprinkler and Flood/Furrow Irrigation methods were adapted from LIS CLSMF2.5 irrigation module + ! (Rodell et al., 2024 (Under Review) + ! Drip irrigation method calculation is similar to that of sprinkler, albeit the drip irrigation + ! method assumes a 10% water loss. (Source FAO) + ! + ! (1) EXPORTS - MODEL OUTPUTS TO THE LAND MODEL (IRRIGATION RATES): + ! 1) IRRG_RATE_SPR [kg m-2 s-1] + ! 2) IRRG_RATE_DRP [kg m-2 s-1] + ! 3) IRRG_RATE_FRW [kg m-2 s-1] + ! 4) IRRG_RATE_PDY [kg m-2 s-1] + ! 5) IRRG_RATE_TOT [kg m-2 s-1] (diagnostic only) + ! + ! (2) IRRIGATED AND PADDY TILES: + ! During land BC's generation, the fraction of irrigated crops and paddy is set to zero + ! if their sum is below an irrigation threshold (default 1%). + ! Irrigated fractions can be irrigated with sprinkler, drip, and furrow, + ! while paddy fractions can only be irrigated using the flood irrigation method. + ! Vegetation characteristics and vegetation dynamic parameters for irrigated + ! crops and paddy tiles were taken from the nearest grass or cropland tile. + ! + ! (3) MODEL INPUTS: + ! SMWP : rootzone soil moisture content at wilting point [mm] + ! SMSAT : rootzone soil moisture content at saturation [mm] + ! SMREF : rootzone soil moisture is at lower tercile of RZ soil moisture range [mm] + ! SMCNT : currrent rootzone soil moisture content [mm] + ! RZDEF : rootzone moisture deficit to reach complete soil saturation for paddy [mm] + ! LOCAL_HOUR to set irrigation switch. + ! + ! (4) SEASONAL CYLCE OF CROP WATER DEMAND: + ! The module provides 2 options to determine the seasonal cycle of crop water demand: + ! 4.1) IRRG_TRIGGER: 0 - SUBROUTINE irrigrate_lai_trigger + ! The LAI-based trigger (Default and the current LIS implementation) + ! uses precomputed minimum and maximum LAI on irrigateed pixels to determine + ! beginning and end of crop growing seasons. + ! + ! This LAI-based trigger is also equipped with an additional control parameter, IRRG_METHOD, + ! which is good to choose the method of irrigation that would run on corresponding fraction: + ! 0: (Default) All 4 methods (sprinkler/drip/furrow/paddy) concurrently, according to specified area fracs + ! 1: Only sprinkler irrigation on *entire* tile (regardless of IRRIGFRAC or PADDYFRAC) + ! 2: Only drip irrigation on *entire* tile (regardless of IRRIGFRAC or PADDYFRAC) + ! 3: Only furrow/flood irrigation on *entire* tile (regardless of IRRIGFRAC or PADDYFRAC) + ! + ! IRRG_TRIGGER: 0 SPECIFIC INPUTS: + ! IRRG_IRRIGFRAC : fraction of tile covered by sprinkler/drip/furrow-irrigated crops; + ! ranges between 0 and 1 (if IRRG_IRRIGFRAC + IRRG_PADDYFRAC > Irrigation Threshold) + ! IRRG_PADDYFRAC : fraction of tile covered by paddy; + ! ranges between 0 and 1 (if IRRG_IRRIGFRAC + IRRG_PADDYFRAC > Irrigation Threshold) + ! + ! Within IRRIGFRAC, allocation by irrigation method: + ! + ! Note: IRRG_IRRIGFRAC_SPR + IRRG_IRRIGFRAC_DRP + IRRG_IRRIGFRAC_FRW = 1. + ! + ! IRRG_IRRIGFRAC_SPR : fraction of IRRG_IRRIGFRAC equipped for sprinkler irrigation + ! IRRG_IRRIGFRAC_DRP : fraction of IRRG_IRRIGFRAC equipped for drip irrigation + ! IRRG_IRRIGFRAC_FRW : fraction of IRRG_IRRIGFRAC equipped for flood/furrow irrigation + ! + ! LAI : time varying Leaf Area Index from the model + ! IRRG_LAIMIN : Minimum LAI spatially averaged over the irrigated tile fraction + ! IRRG_LAIMAX : Maximum LAI spatially averaged over the irrigated tile fraction + ! + ! 4.2) IRRG_TRIGGER: 1 - SUBROUTINE irrigrate_crop_calendar + ! Uses 26 crop calendars based on monthly crop growing areas of below crops. + ! 1 Wheat 14 Oil palm + ! 2 Maize 15 Rape seed / Canola + ! 3 Rice 16 Groundnuts / Peanuts + ! 4 Barley 17 Pulses + ! 5 Rye 18 Citrus + ! 6 Millet 19 Date palm + ! 7 Sorghum 20 Grapes / Vine + ! 8 Soybeans 21 Cotton + ! 9 Sunflower 22 Cocoa + ! 10 Potatoes 23 Coffee + ! 11 Cassava 24 Others perennial + ! 12 Sugar cane 25 Fodder grasses + ! 13 Sugar beet 26 Others annual + ! + ! IRRG_TRIGGER: 1 SPECIFIC INPUTS: + ! DOFYR : day of year + ! IRRG_TYPE : Preferred Irrig method (NTILES, 26) - + ! 0 CONCURRENT (default), + ! 1 SPRINKLER ONLY + ! 2 DRIP ONLY + ! 3 FLOOD/FURROW ONLY, and + ! <0 AVOID this method + ! IRRG_CROPIRRIGFRAC: Crop irrigated fraction (NTILES, 26) (per Section 2, fractions have been + ! adjusted such that IRRG_CROPIRRIGFRAC=1. on paddy tiles; the sum of available + ! crop fractions equals 1. on irrigated crop tiles; + ! and is zero on non-irrigated tiles. + ! IRRG_DOY_PLANT : DOY start planting (NTILES, 2, 26) - up to two seasons + ! IRRG_DOY_HARVEST : DOY end harvesting (NTILES, 2, 26) - up to two seasons + ! If IRRG_DOY_PLANT = IRRG_DOY_HARVEST = 998, the crop is not grown on that tile + ! + ! (5) MODEL UPDATES (OPTIONAL INTERNALS): + ! SRATE, DRATE, and FRATE contain irrigation rates applied on individual fractions at any given time. + ! The second dimensions of 2D arrays is for different crop fractions i.e. the second dimension is 2 for above + ! IRRG_TRIGGER: 0 to separately store irrigation rates in irrigated crop and paddy fractions. + ! It would be 26 for IRRG_TRIGGER: 1. + ! The crop calendar implemetation (IRRG_TRIGGER: 1) computes IRRG_RATE_SPR, IRRG_RATE_DRP, + ! IRRG_RATE_FRW, and IRRG_RATE_PDY as weighted averages of irrigation rates from + ! all active crops in SRATE, DRATE and FRATE arrays. + + PRIVATE + + INTEGER, PARAMETER, PUBLIC :: IRRG_NCROPS = 26, IRRG_NSEASONS = 2 + + type, public :: irrig_params + + ! Below parameters can be set via RC file. + + ! IRRGRR: Do we really want to hardwire defaults here *and* in GEOS_SurfaceGridComp.rc ??? + + REAL :: irrig_thres = 0.01 ! threshold of tile fraction to turn the irrigation model on. + REAL :: lai_thres = 0.6 ! threshold of LAI range to turn irrigation on + REAL :: efcor = 25.0 ! efficiency correction (% water loss: efcor = 0% denotes 100% efficient water use) + REAL :: MIDS_LNGTH = 0.6 ! Mid-season length as a fraction of crop growing season length (to be used with IRRG_TRIGGER: 1) + + ! Sprinkler parameters + ! -------------------- + REAL :: sprinkler_stime = 6.0 ! sprinkler irrigation start time (local) [hours] + REAL :: sprinkler_dur = 4.0 ! sprinkler irrigation duration [hours] + REAL :: sprinkler_thres = 0.7 ! threshold for soil moisture availability ("MA") to trigger sprinkler irrigation [dim-less] + + ! Drip parameters + ! --------------- + REAL :: drip_stime = 8.0 ! drip irrigation start time (local) [hours] + REAL :: drip_dur = 8.0 ! drip irrigation duration [hours] + REAL :: drip_thres = 0.7 ! threshold for soil moisture availability ("MA") to trigger drip irrigation [dim-less] + + ! Flood parameters + ! ---------------- + REAL :: flood_stime = 6.0 ! flood irrigation start time (local) [hours] + REAL :: flood_dur = 8.0 ! flood irrigation duration [hours] + REAL :: flood_thres = 0.6 ! threshold for soil moisture availability ("MA") to trigger flood irrigation [dim-less] + + + + end type irrig_params + + type, public, extends (irrig_params) :: irrigation_model + + contains + + ! public + procedure, public :: init_model + generic, public :: run_model => irrigrate_lai_trigger, irrigrate_crop_calendar + generic, public :: update_irates => update_irates_lai, update_irates_ccalendar + + ! private + procedure, private :: irrigrate_lai_trigger + procedure, private :: irrigrate_crop_calendar + procedure, private :: cwd => crop_water_deficit + procedure, private :: irrig_by_method + procedure, private :: update_irates_lai + procedure, private :: update_irates_ccalendar + + end type irrigation_model + +contains + + ! ---------------------------------------------------------------------------- + + SUBROUTINE init_model (IP, SURFRC) + + implicit none + + class (irrigation_model), intent(inout) :: IP + CHARACTER(*), INTENT(IN) :: SURFRC + + type(irrig_params) :: DP + type(ESMF_Config) :: SCF + + integer :: status, RC + character(len=ESMF_MAXSTR) :: Iam + + Iam='IRRIGATION_MODULE: init_model' + + SCF = ESMF_ConfigCreate(__RC__) + + CALL ESMF_ConfigLoadFile (SCF,SURFRC,rc=status) ; VERIFY_(STATUS) + + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_SPR_STIME:' , VALUE=IP%sprinkler_stime, DEFAULT=DP%sprinkler_stime, __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_SPR_DUR:' , VALUE=IP%sprinkler_dur, DEFAULT=DP%sprinkler_dur , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_SPR_THRES:' , VALUE=IP%sprinkler_thres, DEFAULT=DP%sprinkler_thres, __RC__ ) + + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_DRP_STIME:' , VALUE=IP%drip_stime, DEFAULT=DP%drip_stime , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_DRP_DUR:' , VALUE=IP%drip_dur, DEFAULT=DP%drip_dur , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_DRP_THRES:' , VALUE=IP%drip_thres, DEFAULT=DP%drip_thres , __RC__ ) + + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_FLD_STIME:' , VALUE=IP%flood_stime, DEFAULT=DP%flood_stime , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_FLD_DUR:' , VALUE=IP%flood_dur, DEFAULT=DP%flood_dur , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_FLD_THRES:' , VALUE=IP%flood_thres, DEFAULT=DP%flood_thres , __RC__ ) + + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_EFCOR:' , VALUE=IP%efcor, DEFAULT=DP%efcor , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_LAI_THRES:' , VALUE=IP%lai_thres, DEFAULT=DP%lai_thres , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_MIDS_LNGTH:', VALUE=IP%MIDS_LNGTH, DEFAULT=DP%MIDS_LNGTH , __RC__ ) + CALL ESMF_ConfigGetAttribute (SCF, label='IRRG_FRAC_THRES:', VALUE=IP%irrig_thres, DEFAULT=DP%irrig_thres , __RC__ ) + CALL ESMF_ConfigDestroy (SCF, __RC__) + + END SUBROUTINE init_model + + ! ---------------------------------------------------------------------------- + + SUBROUTINE irrigrate_lai_trigger (this,IRRG_METHOD, local_hour, & + IRRG_IRRIGFRAC, IRRG_PADDYFRAC, IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW, & + SMWP, SMSAT, SMREF, SMCNT, LAI, IRRG_LAIMIN,IRRG_LAIMAX, RZDEF, & + IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, & + SRATE, DRATE, FRATE) + + implicit none + + class (irrigation_model), intent(inout) :: this + + integer, intent (in) :: IRRG_METHOD + + real, dimension (:), intent (in) :: local_hour + real, dimension (:), intent (in) :: IRRG_IRRIGFRAC, IRRG_PADDYFRAC, IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW + real, dimension (:), intent (in) :: SMWP, SMSAT, SMREF, SMCNT, LAI, IRRG_LAIMIN, IRRG_LAIMAX, RZDEF + + real, dimension (:), intent (inout) :: IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY + real, dimension (:,:), intent (inout) :: SRATE, DRATE, FRATE + + ! local variables + + INTEGER :: NTILES, N, crop + REAL :: ma, H1, H2, HC, IT, ROOTFRAC, LAITHRES + logical :: season_end + + integer :: status, RC ! required for _ASSERT macro + character(len=ESMF_MAXSTR) :: Iam + + ! ------------------------------------------------------ + + Iam='IRRIGATION_MODULE: irrigrate_lai_trigger' + + NTILES = SIZE (IRRG_IRRIGFRAC) + + TILE_LOOP : DO N = 1, NTILES + IF(IRRG_LAIMAX (N) > IRRG_LAIMIN (N)) THEN + LAITHRES = IRRG_LAIMIN (N) + this%lai_thres * (IRRG_LAIMAX (N) - IRRG_LAIMIN (N)) + ROOTFRAC = MIN((LAI(N) - IRRG_LAIMIN (N)) / (IRRG_LAIMAX(N) - IRRG_LAIMIN(N)) ,1.0) + ELSE + ROOTFRAC = 0. + ENDIF + HC = local_hour(n) + + season_end = .true. + + CHECK_LAITHRES : IF (LAI(N) >= LAITHRES) THEN + season_end = .false. + CHECK_IRRIGFRACS: IF ((IRRG_IRRIGFRAC(N) > 0.).OR.(IRRG_PADDYFRAC(N)>0.)) THEN + + !----------------------------------------------------------------------------- + ! Get the rootzone moisture availability to the plant + !----------------------------------------------------------------------------- + + if (IRRG_IRRIGFRAC(N) > 0.) then + + if(SMREF(N) > SMWP(N))then + ma = (SMCNT(N) - SMWP(N)) /(SMREF(N) - SMWP(N)) + else + ma = -1. + endif + + if(ma >= 0) then + + SELECT CASE (IRRG_METHOD) + + CASE (0) ! CONCURRENTLY SPRINKER + DRIP + FURROW + PADDY on corresponding fractions + + call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), & + SRATE = SRATE (N,1), & + DRATE = DRATE (N,1), & + FRATE = FRATE (N,1)) + + SRATE (N,1) = SRATE (N,1)*IRRG_IRRIGFRAC_SPR(N) + DRATE (N,1) = DRATE (N,1)*IRRG_IRRIGFRAC_DRP(N) + FRATE (N,1) = FRATE (N,1)*IRRG_IRRIGFRAC_FRW(N) + + CASE (1) ! SPRINKLER only + + call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), & + SRATE = SRATE (N,1)) + + DRATE (N,1) = 0. + FRATE (N,1) = 0. + + CASE (2) ! DRIP only + + call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), & + DRATE = DRATE (N,1)) + + SRATE (N,1) = 0. + FRATE (N,1) = 0. + + CASE (3) ! FURROW only + + call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), & + FRATE = FRATE (N,1)) + + SRATE (N,1) = 0. + DRATE (N,1) = 0. + + CASE DEFAULT + + _ASSERT( .FALSE., 'irrigrate_lai_trigger(): IRRG_METHOD can only be 0, 1, 2, or 3') + + END SELECT + endif + endif + + if (IRRG_PADDYFRAC (N) > 0.) then + + H1 = this%flood_stime + H2 = this%flood_stime + this%flood_dur + if ((HC >= H1).AND.(HC < H2)) then + ! use RZDEF at H1 during H1 <= HC < H2 to compute irrigrate for paddy. + if(abs(H1 - HC) < 1./3600.) FRATE (N,2) = RZDEF(N) *0.1/(H2 - H1)/ 3600. + else + FRATE (N,2) = 0. + endif + SRATE (N,2) = 0. + DRATE (N,2) = 0. + endif + + ELSE + + SRATE (N,:) = 0. + DRATE (N,:) = 0. + FRATE (N,:) = 0. + + ENDIF CHECK_IRRIGFRACS + ENDIF CHECK_LAITHRES + + ! turn off irrigation if LAI is smaller than the LAI trigger marking end of the season + if(season_end) then + DO crop = 1, 2 ! With LAI trigger Crop 1 = Irrigated Fraction, Crop 2 = Paddy Fraction + SRATE (N,crop) = 0. + DRATE (N,crop) = 0. + FRATE (N,crop) = 0. + END DO + endif + + END DO TILE_LOOP + + ! Update IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY EXPORTS to be sent to land models. + + call this%update_irates( IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, & + IRRG_IRRIGFRAC, IRRG_PADDYFRAC, SRATE, DRATE, FRATE ) + + END SUBROUTINE irrigrate_lai_trigger + + ! ---------------------------------------------------------------------------- + + SUBROUTINE irrigrate_crop_calendar(this,dofyr,local_hour, & + IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW, & + IRRG_CROPIRRIGFRAC,IRRG_DOY_PLANT, IRRG_DOY_HARVEST, IRRG_TYPE , & + SMWP,SMSAT,SMREF,SMCNT, RZDEF, & + IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, SRATE, DRATE, FRATE) + + implicit none + + class(irrigation_model), intent(inout) :: this + integer, intent(in) :: dofyr + real, dimension(:), intent(in) :: local_hour, IRRG_IRRIGFRAC_SPR, IRRG_IRRIGFRAC_DRP, IRRG_IRRIGFRAC_FRW + real, dimension(:), intent(in) :: SMWP, SMSAT, SMREF, SMCNT, RZDEF + real, dimension(:,:), intent(in) :: IRRG_CROPIRRIGFRAC ! IRRG_NCROPS + real, dimension(:,:), intent(in) :: IRRG_TYPE ! IRRG_NCROPS + real, dimension(:,:,:), intent(in) :: IRRG_DOY_PLANT ! IRRG_NSEASONS, IRRG_NCROPS + real, dimension(:,:,:), intent(in) :: IRRG_DOY_HARVEST ! IRRG_NSEASONS, IRRG_NCROPS + real, dimension(:), intent(inout) :: IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY + real, dimension(:,:), intent(inout) :: SRATE, DRATE, FRATE + + ! local variables + + INTEGER :: NTILES, N, crop, sea, ITYPE, I + REAL :: ma, H1, H2, HC, IT, ROOTFRAC, void_frac + logical :: season_end (IRRG_NCROPS) + + integer :: status, RC ! required for _ASSERT macro + character(len=ESMF_MAXSTR) :: Iam + + ! ---------------------------------------------------------- + + Iam='IRRIGATION_MODULE: irrigrate_crop_calendar' + + NTILES = SIZE (local_hour) + + TILE_LOOP : DO N = 1, NTILES + HC = local_hour(n) + IF_IRR: if(SUM(IRRG_CROPIRRIGFRAC(N,:)) > 0.) then + ! the tile is irrigated crop or paddy + season_end = .true. + CROP_LOOP: DO crop = 1, IRRG_NCROPS + CROP_IN_TILE: if(IRRG_CROPIRRIGFRAC(N,crop) > 0.) then + ! crop is grown in this tile + TWO_SEASONS: do sea = 1, IRRG_NSEASONS + IS_CROP: IF(IRRG_DOY_PLANT(N, sea, crop) /= 998) THEN + ! crop is grown in sea + IS_SEASON: IF(IS_WITHIN_SEASON(dofyr,NINT(IRRG_DOY_PLANT(N, sea, crop)),NINT(IRRG_DOY_HARVEST(N, sea, crop)))) THEN + ! dofyr falls within the crop season + season_end(crop) = .false. + PADDY_OR_CROP: if (crop == 3) then + ! PADDY TILE + H1 = this%flood_stime + H2 = this%flood_stime + this%flood_dur + if ((HC >= H1).AND.(HC < H2)) then + ! use RZDEF at H1 during H1 <= HC < H2 to compute irrigrate. + if(abs(H1 - HC) < 1./3600.) FRATE (N,crop) = RZDEF(N) /(H2 - H1)/ 3600. + else + FRATE (N,crop) = 0. + endif + SRATE (N,crop) = 0. + DRATE (N,crop) = 0. + + else + + ! IRRIGATED CROP: compute sum of irrigrates from 25 crops. + + ROOTFRAC = CROP_SEASON_STAGE (this%MIDS_LNGTH, dofyr,NINT(IRRG_DOY_PLANT(N, sea, crop)),NINT(IRRG_DOY_HARVEST(N, sea, crop))) + if(SMREF(N) > SMWP(N))then + ma = (SMCNT(N) - SMWP(N)) /(SMREF(N) - SMWP(N)) + else + ma = -1. + endif + + SOILM: if(ma >= 0) then + + ITYPE = NINT(IRRG_TYPE(N,crop)) + + CROP_IMETHOD: if (ITYPE == 0) then + + ! concurrently on sprinkler, drip and flood fractions + call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), & + SRATE = SRATE (N,crop), & + DRATE = DRATE (N,crop), & + FRATE = FRATE (N,crop)) + + SRATE (N,crop) = SRATE (N,crop)*IRRG_IRRIGFRAC_SPR(N) + DRATE (N,crop) = DRATE (N,crop)*IRRG_IRRIGFRAC_DRP (N) + FRATE (N,crop) = FRATE (N,crop)*IRRG_IRRIGFRAC_FRW (N) + + elseif (ITYPE > 0) then + + ! only this method + if (ITYPE == 1) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), SRATE = SRATE (N,crop)) + if (ITYPE == 2) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), DRATE = DRATE (N,crop)) + if (ITYPE == 3) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), FRATE = FRATE (N,crop)) + + elseif (ITYPE < 0) then + + ! crop does not use IRRG_METHOD -(ITYPE) + void_frac = 0. + DO I = 1,3 + if(I == ABS(ITYPE))then + ! this itype isn't used by this crop other 2 fractions equally share this fraction + if (I == 1) then + void_frac = IRRG_IRRIGFRAC_SPR(N)/2. + SRATE(N,crop) = 0. + endif + if (I == 2) then + void_frac = IRRG_IRRIGFRAC_DRP (N)/2. + DRATE(N,crop) = 0. + endif + if (I == 3)then + void_frac = IRRG_IRRIGFRAC_FRW (N)/2. + FRATE(N,crop) = 0. + endif + else + if (I == 1) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), SRATE = SRATE (N,crop)) + if (I == 2) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), DRATE = DRATE (N,crop)) + if (I == 3) call this%irrig_by_method (HC, ma, ROOTFRAC, SMCNT(N), SMREF(N), FRATE = FRATE (N,crop)) + endif + END DO + DO I = 1,3 + if(I /= ABS(ITYPE))then + if (I == 1) SRATE (N,crop) = SRATE (N,crop)*(IRRG_IRRIGFRAC_SPR(N) + void_frac) + if (I == 2) DRATE (N,crop) = DRATE (N,crop)*(IRRG_IRRIGFRAC_DRP (N) + void_frac) + if (I == 3) FRATE (N,crop) = FRATE (N,crop)*(IRRG_IRRIGFRAC_FRW (N) + void_frac) + endif + ENDDO + + endif CROP_IMETHOD + endif SOILM + ENDIF PADDY_OR_CROP + ENDIF IS_SEASON + end IF IS_CROP + end do TWO_SEASONS + endif CROP_IN_TILE + END DO CROP_LOOP + + ! turn off irrigation for crops that ended the season + DO crop = 1, IRRG_NCROPS + if(season_end(crop)) then + SRATE (N,crop) = 0. + DRATE (N,crop) = 0. + FRATE (N,crop) = 0. + endif + END DO + + endif IF_IRR + END DO TILE_LOOP + + ! Update IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY EXPORTS to be sent to land models + ! They are weighted averaged over 26 crop fractions. + + call this%update_irates( IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY, & + IRRG_CROPIRRIGFRAC, SRATE, DRATE, FRATE ) + + END SUBROUTINE irrigrate_crop_calendar + + ! ---------------------------------------------------------------------------- + + SUBROUTINE update_irates_lai (this,IRRG_RATE_SPR,IRRG_RATE_DRP,IRRG_RATE_FRW,IRRG_RATE_PDY, & + IRRG_IRRIGFRAC,IRRG_PADDYFRAC,SRATE,DRATE,FRATE) + + implicit none + + class(irrigation_model), intent(inout) :: this + + real, dimension(:), intent(in) :: IRRG_IRRIGFRAC, IRRG_PADDYFRAC + real, dimension(:,:), intent(in) :: SRATE, DRATE, FRATE + real, dimension(:), intent(inout) :: IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY + + ! local variables + + integer :: N, NT + + integer :: status, RC ! required for _ASSERT macro + character(len=ESMF_MAXSTR) :: Iam + + Iam='IRRIGATION_MODULE: update_irates_lai' + + ! INITIALIZE EXPORTS + IRRG_RATE_SPR = 0. + IRRG_RATE_DRP = 0. + IRRG_RATE_FRW = 0. + IRRG_RATE_PDY = 0. + + NT = size (IRRG_IRRIGFRAC) + + _ASSERT( size(SRATE,2)==IRRG_NCROPS, 'SRATE dimension mismatch') + + DO N = 1, NT + IF ((IRRG_IRRIGFRAC(N) + IRRG_PADDYFRAC(N)) > 0.) THEN + IRRG_RATE_SPR (N) = SRATE (N,1) + IRRG_RATE_DRP (N) = DRATE (N,1) + IRRG_RATE_FRW (N) = FRATE (N,1) + IRRG_RATE_PDY (N) = FRATE (N,2) + ENDIF + END DO + + END SUBROUTINE update_irates_lai + + !............................................................................... + + SUBROUTINE update_irates_ccalendar(this,IRRG_RATE_SPR,IRRG_RATE_DRP,IRRG_RATE_FRW,IRRG_RATE_PDY, & + IRRG_CROPIRRIGFRAC,SRATE,DRATE,FRATE) + + implicit none + + class(irrigation_model), intent(inout) :: this + + real, dimension(:,:), intent (in) :: IRRG_CROPIRRIGFRAC ! IRRG_NCROPS + real, dimension(:,:), intent (in) :: SRATE, DRATE, FRATE + real, dimension(:), intent (inout) :: IRRG_RATE_SPR, IRRG_RATE_DRP, IRRG_RATE_FRW, IRRG_RATE_PDY + + ! local variables + + integer :: N, NT, crop + + integer :: status, RC ! required for _ASSERT macro + character(len=ESMF_MAXSTR) :: Iam + + Iam='IRRIGATION_MODULE: update_irates_ccalendar' + + ! INITIALIZE EXPORTS + IRRG_RATE_SPR = 0. + IRRG_RATE_DRP = 0. + IRRG_RATE_FRW = 0. + IRRG_RATE_PDY = 0. + + _ASSERT( size (SRATE,2)==IRRG_NCROPS, 'SRATE dimension mismatch') + + NT = size (IRRG_RATE_SPR) + DO N = 1, NT + if(SUM(IRRG_CROPIRRIGFRAC(N,:)) > 0.) then + DO crop = 1, IRRG_NCROPS + IRRG_RATE_SPR(N) = IRRG_RATE_SPR(N) + SRATE (N,crop)*IRRG_CROPIRRIGFRAC(N,crop)/SUM(IRRG_CROPIRRIGFRAC(N,:)) + IRRG_RATE_DRP(N) = IRRG_RATE_DRP(N) + DRATE (N,crop)*IRRG_CROPIRRIGFRAC(N,crop)/SUM(IRRG_CROPIRRIGFRAC(N,:)) + if (crop==3) then + ! If crop is rice (crop ==3) then use flood irrigation. Otherwise use furrow irrigation. + IRRG_RATE_PDY(N) = IRRG_RATE_PDY(N) + FRATE (N,crop)*IRRG_CROPIRRIGFRAC(N,crop)/SUM(IRRG_CROPIRRIGFRAC(N,:)) + else + IRRG_RATE_FRW(N) = IRRG_RATE_FRW(N) + FRATE (N,crop)*IRRG_CROPIRRIGFRAC(N,crop)/SUM(IRRG_CROPIRRIGFRAC(N,:)) + endif + END DO + endif + END DO + + END SUBROUTINE update_irates_ccalendar + + ! ---------------------------------------------------------------------------- + + SUBROUTINE irrig_by_method (this, HC, ma, ROOTFRAC, SMCNT, SMREF, SRATE, DRATE, FRATE) + + implicit none + class (irrigation_model), intent(inout) :: this + REAL, intent(in) :: HC, ma, ROOTFRAC,SMCNT, SMREF + REAL, optional, intent(inout) :: SRATE, DRATE, FRATE + REAL :: H1, H2, IT + + if(present (SRATE)) then + ! SPRINKLER IRRIGATION + H1 = this%sprinkler_stime + H2 = this%sprinkler_stime + this%sprinkler_dur + IT = this%sprinkler_thres + if ((HC >= H1).AND.(HC < H2)) then + ! The model uses rootzone soil moisture state at H1 to compute irrigation + ! rates for the day and maintains the same rate through out the irrigation + ! duration (H1 <= HC < H2). + if((ma <= IT).AND.(abs(H1 - HC) < 1./3600.)) & + SRATE = this%cwd (ROOTFRAC,SMCNT,SMREF,this%efcor)/(H2 - H1)/3600. + else + SRATE = 0. + endif + endif + + if(present (DRATE)) then + ! DRIP IRRIGATION + H1 = this%drip_stime + H2 = this%drip_stime + this%drip_dur + IT = this%drip_thres + if ((HC >= H1).AND.(HC < H2)) then + ! use SMCNT at H1 during H1 <= HC < H2 to compute irrigrate. + ! Notice drip uses the same soil moisture threshold of sprinkler but with 10.% efficiency correction. + if((ma <= IT).AND.(abs(H1 - HC) < 1./3600.)) & + DRATE = this%cwd(ROOTFRAC,SMCNT,SMREF,10.)/(H2 - H1)/3600. + else + DRATE = 0. + endif + endif + + if(present (FRATE)) then + ! FURROW IRRIGATION + H1 = this%flood_stime + H2 = this%flood_stime + this%flood_dur + IT = this%flood_thres + if ((HC >= H1).AND.(HC < H2)) then + ! use SMCNT at H1 during H1 <= HC < H2 to compute irrigrate. + ! Notice Furrow irrigation uses the same soil moisture threshold of sprinkler but with + ! the efficiency correction increased by 15 (e.g., Field application efficiency Sprinkler 75%, Surface Irrigation 60%. + ! Source FAO) + if((ma <= IT).AND.(abs(H1 - HC) < 1./3600.)) & + FRATE = this%cwd (ROOTFRAC,SMCNT,SMREF,this%efcor+15.)/(H2 - H1)/3600. + else + FRATE = 0. + endif + endif + + END SUBROUTINE irrig_by_method + + ! ---------------------------------------------------------------------------- + + REAL FUNCTION crop_water_deficit (this, rootfrac, asmc, smcref, efcor) + + implicit none + class(irrigation_model),intent(inout):: this + real, intent (in) :: rootfrac, asmc, smcref, efcor + + crop_water_deficit = 0. + if(smcref > asmc) crop_water_deficit = rootfrac*(smcref - asmc)*100.0/(100.0-efcor) + + END FUNCTION crop_water_deficit + + ! ---------------------------------------------------------------------------- + + logical FUNCTION IS_WITHIN_SEASON (DOY,DP, DH) + + implicit none + integer, intent (in) :: DOY,DP, DH + + IS_WITHIN_SEASON = .false. + if(DH > DP) then + if((DOY >= DP).AND.(DOY <= DH)) IS_WITHIN_SEASON = .true. + elseif (DH < DP) then + if((DOY >= DP).AND.(DOY <= 366)) IS_WITHIN_SEASON = .true. + if((DOY >= 1).AND.(DOY <= DH)) IS_WITHIN_SEASON = .true. + endif + + end FUNCTION IS_WITHIN_SEASON + + ! ---------------------------------------------------------------------------- + + real FUNCTION CROP_SEASON_STAGE (MSL, DOY,DP, DH) + + ! MSL : mid season length [-] as a fraction of the length of the season + ! DOY : doy of year + ! DP : plant date + ! DH : harvest date + ! MSL + ! 1.0 ___________________________ + ! /| |\ + ! / | | \ + ! / | | \ + ! / | | \ + ! --------------------- DOY ----------------------------> + ! t0 t1 t2 t3 + ! DP DH + ! |<---- SEAL -->| + + implicit none + real, intent (in) :: MSL + integer, intent (in) :: DOY,DP, DH + real :: seal, t0, t1, t2, t3, CTime + + CROP_SEASON_STAGE = 0. + + if(DH > DP) then + seal = real (DH - DP + 1) + CTime = real (DOY) + else + ! the crop season is fall-to-spring + seal = real(366 - DP + 1 + DH) + if((DOY >= DP).AND.(DOY <= 366)) CTIME = real (DOY) + if((DOY >= 1).AND.(DOY <= DH)) CTIME = real (DOY) + 365. + endif + + t0 = real (DP) + t1 = t0 + seal * (1. - MSL)/2. ! assumes equal development and late periods + t2 = t1 + seal*MSL + t3 = t0 + seal + + if (ctime < t1) CROP_SEASON_STAGE = (CTIME -t0)/(t1 - t0) + if ((t1 <= ctime).and.(ctime < t2)) CROP_SEASON_STAGE = 1. + if (ctime >= t2) CROP_SEASON_STAGE = (t3 - ctime)/(t3 - t2) + + if(CROP_SEASON_STAGE > 1.) CROP_SEASON_STAGE = 1. + + end FUNCTION CROP_SEASON_STAGE + +END MODULE IRRIGATION_MODULE diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 index 4c50fedcd..5728fe015 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_wrap_state.F90 @@ -31,7 +31,7 @@ module catch_wrap_stateMod real :: SURFLAY real :: FWETC, FWETL logical :: USE_FWET_FOR_RUNOFF - integer :: RUN_IRRIG, IRRIG_METHOD + integer :: RUN_IRRIG, IRRG_METHOD end type T_CATCH_STATE type CATCH_WRAP @@ -142,7 +142,7 @@ subroutine surface_params_to_wrap_state(statePtr, scf, rc) call MAPL_GetResource( SCF, statePtr%N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) call MAPL_GetResource( SCF, statePtr%AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0, __RC__ ) call MAPL_GetResource( SCF, statePtr%RUN_IRRIG, label='RUN_IRRIG:', DEFAULT=0, __RC__ ) - call MAPL_GetResource( SCF, statePtr%IRRIG_METHOD, label='IRRIG_METHOD:', DEFAULT=0, __RC__ ) + call MAPL_GetResource( SCF, statePtr%IRRG_METHOD, label='IRRG_METHOD:', DEFAULT=0, __RC__ ) select type (statePtr) type is (T_CATCHCN_STATE) ! CATCHCN 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 aab31393d..8bbd75581 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 @@ -56,7 +56,7 @@ MODULE lsm_routines PUBLIC :: SIBALB, catch_calc_soil_moist, catch_calc_zbar, catch_calc_peatclsm_waterlevel PUBLIC :: catch_calc_subtile2tile PUBLIC :: gndtmp, catch_calc_tp, catch_calc_wtotl, catch_calc_ght, catch_calc_FT - PUBLIC :: dampen_tc_oscillations, irrigation_rate + PUBLIC :: dampen_tc_oscillations INTERFACE catch_calc_zbar MODULE PROCEDURE catch_calc_zbar_scalar @@ -1744,8 +1744,9 @@ subroutine catch_calc_soil_moist( & ! 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. + ! to inconsistency. REQUIRES presence of optional "sfmc", "rzmc", and "prmc". ! ! Optional outputs "smfcun", "rzmcun", "prmcun" are surface, ! root zone, and profile moisture content for unsaturated areas only, @@ -2729,239 +2730,8 @@ subroutine dampen_tc_oscillations( dtstep, tair, tc_old, tc_new_in, & end subroutine dampen_tc_oscillations - ! ******************************************************************** - - SUBROUTINE irrigation_rate (IRRIG_METHOD, & - NTILES, AGCM_HH, AGCM_MI, AGCM_S, lons, IRRIGFRAC, PADDYFRAC, & - CLMPT,CLMST, CLMPF, CLMSF, LAIMAX, LAIMIN, LAI, & - POROS, WPWET, VGWMAX, RZMC, IRRIGRATE) - - ! !DESCRIPTION: - ! - ! NOTE: This is an experimental feature under development. - ! - ! Calculate water requirement and apply the amount to precipitation. - ! - ! Irrigate when available root zone soil moisture falls below tunable - ! irrigation threshold parameter. - ! Below GRIPC irrigated data provide fractions of croplands and paddy croplands. - ! The irrigation model is applied on a tile if: - ! (1) the irrigated fraction of the tile is greater than 0. AND - ! (2) primary or secondary type in the tile is CLM4 type 16 (cropland) AND - ! (3) LAI exceeds the LAI theshhold (60% of LAI range) - ! - ! GRIPC croplands and paddy croplands fractions determine whether to apply - ! either sprinkler or flood OR both irrigation methods. Each method has - ! its own local start times, durations and irrigation threshold parameters. - ! - ! We assume plants need available soil moisture stay above 1/3 of soil moisture range - ! [ wilting - saturation] - ! Irrigation amount is scaled to grid total crop fraction when intensity - ! is less than the fraction. Irrigation is expanded to non-crop, non-forest, - ! non-baresoil/urban tiles if intensity exceeds grid total crop fraction. - ! In latter case, scaled irrigation is applied to grassland first, - ! then further applied over the rest of tiles equally if the intensity - ! exceeds grassland fraction as well. - ! - ! Optionally efficiency correction is applied to account for field loss. - ! - ! REVISION HISTORY: - ! - ! Aug 2018: Sarith Mahanama ; Version 1 adapted from LIS subroutine clsmf25_getirrigationstates.F90 - - implicit none - - ! INPUTS - ! ------ - integer, intent (in) :: IRRIG_METHOD, NTILES, AGCM_HH, AGCM_MI, AGCM_S - real , intent (in), dimension (ntiles) :: lons, IRRIGFRAC, PADDYFRAC, LAIMAX, & - LAIMIN, LAI, CLMPT,CLMST, CLMPF, CLMSF, POROS, WPWET, VGWMAX, RZMC - ! IRRIG_METHOD : 0 sprinkler and flood combined; 1 sprinkler irrigation ; 2 flood irrigation - ! AGCM_HH / AGCM_MI / AGCM_S/ lons : Current hour, minute, second (UTC) and longitude - - ! Irrigation hotspots : Using the Global Rain-Fed, Irrigated, and Paddy Croplands (GRIPC) Dataset (Salmon et al., 2015) - ! Salmon JM, Friedl MA, Frolking S, Wisser D and Douglas EM: Global rain-fed, irrigated, - ! and paddy croplands: A new high resolution map derived from remote sensing, crop - ! inventories and climate data, Int. J. Appl. Earth Obs. Geoinf, 38, 321–334, - ! doi:10.1016/j.jag.2015.01.014, 2015. - - ! IRRIGFRAC : Fraction of irrigated croplands [-] = total number of 500m irrigated croplands pixels in the tile / - ! total number of 500m pixels in the tile - ! PADDYFRAC : Fraction of paddy croplands [-] = total number of 500m paddy croplands pixels in the tile / - ! total number of 500m pixels in the tile - ! LAIMAX / LAIMIN / LAI : Maximum, minimum and current Leaf Area Index - ! CLMPT / CLMST : CLM4 primary and secondary types (Note type 16 is cropland) - ! CLMPF / CLMSF : CLM4 fractions of primary and secondary types - ! POROS / WPWET / VGWMAX / RZMC : porosity [m3/m3], wilting point wetness [-], maximum and current root zone soil moisture content [m3/m3] - - ! ONLY output - ! ----------- - real , intent (out), dimension (ntiles) :: IRRIGRATE - - real, parameter :: efcor = 76.0 ! Efficiency Correction (%) - - ! Sprinkler parameters - ! -------------------- - real, parameter :: otimess = 6.0 ! local trigger check start time [hour] - real, parameter :: irrhrs = 4.0 ! duration of irrigation hours - real, parameter :: sprinkler_thersh = 0.5 ! soil moisture threshhold to trigger sprinkler irrigation - - ! Drip parameters (not currently implemented) - ! ------------------------------------------- - real, parameter :: otimeds = 6.0 ! local trigger check start time [hour] - real, parameter :: irrhrd = 12.0 ! duration of irrigation hours - - ! Flood parameters - ! ---------------- - real, parameter :: otimefs = 6.0 ! local trigger check start time [hour] - real, parameter :: irrhrf = 1.0 ! duration of irrigation hours - real, parameter :: flood_thersh = 0.25 ! soil moisture threshhold to trigger flood irrigation - - ! local vars - ! ---------- - real :: smcwlt, smcref, smcmax, asmc, laithresh, laifac, RZDEP, vfrac, ma, & - otimee, irrig_thresh, IrrigScale, s_irate, f_irate, local_long, local_hour - integer :: n, t, vtyp - - IRRIGRATE (:) = 0. - - TILE_LOOP : DO N = 1, NTILES - - local_long = 180. * lons(n) / PIE ! local logitude [degrees] - local_hour = AGCM_HH + AGCM_MI / 60. + AGCM_S / 3600. + 12.* local_long /180. ! local time [hours] - if (local_hour >= 24.) local_hour = local_hour - 24. - if (local_hour < 0.) local_hour = local_hour + 24. - - laithresh = laimin (n) + 0.60 * (laimax (n) - laimin (n)) - if(laimax (n) /= laimin (n)) then - laifac = (lai(n) - laimin (n)) / (laimax(n) - laimin(n)) - else - laifac = 0. - endif - - RZDEP = laifac * VGWMAX (n) / poros (n) ! root zone depth [mm] - smcwlt = RZDEP * wpwet (n) * poros (n) ! RZ soil moisture content at wilting point [mm] - smcref = RZDEP * (wpwet (n) + 0.333 * (1. - wpwet (n))) * poros(n) ! RZ reference soil moisture content [mm] - smcmax = RZDEP * poros (n) ! RZ soil moisture at saturatopm [mm] - asmc = RZDEP * rzmc (n) ! actual RZ soil moisture content [mm] - - CHECK_IRRIG_INTENSITY : IF ((IRRIGFRAC(N) + PADDYFRAC(N)) > 0.) THEN - - s_irate = 0. - f_irate = 0. - - TWO_CLMTYPS : DO t = 1, 2 - - if (t == 1) then - ! Primary CLM fraction - vtyp = NINT (CLMPT (n)) - vfrac = CLMPF (n) - endif - - if (t == 2) then - ! Secondary CLM fraction - vtyp = NINT (CLMST (n)) - vfrac = CLMSF (n) - endif - - CHECK_CROP_LAITHRESH : IF ((vtyp == 16).and.(vfrac > 0.).and.(lai(n) >= laithresh).and.(laifac > 0.)) THEN - - !----------------------------------------------------------------------------- - ! Compute irrigation scale parameter : - ! Scale the irrigation intensity to the crop % when intensity < crop%. - ! Expand irrigation for non-crop, non-forest when intensity > crop % - ! in preference order of grassland first then rest. - !----------------------------------------------------------------------------- - - IF ((IRRIGFRAC(N) + PADDYFRAC(N)) < vfrac) THEN - IrrigScale = vfrac / (IRRIGFRAC(N) + PADDYFRAC(N)) - ELSE - IrrigScale = 1. - ENDIF - - !----------------------------------------------------------------------------- - ! Get the root zone moisture availability to the plant - !----------------------------------------------------------------------------- - - if(smcref.ge.smcwlt) then - ma = (asmc - smcwlt) /(smcref - smcwlt) - else - ma = -1 - endif - - SELECT CASE (IRRIG_METHOD) - - !-------------------------------------------------------------------------------------------------------------------------- - ! IRRIGRATE : irrigation rate required to fill up water deficit before END OF IRRIGATION PERIOD (otimee - local_hour) - !-------------------------------------------------------------------------------------------------------------------------- - - CASE (0) - ! SPRINKLER AND FLOOD IRRIGATION COMBINED - ! --------------------------------------- - C_SPRINKLER : IF((IRRIGFRAC (N) > 0.).and.(local_hour >= otimess).and. (local_hour < otimess + irrhrs)) THEN - otimee = otimess + irrhrs ; irrig_thresh = sprinkler_thersh - IF ((ma <= irrig_thresh).and.(ma.ge.0)) THEN - s_irate = crop_water_deficit (IRRIGFRAC (N) * irrigScale, asmc, smcref, efcor) / (otimee - local_hour) /3600.0 - ENDIF - ENDIF C_SPRINKLER - - C_FLOOD : IF((PADDYFRAC (N) > 0.).and.(local_hour >= otimefs).and. (local_hour <= otimefs + irrhrf)) THEN - otimee = otimefs + irrhrf ; irrig_thresh = flood_thersh - IF ((ma <= irrig_thresh).and.(ma.ge.0)) THEN - f_irate = crop_water_deficit (PADDYFRAC (N) * irrigScale, asmc, smcref, efcor) / (otimee - local_hour) /3600.0 - ENDIF - ENDIF C_FLOOD - - IRRIGRATE (N) = (s_irate * IRRIGFRAC (N) + f_irate * PADDYFRAC (N)) / (IRRIGFRAC(N) + PADDYFRAC(N)) ! weighted averaged sprinkler + flood - - CASE (1) - ! SPRINKLER IRRIGATION ONLY - ! ------------------------- - SPRINKLER : IF(((IRRIGFRAC (N) + PADDYFRAC (N)) > 0.).and.(local_hour >= otimess).and. (local_hour < otimess + irrhrs)) THEN - otimee = otimess + irrhrs ; irrig_thresh = sprinkler_thersh - IF ((ma <= irrig_thresh).and.(ma.ge.0)) THEN - IRRIGRATE (N) = crop_water_deficit ((IRRIGFRAC (N) + PADDYFRAC (N)) * irrigScale, asmc, smcref, efcor) / & - (otimee - local_hour) /3600.0 - ENDIF - ENDIF SPRINKLER - - CASE (2) - ! FLOOD IRRIGATION ONLY - ! --------------------- - FLOOD : IF(((IRRIGFRAC (N) + PADDYFRAC (N)) > 0.).and.(local_hour >= otimefs).and. (local_hour <= otimefs + irrhrf)) THEN - otimee = otimefs + irrhrf ; irrig_thresh = flood_thersh - IF ((ma <= irrig_thresh).and.(ma.ge.0)) THEN - IRRIGRATE (N) = crop_water_deficit ((IRRIGFRAC (N) +PADDYFRAC (N)) * irrigScale, asmc, smcref, efcor) / & - (otimee - local_hour) /3600.0 - ENDIF - ENDIF FLOOD - - CASE DEFAULT - print *, 'IN IRRIGATION_RATE : IRRIGATION_METHOD can be 0, 1, or 2' - call exit(1) - END SELECT - END IF CHECK_CROP_LAITHRESH - END DO TWO_CLMTYPS - END IF CHECK_IRRIG_INTENSITY - END DO TILE_LOOP - - END SUBROUTINE irrigation_rate - - ! ******************************************************************** - - REAL FUNCTION crop_water_deficit (IrrigScale, asmc, smcref, efcor) - - implicit none - - real, intent (in) :: IrrigScale, asmc, smcref, efcor - real :: twater - - twater = smcref - asmc - twater = twater * IrrigScale ! Scale the irrigation intensity - crop_water_deficit = twater*(100.0/(100.0-efcor)) ! Apply efficiency correction - - END FUNCTION crop_water_deficit + ! ******************************************************************* + - ! ******************************************************************** END MODULE lsm_routines diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index 1ba28b70b..817ce7392 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -164,20 +164,54 @@ # ---- Run irrigation module # +# NOTE: The irrigation model needs the irrigation parameter file ('irrigation_IMxJM_DL.dat') in BCSDIR to run the model. # 0 : Do NOT run irrigation module (default) -# 1 : YES, run irrigation module +# 1 : YES, run irrigation module # # GEOSagcm=>RUN_IRRIG: 0 # GEOSldas=>RUN_IRRIG: 0 - -# ---- Irrigation model method # -# 0 : Sprinkler and Flood irrigation combined (default) -# 1 : Sprinkler irrigation only -# 2 : Flood irrigation only +# ---- Irrigation trigger +# +# 0 : (Default) LAI-based trigger turns irrigation on if LAI >= (LAImin + IRRG_LAI_THRES * (LAImax - LAImin)) +# 1 : Use planting and harvesting times from 26 crop calendars +# +# GEOSagcm=>IRRG_TRIGGER: 0 +# GEOSldas=>IRRG_TRIGGER: 0 +# +# ---- Irrigation method +# +# The LAI-based trigger (IRRG_TRIGGER: 0) offers 4 choices of irrigation method: +# +# 0 : CONCURRENT sprinkler, drip, and furrow irrigation for specified tile fractions (default) +# 1 : Sprinkler irrigation on entire tile +# 2 : Drip irrigation on entire tile +# 3 : Furrow irrigation on entire tile +# +# GEOSagcm=>IRRG_METHOD: 0 +# GEOSldas=>IRRG_METHOD: 0 +# +# ----- Miscellaneous irrigation parameters +# +# GEOSldas=>IRRG_FRAC_THRES: 0.01 # threshold of tile fraction to turn the irrigation model on +# GEOSldas=>IRRG_LAI_THRES: 0.6 # threshold of LAI range to turn irrigation on (for IRRG_TRIGGER: 0) +# +# GEOSldas=>IRRG_SPR_STIME: 6.0 # sprinkler irrigatrion start time [hours] +# GEOSldas=>IRRG_SPR_DUR: 4.0 # sprinkler irrigation duration [hours] +# GEOSldas=>IRRG_SPR_THRES: 0.7 # threshold for soil moisture availability ("MA") to trigger sprinkler irrigation [dim-less] +# +# GEOSldas=>IRRG_DRP_STIME: 8.0 # drip irrigatrion start time [hours] +# GEOSldas=>IRRG_DRP_DUR: 8.0 # drip irrigation duration [hours] +# GEOSldas=>IRRG_DRP_THRES: 0.7 # threshold for soil moisture availability ("MA") to trigger drip irrigation [dim-less] +# +# GEOSldas=>IRRG_FLD_STIME: 6.0 # flood irrigatrion start time [hours] +# GEOSldas=>IRRG_FLD_DUR: 8.0 # flood irrigation duration [hours] +# GEOSldas=>IRRG_FLD_THRES: 0.7 # threshold for soil moisture availability ("MA") to trigger flood irrigation [dim-less] # -# GEOSagcm=>IRRIG_METHOD: 0 -# GEOSldas=>IRRIG_METHOD: 0 +# GEOSldas=>IRRG_EFCOR: 25.0 # efficiency correction (% water loss: efcor = 0% denotes 100% efficient water use) +# GEOSldas=>IRRG_MIDS_LNGTH: 0.6 # mid-season length as fraction of crop growing season length (for IRRG_TRIGGER: 1); +# # length of development season = length of end season = (1 - MIDS_LENGTH)/2; +# # see function CROP_SEASON_STAGE() #--------------------------------------------------------# diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt index c52656386..ed94cbb4e 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt @@ -7,6 +7,7 @@ rasterize.F90 read_riveroutlet.F90 CubedSphere_GridMod.F90 rmTinyCatchParaMod.F90 +module_irrig_params.F90 zip.c util.c ) @@ -17,7 +18,7 @@ endif () set_source_files_properties(mkMITAquaRaster.F90 PROPERTIES COMPILE_FLAGS "${BYTERECLEN}") -esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared ESMF::ESMF NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran) +esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared GEOSirrigation_GridComp ESMF::ESMF NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran) if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) target_compile_definitions(${this} PRIVATE USE_EXTERNAL_FINDLOC) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro index d5584b3dd..e827a676d 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro @@ -3102,9 +3102,9 @@ if file_test ('limits.idl') then restore,'limits.idl' id = NCDF_OPEN (filename, /NOWRITE) -NCDF_VARGET, id,'SPRINKLERFR',SPRINKLERV -NCDF_VARGET, id,'DRIPFR',DRIPV -NCDF_VARGET, id,'FLOODFR',FLOODV +NCDF_VARGET, id,'IRRG_IRRIGFRAC_SPR',SPRINKLERV +NCDF_VARGET, id,'IRRG_IRRIGFRAC_DRP',DRIPV +NCDF_VARGET, id,'IRRG_IRRIGFRAC_FRW',FLOODV NCDF_CLOSE, id @@ -3197,8 +3197,8 @@ if file_test ('limits.idl') then restore,'limits.idl' id = NCDF_OPEN (filename, /NOWRITE) -NCDF_VARGET, id,'LAIMIN',LAI_MNV -NCDF_VARGET, id,'LAIMAX',LAI_MXV +NCDF_VARGET, id,'IRRG_LAIMIN',LAI_MNV +NCDF_VARGET, id,'IRRG_LAIMAX',LAI_MXV NCDF_CLOSE, id @@ -3313,14 +3313,14 @@ if file_test ('limits.idl') then restore,'limits.idl' id = NCDF_OPEN (filename, /NOWRITE) -NCDF_VARGET, id,'IRRIGFRAC',IRRIGFRACV -NCDF_VARGET, id,'PADDYFRAC',PADDYFRACV +NCDF_VARGET, id,'IRRG_IRRIGFRAC',IRRIGFRAC_V +NCDF_VARGET, id,'IRRG_PADDYFRAC',IRRG_PADDYFRACV NCDF_VARGET, id,'RAINFEDFRAC',RAINFEDFRACV NCDF_CLOSE, id -IRRIGFRACV (where (IRRIGFRACV gt 1.)) = !VALUES.F_NAN -PADDYFRACV (where (PADDYFRACV gt 1.)) = !VALUES.F_NAN +IRRIGFRAC_V (where (IRRIGFRAC_V gt 1.)) = !VALUES.F_NAN +IRRG_PADDYFRACV (where (IRRG_PADDYFRACV gt 1.)) = !VALUES.F_NAN RAINFEDFRACV(where (RAINFEDFRACV gt 1.)) = !VALUES.F_NAN im = n_elements(tile_id[*,0]) @@ -3332,22 +3332,22 @@ dy = 180. / jm x = indgen(im)*dx -180. + dx/2. y = indgen(jm)*dy -90. + dy/2. -IRRIGFRAC = REPLICATE (!VALUES.F_NAN,IM, JM) -PADDYFRAC = REPLICATE (!VALUES.F_NAN,IM, JM) +IRRG_IRRIGFRAC = REPLICATE (!VALUES.F_NAN,IM, JM) +IRRG_PADDYFRAC = REPLICATE (!VALUES.F_NAN,IM, JM) RAINFEDFRAC = REPLICATE (!VALUES.F_NAN,IM, JM) for j = 0l, jm -1l do begin for i = 0l, im -1 do begin if(tile_id[i,j] gt 0) then begin - IRRIGFRAC (i,j) = IRRIGFRACV (tile_id[i,j] -1) - PADDYFRAC (i,j) = PADDYFRACV (tile_id[i,j] -1) + IRRG_IRRIGFRAC (i,j) = IRRIGFRAC_V (tile_id[i,j] -1) + IRRG_PADDYFRAC (i,j) = IRRG_PADDYFRACV (tile_id[i,j] -1) RAINFEDFRAC(i,j) = RAINFEDFRACV (tile_id[i,j] -1) endif endfor endfor -IRRIGFRAC (where (IRRIGFRAC eq 0.)) = !VALUES.F_NAN -PADDYFRAC (where (PADDYFRAC eq 0.)) = !VALUES.F_NAN +IRRG_IRRIGFRAC (where (IRRG_IRRIGFRAC eq 0.)) = !VALUES.F_NAN +IRRG_PADDYFRAC (where (IRRG_PADDYFRAC eq 0.)) = !VALUES.F_NAN RAINFEDFRAC(where (RAINFEDFRAC eq 0.)) = !VALUES.F_NAN colors = indgen (21) + 140 @@ -3404,10 +3404,10 @@ if file_test ('limits.idl') then restore,'limits.idl' id = NCDF_OPEN (filename, /NOWRITE) -NCDF_VARGET, id,'IRRIGPLANT',plantv -NCDF_VARGET, id,'IRRIGHARVEST',harvestv -NCDF_VARGET, id,'CROPIRRIGFRAC',Fracv -NCDF_VARGET, id,'IRRIGTYPE',irrigtypev +NCDF_VARGET, id,'IRRG_DOY_PLANT',plantv +NCDF_VARGET, id,'IRRG_DOY_HARVEST',harvestv +NCDF_VARGET, id,'IRRG_CROPIRRIGFRAC',Fracv +NCDF_VARGET, id,'IRRG_TYPE',IRRG_TYPE_V NCDF_VARGET, id,'CROPCLASSNAME',cropname NCDF_CLOSE, id @@ -3426,7 +3426,7 @@ y = fltarr (IM,JM) PLANT = REPLICATE (!VALUES.F_NAN,IM, JM, 2, 26) HARVEST = REPLICATE (!VALUES.F_NAN,IM, JM, 2, 26) FRAC = REPLICATE (!VALUES.F_NAN,IM, JM, 26) -IRRIGTYPE = REPLICATE (!VALUES.F_NAN,IM, JM, 26) +IRRG_TYPE = REPLICATE (!VALUES.F_NAN,IM, JM, 26) for j = 0l, jm -1l do begin for i = 0l, im -1 do begin @@ -3436,7 +3436,7 @@ for j = 0l, jm -1l do begin PLANT (i,j,*,*) = PLANTV (tile_id[i,j] -1,*,*) HARVEST (i,j,*,*) = HARVESTV (tile_id[i,j] -1,*,*) FRAC (i,j,*) = FRACV (tile_id[i,j] -1,*) - IRRIGTYPE(i,j,*) = IRRIGTYPEV (tile_id[i,j] -1,*) + IRRG_TYPE(i,j,*) = IRRG_TYPE_V (tile_id[i,j] -1,*) endif endfor endfor @@ -3544,7 +3544,7 @@ for n = 0, 25 do begin endif if (col eq 3) then begin - ptitle = ' : IRRIGTYPE' + ptitle = ' : IRRG_TYPE' data_grid = data4 endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py index 386d834c4..632fbeb1d 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py @@ -199,6 +199,7 @@ def ask_questions(default_grid="Cubed-Sphere"): "v11 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2", \ "v12 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix", \ "v13 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix + mean land elevation fix", \ + "v14 : v13 + Irrigation", \ "ICA : Icarus (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus/)", \ "GM4 : Ganymed-4_0 (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Ganymed-4_0/)", \ "F25 : Fortuna-2_5 (archived*: n/a)"], diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 index fcafbf333..649e5009e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 @@ -21,12 +21,12 @@ PROGRAM mkCatchParam ! ! Sarith Mahanama - March 23, 2012 ! Email: sarith.p.mahanama@nasa.gov - use MAPL, only: MAPL_ease_extent, MAPL_ReadTilingNC4 + + use MAPL, ONLY: MAPL_ease_extent, MAPL_ReadTilingNC4 use rmTinyCatchParaMod use process_hres_data use MAPL_ExceptionHandling - - ! use module_irrig_params, ONLY : create_irrig_params + use module_irrig_params, ONLY : create_irrig_params implicit none @@ -809,10 +809,17 @@ PROGRAM mkCatchParam write (log_file,'(a)')' ' endif - ! inquire(file='clsm/irrig.dat', exist=file_exists) - ! if (.not.file_exists) call create_irrig_params (nc,nr,fnameRst) - ! write (log_file,'(a)')'Done computing irrigation model parameters ...............13' - + if(IRRIGBCS) then + tmpstring = 'Step 15: Irrigation' + inquire(file='clsm/irrig.dat', exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)') trim(tmpstring) + write (log_file,'(a)')' Creating file...' + call create_irrig_params (nc,nr,fnameRst) + write (log_file,'(a)') ' Done computing irrigation model parameters........' + endif + endif + write (log_file,'(a)')'============================================================' write (log_file,'(a)')'DONE creating CLSM data files...............................' write (log_file,'(a)')'============================================================' diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/module_irrig_params.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/module_irrig_params.F90 new file mode 100755 index 000000000..e3b361a65 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/module_irrig_params.F90 @@ -0,0 +1,1612 @@ +#define VERIFY_(A) IF(A/=0)THEN;PRINT *,'ERROR AT LINE ', __LINE__;STOP;ENDIF +#define ASSERT_(A) if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif + +module module_irrig_params +! Version 1 - Sarith Mahanama +! Version 2 - Stefano Casirati - LAI min-max obtained from LAI climatology boundary conditions + + use rmTinyCatchParaMod, ONLY : RegridRaster,regridrasterreal + use process_hres_data, ONLY : get_country_codes + use MAPL + use irrigation_module, ONLY : NCROPS => IRRG_NCROPS + + implicit none + + INCLUDE 'netcdf.inc' + + private + + public :: create_irrig_params + +contains + + subroutine create_irrig_params (nc, nr, gfile) + + implicit none + + integer , intent (in) :: nc, nr + character(*) , intent (in) :: gfile + REAL, PARAMETER :: UNDEF = -9999., UNDEFG = 1.e15 + + ! GRIPC data + ! ---------- + + integer, parameter :: NX_gripc = 86400 + integer, parameter :: NY_gripc = 43200, NY_GripcData = 36000 + character*300, parameter :: GRIPC_file = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/crop_fraction_data/v1/irrigtype_salmon2013.flt' + real, allocatable, dimension (:) :: IGRIPC, RGRIPC, PGRIPC, NGRIPC + + ! MIRCA2000 data + ! -------------- + + integer, parameter :: NX_mirca = 4320 + integer, parameter :: NY_mirca = 2160 + integer, parameter :: NMON = 12, STRLEN = 20 + real, parameter :: DXY_mirca= 360./REAL(NX_mirca) + real, parameter :: lat1_mirca = 90.0 - DXY_mirca / 2.0 !1st grid center lat + real, parameter :: lon1_mirca = -180.0 + DXY_mirca / 2.0 !1st grid center lon + character*300, parameter :: MIRCA_pathIrr = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/crop_fraction_data/v1/irrigated/crop_' + character*300, parameter :: MIRCA_pathRain = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/crop_fraction_data/v1/rainfed/crop_' + real, allocatable, dimension(:,:,:) :: MIFRAC, MRFRAC + + ! Global Irrigated Area data (GIA) + ! -------------------------------- + + integer, parameter :: NX_GIA = 43200 + integer, parameter :: NY_GIA = 21600, NY_GIAData = 18000 + character*300, parameter :: GIA_file = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/irrigation_class/v1/global_irrigated_areas.nc4' + real, allocatable, dimension (:) :: GIAFRAC + + ! LAI data + ! -------- + + integer, parameter :: NX_LAI = 86400 + integer, parameter :: NY_LAI = 43200 + character*300, parameter :: LAI_file = '/discover/nobackup/projects/lis/LS_PARAMETERS/MODIS/MCD15A2H.006/MCD15A2H.006_LAI_YYYY' + + ! Irrigation Method + ! ----------------- + character*300, parameter :: IM_path = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/country_code_IMethod/v1/' + + ! Global/Local variables + ! ---------------------- + + integer, allocatable, dimension (:,:) :: tile_id + integer :: i,j, NTILES, STATUS, tindex1,pfaf1, NCOutID + real, allocatable, dimension (:) :: tile_lon, tile_lat + real :: minlon,maxlon,minlat,maxlat + + !The codes for the 26 crop classes are as follows: + !1 Wheat + !2 Maize + !3 Rice + !4 Barley + !5 Rye + !6 Millet + !7 Sorghum + !8 Soybeans + !9 Sunflower + !10 Potatoes + !11 Cassava + !12 Sugar cane + !13 Sugar beet + !14 Oil palm + !15 Rape seed / Canola + !16 Groundnuts / Peanuts + !17 Pulses + !18 Citrus + !19 Date palm + !20 Grapes / Vine + !21 Cotton + !22 Cocoa + !23 Coffee + !24 Others perennial + !25 Fodder grasses + !26 Others annual + + character(len=STRLEN),dimension(ncrops) :: cname = (/"Wheat ", & + "Maize ", "Rice ", "Barley ", & + "Rye ", "Millet ", "Sorghum ", & + "Soybeans ", "Sunflower ", "Potatoes ", & + "Cassava ", "Sugar cane ", "Sugar beet ", & + "Oil palm ", "Rape seed /Canola ", "Groundnuts/ Peanuts ", & + "Pulses ", "Citrus ", "Date palm ", & + "Grapes / Vine ", "Cotton ", "Cocoa ", & + "Coffee ", "Other sperennial ", "Fodder grasses ", & + "Others annual "/) + + ! (1) Reading rst file and NTILES + ! ------------------------------- + + open (10,file= trim(gfile)//'.rst',status='old',action='read', & + form='unformatted',convert='little_endian') + allocate (tile_id (1:nc,1:nr)) + + do j=1,nr + read(10)tile_id(:,j) + end do + close (10,status='keep') + + open (10,file='clsm/catchment.def',status='old',action='read', form='formatted') + read (10, *) ntiles + allocate (tile_lon (1:NTILES)) + allocate (tile_lat (1:NTILES)) + do i = 1, NTILES + read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat + tile_lon(i) = (minlon + maxlon)/2. + tile_lat(i) = (minlat + maxlat)/2. + end do + close (10, status = 'keep') + + call OpenFile + + ! (3) Process GIA + ! ----------------- + + allocate (GIAFRAC (NTILES)) + call ReadProcess_GIA (NC, NR, NTILES, tile_id, GIAFRAC) + + ! (4) Process GRIPC and LAI (Min-Max) + ! ---------------------------------- + + allocate (IGRIPC (NTILES)) + allocate (RGRIPC (NTILES)) + allocate (PGRIPC (NTILES)) + allocate (NGRIPC (NTILES)) + call ReadProcess_GRIPC (NC, NR, NTILES, tile_id, IGRIPC, RGRIPC, PGRIPC, NGRIPC) + + ! (1) Process MIRCA2000 + ! --------------------- + + allocate (MIFRAC (NTILES, 12, NCROPS)) + allocate (MRFRAC (NTILES, 12, NCROPS)) + call ReadProcess_MIRCA (NC, NR, NTILES, tile_id, MIFRAC, MRFRAC) + + ! (5) Create parameter file for the model + ! --------------------------------------- + + call MergeData (NTILES) + + return + + contains + + ! =========================================================================================== + + SUBROUTINE OpenFile + + implicit none + integer :: i, m, n,l, lid, mid, cid, vid, sid + integer, dimension(8) :: date_time_values + character (22) :: time_stamp + character (len=STRLEN) :: ThisCrop + real :: abm_int, peatf_r, gdp_r, hdm_r + real, dimension (:), allocatable :: field_cap + + status = NF_CREATE ('clsm/irrig.dat' , NF_NETCDF4, NCOutID) ; VERIFY_(STATUS) + status = NF_DEF_DIM(NCOutID, 'tile' , NTILES, lid) ; VERIFY_(STATUS) + status = NF_DEF_DIM(NCOutID, 'unknown_dim1', NCROPS, cid) ; VERIFY_(STATUS) + status = NF_DEF_DIM(NCOutID, 'unknown_dim2', 2, mid) ; VERIFY_(STATUS) + status = NF_DEF_DIM(NCOutID, 'strlen' , STRLEN, sid) ; VERIFY_(STATUS) + + ! GIA -> GRIPC MERGE + ! ------------------ + + status = NF_DEF_VAR(NCOutID, 'RAINFEDFRAC' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of rainfed cropland'), & + 'fraction of rainfed cropland') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_IRRIGFRAC' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of irrigated cropland'), & + 'fraction of irrigated cropland') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_PADDYFRAC' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of paddy cropland'), & + 'fraction of paddy cropland') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'FIELDCAP' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('soil field capacity'), & + 'soil field capacity') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 5,'m3/m3') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + ! GIA- GRIPC -> MIRCA + ! ------------------- + + status = NF_DEF_VAR(NCOutID, 'CROPCLASSNAME' , NF_CHAR, 2 ,(/sid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('Crop Class Name'), & + 'Crop Class Name') ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_CROPIRRIGFRAC' , NF_FLOAT, 2 ,(/lid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('Crop irrigated fraction'), & + 'Crop irrigated fraction') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'CROPRAINFEDFRAC' , NF_FLOAT, 2 ,(/lid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('Crop rainfed fraction'), & + 'Crop rainfed fraction') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + ! Crop calendar + ! ------------- + + status = NF_DEF_VAR(NCOutID, 'IRRG_DOY_PLANT' , NF_FLOAT, 3 ,(/lid, mid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('DOY start planting'), & + 'DOY start planting') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 4,'days') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_DOY_HARVEST' , NF_FLOAT, 3 ,(/lid, mid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('DOY end harvesting'), & + 'DOY end harvesting') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 4,'days') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'RAINFEDPLANT' , NF_FLOAT, 3 ,(/lid, mid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('DOY start planting'), & + 'DOY start planting') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 4,'days') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'RAINFEDHARVEST' , NF_FLOAT, 3,(/lid, mid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('DOY end harvesting'), & + 'DOY end harvesting') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 4,'days') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + + ! IRRIG TYPE + ! ---------- + + status = NF_DEF_VAR(NCOutID, 'IRRG_TYPE' , NF_FLOAT, 2 ,(/lid, cid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', & + LEN_TRIM('Preferred Irrig Type : Concurrent (0) SPRINKLER(1) DRIP(2) FLOOD(3) AVOID (negative)'), & + 'Preferred Irrig Type : Concurrent (0) SPRINKLER(1) DRIP(2) FLOOD(3) AVOID (negative)') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_IRRIGFRAC_SPR' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of sprinkler irrigation'), & + 'fraction of sprinkler irrigation') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_IRRIGFRAC_DRP' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of drip irrigation'), & + 'fraction of drip irrigation') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_IRRIGFRAC_FRW' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('fraction of flood irrigation'), & + 'fraction of flood irrigation') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + ! LAI + ! --- + status = NF_DEF_VAR(NCOutID, 'IRRG_LAIMIN' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('Minimum LAI irrigated crops'), & + 'Minimum LAI irrigated crops') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + status = NF_DEF_VAR(NCOutID, 'IRRG_LAIMAX' , NF_FLOAT, 1 ,(/lid/), vid) ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'long_name', LEN_TRIM('Maximum LAI irrigated crops'), & + 'Maximum LAI irrigated crops') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, vid, 'units', 1,'-') ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'missing_value', NF_REAL,1, UNDEFG) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'add_offset', NF_REAL,1, 0.) ; VERIFY_(STATUS) + status = NF_PUT_ATT_REAL(NCOutID, vid, 'scale_factor', NF_REAL,1, 1.) ; VERIFY_(STATUS) + + call date_and_time(VALUES=date_time_values) + + write (time_stamp,'(i4.4,a1,i2.2,a1,i2.2,1x,a2,1x,i2.2,a1,i2.2,a1,i2.2)') & + date_time_values(1),'-',date_time_values(2),'-',date_time_values(3),'at', & + date_time_values(5),':',date_time_values(6),':',date_time_values(7) + + status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'CreatedBy', LEN_TRIM('Sarith Mahanama, Stefano Casirati'), & + 'Sarith Mahanama, Stefano Casirati') ; VERIFY_(STATUS) + status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'Contact' , LEN_TRIM('sarith.p.mahanama@nasa.gov'), & + 'sarith.p.mahanama@nasa.gov') + status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'Date' , LEN_TRIM(time_stamp),trim(time_stamp)) + + status = NF_ENDDEF(NCOutID) + + DO n = 1, NCROPS + i = LEN_TRIM(cname(n)) + ThisCrop = trim(cname(n)) + status = NF_PUT_VARA_text(NCOutID,VarID(NCOutID,'CROPCLASSNAME') ,(/1,n/),(/i,1/),ThisCrop (1:i)) ; VERIFY_(STATUS) + END DO + + ! Put field capacity + + open (10,file='clsm/CLM4.5_abm_peatf_gdp_hdm_fc', & + form='formatted',status='unknown',action = 'read') + allocate (field_cap(1:NTILES)) + + do n = 1, NTILES + read (10,'(2I10, i3, f8.4, f8.2, f10.2, f8.4)' ) i, vid, abm_int, peatf_r, gdp_r, hdm_r, field_cap(n) + end do + + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'FIELDCAP' ) ,(/1/),(/NTILES/),field_cap ) ; VERIFY_(STATUS) + close (10, status = 'keep') + deallocate (field_cap) + + + END SUBROUTINE OpenFile + + ! ----------------------------------------------------------------------------------------- + + SUBROUTINE MergeData (NTILES) + + implicit none + integer, intent (in) :: NTILES + real, allocatable, dimension (:,:) :: MI, MR + real :: MICROP, MRCROP, MIRICEA, MRRICEA, MICROPA, MRCROPA, DF, SF, FF, ITYPE(3) + integer :: i, j, m, n,l, t + real, dimension (:), ALLOCATABLE :: sprinkler, drip, flood + integer :: nc, day1, dayL, day1_2, dayL_2 + integer, dimension (12) :: fmonth, fmonth2, fmonth3 + integer, dimension (12) :: DOY_MidMonth, DOY_BegMonth, DOY_EndMonth + logical, dimension (4) :: found = .false. + integer, allocatable , dimension (:) :: crop_mons + real, allocatable, dimension (:,:,:) :: IRRG_DOY_PLANT, IRRG_DOY_HARVEST, RAINFEDPLANT, RAINFEDHARVEST + real, allocatable, dimension (:,:) :: IRRG_TYPE + + data DOY_BegMonth / 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335/ + data DOY_MidMonth /15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349/ + data DOY_EndMonth /31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 366/ + + ALLOCATE (FLOOD (1:NTILES)) + ALLOCATE (SPRINKLER(1:NTILES)) + ALLOCATE (DRIP (1:NTILES)) + CALL ReadProcess_IMethod (NTILES, sprinkler, drip, flood) + + ! MERGING PROCEDURE + ! ================= + ! CROP CALENDAR + ! GIA GIA-GRIPC GIA-GRIPC-MIRCA ------------- + ! --- --------- --------------- + ! CAL_STEP 1 : USE MIRCA monthly crop fractions : + ! -> NO -> 100% iCrop -> YES scale to match GIA-GRIPC 1) Year around : 366 + ! | | 1) MIRCA type 3 is paddy 2) if not year around check whether 1 or 2 seasons + ! | | 2) Scale fractions to match + ! | | the sum of individual crops CAL_STEP 2 : Look for gaps in calendar for GIA-GRIPC-MIRCA and + ! GIA -> GRIPC_| -> MIRCA -| to GIA-GRIPC use nearest neighbor with similar crop type + ! HAVE | | iCrop | + ! | | | -> YES Scale to match GIA-GRIPC + ! | |-> iCrop | | | 1) MIRCA type 3 is paddy + ! | | | | | 2) Scale fractions to match + ! | | | | | the sum of individual crops + ! |-> YES|-> Paddy | |-> NO ->MIRCA -| to GIA-GRIPC + ! | | rCrop | 3) set MIRCA rCrop to zero. + ! | | | + ! |-> rCrop | |-> NO Plant Wheat GIA-GRIPC frac + ! | + ! | + ! | -> YES Max (GIA-GRIPC, MIRCA rCrop) + ! | | 1) MIRCA type 3 is paddy + ! | | 2) Scale fractions to match + ! | | the sum of individual crops + ! | ->MIRCA -| to GIA-GRIPC-MIRCA + ! rCrop | + ! |-> NO Plant Wheat + + allocate (MI (1 : NTILES, 1 : NCROPS)) + allocate (MR (1 : NTILES, 1 : NCROPS)) + allocate (IRRG_DOY_PLANT (1 : NTILES, 1 : 2, 1 : NCROPS)) + allocate (IRRG_DOY_HARVEST (1 : NTILES, 1 : 2, 1 : NCROPS)) + allocate (RAINFEDPLANT (1 : NTILES, 1 : 2, 1 : NCROPS)) + allocate (RAINFEDHARVEST (1 : NTILES, 1 : 2, 1 : NCROPS)) + allocate (IRRG_TYPE (1 : NTILES, 1 : NCROPS)) + + MI = 0. + MR = 0. + IRRG_TYPE = 0 + + IRRG_DOY_PLANT = 998 + IRRG_DOY_HARVEST = 998 + RAINFEDPLANT = 998 + RAINFEDHARVEST = 998 + + ! Compute annual maximum fractions from MIRCA monthly fractions : + ! (1) crop specific and (2) paddy and irrigated crops, seperately for rainfed and irrigated + + DO m = 1,12 + DO I = 1, NTILES + IF((maxval(MIFRAC (I,m,:)) > 0.).OR.(maxval(MRFRAC (I,m,:)) > 0.)) then + DO N = 1, NCROPS + ! Maximum Crop Fraction over 12 months + IF(MI (I,n) < MIFRAC (I,m,n)) MI (I,N) = MIFRAC (I,m,n) + IF(MR (I,n) < MRFRAC (I,m,n)) MR (I,N) = MRFRAC (I,m,n) + END DO + ENDIF + END DO + END DO + + ! CROP CALENDAR CAL_STEP1 : Compute plant/harvest dates and create HYBRID of GRIPC and MIRCA fractions + ! ---------------------------------------------------------------------------------------------------- + + DO I = 1, NTILES + IF((maxval(MI (I,:)) > 0.).OR.(maxval(MR (I,:)) > 0.).OR.(IGRIPC (I) > 0.).OR.(PGRIPC (I) > 0.).OR.(RGRIPC (I) > 0.)) THEN + + ! Crop planting/Harvesting days + ! ----------------------------- + ! OUTPPUTS IRRG_DOY_PLANT, IRRG_DOY_HARVEST, RAINFEDPLANT, RAINFEDHARVEST + + forall (m=1:12) fmonth3(m) = m + + DO N = 1, NCROPS + + fmonth = 0. + fmonth2 = 0. + + DO t = 1,2 + + if (t == 1) nc = count (MIFRAC (i,:,n) > 0.) + if (t == 2) nc = count (MRFRAC (i,:,n) > 0.) + + if(nc > 0) then + + if(nc == 12) then + ! year around + + day1 = 1 + dayL = 366 + day1_2 =998 + dayL_2 =998 + + else + + fmonth = 0. + fmonth2 = 0. + day1 = 998 + dayL = 998 + day1_2 = 998 + dayL_2 = 998 + + if (t == 1) forall (m=1:12) fmonth(m) = ceiling (MIFRAC (i,m,n)) + if (t == 2) forall (m=1:12) fmonth(m) = ceiling (MRFRAC (i,m,n)) + + fmonth2(1) = 1 + do m = 2,12 + if(fmonth(m) == fmonth(m-1)) then + fmonth2 (m) = fmonth2(m-1) + else + fmonth2 (m) = fmonth2(m-1) + 1 + endif + end do + + if(maxval (fmonth2) > 3) then + + ! This crop grows in 2 seasons + ! ............................ + + allocate (crop_mons (1:NC)) + crop_mons = pack(fmonth3, mask = (fmonth > 0.)) + found = .false. + + if(fmonth(1) == 1) then + if(fmonth(12) == 0) then + ! Season begins on Jan 1 + day1 = DOY_BegMonth(crop_mons(1)) + found(1) = .true. + do m = 1, nc-1 + if((crop_mons(m+1) - crop_mons(m)) > 1) then + dayL = DOY_EndMonth(crop_mons(m)) + day1_2 = DOY_BegMonth(crop_mons(m+1)) + dayL_2 = DOY_EndMonth(crop_mons(nc)) + found(2) = .true. + found(3) = .true. + found(4) = .true. + exit + endif + enddo + else + ! season one begins in the fall + do m = 1, nc-1 + if((crop_mons(m+1) - crop_mons(m)) > 1) then + if(.not.found(2)) then + dayL = DOY_EndMonth(crop_mons(m)) + day1_2 = DOY_BegMonth(crop_mons(m+1)) + found(2) = .true. + found(3) = .true. + elseif (.not.found(4)) then + found(4) = .true. + found(1) = .true. + dayL_2 = DOY_EndMonth(crop_mons(m)) + day1 = DOY_BegMonth(crop_mons(m+1)) + endif + endif + end do + endif + else + + ! season 1 brings in the spring + day1 = DOY_BegMonth(crop_mons(1)) + found(1) = .true. + do m = 1, nc-1 + if((crop_mons(m+1) - crop_mons(m)) > 1) then + dayL = DOY_EndMonth(crop_mons(m)) + day1_2 = DOY_BegMonth(crop_mons(m+1)) + dayL_2 = DOY_EndMonth(crop_mons(nc)) + found(2) = .true. + found(3) = .true. + found(4) = .true. + exit + endif + enddo + endif + deallocate (crop_mons) + + else + + ! Single crop season + ! .................. + if((fmonth(1) == 0).and.(fmonth(12) == 0)) then + day1 = DOY_BegMonth (maxloc(fmonth, 1)) + dayL = DOY_EndMonth (maxloc(fmonth2, 1)-1) + else + if((fmonth(1) == 1).and.(fmonth(12) == 1)) then + day1 = DOY_BegMonth (maxloc(fmonth2, 1,mask=(fmonth2 > 2))) + dayL = DOY_EndMonth (maxloc(fmonth2, 1,mask=(fmonth2 == 2))-1) + endif + if((fmonth(1) == 0).and.(fmonth(12) == 1)) then + day1 = DOY_BegMonth (maxloc(fmonth2, 1)) + dayL = DOY_EndMonth (12) + endif + if((fmonth(1) == 1).and.(fmonth(12) == 0)) then + day1 = DOY_BegMonth (1) + dayL = DOY_EndMonth (maxloc(fmonth2, 1)-1) + endif + endif + endif + + if (t == 1) then + IRRG_DOY_PLANT (I,1,N) = day1 + IRRG_DOY_PLANT (I,2,N) = day1_2 + IRRG_DOY_HARVEST (I,1,N) = dayL + IRRG_DOY_HARVEST (I,2,N) = dayL_2 + else + RAINFEDPLANT (I,1,N) = day1 + RAINFEDPLANT (I,2,N) = day1_2 + RAINFEDHARVEST (I,1,N) = dayL + RAINFEDHARVEST (I,2,N) = dayL_2 + endif + endif + endif + END DO + END DO + + ! 1. Main Fractions (OUTPUT) : + ! 1.1 IRRG_IRRIGFRAC : The maximum value between (1) GRIPC irrigfrac, and (2) sum of MIRCA monthly crop frations without rice + ! 1.2 IRRG_PADDYFRAC : The maximum value between (1) GRIPC paddyfrac, and (2) monthly rice fractions from MIRCA + ! 1.3 RAINFEDFRAC : The maximum value between (1) GRIPC rainfedfrac, and (2) sum of MIRCA monthly crop frations + ! 1.4 MI (I,CROPS) : Irrigated crop fractions with rice is the 3rd slice crops = 3 + ! 1.5 MR (I,CROPS) : rainfed crop fractions with rice is the 3rd slice crops = 3 + + MICROP = 0. + MRCROP = 0. + MIRICEA= 0. + MRRICEA= 0. + MICROPA= 0. + MRCROPA= 0. + + DO N = 1, NCROPS + IF (n == 3) THEN + IF(MIRICEA < MI (I,n)) MIRICEA = MI (I,n) + IF(MRRICEA < MR (I,n)) MRRICEA = MR (I,n) + ELSE + MICROP = MICROP + MI (I,n) + MRCROP = MRCROP + MR (I,n) + ENDIF + END DO + + IF(MICROPA < MICROP) MICROPA = MICROP + IF(MRCROPA < MRCROP) MRCROPA = MRCROP + + ! GIA-GRIPC + ! ......... + + IF(GIAFRAC (I) <= 0 ) THEN + ! MASK OUT non-irrigated per GIA + RGRIPC (I) = RGRIPC (I) + IGRIPC (I) + PGRIPC (I) + IGRIPC (I) = 0. + PGRIPC (I) = 0. + ELSE + IF ((IGRIPC (I) + PGRIPC (I)) < 0.) THEN + ! GRIPC does not have data + PGRIPC (I) = 0. + IGRIPC (I) = GIAFRAC (I) + ELSE + ! GRIPC HAVE DATA + MICROP = PGRIPC (I) + IGRIPC (I) + 1.e-15 ! RGRIPC (I) + IF (GIAFRAC (I) > MICROP) THEN + PGRIPC (I) = PGRIPC (I) * GIAFRAC (I) / MICROP + IGRIPC (I) = IGRIPC (I) * GIAFRAC (I) / MICROP + ENDIF + ENDIF + ENDIF + + if ((RGRIPC(I) + IGRIPC(I) + PGRIPC(I) + NGRIPC(I)) > 0.) then + RGRIPC (I) = RGRIPC (I) /(RGRIPC(I) + IGRIPC(I) + PGRIPC(I) + NGRIPC(I)) + NGRIPC(I) = NGRIPC(I) /(RGRIPC(I) + IGRIPC(I) + PGRIPC(I) + NGRIPC(I)) + endif + + ! GIA-GRIPC-MIRCA IRRIGATED CROPS + ! ............................... + + IF (IGRIPC(I) == 0) THEN + DO N = 1, NCROPS + IF(N /= 3) MI (I,N) = 0. + END DO + ELSE + ! IF(MICROPA > IGRIPC (I)) THEN + ! + ! ! MIRCA is the larger fraction + ! IGRIPC (I) = MICROPA + ! ! CALL STOPIT (1, IGRIPC (I), MICROPA, MI (I,:)) + ! ELSE + + IF(MICROPA > 0.) THEN + + ! MIRCA has data too, thus scale crop fractions to match GIA-GRIPC (i.e. GIA) + DO N = 1, NCROPS + IF(N /= 3) MI (I,N) = MI (I,N) * IGRIPC (I) / MICROPA + END DO + ! CALL STOPIT (2, IGRIPC (I), MICROPA, MI (I,:)) + + ELSE + + ! MIRCA does not have data but GRIPC has + IF(MRCROPA > 0.) THEN + + ! Looks like MIRCA rainfed frac has data + DO N = 1, NCROPS + IF(N /= 3) THEN + MI (I,N) = MR (I,N) * IGRIPC (I) / MRCROPA + MR (I,N) = 0. + IRRG_DOY_PLANT (I,1,N) = RAINFEDPLANT (I,1,N) + IRRG_DOY_PLANT (I,2,N) = RAINFEDPLANT (I,2,N) + IRRG_DOY_HARVEST (I,1,N) = RAINFEDHARVEST (I,1,N) + IRRG_DOY_HARVEST (I,2,N) = RAINFEDHARVEST (I,2,N) + ENDIF + END DO + ! CALL STOPIT (3, IGRIPC (I), MICROPA, MI (I,:)) + MRCROPA = 0. + + ELSE + + ! MIRCA irrigated and rainfed do not have data plant some wheat + MI (I,1) = IGRIPC (I) + IRRG_DOY_PLANT (I,1,1) = 999 + IRRG_DOY_PLANT (I,2,1) = 0 + IRRG_DOY_HARVEST (I,1,1) = 999 + IRRG_DOY_HARVEST (I,2,1) = 0 + ! CALL STOPIT (4, IGRIPC (I), MICROPA, MI (I,:)) + ENDIF + ! ENDIF + ENDIF + ENDIF + + ! CALL STOPIT (5, IGRIPC (I), MICROPA, MI (I,:)) + + ! GIA-GRIPC-MIRCA PADDY + ! ..................... + IF(PGRIPC (I) == 0.) THEN + MI (I,3) = 0. + ELSE + ! IF(MIRICEA > PGRIPC (I)) THEN + ! + ! ! MIRCA is the larger fraction + ! PGRIPC (I) = MIRICEA + ! + ! ELSE + + IF(MIRICEA > 0.) THEN + + ! MIRCA has data too, thus scale crop fractions to match GRIPC + MI (I,3) = MI (I,3) * PGRIPC (I) / MIRICEA + + ELSE + + ! MIRCA does not have data but GRIPC has + IF(MRRICEA > 0.) THEN + + ! Looks like MIRCA rainfed frac has rice + MI (I,3) = MR (I,3) * PGRIPC (I) / MRRICEA + MR (I,3) = 0. + MRRICEA = 0. + IRRG_DOY_PLANT (I,1,3) = RAINFEDPLANT (I,1,3) + IRRG_DOY_PLANT (I,2,3) = RAINFEDPLANT (I,2,3) + IRRG_DOY_HARVEST (I,1,3) = RAINFEDHARVEST (I,1,3) + IRRG_DOY_HARVEST (I,2,3) = RAINFEDHARVEST (I,2,3) + ELSE + + ! MIRCA irrigated and rainfed do not have data plant rice to PGRIPC + MI (I,3) = PGRIPC (I) + + ! Get crop planting days for the nearest neighbor later + + IRRG_DOY_PLANT (I,1,3) = 999 + IRRG_DOY_PLANT (I,2,3) = 0 + IRRG_DOY_HARVEST (I,1,3) = 999 + IRRG_DOY_HARVEST (I,2,3) = 0 + + ENDIF + ENDIF + ! ENDIF + ENDIF + + ! GIA-GRIPC-MIRCA Rainfed CROPS + ! ............................. + + IF( RGRIPC (I) == 0.) THEN + MR (I,:) = 0. + ELSE + + ! IF(MRCROPA + MRRICEA > RGRIPC (I)) THEN + ! + ! ! MIRCA is the larger fraction + ! RGRIPC (I) = MRCROPA + MRRICEA + ! + ! ELSE + IF(MRCROPA + MRRICEA > 0.) THEN + + ! MIRCA has data too, thus scale crop fractions to match GRIPC + DO N = 1, NCROPS + IF(N /= 3) MR (I,N) = MR (I,N) * RGRIPC (I) / (MRCROPA + MRRICEA) + END DO + + IF (MRRICEA > 0.) MR (I,3) = MR (I,3) * RGRIPC (I) / (MRCROPA + MRRICEA) + + ELSE + + ! MIRCA does not have data but GRIPC has + ! MIRCA rainfed do not have data plant some wheat + MR (I,1) = RGRIPC (I) + RAINFEDPLANT (I,1,1) = 999 + RAINFEDPLANT (I,2,1) = 0 + RAINFEDHARVEST (I,1,1) = 999 + RAINFEDHARVEST (I,2,1) = 0 + ENDIF + ENDIF + ! ENDIF + + ! IRRG_TYPE + + DO N = 1, NCROPS + FF = FLOOD (I) + SF = SPRINKLER (I) + DF = DRIP (I) + IF(N == 3) THEN ! Rice + SF = 0. + DF = 0. + FF = 1. + IRRG_TYPE (I, N) = 3 ! Always flood + ENDIF + IF(N == 10) IRRG_TYPE (I, N) = -1 ! never sprinkler + IF(N == 22) IRRG_TYPE (I, N) = -1 ! never sprinkler + !IF(N == 10) SF = 0. ! Date palm + !IF(N == 22) SF = 0. ! Cocoa + !ITYPE = (/SF, DF, FF/) + !IRRG_TYPE (I, N) = maxloc(ITYPE, 1) + END DO + ENDIF + END DO + + ! CROP CALENDAR CAL_STEP2 Update missing plant harvest dates + ! ---------------------------------------------------------- + + DO I = 1, NTILES + DO N = 1,3,2 + + ! fill missing crop plant/harvest DOYs in irrigated crops + ! ....................................................... + + IF((IRRG_DOY_PLANT(I,1,N) == 999).AND.(MI (I,N) > 0.)) THEN + l = getNeighbor (I,day_in = IRRG_DOY_PLANT (:,1,N)) + IRRG_DOY_PLANT (I,1,N) = IRRG_DOY_PLANT (l,1,N) + IRRG_DOY_HARVEST(I,1,N) = IRRG_DOY_HARVEST(l,1,N) + if(N == 1) then + IF(RAINFEDPLANT(I,1,N) == 999) THEN + RAINFEDPLANT (I,1,N) = IRRG_DOY_PLANT (l,1,N) + RAINFEDHARVEST(I,1,N) = IRRG_DOY_HARVEST(l,1,N) + ENDIF + endif + endif + + IF(N == 1) THEN + ! fill missing crop plant/harvest DOYs in rainfed crops + ! ..................................................... + ! temperorily commented out to save time, because we don't irrigate here anyway + ! IF((RAINFEDPLANT(I,1,N) == 999).AND.(MR (I,N) > 0.)) THEN + ! print *,'RAINFEDPLANT(I,1,N)',I,N, MR (I,N) + ! l = getNeighbor (I,day_in = IRRG_DOY_PLANT (:,1,N)) + ! RAINFEDPLANT (I,1,N) = IRRG_DOY_PLANT (l,1,N) + ! RAINFEDHARVEST(I,1,N) = IRRG_DOY_HARVEST(l,1,N) + ! endif + ENDIF + END DO + if( ((IRRG_DOY_PLANT (I,1,1) == 999).and.(MI (i,1) > 0.)) .OR. & + ((IRRG_DOY_PLANT (I,1,3) == 999).and.(MI (i,3) > 0.)) ) then + print *, i, IRRG_DOY_PLANT (I,1,1:3), MI (i,1:3) + stop + endif + END DO + + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_IRRIGFRAC' ) ,(/1/),(/NTILES/),IGRIPC ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_PADDYFRAC' ) ,(/1/),(/NTILES/),PGRIPC ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'RAINFEDFRAC') ,(/1/),(/NTILES/),RGRIPC ) ; VERIFY_(STATUS) + + do n = 1,NCROPS + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_CROPIRRIGFRAC' ) ,(/1,n /),(/NTILES,1 /), MI (:,n) ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'CROPRAINFEDFRAC' ) ,(/1,n /),(/NTILES,1 /), MR (:,n) ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_TYPE' ) ,(/1,n /),(/NTILES,1 /), IRRG_TYPE (:,n) ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_DOY_PLANT' ) ,(/1,1,n/),(/NTILES,1,1/), IRRG_DOY_PLANT (:,1,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_DOY_PLANT' ) ,(/1,2,n/),(/NTILES,1,1/), IRRG_DOY_PLANT (:,2,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_DOY_HARVEST' ) ,(/1,1,n/),(/NTILES,1,1/), IRRG_DOY_HARVEST(:,1,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_DOY_HARVEST' ) ,(/1,2,n/),(/NTILES,1,1/), IRRG_DOY_HARVEST(:,2,n)) ; VERIFY_(STATUS) + + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'RAINFEDPLANT' ) ,(/1,1,n/),(/NTILES,1,1/), RAINFEDPLANT (:,1,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'RAINFEDPLANT' ) ,(/1,2,n/),(/NTILES,1,1/), RAINFEDPLANT (:,2,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'RAINFEDHARVEST' ) ,(/1,1,n/),(/NTILES,1,1/), RAINFEDHARVEST (:,1,n)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'RAINFEDHARVEST' ) ,(/1,2,n/),(/NTILES,1,1/), RAINFEDHARVEST (:,2,n)) ; VERIFY_(STATUS) + end do + + status = NF_CLOSE(NCOutID) + + END SUBROUTINE MergeData + + !---------------------------------------------------------------------------------------- + + SUBROUTINE ReadProcess_IMethod (NTILES, f_sprink, f_drip, f_flood) + + implicit none + + integer, INTENT (IN) :: NTILES + real, dimension (:), INTENT(INOUT) :: f_sprink, f_drip, f_flood + integer :: i,j,n, k, N_METHOD, cnt_code, st_code + character*2 :: ST_NAME + character*3 :: CNT_ABR + integer, parameter :: N_STATES = 50, N_COUNTRY = 256 + real :: s_dum, d_dum, f_dum + real, dimension (:), allocatable :: us_sprink, us_drip, us_flood, us_tarea + real, dimension (:), allocatable :: sprink, drip, flood, tarea + integer, dimension (:),pointer :: index_range + character*2, dimension (:),pointer :: ST_NAME_ABR + character*3, dimension (:),pointer :: CNT_NAME_ABR + + ! Read state fractions + print *, ' ' + print *, '.........................................................................' + print *, 'PROCESSING IRRIGATION METHOD DATA ' + + open (10, file = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/country_code_IMethod/v1/US_IMethod.2015' , form = 'formatted', status ='old', action = 'read') + open (11, file = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/country_code_IMethod/v1/Global_IMethod.data', form = 'formatted', status ='old', action = 'read') + + READ (11, *) N_METHOD + + call get_country_codes (index_range=index_range, ST_NAME_ABR = ST_NAME_ABR, CNT_NAME_ABR = CNT_NAME_ABR) + + allocate (sprink (0:N_COUNTRY)) + allocate (drip (0:N_COUNTRY)) + allocate (flood (0:N_COUNTRY)) + allocate (tarea (0:N_COUNTRY)) + allocate (us_sprink(1:N_STATES )) + allocate (us_drip (1:N_STATES )) + allocate (us_flood (1:N_STATES )) + allocate (us_tarea (1:N_STATES )) + + sprink = 0. + drip = 0. + flood = 0. + tarea = 0. + us_sprink = 0. + us_drip = 0. + us_flood = 0. + us_tarea = 0. + + do i = 1, N_METHOD + read (11, *) CNT_ABR,s_dum, d_dum, f_dum + do k = 1, N_COUNTRY + if(CNT_ABR == CNT_NAME_ABR(k)) then + sprink(index_range(k)) = s_dum + drip (index_range(k)) = d_dum + flood (index_range(k)) = f_dum + endif + end do + end do + + tarea = sprink + drip + flood + where (tarea > 0.) + sprink = sprink / tarea + drip = drip / tarea + flood = flood / tarea + endwhere + + do i = 1, N_STATES + read (10, *) ST_NAME,s_dum, d_dum, f_dum + do k = 1, N_STATES + if(ST_NAME == ST_NAME_ABR(k)) then + us_sprink(k) = s_dum + us_drip (k) = d_dum + us_flood (k) = f_dum + endif + end do + end do + + us_tarea = us_sprink + us_drip + us_flood + where (us_tarea > 0.) + us_sprink = us_sprink / us_tarea + us_drip = us_drip / us_tarea + us_flood = us_flood / us_tarea + endwhere + + close (10, status = 'keep') + close (11, status = 'keep') + + ! map irrig method fractions + + open (10,file='clsm/country_and_state_code.data', & + form='formatted',status='old', action = 'read') + + ! allocate and initialize + + f_flood = 1. + f_sprink = 0. + f_drip = 0. + + tile_loop : do n = 1, NTILES + read (10, '(i10, 2I4)') j, cnt_code, st_code + if (cnt_code < 257) then + if(tarea (cnt_code) > 0.) then + f_flood (n) = flood (cnt_code) + f_sprink(n) = sprink(cnt_code) + f_drip (n) = drip (cnt_code) + endif + endif + + ! overwrite with US state methods + if(st_code /= 999) then + f_flood (n) = us_flood (st_code) + f_sprink(n) = us_sprink(st_code) + f_drip (n) = us_drip (st_code) + endif + end do tile_loop + + close (10, status = 'keep') + + CALL update_IMethod_bycounty (NTILES, f_sprink, f_drip, f_flood) + + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_IRRIGFRAC_SPR') ,(/1/),(/NTILES/), f_sprink) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_IRRIGFRAC_DRP' ) ,(/1/),(/NTILES/), f_drip ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCOutID,VarID(NCOutID,'IRRG_IRRIGFRAC_FRW' ) ,(/1/),(/NTILES/), f_flood ) ; VERIFY_(STATUS) + + print *, 'DONE PROCESSING IRRIGATION METHOD DATA ' + + deallocate (us_sprink, us_drip, us_flood, us_tarea, sprink, drip, flood, tarea) + + END SUBROUTINE ReadProcess_IMethod + + !---------------------------------------------------------------------------------------- + + SUBROUTINE update_IMethod_bycounty (NTILES, f_sprink, f_drip, f_flood) + + implicit none + integer, INTENT (IN) :: NTILES + real, dimension (:), INTENT(INOUT) :: f_sprink, f_drip, f_flood + integer, parameter :: NX_cb = 43200, NY_cb = 21600, NY_cbData = 10800 + integer, parameter :: cb_states = 72, cb_county = 900, cb_countyUS = 3220 + integer :: i,j, n, status, ncid, I0(1), j0(1),SS, CCC + real, dimension(:,:),allocatable :: SFR, DFR, FFR + integer, dimension (:),allocatable :: GEOID + integer, dimension(:,:),allocatable :: POLYID + real (kind =8) :: XG(NX_cb),YG(NY_cb), y0, x0, dxy + integer :: ii(NX_cb),jj(NY_cb) + + allocate (SFR (1:cb_county,1:cb_states)) + allocate (DFR (1:cb_county,1:cb_states)) + allocate (FFR (1:cb_county,1:cb_states)) + allocate (GEOID (1:cb_countyUS)) + allocate (POLYID(1:NX_cb,1:NY_cb)) + + POLYID = -9999 + + status = NF_OPEN ('/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/fraction_drip_flood_sprinkler/v1/cb_2015_us_county_30arcsec.nc4',NF_NOWRITE, ncid) ; VERIFY_(STATUS) + do j = 1, NY_cbData + status = NF_GET_VARA_INT(NCID,VarID(NCID,'POLYID') ,(/1,j/),(/NX_cb, 1/), POLYID (:,NY_cb - j + 1)) ; VERIFY_(STATUS) ! reading north to south + end do + do j = 1, cb_states + status = NF_GET_VARA_REAL(NCID,VarID(NCID,'IRRG_IRRIGFRAC_SPR') ,(/1,j/),(/cb_county, 1/), SFR (:,j)) ; VERIFY_(STATUS) + status = NF_GET_VARA_REAL(NCID,VarID(NCID,'IRRG_IRRIGFRAC_DRP' ) ,(/1,j/),(/cb_county, 1/), DFR (:,j)) ; VERIFY_(STATUS) + status = NF_GET_VARA_REAL(NCID,VarID(NCID,'IRRG_IRRIGFRAC_FRW' ) ,(/1,j/),(/cb_county, 1/), FFR (:,j)) ; VERIFY_(STATUS) + end do + status = NF_GET_VARA_INT(NCID,VarID(NCID,'GEOID' ) ,(/1/),(/cb_countyUS/), GEOID) ; VERIFY_(STATUS) + status = NF_CLOSE(NCID) ; VERIFY_(STATUS) + + dxy = 360.d0/NX_cb + do i = 1, NX_cb + xg(i) = (i-1)*dxy -180.d0 + dxy/2.d0 + end do + do i = 1, NY_cb + yg(i) = (i-1)*dxy -90.d0 + dxy/2.d0 + end do + + do n = 1, NTILES + + x0 = dble (tile_lon(n)) + y0 = dble (tile_lat(n)) + II = 0 + JJ = 0 + WHERE ((xg >= x0).and.(xg < x0 + dxy)) II = 1 + WHERE ((yg >= y0).and.(yg < y0 + dxy)) JJ = 1 + + I0 = FINDLOC(II,1) + J0 = FINDLOC(JJ,1) + + if((POLYID(I0(1), J0(1)) >= 1).AND.(POLYID(I0(1), J0(1)) <= cb_countyUS)) then + SS = GEOID(POLYID(I0(1),J0(1))) / 1000 + CCC= GEOID(POLYID(I0(1),J0(1))) - SS*1000 + f_sprink (n) = SFR (CCC,SS) + f_drip (n) = DFR (CCC,SS) + f_flood (n) = FFR (CCC,SS) + endif + + END DO + + deallocate (SFR, DFR, FFR, GEOID, POLYID) + + END SUBROUTINE update_IMethod_bycounty + + !---------------------------------------------------------------------------------------- + + SUBROUTINE ReadProcess_MIRCA (NC, NR, NTILES, tile_id, MIFRAC, MRFRAC) + + implicit none + + INTEGER, INTENT (IN) :: NTILES, NC, NR + INTEGER, INTENT (IN), DIMENSION(:,:):: tile_id + REAL,DIMENSION(:,:,:),INTENT(INOUT) :: MIFRAC, MRFRAC + character*2 :: TT + integer :: i, j, n, m + real,dimension (:), allocatable :: read_ir , read_rn, cnt_pix1, cnt_pix2 + real,dimension (:,:,:), allocatable :: mon_ir , mon_rn + real, pointer, dimension (:,:) :: var_raster1, var_raster2 + real,parameter :: radius = MAPL_radius, pi = MAPL_PI + real :: D2R, latc, lonc, area + + ! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ + ! + ! NOTE: "!$" is for conditional compilation + ! + logical :: running_omp = .false. + ! + !$ integer :: omp_get_thread_num, omp_get_num_threads + ! + integer :: n_threads=1, li, ui, t_count, n_used + ! + integer, dimension(:), allocatable :: low_ind, upp_ind + ! + ! ------------------------------------------------------------------ + + ! ----------- OpenMP PARALLEL ENVIRONMENT ---------------------------- + ! + ! FIND OUT WHETHER -omp FLAG HAS BEEN SET DURING COMPILATION + ! + !$ running_omp = .true. ! conditional compilation + ! + ! ECHO BASIC OMP VARIABLES + ! + !$OMP PARALLEL DEFAULT(NONE) SHARED(running_omp,n_threads) + ! + !$OMP SINGLE + ! + !$ n_threads = omp_get_num_threads() + ! + !$ write (*,*) 'running_omp = ', running_omp + !$ write (*,*) + !$ write (*,*) 'parallel OpenMP with ', n_threads, 'threads' + !$ write (*,*) + !$OMP ENDSINGLE + ! + !$OMP CRITICAL + !$ write (*,*) 'thread ', omp_get_thread_num(), ' alive' + !$OMP ENDCRITICAL + ! + !$OMP BARRIER + ! + !$OMP ENDPARALLEL + + MIFRAC = 0. + MRFRAC = 0. + + n_used = MIN(n_threads, NCROPS/2) + if(NTILES > 6000000) n_used = 1 ! otherwise multiple threads will run out of virtual memory + allocate(low_ind(n_used)) + allocate(upp_ind(n_used)) + low_ind(1) = 1 + upp_ind(n_used) = NCROPS + + if (running_omp) then + do i=1,n_used-1 + upp_ind(i) = low_ind(i) + (NCROPS/n_used) - 1 + low_ind(i+1) = upp_ind(i) + 1 + end do + end if + + D2R= PI/180. + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP COPYIN (read_ir, read_rn,var_raster1, var_raster2, mon_ir, mon_rn, cnt_pix1, cnt_pix2) & + !$OMP SHARED( n_used, low_ind, upp_ind, tile_id, NTILES, MIFRAC,MRFRAC, NC, NR, D2R) & + !$OMP PRIVATE(n,i,j,m,tt,t_count, read_ir, read_rn,var_raster1, var_raster2, mon_ir, mon_rn, cnt_pix1, cnt_pix2, latc, area) + + allocate (read_ir (1:NX_mirca*12)) + allocate (read_rn (1:NX_mirca*12)) + allocate (mon_ir (1:NX_mirca,1:NY_mirca,1:12)) + allocate (mon_rn (1:NX_mirca,1:NY_mirca,1:12)) + allocate (var_raster1 (1:nc,1:nr)) + allocate (var_raster2 (1:nc,1:nr)) + allocate (cnt_pix1 (1:NTILES)) + allocate (cnt_pix2 (1:NTILES)) + + !$OMP DO + + DO t_count = 1,n_used + + CROP_TYPE : DO n = low_ind(t_count),upp_ind(t_count) + + write (tt, '(i2.2)') n + + print *, ' ' + print *, '.........................................................................' + print *, 'PROCESSING MIRCA : crop_', tt,'_irrigated_12.flt' + print *, 'PROCESSING MIRCA : crop_', tt,'_rainfed_12.flt' + + mon_ir = UNDEF + mon_rn = UNDEF + + open (50 + t_count, file = trim(MIRCA_pathIrr)//tt//'_irrigated_12.flt', action = 'read', & + form = 'unformatted', access='direct', recl=NX_mirca*12) + open (100 + t_count, file = trim(MIRCA_pathRain)//tt//'_rainfed_12.flt', action = 'read', & + form = 'unformatted', access='direct', recl=NX_mirca*12) + + mirca_rows : do j = 1, NY_mirca + read(50 + t_count,rec= NY_mirca - J + 1) read_ir + read(100 + t_count,rec= NY_mirca - J + 1) read_rn + latc = lat1_mirca - (j-1) * DXY_MIRCA + area = (sin(d2r*(latc+0.5*dxy_mirca)) - sin(d2r*(latc-0.5*dxy_mirca)))*(dxy_mirca*d2r) + area = area * radius * radius / 10000. ! in ha + + mirca_COLS : do i = 1, NX_MIRCA + mon_ir (i,j,:) = read_ir ((i-1)*12 + 1: (i-1)*12 + 12)/area + mon_rn (i,j,:) = read_rn ((i-1)*12 + 1: (i-1)*12 + 12)/area + end do mirca_cols + + end do mirca_rows + + close (50 + t_count, status ='keep') + close (100 + t_count, status ='keep') + + ! Grid 2 tile + ! ----------- + do m = 1,12 + + var_raster1 = 0. + var_raster2 = 0. + + call RegridRasterReal(mon_ir(:,:,m),var_raster1) + call RegridRasterReal(mon_rn(:,:,m),var_raster2) + + cnt_pix1 = 0. + cnt_pix2 = 0. + + do j = 1,nr + do i = 1,nc + if((var_raster1 (i,j) > 0.).and.(tile_id (i,j) >= 1).AND.(tile_id (i,j) <= NTILES)) then + MIFRAC (tile_id(i,j),m,n) = MIFRAC (tile_id(i,j),m,n) + var_raster1 (i,j) + cnt_pix1 (tile_id(i,j)) = cnt_pix1 (tile_id(i,j)) + 1. + endif + if((var_raster2 (i,j) > 0.).and.(tile_id (i,j) >= 1).AND.(tile_id (i,j) <= NTILES)) then + MRFRAC (tile_id(i,j),m,n) = MRFRAC (tile_id(i,j),m,n) + var_raster2 (i,j) + cnt_pix2 (tile_id(i,j)) = cnt_pix2 (tile_id(i,j)) + 1. + endif + end do + end do + + do i = 1, NTILES + if(cnt_pix1(i) > 0.) MIFRAC (i,m,n) = MIFRAC (i,m,n)/cnt_pix1(i) + if(cnt_pix2(i) > 0.) MRFRAC (i,m,n) = MRFRAC (i,m,n)/cnt_pix2(i) + end do + + end do + END DO CROP_TYPE + END DO + !$OMP END DO + !$OMP END PARALLEL + + print *,'DONE MIRCA PROCESSING' + END SUBROUTINE ReadProcess_MIRCA + + !---------------------------------------------------------------------------------------- + + SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NGRIPC) + + INTEGER, INTENT (IN) :: NTILES, NC, NR + INTEGER, INTENT (IN), DIMENSION(:,:):: tile_id + REAL,DIMENSION(:),INTENT(INOUT) :: IGRIPC, RGRIPC, PGRIPC, NGRIPC + real, allocatable :: var_in (:,:), tot_cnt (:), min_cnt(:), max_cnt(:) + integer :: i,j, n, r, status, DOY, NCID, NCIDW,xid, yid,vid + integer, pointer :: iraster (:,:) + character*3 :: DDD + integer*2, allocatable, dimension (:,:) :: Lai_clim + real, allocatable, dimension (:,:) :: clim_min, clim_max, clim_lai + real, allocatable, dimension (:) :: LAI_MIN, LAI_MAX, lai + real,allocatable :: dum,yr,mn,dy,nt + logical :: write_lai = .false. + + !V2: Min Max LAI from LAI Climatology for consistence + + allocate( var_in(NX_gripc,NY_gripc)) + var_in = UNDEF + + open ( 10, file = trim(GRIPC_file), form = 'unformatted', access='direct', recl=(NX_gripc)) + + !- Read input file:: + + do j = 1, NY_gripcdata + r = NY_gripc -j + 1 + read(10,rec=j) var_in(:, r) + do i = 1, NX_gripc + if( var_in(i, r) == 0. ) var_in(i, r) = -9999. + if( var_in(i, r) == 4. ) var_in(i, r) = -9999. + end do + end do + close( 10 ) + + allocate(iraster(NX_gripc,NY_gripc),stat=STATUS); VERIFY_(STATUS) + call RegridRaster(tile_id,iraster) + + allocate (tot_cnt (1:ntiles)) + allocate (min_cnt (1:ntiles)) + allocate (max_cnt (1:ntiles)) + + RGRIPC = 0. + IGRIPC = 0. + PGRIPC = 0. + tot_cnt = 0. + min_cnt = 0. + max_cnt = 0. + allocate (LAI_MIN (NTILES)) + allocate (LAI_MAX (NTILES)) + allocate (lai (NTILES)) + LAI_MIN = -9999. + LAI_MAX = -9999. + + do j = 1,NY_gripc + do i = 1,NX_gripc + if((iraster (i,j) >=1).and.(iraster (i,j) <=ntiles)) then + tot_cnt (iraster (i,j)) = tot_cnt (iraster (i,j)) + 1. + if (var_in(i,j) == 1) RGRIPC(iraster (i,j)) = RGRIPC(iraster (i,j)) + 1. + if (var_in(i,j) == 2) IGRIPC(iraster (i,j)) = IGRIPC(iraster (i,j)) + 1. + if (var_in(i,j) == 3) PGRIPC(iraster (i,j)) = PGRIPC(iraster (i,j)) + 1. + if (var_in(i,j) == 4) NGRIPC(iraster (i,j)) = NGRIPC(iraster (i,j)) + 1. + + endif + end do + end do + + RGRIPC = RGRIPC / tot_cnt + IGRIPC = IGRIPC / tot_cnt + PGRIPC = PGRIPC / tot_cnt + NGRIPC = NGRIPC / tot_cnt + + allocate(dum) + allocate(yr) + allocate(mn) + allocate(dy) + allocate(nt) + + open (43,file='clsm/lai.dat', & + form='unformatted',status='unknown',convert='little_endian',action='read') + LAI_MAX=-9999. + LAI_MIN=9999. + do j = 1, 48 + read(43) yr,mn,dy,dum,dum,dum,yr,mn,dy,dum,dum,dum,nt,dum + read(43) lai + do i = 1, NTILES + if (lai(i)>LAI_MAX(i)) LAI_MAX(i)=lai(i) + if (lai(i) 0).and.(tile_id (i,j) <= NTILES)) then + cnt_pix1(tile_id(i,j)) = cnt_pix1(tile_id(i,j)) + 1. + if((irrig (i,j) > 0) .AND. (irrig (i,j) < 4)) then + GIAFRAC(tile_id (i,j)) = GIAFRAC(tile_id (i,j)) + 1. + endif + endif + enddo geos_cols + end do geos_rows + where (cnt_pix1 > 0) GIAFRAC = GIAFRAC / cnt_pix1 + deallocate (irrig, var_in) + deallocate (cnt_pix1) + print *,'DONE PROCESSING GIA' + RETURN + END SUBROUTINE ReadProcess_GIA + + ! ---------------------------------------------------------------------- + + integer function VarID (NCFID, VNAME) + + integer, intent (in) :: NCFID + character(*), intent (in) :: VNAME + integer :: status + + STATUS = NF_INQ_VARID (NCFID, trim(VNAME) ,VarID) + IF (STATUS .NE. NF_NOERR) & + CALL HANDLE_ERR(STATUS, trim(VNAME)) + + end function VarID + + ! ----------------------------------------------------------------------- + + SUBROUTINE HANDLE_ERR(STATUS, Line) + + INTEGER, INTENT (IN) :: STATUS + CHARACTER(*), INTENT (IN) :: Line + + IF (STATUS .NE. NF_NOERR) THEN + PRINT *, trim(Line),': ',NF_STRERROR(STATUS) + STOP 'Stopped' + ENDIF + + END SUBROUTINE HANDLE_ERR + + ! ----------------------------------------------------------------------------- + + integer function getNeighbor (tid_in, lai_in, day_in) + + implicit none + integer, intent (in) :: tid_in + real, optional, dimension (NTILES) :: lai_in, day_in + integer :: i, nplus + logical :: tile_found + logical, allocatable, dimension (:) :: mask + integer, allocatable, dimension (:) :: sub_tid + real , allocatable, dimension (:) :: sub_lon, sub_lat, rev_dist + real :: dw, min_lon, max_lon, min_lat, max_lat + integer, allocatable, dimension (:) :: TILEID + + allocate (mask (1: NTILES)) + allocate (TILEID (1: NTILES)) + forall (i=1:NTILES) TILEID (i) = i + + dw = 0.5 + getNeighbor = -9999 + + ZOOMOUT : do + + tile_found = .false. + + ! Min/Max lon/lat of the working window + ! ------------------------------------- + + min_lon = MAX(tile_lon (tid_in) - dw, -180.) + max_lon = MIN(tile_lon (tid_in) + dw, 180.) + min_lat = MAX(tile_lat (tid_in) - dw, -90.) + max_lat = MIN(tile_lat (tid_in) + dw, 90.) + + mask = .false. + if(present (lai_in)) then + mask = ((tile_lat >= min_lat .and. tile_lat <= max_lat).and.(tile_lon >= min_lon .and. tile_lon <= max_lon).and.(lai_in >= 0.)) + endif + if(present (day_in)) then + mask = ((tile_lat >= min_lat .and. tile_lat <= max_lat).and.(tile_lon >= min_lon .and. tile_lon <= max_lon).and.(day_in < 998)) + endif + nplus = count(mask = mask) + + if(nplus < 0) then + dw = dw + 0.5 + CYCLE + endif + + allocate (sub_tid (1:nplus)) + allocate (sub_lon (1:nplus)) + allocate (sub_lat (1:nplus)) + allocate (rev_dist (1:nplus)) + + sub_tid = PACK (TILEID , mask= mask) + sub_lon = PACK (tile_lon, mask= mask) + sub_lat = PACK (tile_lat, mask= mask) + + ! compute distance from the tile + + sub_lat = sub_lat * MAPL_PI/180. + sub_lon = sub_lon * MAPL_PI/180. + + SEEK : if(getNeighbor < 0) then + + rev_dist = 1.e20 + + do i = 1,nplus + + rev_dist(i) = haversine(to_radian(tile_lat(tid_in)), to_radian(tile_lon(tid_in)), & + sub_lat(i), sub_lon(i)) + + end do + + FOUND : if(minval (rev_dist) < 1.e19) then + if(present (lai_in)) then + if(lai_in(sub_tid(minloc(rev_dist,1))) >= 0.) then + getNeighbor = sub_tid(minloc(rev_dist,1)) + tile_found = .true. + endif + endif + if(present (day_in)) then + if(day_in(sub_tid(minloc(rev_dist,1))) < 998) then + getNeighbor = sub_tid(minloc(rev_dist,1)) + tile_found = .true. + endif + endif + endif FOUND + + endif SEEK + + deallocate (sub_tid, sub_lon, sub_lat, rev_dist) + + if(tile_found) GO TO 100 + + ! if not increase the window size + dw = dw + 0.5 + + end do ZOOMOUT + +100 continue + + deallocate (mask) + + end function getNeighbor + + ! ***************************************************************************** + + ! IRRGRR - duplicates same in Utils/mk_restarts/getids.F90 + + function to_radian(degree) result(rad) + + real,intent(in) :: degree + real :: rad + + rad = degree*MAPL_PI/180. + + end function to_radian + + ! ***************************************************************************** + + ! IRRGRR - duplicates same in Utils/mk_restarts/getids.F90 + + real function haversine(deglat1,deglon1,deglat2,deglon2) + ! great circle distance -- adapted from Matlab + real,intent(in) :: deglat1,deglon1,deglat2,deglon2 + real :: a,c, dlat,dlon,lat1,lat2 + real,parameter :: radius = MAPL_radius + + ! dlat = to_radian(deglat2-deglat1) + ! dlon = to_radian(deglon2-deglon1) + ! lat1 = to_radian(deglat1) + ! lat2 = to_radian(deglat2) + dlat = deglat2-deglat1 + dlon = deglon2-deglon1 + lat1 = deglat1 + lat2 = deglat2 + a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2 + if(a>=0. .and. a<=1.) then + c = 2*atan2(sqrt(a),sqrt(1-a)) + haversine = radius*c / 1000. + else + haversine = 1.e20 + endif + end function + + ! ***************************************************************************** + end subroutine create_irrig_params + + END module module_irrig_params + +! ----------------------------------------------------------------------------- + +!PROGRAM irrig_model +! +! use module_irrig_params +! +! call create_irrig_params (43200, 21600, 'rst/SMAP_EASEv2_M36_964x406') +! +!END PROGRAM irrig_model diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 index 307934f61..d9a31a1f2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 @@ -61,6 +61,7 @@ module rmTinyCatchParaMod logical, public, save :: use_PEATMAP = .false. logical, public, save :: jpl_height = .false. + logical, public, save :: IRRIGBCS = .false. character*8, public, save :: LAIBCS = 'UNDEF' character*6, public, save :: SOILBCS = 'UNDEF' character*6, public, save :: MODALB = 'UNDEF' @@ -242,6 +243,17 @@ SUBROUTINE init_bcs_config(LBCSV) use_PEATMAP = .true. jpl_height = .true. + case ("v14") + LAIBCS = 'MODGEO' + SOILBCS = 'HWSD_b' + MODALB = 'MODIS2' + SNOWALB = 'MODC061v2' + OUTLETV = "v2" + GNU = 1.0 + use_PEATMAP = .true. + jpl_height = .true. + IRRIGBCS = .true. + case default print *,'init_bcs_config(): unknown land boundary conditions version (LBCSV)'