diff --git a/CHANGELOG.md b/CHANGELOG.md index 949a160c..54b7f09e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added functionality to simulate landice tiles. +- Added functionality to read nc4-formatted tile file. + ### Changed ### Fixed diff --git a/CMakeLists.txt b/CMakeLists.txt index ffb1ef47..2de54c45 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ esma_add_library(${this} SRCS GEOS_LdasGridComp.F90 SUBCOMPONENTS ${alldirs} SUBDIRS LDAS_Shared - DEPENDENCIES GEOSland_GridComp makebcs MAPL + DEPENDENCIES GEOSland_GridComp GEOSlandice_GridComp makebcs MAPL INCLUDES ${INC_ESMF}) esma_add_subdirectory(GEOSldas_App) diff --git a/GEOS_LdasGridComp.F90 b/GEOS_LdasGridComp.F90 index 4c1b13ed..cf481aa1 100644 --- a/GEOS_LdasGridComp.F90 +++ b/GEOS_LdasGridComp.F90 @@ -9,27 +9,28 @@ module GEOS_LdasGridCompMod use ESMF use MAPL_Mod - use GEOS_MetforceGridCompMod, only: MetforceSetServices => SetServices - use GEOS_LandGridCompMod, only: LandSetServices => SetServices - use GEOS_LandPertGridCompMod, only: LandPertSetServices => SetServices - use GEOS_EnsGridCompMod, only: EnsSetServices => SetServices + use GEOS_MetforceGridCompMod, only: MetforceSetServices => SetServices + use GEOS_LandGridCompMod, only: LandSetServices => SetServices + use GEOS_LandPertGridCompMod, only: LandPertSetServices => SetServices + use GEOS_EnsGridCompMod, only: EnsSetServices => SetServices use GEOS_LandAssimGridCompMod, only: LandAssimSetServices => SetServices + use GEOS_LandiceGridCompMod, only: LandiceSetServices => SetServices - use EASE_conv, only: ease_inverse - use LDAS_TileCoordType, only: tile_coord_type , T_TILECOORD_STATE, TILECOORD_WRAP - use LDAS_TileCoordType, only: grid_def_type, io_grid_def_type, operator (==) + use EASE_conv, only: ease_inverse + use LDAS_TileCoordType, only: tile_coord_type , T_TILECOORD_STATE, TILECOORD_WRAP + use LDAS_TileCoordType, only: grid_def_type, io_grid_def_type, operator (==) use LDAS_TileCoordRoutines, only: get_minExtent_grid, get_ij_ind_from_latlon, io_domain_files - use LDAS_ConvertMod, only: esmf2ldas - use LDAS_PertRoutinesMod, only: get_pert_grid - use LDAS_ensdrv_functions,ONLY: get_io_filename - use LDAS_DateTimeMod,ONLY: date_time_type - use LDAS_ensdrv_mpi, only: MPI_tile_coord_type, MPI_grid_def_type - use LDAS_ensdrv_mpi, only: init_MPI_types,mpicomm,numprocs,myid - use LDAS_ensdrv_mpi, only: root_proc - use LDAS_ensdrv_Globals, only: logunit,logit,root_logit,echo_clsm_ensdrv_glob_param, get_ensid_string - use catch_constants, only: echo_catch_constants - use StieglitzSnow, only: StieglitzSnow_echo_constants - use SurfParams, only: SurfParams_init + use LDAS_ConvertMod, only: esmf2ldas + use LDAS_PertRoutinesMod, only: get_pert_grid + use LDAS_ensdrv_functions, only: get_io_filename + use LDAS_DateTimeMod, only: date_time_type + use LDAS_ensdrv_mpi, only: MPI_tile_coord_type, MPI_grid_def_type + use LDAS_ensdrv_mpi, only: init_MPI_types,mpicomm,numprocs,myid + use LDAS_ensdrv_mpi, only: root_proc + use LDAS_ensdrv_Globals, only: logunit,logit,root_logit,echo_clsm_ensdrv_glob_param, get_ensid_string + use catch_constants, only: echo_catch_constants + use StieglitzSnow, only: StieglitzSnow_echo_constants + use SurfParams, only: SurfParams_init implicit none @@ -50,6 +51,7 @@ module GEOS_LdasGridCompMod ! All children integer,allocatable :: LAND(:) + integer,allocatable :: LANDICE(:) integer,allocatable :: LANDPERT(:) integer,allocatable :: METFORCE(:) integer :: ENSAVG, LANDASSIM @@ -59,6 +61,8 @@ module GEOS_LdasGridCompMod logical :: land_assim logical :: mwRTM logical :: ensemble_forcing ! switch between deterministic and ensemble forcing + logical :: with_landice ! true if landice tiles requested by config + logical :: with_land ! true if land tiles requested by config contains @@ -87,14 +91,16 @@ subroutine SetServices(gc, rc) character(len=ESMF_MAXSTR) :: ensid_string,childname character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file, ENS_FORCING_STR integer :: ens_id_width + integer, allocatable :: tile_types(:) ! Local variables type(T_TILECOORD_STATE), pointer :: tcinternal type(TILECOORD_WRAP) :: tcwrap type(ESMF_Config) :: CF - integer :: LSM_CHOICE + integer :: LSM_CHOICE, nra integer :: FIRST_ENS_ID - + logical :: isPresent + ! Begin... ! Get my name and setup traceback handle @@ -150,15 +156,42 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) ensemble_forcing = (trim(ENS_FORCING_STR) == 'YES') + call ESMF_ConfigFindLabel( CF, LABEL="TILE_TYPES:", isPresent=isPresent, _RC) + if (isPresent) then + nra = ESMF_ConfigGetLen( CF, _RC) + allocate(tile_types(nra)) + call ESMF_ConfigFindLabel( CF, LABEL="TILE_TYPES:", _RC) + call ESMF_ConfigGetAttribute( CF, valueList=tile_types, count=nra, _RC) + else + ! default + tile_types = [MAPL_LAND] + endif + + with_landice = .false. + with_land = .false. +! with_lake = .false. + + if (any(tile_types == MAPL_LANDICE)) with_landice = .true. + if (any(tile_types == MAPL_LAND )) with_land = .true. +! if (any(tile_types == MAPL_LAKE )) with_lake = .true. + call MAPL_GetResource ( MAPL, LAND_ASSIM_STR, Label="LAND_ASSIM:", DEFAULT="NO", RC=STATUS) VERIFY_(STATUS) LAND_ASSIM_STR = ESMF_UtilStringUpperCase(LAND_ASSIM_STR, rc=STATUS) VERIFY_(STATUS) land_assim = (trim(LAND_ASSIM_STR) /= 'NO') + if (land_assim .and. .not. with_land) then + _ASSERT( .false., "Inconsistent configuration: Land assimilation requested but no land tiles requested.") + endif + call MAPL_GetResource ( MAPL, mwRTM_file, Label="LANDASSIM_INTERNAL_RESTART_FILE:", DEFAULT='', RC=STATUS) VERIFY_(STATUS) mwRTM = ( len_trim(mwRTM_file) /= 0 ) + if (mwRTM .and. .not. with_land) then + print*, "Warning: Requested mwRTM file is not necessary because no land tiles requested." + mwRTM = .false. + endif call MAPL_GetResource ( MAPL, LSM_CHOICE, Label="LSM_CHOICE:", DEFAULT=1, RC=STATUS) if (LSM_CHOICE /=1 ) then @@ -171,7 +204,8 @@ subroutine SetServices(gc, rc) allocate(METFORCE(1)) endif - allocate(LAND(NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE)) + if (with_land) allocate(LAND( NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE)) + if (with_landice) allocate(LANDICE(NUM_ENSEMBLE)) ! ens_id_with = 2 + number of digits = total number of chars in ensid_string ("_eXXXX") ! @@ -205,105 +239,89 @@ subroutine SetServices(gc, rc) call get_ensid_string(ensid_string, ens_id, ens_id_width, NUM_ENSEMBLE) - childname='LANDPERT'//trim(ensid_string) - LANDPERT(i) = MAPL_AddChild(gc, name=childname, ss=LandPertSetServices, rc=status) - VERIFY_(status) + if (with_land) then + childname='LANDPERT'//trim(ensid_string) + LANDPERT(i) = MAPL_AddChild(gc, name=childname, ss=LandPertSetServices, rc=status) + VERIFY_(status) - childname='LAND'//trim(ensid_string) - LAND(i) = MAPL_AddChild(gc, name=childname, ss=LandSetServices, rc=status) - VERIFY_(status) + childname='LAND'//trim(ensid_string) + LAND(i) = MAPL_AddChild(gc, name=childname, ss=LandSetServices, rc=status) + VERIFY_(status) + endif + + if (with_landice) then + childname='LANDICE'//trim(ensid_string) + LANDICE(i) = MAPL_AddChild(gc, name=childname, ss=LandiceSetServices, rc=status) + VERIFY_(status) + endif enddo - ENSAVG = MAPL_AddChild(gc, name='ENSAVG', ss=EnsSetServices, rc=status) - VERIFY_(status) - - if(land_assim .or. mwRTM ) then - LANDASSIM = MAPL_AddChild(gc, name='LANDASSIM', ss=LandAssimSetServices, rc=status) + if (with_land) then + ENSAVG = MAPL_AddChild(gc, name='ENSAVG', ss=EnsSetServices, rc=status) VERIFY_(status) - endif + + if(land_assim .or. mwRTM ) then + LANDASSIM = MAPL_AddChild(gc, name='LANDASSIM', ss=LandAssimSetServices, rc=status) + VERIFY_(status) + endif - ! Connections - do i=1,NUM_ENSEMBLE - ! -METFORCE-feeds-LANDPERT's-imports- - k = 1 - if ( ensemble_forcing ) k = i - call MAPL_AddConnectivity( & - gc, & - SHORT_NAME = ['Tair ', 'Qair ', 'Psurf ', 'Rainf_C', 'Rainf ', & - 'Snowf ', 'LWdown ', 'SWdown ', 'PARdrct', 'PARdffs', & - 'Wind ', 'RefH '], & - SRC_ID = METFORCE(k), & - DST_ID = LANDPERT(i), & - rc = status & - ) - VERIFY_(status) - ! -LANDPERT-feeds-LAND's-imports- - call MAPL_AddConnectivity( & - gc, & - SRC_NAME = ['TApert ', 'QApert ', 'UUpert ', & - 'UWINDLMTILEpert', 'VWINDLMTILEpert', 'PCUpert ', & - 'PLSpert ', 'SNOpert ', 'DRPARpert ', & - 'DFPARpert ', 'DRNIRpert ', 'DFNIRpert ', & - 'DRUVRpert ', 'DFUVRpert ', 'LWDNSRFpert '], & - SRC_ID = LANDPERT(i), & - DST_NAME = ['TA ', 'QA ', 'UU ', 'UWINDLMTILE',& - 'VWINDLMTILE', 'PCU ', 'PLS ', 'SNO ',& - 'DRPAR ', 'DFPAR ', 'DRNIR ', 'DFNIR ',& - 'DRUVR ', 'DFUVR ', 'LWDNSRF '], & - DST_ID = LAND(i), & - rc = status & - ) + ! Connections + do i=1,NUM_ENSEMBLE + k = 1 + if ( ensemble_forcing ) k = i + ! -LANDPERT-feeds-LAND's-imports- + call MAPL_AddConnectivity( & + gc, & + SRC_NAME = ['TApert ', 'QApert ', 'UUpert ', & + 'UWINDLMTILEpert', 'VWINDLMTILEpert', 'PCUpert ', & + 'PLSpert ', 'SNOpert ', 'DRPARpert ', & + 'DFPARpert ', 'DRNIRpert ', 'DFNIRpert ', & + 'DRUVRpert ', 'DFUVRpert ', 'LWDNSRFpert '], & + SRC_ID = LANDPERT(i), & + DST_NAME = ['TA ', 'QA ', 'UU ', 'UWINDLMTILE',& + 'VWINDLMTILE', 'PCU ', 'PLS ', 'SNO ',& + 'DRPAR ', 'DFPAR ', 'DRNIR ', 'DFNIR ',& + 'DRUVR ', 'DFUVR ', 'LWDNSRF '], & + DST_ID = LAND(i), & + rc = status & + ) + VERIFY_(status) + + ! -LAND-feeds-LANDPERT's-imports- + call MAPL_AddConnectivity( & + gc, & + SRC_NAME = ['TC ','CATDEF ','RZEXC ','SRFEXC ','WESNN1 ','WESNN2 ','WESNN3 ', & + 'GHTCNT1','GHTCNT2','GHTCNT3','GHTCNT4','GHTCNT5','GHTCNT6', & + 'HTSNNN1','HTSNNN2','HTSNNN3','SNDZN1 ','SNDZN2 ','SNDZN3 '], & + SRC_ID = LAND(i), & + DST_NAME = ['TCPert ','CATDEFPert ','RZEXCPert ','SRFEXCPert ','WESNN1Pert ', & + 'WESNN2Pert ','WESNN3Pert ','GHTCNT1Pert','GHTCNT2Pert', & + 'GHTCNT3Pert','GHTCNT4Pert','GHTCNT5Pert','GHTCNT6Pert', & + 'HTSNNN1Pert','HTSNNN2Pert','HTSNNN3Pert','SNDZN1Pert ', & + 'SNDZN2Pert ','SNDZN3Pert '], & + DST_ID = LANDPERT(i), & + rc = status & + ) VERIFY_(status) - ! -METFORCE-feeds-LAND's-imports- - call MAPL_AddConnectivity( & - gc, & - SRC_NAME = ['Psurf', 'RefH ', & - 'DUDP ', 'DUSV ', 'DUWT ', 'DUSD ', 'BCDP ', 'BCSV ', & - 'BCWT ', 'BCSD ', 'OCDP ', 'OCSV ', 'OCWT ', 'OCSD ', & - 'SUDP ', 'SUSV ', 'SUWT ', 'SUSD ', 'SSDP ', 'SSSV ' ], & - SRC_ID = METFORCE(k), & - DST_NAME = ['PS ', 'DZ ', & - 'DUDP', 'DUSV', 'DUWT', 'DUSD', 'BCDP', 'BCSV', & - 'BCWT', 'BCSD', 'OCDP', 'OCSV', 'OCWT', 'OCSD', & - 'SUDP', 'SUSV', 'SUWT', 'SUSD', 'SSDP', 'SSSV' ], & - DST_ID = LAND(i), & - rc = status & - ) - VERIFY_(status) - ! -LAND-feeds-LANDPERT's-imports- - call MAPL_AddConnectivity( & - gc, & - SRC_NAME = ['TC ','CATDEF ','RZEXC ','SRFEXC ','WESNN1 ','WESNN2 ','WESNN3 ', & - 'GHTCNT1','GHTCNT2','GHTCNT3','GHTCNT4','GHTCNT5','GHTCNT6', & - 'HTSNNN1','HTSNNN2','HTSNNN3','SNDZN1 ','SNDZN2 ','SNDZN3 '], & - SRC_ID = LAND(i), & - DST_NAME = ['TCPert ','CATDEFPert ','RZEXCPert ','SRFEXCPert ','WESNN1Pert ', & - 'WESNN2Pert ','WESNN3Pert ','GHTCNT1Pert','GHTCNT2Pert', & - 'GHTCNT3Pert','GHTCNT4Pert','GHTCNT5Pert','GHTCNT6Pert', & - 'HTSNNN1Pert','HTSNNN2Pert','HTSNNN3Pert','SNDZN1Pert ', & - 'SNDZN2Pert ','SNDZN3Pert '], & - DST_ID = LANDPERT(i), & - rc = status & - ) - VERIFY_(status) - enddo + enddo - if(land_assim .or. mwRTM) then - ! -LAND-feeds-LANDASSIM's-imports- - ! Catchment model parameters from first LAND ens member, assumes no parameter perturbations! - call MAPL_AddConnectivity( & - gc, & - SHORT_NAME = ['POROS ', 'COND ','PSIS ','BEE ','WPWET ','GNU ','VGWMAX', & - 'BF1 ', 'BF2 ','BF3 ','CDCR1 ','CDCR2 ','ARS1 ', & - 'ARS2 ', 'ARS3 ','ARA1 ','ARA2 ','ARA3 ','ARA4 ', & - 'ARW1 ', 'ARW2 ','ARW3 ','ARW4 ','TSA1 ','TSA2 ','TSB1 ', & - 'TSB2 ', 'ATAU ','BTAU ','ITY ','Z2CH ' ], & - SRC_ID = LAND(1), & ! Note (1) ! - DST_ID = LANDASSIM, & - rc = status & - ) - VERIFY_(status) - endif + if(land_assim .or. mwRTM) then + ! -LAND-feeds-LANDASSIM's-imports- + ! Catchment model parameters from first LAND ens member, assumes no parameter perturbations! + call MAPL_AddConnectivity( & + gc, & + SHORT_NAME = ['POROS ', 'COND ','PSIS ','BEE ','WPWET ','GNU ','VGWMAX', & + 'BF1 ', 'BF2 ','BF3 ','CDCR1 ','CDCR2 ','ARS1 ', & + 'ARS2 ', 'ARS3 ','ARA1 ','ARA2 ','ARA3 ','ARA4 ', & + 'ARW1 ', 'ARW2 ','ARW3 ','ARW4 ','TSA1 ','TSA2 ','TSB1 ', & + 'TSB2 ', 'ATAU ','BTAU ','ITY ','Z2CH ' ], & + SRC_ID = LAND(1), & ! Note (1) ! + DST_ID = LANDASSIM, & + rc = status & + ) + VERIFY_(status) + endif + end if !with_land call MAPL_TimerAdd(gc, name="Initialize", rc=status) VERIFY_(status) @@ -367,6 +385,8 @@ subroutine Initialize(gc, import, export, clock, rc) ! MAPL variables type(MAPL_LocStream) :: surf_locstream type(MAPL_LocStream) :: land_locstream + type(MAPL_LocStream) :: landice_locstream + type(MAPL_LocStream) :: force_locstream type(MAPL_MetaComp), pointer :: MAPL=>null() ! GC's MAPL obj type(MAPL_MetaComp), pointer :: CHILD_MAPL=>null() ! Child's MAPL obj @@ -382,9 +402,8 @@ subroutine Initialize(gc, import, export, clock, rc) character(len=ESMF_MAXSTR) :: LAND_PARAMS character(len=ESMF_MAXSTR) :: grid_type - integer :: total_nt,land_nt_local,i,j - real, pointer :: LandTileLats(:) - real, pointer :: LandTileLons(:) + integer :: i, j, k, n + integer :: total_nt, local_nt, tile_id integer, pointer :: local_id(:) real(ESMF_KIND_R8), pointer :: centerX(:,:) real(ESMF_KIND_R8), pointer :: centerY(:,:) @@ -400,7 +419,7 @@ subroutine Initialize(gc, import, export, clock, rc) type(tile_coord_type), dimension(:), pointer :: tile_coord_f => null() integer,dimension(:),pointer :: f2g - integer :: N_catf + integer :: N_f ! number of tiles in full domain integer :: LSM_CHOICE type(grid_def_type) :: tile_grid_g, pert_grid_g @@ -416,6 +435,7 @@ subroutine Initialize(gc, import, export, clock, rc) real :: DT, DT_Solar type(ESMF_Alarm) :: SolarAlarm type(ESMF_TimeInterval) :: Solar_DT + integer, allocatable :: mask(:) ! Begin... @@ -551,69 +571,89 @@ subroutine Initialize(gc, import, export, clock, rc) ) VERIFY_(status) - call MAPL_TimerOff(MAPL, "-LocStreamCreate") - ! Get children and their im/ex states from MAPL obj call MAPL_Get(MAPL, GCS=gcs, GCNAMES=gcnames, rc=status) VERIFY_(status) ! Create LAND's locstreams as subset of Surface locstream ! and add it to the children's MAPL objects + allocate(mask(0)) + if (with_land) then + call MAPL_LocStreamCreate( & + land_locstream, & + surf_locstream, & + name=gcnames(LAND(1)), & + mask=[MAPL_LAND], & + rc=status & + ) + VERIFY_(status) + mask =[mask,MAPL_LAND] + endif - call MAPL_TimerOn(MAPL, "-LocStreamCreate") - call MAPL_LocStreamCreate( & - land_locstream, & - surf_locstream, & - name=gcnames(LAND(1)), & - mask=[MAPL_LAND], & - rc=status & + if (with_landice) then + call MAPL_LocStreamCreate( & + landice_locstream, & + surf_locstream, & + name=gcnames(LANDICE(1)), & + mask=[MAPL_LANDICE], & + rc=status & ) + VERIFY_(status) + mask = [mask, MAPL_LANDICE] + endif + + call MAPL_LocStreamCreate( & + force_locstream, & + surf_locstream, & + name=gcnames(METFORCE(1)), & + mask=mask, & + rc=status & + ) VERIFY_(status) + call MAPL_TimerOff(MAPL, "-LocStreamCreate") ! Convert LAND's LocStream to LDAS' tile_coord and save it in the GridComp ! -get-tile-information-from-land's-locstream- - call MAPL_LocStreamGet( & - land_locstream, & - NT_LOCAL=land_nt_local, & - TILELATS=LandTileLats, & - TILELONS=LandTileLons, & - LOCAL_ID=local_id , & - rc=status & - ) + call MAPL_LocStreamGet( & + force_locstream, & + NT_LOCAL=local_nt, & + LOCAL_ID=local_id, & + rc=status & + ) VERIFY_(status) - + ! -get-component's-internal-state- call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) tcinternal => tcwrap%ptr ! -allocate-memory-for-tile-coord- - allocate(tcinternal%tile_coord(land_nt_local), stat=status) + allocate(tcinternal%tile_coord(local_nt), stat=status) VERIFY_(status) - allocate(tcinternal%l2f(land_nt_local)) + allocate(tcinternal%l2f(local_nt)) VERIFY_(status) - call MAPL_GetResource ( MAPL, out_path, Label="OUT_PATH:", DEFAULT="./", RC=STATUS) call MAPL_GetResource ( MAPL, exp_id, Label="EXP_ID:", DEFAULT="./", RC=STATUS) call init_MPI_types() - call MPI_Reduce(land_nt_local,total_nt,1,MPI_INT,MPI_SUM,0,mpicomm,mpierr); + call MPI_AllReduce(local_nt,total_nt,1,MPI_INT,MPI_SUM,mpicomm,mpierr); decomf = get_io_filename(trim(out_path), trim(exp_id), 'ldas_domdecomp', date_time=start_time, & dir_name='rc_out', file_ext='.txt') + if (IamRoot) then - call io_domain_files('r',trim(out_path), trim(exp_id),N_catf,f2g,tile_coord_f,tile_grid_g,tile_grid_f,RC) + call io_domain_files('r',trim(out_path), trim(exp_id),N_f,f2g,tile_coord_f,tile_grid_g,tile_grid_f,RC) ! WY notes: f2g == tile_coord_f%tile_id deallocate(f2g) - print*, "Number of tiles: ", N_catf - if(N_catf /= total_nt) then + print*, "Number of tiles: ", N_f + if(N_f /= total_nt) then print*, "total_nt = ", total_nt stop "tiles number not equal" endif open(10,file= trim(decomf), action='write') - write(10,*) N_catf + write(10,*) N_f close(10) call io_grid_def_type('w', logunit, tile_grid_f, 'tile_grid_f') @@ -635,12 +675,12 @@ subroutine Initialize(gc, import, export, clock, rc) ! arrive here when tile_grid_g is cube-sphere and pert_grid_g is lat/lon after call to get_pert_grid() above !1) get pert_i_indg, pert_j_indg for tiles in (full) domain relative to pert_grid_g - do i = 1, N_catf + do i = 1, N_f call get_ij_ind_from_latlon(pert_grid_g,tile_coord_f(i)%com_lat,tile_coord_f(i)%com_lon, & tile_coord_f(i)%pert_i_indg,tile_coord_f(i)%pert_j_indg) enddo !2) determine pert_grid_f - pert_grid_f = get_minExtent_grid(N_catf, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, & + pert_grid_f = get_minExtent_grid(N_f, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, & tile_coord_f%min_lon, tile_coord_f%min_lat, tile_coord_f%max_lon, tile_coord_f%max_lat, & pert_grid_g) @@ -654,34 +694,37 @@ subroutine Initialize(gc, import, export, clock, rc) endif endif - call MPI_BCAST(N_catf,1,MPI_INTEGER,0,mpicomm,mpierr) - if (.not. IamRoot) allocate(tile_coord_f(N_catf)) - - call MPI_BCAST(tile_coord_f,N_catf, MPI_tile_coord_type,0,mpicomm, mpierr) + if (.not. IamRoot) allocate(tile_coord_f(total_nt)) + call MPI_Barrier(mpicomm, mpierr) + call MPI_BCAST(tile_coord_f, total_nt, MPI_tile_coord_type,0,mpicomm, mpierr) call MPI_BCAST(pert_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) call MPI_BCAST(pert_grid_f, 1, MPI_grid_def_type, 0,mpicomm, mpierr) call MPI_BCAST(tile_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) block - integer, allocatable :: f2tile_id(:), tile_id2f(:) - integer :: max_id - allocate(f2tile_id(N_catf)) - f2tile_id = tile_coord_f%tile_id - - max_id = maxval(f2tile_id) - allocate(tile_id2f(max_id),source = 0) - do i = 1, N_catf - tile_id2f(f2tile_id(i)) = i + i = 1 + j = 0 + do k = 1, local_nt + do n = i, total_nt + tile_id = tile_coord_f(n)%tile_id + if (local_id(k) == tile_id) then + tcinternal%l2f(k) = n + i = n + 1 + j = j + 1 + exit + endif + enddo enddo - tcinternal%l2f = tile_id2f(local_id) + if (j /= local_nt) then + stop "tile distributtion is wrong. cannot find the right tile" + endif tcinternal%tile_coord = tile_coord_f(tcinternal%l2f) - deallocate(f2tile_id, tile_id2f) end block do i = 0, numprocs-1 if( i == myid) then open(10,file= trim(decomf), action='write',position='append') - do j = 1, land_nt_local + do j = 1, local_nt write(10,*) local_id(j), myid enddo close(10) @@ -691,7 +734,7 @@ subroutine Initialize(gc, import, export, clock, rc) allocate(tcinternal%tile_coord_f,source = tile_coord_f) - pert_grid_l = get_minExtent_grid(land_nt_local, & + pert_grid_l = get_minExtent_grid(local_nt, & tcinternal%tile_coord%pert_i_indg, tcinternal%tile_coord%pert_j_indg, & tcinternal%tile_coord%min_lon, tcinternal%tile_coord%min_lat, & tcinternal%tile_coord%max_lon, tcinternal%tile_coord%max_lat, & @@ -705,36 +748,51 @@ subroutine Initialize(gc, import, export, clock, rc) do i = 1, NUM_ENSEMBLE call MAPL_GetObjectFromGC(gcs(METFORCE(i)), CHILD_MAPL, rc=status) VERIFY_(status) ! CHILD = METFORCE - call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) + call MAPL_Set(CHILD_MAPL, LocStream=force_locstream, rc=status) VERIFY_(status) call ESMF_UserCompSetInternalState(gcs(METFORCE(i)), 'TILE_COORD', tcwrap, status) VERIFY_(status) + ! exit after i=1 if using deterministic forcing if (.not. ensemble_forcing) exit enddo - call MAPL_GetObjectFromGC(gcs(ENSAVG), CHILD_MAPL, rc=status) - VERIFY_(status) ! CHILD = ens_avg - call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) - VERIFY_(status) - - do i = 1,NUM_ENSEMBLE - call MAPL_GetObjectFromGC(gcs(LAND(i)), CHILD_MAPL, rc=status) - VERIFY_(status) + if ( with_land) then + call MAPL_GetObjectFromGC(gcs(ENSAVG), CHILD_MAPL, rc=status) + VERIFY_(status) ! CHILD = ens_avg call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) VERIFY_(status) - call MAPL_GetObjectFromGC(gcs(LANDPERT(i)), CHILD_MAPL, rc=status) - VERIFY_(status) ! CHILD = LANDPERT - call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) - VERIFY_(status) - ! Add LAND's tile_coord to children's GridComps - call ESMF_UserCompSetInternalState(gcs(LAND(i)), 'TILE_COORD', tcwrap, status) - VERIFY_(status) - call ESMF_UserCompSetInternalState(gcs(LANDPERT(i)), 'TILE_COORD', tcwrap, status) - VERIFY_(status) + endif + + do i = 1,NUM_ENSEMBLE + if (with_land) then + call MAPL_GetObjectFromGC(gcs(LAND(i)), CHILD_MAPL, rc=status) + VERIFY_(status) + call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) + VERIFY_(status) + + call MAPL_GetObjectFromGC(gcs(LANDPERT(i)), CHILD_MAPL, rc=status) + VERIFY_(status) ! CHILD = LANDPERT + call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) + VERIFY_(status) + + ! Add LAND's tile_coord to children's GridComps + call ESMF_UserCompSetInternalState(gcs(LAND(i)), 'TILE_COORD', tcwrap, status) + VERIFY_(status) + call ESMF_UserCompSetInternalState(gcs(LANDPERT(i)), 'TILE_COORD', tcwrap, status) + VERIFY_(status) + endif + + if (with_landice) then + call MAPL_GetObjectFromGC(gcs(LANDICE(i)), CHILD_MAPL, rc=status) + VERIFY_(status) + call MAPL_Set(CHILD_MAPL, LocStream=landice_locstream, rc=status) + VERIFY_(status) + endif + enddo - if (land_assim .or. mwRTM) then + if (land_assim .or. mwRTM ) then call MAPL_GetObjectFromGC(gcs(LANDASSIM), CHILD_MAPL, rc=status) VERIFY_(status) call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status) @@ -771,7 +829,6 @@ subroutine Initialize(gc, import, export, clock, rc) if ( IamRoot) call echo_catch_constants(logunit) if ( IamRoot) call StieglitzSnow_echo_constants(logunit) - ! Turn timer off call MAPL_TimerOff(MAPL, "Initialize") call MAPL_TimerOff(MAPL, "TOTAL") @@ -822,7 +879,7 @@ subroutine Run(gc, import, export, clock, rc) type(MAPL_MetaComp), pointer :: MAPL ! Misc variables - integer :: igc,i, ens_id, FIRST_ENS_ID, ens_id_width + integer :: igc, i, ens_id, FIRST_ENS_ID, ens_id_width, k logical :: IAmRoot integer :: LSM_CHOICE type (ESMF_Field) :: field @@ -868,99 +925,133 @@ subroutine Run(gc, import, export, clock, rc) call ESMF_TimePrint(ModelTimeCur, options='string', rc=status) VERIFY_(status) end if - + !phase2 initialization ( executed once) !adjust mean of perturbed forcing or Progn - do i = 1,NUM_ENSEMBLE - igc = LANDPERT(i) - call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) - VERIFY_(status) - call MAPL_TimerOff(MAPL, gcnames(igc)) - enddo - - ! Run children GridComps (in order) - ! Generate raw perturbed force and progn - do i = 1,NUM_ENSEMBLE - igc = LANDPERT(i) - call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status) - VERIFY_(status) - call MAPL_TimerOff(MAPL, gcnames(igc)) - enddo + if (with_land) then + do i = 1,NUM_ENSEMBLE + igc = LANDPERT(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) + VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) + enddo + ! Run children GridComps (in order) + ! Generate raw perturbed force and progn + do i = 1,NUM_ENSEMBLE + igc = LANDPERT(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status) + VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) + enddo + endif ! with_land do i = 1, NUM_ENSEMBLE igc = METFORCE(i) call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, userRC=status) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) VERIFY_(status) call MAPL_TimerOff(MAPL, gcnames(igc)) ! exit after i=1 if using deterministic forcing if (.not. ensemble_forcing) exit enddo - - do i = 1,NUM_ENSEMBLE - - !ApplyForcePert - igc = LANDPERT(i) + ! distribute surface met forcing (export forcing to imports of land, landpert, and landice) + do i = 1, NUM_ENSEMBLE + k = 1 + if (ensemble_forcing) k = i + igc = METFORCE(k) call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=3, userRC=status) - VERIFY_(status) - call MAPL_TimerOff(MAPL, gcnames(igc)) - ! Use landpert's output as the input to calculate the ensemble average forcing - ! W.J note: So far it is only for the Catchment model. - ! To make CatchmentCN work with assim, the export from landgrid and catchmentCN grid need to be modified. - if ( LSM_CHOICE == 1 ) then - call ESMF_GridCompRun(gcs(ENSAVG), importState=gex(igc), exportState=gex(ENSAVG), clock=clock,phase=1, userRC=status) + if (with_land) then + call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LAND(i)), clock=clock, phase=2, userRC=status) + VERIFY_(status) + call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LANDPERT(i)), clock=clock, phase=3, userRC=status) VERIFY_(status) endif - ! Run the land model - igc = LAND(i) - call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) - VERIFY_(status) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status) - VERIFY_(status) + if (with_landice) then + call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LANDICE(i)), clock=clock, phase=4, userRC=status) + VERIFY_(status) + endif call MAPL_TimerOff(MAPL, gcnames(igc)) + enddo - ! ApplyPrognPert - moved: now before calculating ensemble average that is picked up by land analysis and HISTORY; reichle 28 May 2020 - igc = LANDPERT(i) - call MAPL_TimerOn(MAPL, gcnames(igc)) - call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=4, userRC=status) - VERIFY_(status) - call MAPL_TimerOff(MAPL, gcnames(igc)) + do i = 1,NUM_ENSEMBLE - ! Use LAND's output as the input to calculate the ensemble average - igc = LAND(i) - if (LSM_CHOICE == 1) then - ! collect cat_param - ens_id = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID - call get_ensid_string(ensid_string, ens_id, ens_id_width, NUM_ENSEMBLE) + if (with_land) then + !ApplyForcePert + igc = LANDPERT(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=3, userRC=status) + VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) - member_name = 'CATCH'//trim(ensid_string)//"_Exports" + ! Use landpert's output as the input to calculate the ensemble average forcing + ! W.J note: So far it is only for the Catchment model. + ! To make CatchmentCN work with assim, the export from landgrid and catchmentCN grid need to be modified. + if ( LSM_CHOICE == 1 ) then + call ESMF_GridCompRun(gcs(ENSAVG), importState=gex(igc), exportState=gex(ENSAVG), clock=clock,phase=1, userRC=status) + VERIFY_(status) + endif - call ESMF_StateGet(gex(igc), trim(member_name), member_export, _RC) - call ESMF_StateGet(gex(igc), "Z2CH", field, _RC) - call ESMF_StateAddReplace(member_export, [field],_RC) - call ESMF_StateGet(gex(igc), "LAI", field, _RC) - call ESMF_StateAddReplace(member_export, [field],_RC) + ! Run the land model + igc = LAND(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) + VERIFY_(status) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status) + VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) + endif ! with_land - call ESMF_GridCompRun(gcs(ENSAVG), importState=member_export, exportState=gex(ENSAVG), clock=clock,phase=3, userRC=status) + if (with_landice) then + igc = LANDICE(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status) VERIFY_(status) - call ESMF_GridCompRun(gcs(ENSAVG), importState=member_export, exportState=gex(ENSAVG), clock=clock,phase=2, userRC=status) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status) VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) + endif ! with_land_ice + + if (with_land) then + ! ApplyPrognPert - moved: now before calculating ensemble average that is picked up by land analysis and HISTORY; reichle 28 May 2020 + igc = LANDPERT(i) + call MAPL_TimerOn(MAPL, gcnames(igc)) + call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=4, userRC=status) + VERIFY_(status) + call MAPL_TimerOff(MAPL, gcnames(igc)) + + ! Use LAND's output as the input to calculate the ensemble average + igc = LAND(i) + if (LSM_CHOICE == 1) then + ! collect cat_param + ens_id = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID + call get_ensid_string(ensid_string, ens_id, ens_id_width, NUM_ENSEMBLE) + + member_name = 'CATCH'//trim(ensid_string)//"_Exports" + + call ESMF_StateGet(gex(igc), trim(member_name), member_export, _RC) + call ESMF_StateGet(gex(igc), "Z2CH", field, _RC) + call ESMF_StateAddReplace(member_export, [field],_RC) + call ESMF_StateGet(gex(igc), "LAI", field, _RC) + call ESMF_StateAddReplace(member_export, [field],_RC) - if( mwRTM ) then - ! Calculate ensemble-average L-band Tb using LAND's output (add up and normalize after last member has been added) - call ESMF_GridCompRun(gcs(LANDASSIM), importState=member_export, exportState=gex(LANDASSIM), clock=clock,phase=3, userRC=status) + call ESMF_GridCompRun(gcs(ENSAVG), importState=member_export, exportState=gex(ENSAVG), clock=clock,phase=3, userRC=status) + VERIFY_(status) + call ESMF_GridCompRun(gcs(ENSAVG), importState=member_export, exportState=gex(ENSAVG), clock=clock,phase=2, userRC=status) VERIFY_(status) - endif - endif + if( mwRTM ) then + ! Calculate ensemble-average L-band Tb using LAND's output (add up and normalize after last member has been added) + call ESMF_GridCompRun(gcs(LANDASSIM), importState=member_export, exportState=gex(LANDASSIM), clock=clock,phase=3, userRC=status) + VERIFY_(status) + endif + endif + endif ! with_land enddo if ( mwRTM .and. LSM_CHOICE == 1 ) then diff --git a/GEOSldas_App/GEOSldas_HIST.rc b/GEOSldas_App/GEOSldas_HIST.rc index ce692ffb..d9d09dd8 100644 --- a/GEOSldas_App/GEOSldas_HIST.rc +++ b/GEOSldas_App/GEOSldas_HIST.rc @@ -25,6 +25,8 @@ COLLECTIONS: # 'inst3_2d_lndfcstana_Nx' # 'const_1d_lnd_Nt' # 'const_2d_lnd_Nx' +# 'tavg24_2d_glc_Nx' +# 'tavg24_1d_glc_Nt' :: #CUBE GRID_LABELS: PC720x361-DC @@ -543,5 +545,72 @@ COLLECTIONS: 'TPSURF_ANA_ENSSTD' , 'LANDASSIM' , 'TSURF_ANA_ENSSTD' , 'TP1_ANA_ENSSTD' , 'LANDASSIM' , 'TSOIL1_ANA_ENSSTD' :: + + tavg24_2d_glc_Nx.descr: '2d,Daily,Time-Averaged,Single-Level,Land Ice Diagnostics', + tavg24_2d_glc_Nx.nbits: 12, + tavg24_2d_glc_Nx.template: '%y4%m2%d2_%h2%n2z.nc4' , + tavg24_2d_glc_Nx.mode: 'time-averaged' , + tavg24_2d_glc_Nx.frequency: 240000 , + tavg24_2d_glc_Nx.ref_time: 000000 , + tavg24_2d_glc_Nx.format: 'CFIO' , + tavg24_2d_glc_Nx.regrid_exch: '../input/tile.data' , + tavg24_2d_glc_Nx.regrid_name: 'GRIDNAME' , + tavg24_2d_glc_Nx.grid_label: PC720x361-DC , #comment this line out for cube face output + tavg24_2d_glc_Nx.deflate: 1, + tavg24_2d_glc_Nx.fields: 'ACCUM' , 'LANDICE' , + 'ALBVR' , 'LANDICE' , 'ALBVR_GL' , + 'ALBVF' , 'LANDICE' , 'ALBVF_GL' , + 'ALBNR' , 'LANDICE' , 'ALBNR_GL' , + 'ALBNF' , 'LANDICE' , 'ALBNF_GL' , + 'ASNOW_GL' , 'LANDICE' , + 'DELTS' , 'LANDICE' , + 'DNICFLX' , 'LANDICE' , + 'EVAPOUT' , 'LANDICE' , + 'QH' , 'LANDICE' , + 'GHTSKIN' , 'LANDICE' , 'GHTSKIN_GL' , + 'HLATN' , 'LANDICE' , + 'HLWUP' , 'LANDICE' , 'HLWUP_GL' , + 'IMELT' , 'LANDICE' , + 'LWNDSRF' , 'LANDICE' , + 'MELTWTR' , 'LANDICE' , + 'MELTWTRCONT' , 'LANDICE' , + 'RUNOFF' , 'LANDICE' , 'RUNOFF_GL' , + 'SHOUT' , 'LANDICE' , + 'SMELT' , 'LANDICE' , + 'SNICEALB' , 'LANDICE' , + 'SNOMAS_GL' , 'LANDICE' , + 'SNOWALB' , 'LANDICE' , + 'SNOWDP_GL' , 'LANDICE' , + 'SWNDSRF' , 'LANDICE' , + 'TST' , 'LANDICE' , + 'WESNBOT' , 'LANDICE' , + 'WESNEXT' , 'LANDICE' , + 'WESNPERC' , 'LANDICE' , + 'WESNPREC' , 'LANDICE' , + :: + tavg24_1d_glc_Nt.descr: 'Tile-space,Daily,Time-Averaged,Single-level,Land Ice Diagnostics', + tavg24_1d_glc_Nt.nbits: 12, + tavg24_1d_glc_Nt.template: '%y4%m2%d2_%h2%n2z.bin' , + tavg24_1d_glc_Nt.mode: 'time-averaged' , + tavg24_1d_glc_Nt.frequency: 240000 , + tavg24_1d_glc_Nt.ref_time: 000000 , + tavg24_1d_glc_Nt.fields: 'ASNOW_GL' , 'LANDICE' , + 'DELTS' , 'LANDICE' , + 'EVAPOUT' , 'LANDICE' , + 'GHTSKIN' , 'LANDICE' , 'GHTSKIN_GL' , + 'HLATN' , 'LANDICE' , + 'HLWUP' , 'LANDICE' , 'HLWUP_GL' , + 'LWNDSRF' , 'LANDICE' , + 'MELTWTR' , 'LANDICE' , + 'MELTWTRCONT' , 'LANDICE' , + 'RUNOFF' , 'LANDICE' , 'RUNOFF_GL' , + 'SHOUT' , 'LANDICE' , + 'SNOMAS_GL' , 'LANDICE' , + 'SNOWDP_GL' , 'LANDICE' , + 'SWNDSRF' , 'LANDICE' , + 'TST' , 'LANDICE' , + 'WESNBOT' , 'LANDICE' , + 'WESNEXT' , 'LANDICE' , + :: # ========================== EOF ============================================================== diff --git a/GEOSldas_App/GEOSldas_LDAS.rc b/GEOSldas_App/GEOSldas_LDAS.rc index df92fdb0..c77ae6bd 100644 --- a/GEOSldas_App/GEOSldas_LDAS.rc +++ b/GEOSldas_App/GEOSldas_LDAS.rc @@ -29,7 +29,7 @@ CATCHMENT_OFFLINE: 1 CATCHMENT_SPINUP: 0 -# ---- Choice of land surface model +# ---- Choice of land surface model (for LAND tiles) # # 1 : Catchment model (default) # 2 : CatchmentCN-CLM4.0 @@ -37,6 +37,28 @@ CATCHMENT_SPINUP: 0 LSM_CHOICE: 1 +# ---- Choice of tile type(s) +# +# List of tile types to be included in simulation. +# Use blank space as separator if there is more than one type. +# +# land : 100 (non-glaciated land) +# landice : 20 ( glaciated land) +# lake : 19 [not yet implemented] +# +# For example, include land and landice tiles as follows: +# TILE_TYPES: 100 20 +# +TILE_TYPES: 100 + +# ---- Format of tile file (from bcs directory) +# +# DEFAULT : Use nc4 tile file if it exists, txt tile file otherwise +# TXT : Use txt tile file (e.g., for 0-diff testing) +# +TILE_FILE_FORMAT: DEFAULT + + # ---- Domain definition # # The domain is determined by specifying a lat/lon rectangle in conjunction diff --git a/GEOSldas_App/ldas_setup b/GEOSldas_App/ldas_setup index aeb152d4..e1dbe6be 100755 --- a/GEOSldas_App/ldas_setup +++ b/GEOSldas_App/ldas_setup @@ -13,13 +13,15 @@ import resource import subprocess as sp import shlex import tempfile +import netCDF4 from dateutil import rrule from datetime import datetime from datetime import timedelta from collections import OrderedDict from dateutil.relativedelta import relativedelta from remap_utils import * -from remap_catchANDcn import * +from remap_catchANDcn import * +from remap_lake_landice_saltwater import * from lenkf_j_template import * """ This script is intended to be run from any installed directory with GEOSldas.x and ldas_setup @@ -49,7 +51,7 @@ class LDASsetup: 'VISDF_FILE','CATCH_DEF_FILE','NDVI_FILE', 'NML_INPUT_PATH','HISTRC_FILE','RST_FROM_GLOBAL','JOB_SGMT','NUM_SGMT','POSTPROC_HIST', 'MINLON','MAXLON','MINLAT','MAXLAT','EXCLUDE_FILE','INCLUDE_FILE','MWRTM_PATH','GRIDNAME', - 'ADAS_EXPDIR', 'BCS_RESOLUTION' ] + 'ADAS_EXPDIR', 'BCS_RESOLUTION', 'TILE_FILE_FORMAT' ] self.GEOS_SITE = "@GEOS_SITE@" @@ -106,7 +108,7 @@ class LDASsetup: self.out_path = None self.inpdir = None self.exefyl = None - self.islocal = False + self.isZoomIn = False self.catch = '' self.has_mwrtm = False self.has_vegopacity = False @@ -123,6 +125,9 @@ class LDASsetup: self.bcs_land = '' self.bcs_geom = '' self.bcs_landshared = '' + self.tile_types = '' + self.with_land = False + self.with_landice = False # ------ # Read exe input file which is required to set up the dir @@ -158,6 +163,12 @@ class LDASsetup: print ('\nInputs from execfile:\n') _printdict(self.rqdExeInp) + self.tile_types = self.rqdExeInp.get('TILE_TYPES',"100").split() + if "100" in self.tile_types : + self.with_land = True + if "20" in self.tile_types : + self.with_landice = True + # nens is an integer and =1 for model run self.nens = int(self.rqdExeInp['NUM_LDAS_ENSEMBLE']) # fail if Nens's val is not int assert self.nens>0, 'NUM_LDAS_ENSEMBLE [%d] <= 0' % self.nens @@ -275,46 +286,53 @@ class LDASsetup: self.rqdExeInp['RESTART_PATH'] = self.rqdExeInp['RESTART_PATH']+'/' # make sure catchment and vegdyn restart files ( at least one for each) exist - if 'CATCH_DEF_FILE' not in self.rqdExeInp: - self.rqdExeInp['CATCH_DEF_FILE']= self.bcs_land + 'clsm/catchment.def' - assert os.path.isfile(self.rqdExeInp['CATCH_DEF_FILE']),"[%s] file does not exist " % self.rqdExeInp['CATCH_DEF_FILE'] + if 'CATCH_DEF_FILE' not in self.rqdExeInp : + self.rqdExeInp['CATCH_DEF_FILE']= self.bcs_land + 'clsm/catchment.def' + if (self.with_land) : + assert os.path.isfile(self.rqdExeInp['CATCH_DEF_FILE']),"[%s] file does not exist " % self.rqdExeInp['CATCH_DEF_FILE'] self.rqdExeInp['RST_FROM_GLOBAL'] = 1 - if self.rqdExeInp['RESTART'].isdigit() : - if int(self.rqdExeInp['RESTART']) == 1 : - _numg = int(linecache.getline(self.rqdExeInp['CATCH_DEF_FILE'], 1).strip()) - _numd = _numg - ldas_domain = self.rqdExeInp['RESTART_PATH']+ \ - self.rqdExeInp['RESTART_ID'] + \ - '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'+self.rqdExeInp['RESTART_ID']+'.ldas_domain.txt' - if os.path.isfile(ldas_domain) : - _numd = int(linecache.getline(ldas_domain, 1).strip()) - - if _numg != _numd : - self.rqdExeInp['RST_FROM_GLOBAL'] = 0 + # skip checking. It is users' reponsibility to make it right! + #if self.rqdExeInp['RESTART'].isdigit() : + # if int(self.rqdExeInp['RESTART']) == 1 : + # _numg = int(linecache.getline(self.rqdExeInp['CATCH_DEF_FILE'], 1).strip()) + # _numd = _numg + # ldas_domain = self.rqdExeInp['RESTART_PATH']+ \ + # self.rqdExeInp['RESTART_ID'] + \ + # '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'+self.rqdExeInp['RESTART_ID']+'.ldas_domain.txt' + # if os.path.isfile(ldas_domain) : + # _numd = int(linecache.getline(ldas_domain, 1).strip()) + # + # if _numg != _numd : + # self.rqdExeInp['RST_FROM_GLOBAL'] = 0 self.rqdExeInp['LNFM_FILE'] = '' + tile_file_format = self.rqdExeInp.get('TILE_FILE_FORMAT', 'DEFAULT') if int(self.rqdExeInp['RST_FROM_GLOBAL']) == 1 : - self.rqdExeInp['TILING_FILE'] =glob.glob(self.bcs_geom + '*.til')[0] - self.rqdExeInp['GRN_FILE']= glob.glob(self.bcs_land + 'green_clim_*.data')[0] - self.rqdExeInp['LAI_FILE']= glob.glob(self.bcs_land + 'lai_clim_*.data')[0] - tmp_ = glob.glob(self.bcs_land + 'lnfm_clim_*.data') - if (len(tmp_) ==1) : - self.rqdExeInp['LNFM_FILE'] = tmp_[0] - self.rqdExeInp['NDVI_FILE'] = glob.glob(self.bcs_land + 'ndvi_clim_*.data')[0] - self.rqdExeInp['NIRDF_FILE']= glob.glob(self.bcs_land + 'nirdf_*.dat')[0] - self.rqdExeInp['VISDF_FILE']= glob.glob(self.bcs_land + 'visdf_*.dat')[0] + txt_tile = glob.glob(self.bcs_geom + '*.til') + nc4_tile = glob.glob(self.bcs_geom + '*.nc4') + if tile_file_format.upper() == 'TXT' : self.rqdExeInp['TILING_FILE'] = txt_tile[0] + if tile_file_format.upper() == 'DEFAULT' : self.rqdExeInp['TILING_FILE'] = (txt_tile+nc4_tile)[-1] + + self.rqdExeInp['GRN_FILE']= glob.glob(self.bcs_land + 'green_clim_*.data')[0] + self.rqdExeInp['LAI_FILE']= glob.glob(self.bcs_land + 'lai_clim_*.data')[0] + tmp_ = glob.glob(self.bcs_land + 'lnfm_clim_*.data') + if (len(tmp_) ==1) : + self.rqdExeInp['LNFM_FILE'] = tmp_[0] + self.rqdExeInp['NDVI_FILE'] = glob.glob(self.bcs_land + 'ndvi_clim_*.data')[0] + self.rqdExeInp['NIRDF_FILE']= glob.glob(self.bcs_land + 'nirdf_*.dat')[0] + self.rqdExeInp['VISDF_FILE']= glob.glob(self.bcs_land + 'visdf_*.dat')[0] else : - inpdir=self.rqdExeInp['RESTART_PATH']+self.rqdExeInp['RESTART_ID']+'/input/' - self.rqdExeInp['TILING_FILE'] =os.path.realpath(glob.glob(inpdir+'*tile.data')[0]) - self.rqdExeInp['GRN_FILE']= os.path.realpath(glob.glob(inpdir+'green*data')[0]) - self.rqdExeInp['LAI_FILE']= os.path.realpath(glob.glob(inpdir+'lai*data')[0]) - tmp_ = glob.glob(self.bcs_land + 'lnfm_clim_*.data') - if (len(tmp_) == 1) : - self.rqdExeInp['LNFM_FILE'] = tmp_[0] - self.rqdExeInp['NDVI_FILE']= os.path.realpath(glob.glob(inpdir+'ndvi*data')[0]) - self.rqdExeInp['NIRDF_FILE']= os.path.realpath(glob.glob(inpdir+'nirdf*data')[0]) - self.rqdExeInp['VISDF_FILE']= os.path.realpath(glob.glob(inpdir+'visdf*data')[0]) + inpdir=self.rqdExeInp['RESTART_PATH']+self.rqdExeInp['RESTART_ID']+'/input/' + self.rqdExeInp['TILING_FILE'] =os.path.realpath(glob.glob(inpdir+'*tile.data')[0]) + self.rqdExeInp['GRN_FILE']= os.path.realpath(glob.glob(inpdir+'green*data')[0]) + self.rqdExeInp['LAI_FILE']= os.path.realpath(glob.glob(inpdir+'lai*data')[0]) + tmp_ = glob.glob(self.bcs_land + 'lnfm_clim_*.data') + if (len(tmp_) == 1) : + self.rqdExeInp['LNFM_FILE'] = tmp_[0] + self.rqdExeInp['NDVI_FILE']= os.path.realpath(glob.glob(inpdir+'ndvi*data')[0]) + self.rqdExeInp['NIRDF_FILE']= os.path.realpath(glob.glob(inpdir+'nirdf*data')[0]) + self.rqdExeInp['VISDF_FILE']= os.path.realpath(glob.glob(inpdir+'visdf*data')[0]) if self.rqdExeInp['RESTART'].isdigit() : if int(self.rqdExeInp['RESTART']) == 2 : @@ -329,11 +347,18 @@ class LDASsetup: in_tilefiles_ = glob.glob(inpdir+'MAPL_*.til') if len(in_tilefiles_) == 0 : in_tilefiles_ = glob.glob(inpdir+'/*.til') + if len(in_tilefiles_) == 0 : + in_tilefiles_ = glob.glob(inpdir+'/*.nc4') + self.in_tilefile =os.path.realpath(in_tilefiles_[0]) - if os.path.isfile(ldas_domain) : + if os.path.isfile(ldas_domain): _numd = int(linecache.getline(ldas_domain, 1).strip()) - self.rqdExeInp['TILING_FILE'] =glob.glob(self.bcs_geom + '*.til')[0] + txt_tile = glob.glob(self.bcs_geom + '*.til') + nc4_tile = glob.glob(self.bcs_geom + '*.nc4') + if tile_file_format.upper() == 'TXT' : self.rqdExeInp['TILING_FILE'] = txt_tile[0] + if tile_file_format.upper() == 'DEFAULT' : self.rqdExeInp['TILING_FILE'] = (txt_tile+nc4_tile)[-1] + self.rqdExeInp['GRN_FILE']= glob.glob(self.bcs_land + 'green_clim_*.data')[0] self.rqdExeInp['LAI_FILE']= glob.glob(self.bcs_land + 'lai_clim_*.data')[0] tmp_ = glob.glob(self.bcs_land + 'lnfm_clim_*.data') @@ -345,9 +370,17 @@ class LDASsetup: self.rqdExeInp['VISDF_FILE']= glob.glob(self.bcs_land + 'visdf_*.dat')[0] if 'GRIDNAME' not in self.rqdExeInp : - tmptile =self.rqdExeInp['TILING_FILE'] - self.rqdExeInp['GRIDNAME'] = linecache.getline(tmptile, 3).strip() - + tmptile = os.path.realpath(self.rqdExeInp['TILING_FILE']) + extension = os.path.splitext(tmptile)[1] + if extension == '.domain': + extension = os.path.splitext(tmptile)[0] + + if extension == '.til': + self.rqdExeInp['GRIDNAME'] = linecache.getline(tmptile, 3).strip() + else: + nc_file = netCDF4.Dataset(tmptile,'r') + self.rqdExeInp['GRIDNAME'] = nc_file.getncattr('Grid_Name') + if 'LSM_CHOICE' not in self.rqdExeInp: self.rqdExeInp['LSM_CHOICE'] = 1 @@ -356,7 +389,8 @@ class LDASsetup: if int(self.rqdExeInp['LSM_CHOICE']) == 2 : self.catch = 'catchcnclm40' - assert int(self.rqdExeInp['LSM_CHOICE']) <= 2, "\nLSM_CHOICE=3 (Catchment-CN4.5) is no longer supported. Please set LSM_CHOICE to 1 (Catchment) or 2 (Catchment-CN4.0)" + if self.with_land: + assert int(self.rqdExeInp['LSM_CHOICE']) <= 2, "\nLSM_CHOICE=3 (Catchment-CN4.5) is no longer supported. Please set LSM_CHOICE to 1 (Catchment) or 2 (Catchment-CN4.0)" if 'POSTPROC_HIST' not in self.rqdExeInp: self.rqdExeInp['POSTPROC_HIST'] = 0 @@ -394,7 +428,7 @@ class LDASsetup: self.domain_def.close() # make sure bcs files exist - if self.rqdExeInp['RESTART'].isdigit() : + if self.rqdExeInp['RESTART'].isdigit() and self.with_land : if int(self.rqdExeInp['RESTART']) >= 1 : y4m2='Y%4d/M%02d' % (self.begDates[0].year, self.begDates[0].month) y4m2d2_h2m2='%4d%02d%02d_%02d%02d' % (self.begDates[0].year, self.begDates[0].month, @@ -442,9 +476,9 @@ class LDASsetup: self.in_tilefile = None # DEAL WITH mwRTM input from exec - self.assim = True if self.rqdExeInp.get('LAND_ASSIM', 'NO').upper() == 'YES' else False + self.assim = True if self.rqdExeInp.get('LAND_ASSIM', 'NO').upper() == 'YES' and self.with_land else False # verify mwrtm file - if 'MWRTM_PATH' in self.rqdExeInp : + if 'MWRTM_PATH' in self.rqdExeInp and self.with_land : self.rqdExeInp['MWRTM_PATH'] = self.rqdExeInp['MWRTM_PATH']+'/'+ self.rqdExeInp['BCS_RESOLUTION']+'/' mwrtm_param_file_ = self.rqdExeInp['MWRTM_PATH']+'mwRTM_param.nc4' vegopacity_file_ = self.rqdExeInp['MWRTM_PATH']+'vegopacity.bin' @@ -740,75 +774,72 @@ class LDASsetup: if 'WEMIN_OUT' in self.rqdExeInp : wemin_out = self.rqdExeInp['WEMIN_OUT'] + tmp_f2g_file = tempfile.NamedTemporaryFile(delete=False) + cmd = self.bindir +'/preprocess_ldas.x c_f2g ' + tile + ' ' + self.domain_def.name + ' '+ self.out_path + ' ' + catchment_def + ' ' + exp_id + ' ' + _y4m2d2h2m2 + ' '+ dzsf + ' ' + tmp_f2g_file.name + ' ' + '_'.join(self.tile_types) - tmp_f2g_file = tempfile.NamedTemporaryFile(delete=False) - cmd = self.bindir +'/preprocess_ldas.x c_f2g ' + tile + ' ' + self.domain_def.name + ' '+ self.out_path + ' ' + catchment_def + ' ' + exp_id + ' ' + _y4m2d2h2m2 + ' '+ dzsf + ' ' + tmp_f2g_file.name - - print ('Creating f2g file: '+ tmp_f2g_file.name +'....\n') + print ('Creating f2g file if necessary: '+ tmp_f2g_file.name +'....\n') print ("cmd: " + cmd) - sp.call(shlex.split(cmd)) + sp.call(shlex.split(cmd)) # check if it is local or global - with open(tmp_f2g_file.name) as f2gfile : - head=[next(f2gfile) for x in range(2)] - if(head[0].strip() != head[1].strip()) : - self.islocal= True + if os.path.getsize(tmp_f2g_file.name) !=0 : + self.isZoomIn= True #os.remove(self.domain_def.name) # update tile domain - if self.islocal: - newlocalTile = tile+'.domain' - print ("\nCreating local tile file :"+ newlocalTile) - print ("\n by excluding land type MAPL_Land_ExcludeFromDomain=1100...\n") - cmd = self.bindir +'/preprocess_ldas.x c_localtile ' + tile + ' ' + newlocalTile + ' '+ tmp_f2g_file.name + if self.isZoomIn: + newZoominTile = tile+'.domain' + print ("\nCreating local tile file :"+ newZoominTile) + print ("\nAdding 1000 to type of tiles to be excluded from domain...\n") + cmd = self.bindir +'/preprocess_ldas.x zoomin_tile ' + tile + ' ' + newZoominTile + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) short_tile=short_tile +'.domain' - tile = newlocalTile + tile = newZoominTile myTile=self.inpdir+'/tile.data' os.symlink(tile,myTile) - - bcs=[self.rqdExeInp['GRN_FILE'], - self.rqdExeInp['LAI_FILE'], - self.rqdExeInp['NDVI_FILE'], - self.rqdExeInp['NIRDF_FILE'], - self.rqdExeInp['VISDF_FILE'] ] - if (self.rqdExeInp['LNFM_FILE'] != ''): - bcs += [self.rqdExeInp['LNFM_FILE']] - if (self.has_vegopacity): - bcs += [self.rqdExeInp['VEGOPACITY_FILE']] - bcstmp=[] - for bcf in bcs : - shutil.copy(bcf, self.bcsdir+'/') - bcstmp=bcstmp+[self.bcsdir+'/'+os.path.basename(bcf)] - bcs=bcstmp - - if self.islocal: - print ("Creating the boundary files for the simulation domain...\n") - bcs_tmp=[] + if self.with_land: + bcs=[self.rqdExeInp['GRN_FILE'], + self.rqdExeInp['LAI_FILE'], + self.rqdExeInp['NDVI_FILE'], + self.rqdExeInp['NIRDF_FILE'], + self.rqdExeInp['VISDF_FILE'] ] + if (self.rqdExeInp['LNFM_FILE'] != ''): + bcs += [self.rqdExeInp['LNFM_FILE']] + if (self.has_vegopacity): + bcs += [self.rqdExeInp['VEGOPACITY_FILE']] + bcstmp=[] for bcf in bcs : - cmd = self.bindir +'/preprocess_ldas.x c_localbc ' + bcf + ' '+ bcf+'.domain' + ' '+ tmp_f2g_file.name - print ("cmd: " + cmd) - sp.call(shlex.split(cmd)) - bcs_tmp=bcs_tmp+[bcf+'.domain'] - bcs=bcs_tmp - - - # link BC - print ("linking bcs...") - bcnames=['green','lai','ndvi','nirdf','visdf'] - if (self.rqdExeInp['LNFM_FILE'] != ''): - bcnames += ['lnfm'] - if (self.has_vegopacity): - bcnames += ['vegopacity'] - for bcln,bc in zip(bcnames,bcs) : - myBC=self.inpdir+'/'+bcln+'.data' - os.symlink(bc,myBC) - - if ("catchcn" in self.catch): - os.symlink(self.bcs_landshared + 'CO2_MonthlyMean_DiurnalCycle.nc4', \ - self.inpdir+'/CO2_MonthlyMean_DiurnalCycle.nc4') + shutil.copy(bcf, self.bcsdir+'/') + bcstmp=bcstmp+[self.bcsdir+'/'+os.path.basename(bcf)] + bcs=bcstmp + + if self.isZoomIn: + print ("Creating the boundary files for the simulation domain...\n") + bcs_tmp=[] + for bcf in bcs : + cmd = self.bindir +'/preprocess_ldas.x zoomin_bc ' + bcf + ' '+ bcf+'.domain' + ' '+ tmp_f2g_file.name + print ("cmd: " + cmd) + sp.call(shlex.split(cmd)) + bcs_tmp=bcs_tmp+[bcf+'.domain'] + bcs=bcs_tmp + + + # link BC + print ("linking bcs...") + bcnames=['green','lai','ndvi','nirdf','visdf'] + if (self.rqdExeInp['LNFM_FILE'] != ''): + bcnames += ['lnfm'] + if (self.has_vegopacity): + bcnames += ['vegopacity'] + for bcln,bc in zip(bcnames,bcs) : + myBC=self.inpdir+'/'+bcln+'.data' + os.symlink(bc,myBC) + + if ("catchcn" in self.catch): + os.symlink(self.bcs_landshared + 'CO2_MonthlyMean_DiurnalCycle.nc4', \ + self.inpdir+'/CO2_MonthlyMean_DiurnalCycle.nc4') # create and link restart print ("Creating and linking restart...") @@ -827,11 +858,7 @@ class LDASsetup: '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/' # pass into remap_config_ldas - sponsorid = self.rqdRmInp['account'] exp_id = self.rqdExeInp['EXP_ID'] - exp_dir = self.exphome - out_bcdir = self.rqdExeInp['BCS_PATH'] - out_tilefile = self.rqdExeInp['TILING_FILE'] RESTART_str = str(self.rqdExeInp['RESTART']) YYYYMMDD = '%4d%02d%02d' % (_start.year, _start.month,_start.day) YYYYMMDDHH= '%4d%02d%02d%02d' % (_start.year, _start.month,_start.day, _start.hour) @@ -851,7 +878,7 @@ class LDASsetup: self.has_landassim_seed = True mk_outdir = self.exphome+'/'+exp_id+'/mk_restarts/' - if (RESTART_str != '1'): + if (RESTART_str != '1' and (self.with_land or self.with_landice)): bcs_path = self.rqdExeInp['BCS_PATH'] while bcs_path[-1] == '/' : bcs_path = bcs_path[0:-1] bc_base = os.path.dirname(bcs_path) @@ -866,7 +893,8 @@ class LDASsetup: config['input']['surface']['catch_tilefile'] = self.in_tilefile config['input']['shared']['expid'] = self.rqdExeInp['RESTART_ID'] config['input']['shared']['yyyymmddhh'] = YYYYMMDDHH - config['input']['shared']['rst_dir'] = os.path.dirname(self.in_rstfile)+'/' + if RESTART_str != 'M': + config['input']['shared']['rst_dir'] = os.path.dirname(self.in_rstfile)+'/' config['input']['surface']['wemin'] = wemin_in config['input']['surface']['catch_model'] = self.catch @@ -905,17 +933,27 @@ class LDASsetup: config['input']['shared']['agrid'] = 'C180' config['input']['shared']['ogrid'] = '1440x720' config['input']['shared']['omodel'] = 'data' - - catch_obj = catchANDcn(config_obj = config) - catch_obj.remap() + config['input']['shared']['MERRA-2'] = True + config['input']['surface']['catch_tilefile'] = '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/GM4/geometry/CF0180x6C_DE1440xPE0720/CF0180x6C_DE1440xPE0720-Pfafstetter.til' + + if self.with_land: + catch_obj = catchANDcn(config_obj = config) + catch_obj.remap() + if self.with_landice: + config['output']['surface']['remap_water'] = True + config['input']['surface']['zoom'] = '2' + landice_obj = lake_landice_saltwater(config_obj = config) + landice_obj.remap() #for ens in self.ensdirs : catchRstFile0 = '' vegdynRstFile0 = '' + landiceRstFile0 = '' for iens in range(self.nens) : ensdir = self.ensdirs[iens] ensid = self.ensids[iens] myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst' + myLandiceRst = myRstDir+'/'+ 'landice' +ensid +'_internal_rst' myVegRst = myRstDir+'/'+'vegdyn'+ensid +'_internal_rst' myPertRst = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst' @@ -935,16 +973,17 @@ class LDASsetup: if not os.path.isfile(vegdynRstFile): # no vegdyn restart from LDASsa if not os.path.isfile(vegdynRstFile0): vegdynRstFile = glob.glob(self.bcs_land + 'vegdyn_*.dat')[0] - else : + else: vegdynRstFile = glob.glob(self.bcs_land + 'vegdyn_*.dat')[0] - catchRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0] + if self.with_land: + catchRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0] # catchment restart file - if os.path.isfile(catchRstFile) : + if os.path.isfile(catchRstFile) and self.with_land : catchLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2 - if self.islocal : + if self.isZoomIn : print( "Creating local catchment restart file... \n") - cmd=self.bindir +'/preprocess_ldas.x c_localcatchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name + cmd=self.bindir +'/preprocess_ldas.x zoomin_catchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name print ("cmd: "+cmd) sp.call(shlex.split(cmd)) else : @@ -958,11 +997,11 @@ class LDASsetup: catchRstFile = catchRstFile0 # vegdyn restart file - if os.path.isfile(vegdynRstFile) : + if os.path.isfile(vegdynRstFile) and self.with_land : vegdynLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst' - if self.islocal : + if self.isZoomIn : print ("Creating the local veg restart file... \n") - cmd=self.bindir + '/preprocess_ldas.x c_localvegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name + cmd=self.bindir + '/preprocess_ldas.x zoomin_vegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) else : @@ -975,29 +1014,60 @@ class LDASsetup: else : vegdynRstFile = vegdynRstFile0 + landiceRstFile = '' + if self.with_landice : + if self.rqdExeInp['RESTART'].isdigit(): + if int(self.rqdExeInp['RESTART']) == 0 or int(self.rqdExeInp['RESTART']) == 2 : + print("RESTART=0 and RESTART=2 not supported for landice tiles. Please use RESTART=M (MERRA-2).") + landiceRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+'landice_internal_rst.'+y4m2d2_h2m2 + else: + landiceRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+'landice_internal_rst.'+YYYYMMDD+'*')[0] + + if os.path.isfile(landiceRstFile) : + landiceLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landice_internal_rst.'+y4m2d2_h2m2 + if self.isZoomIn : + print ("Creating zoom-in of landice restart file... \n") + cmd=self.bindir + '/preprocess_ldas.x zoomin_landicerst '+ landiceRstFile +' ' + landiceLocal + ' '+ tmp_f2g_file.name + print ("cmd: " + cmd) + sp.call(shlex.split(cmd)) + else : + shutil.copy(landiceRstFile,landiceLocal) + + landiceRstFile = landiceLocal + + if '0000' in ensdir : + landiceRstFile0 = landiceRstFile + else : + landiceRstFile = landiceRstFile0 + if (self.has_geos_pert and self.perturb == 1) : pertRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 pertLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 shutil.copy(pertRstFile,pertLocal) pertRstFile = pertLocal - print ('catchRstFile: ' + catchRstFile) - print ('vegdynRstFile: ' + vegdynRstFile) - os.symlink(catchRstFile, myCatchRst) - os.symlink(vegdynRstFile, myVegRst) + if self.with_land : + print ('catchRstFile: ' + catchRstFile) + print ('vegdynRstFile: ' + vegdynRstFile) + os.symlink(catchRstFile, myCatchRst) + os.symlink(vegdynRstFile, myVegRst) + if self.with_landice : + print("link landice restart: " + myLandiceRst) + os.symlink(landiceRstFile, myLandiceRst) if ( self.has_geos_pert and self.perturb == 1 ): - os.symlink(pertRstFile, myPertRst) + os.symlink(pertRstFile, myPertRst) # catch_param restar file catch_param_file = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_catparam.'+y4m2d2_h2m2+'z.bin' - assert os.path.isfile(catch_param_file), "need catch_param file %s" % catch_param_file + if self.with_land: + assert os.path.isfile(catch_param_file), "need catch_param file %s" % catch_param_file if self.has_mwrtm : mwRTMRstFile = self.mwrtm_file mwRTMLocal = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_mwRTMparam.'+y4m2d2_h2m2+'z.nc4' - if self.islocal : + if self.isZoomIn : print ("Creating the local mwRTM restart file... \n") - cmd= self.bindir +'/preprocess_ldas.x c_localmwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name + cmd= self.bindir +'/preprocess_ldas.x zoomin_mwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) @@ -1068,7 +1138,7 @@ class LDASsetup: # get optimzed NX and IMS optimized_distribution_file = tempfile.NamedTemporaryFile(delete=False) print ("Optimizing... decomposition of processes.... \n") - cmd = self.bindir + '/preprocess_ldas.x optimize '+ self.inpdir+'/tile.data '+ str(self.rqdRmInp['ntasks_model']) + ' ' + optimized_distribution_file.name + ' ' + self.rundir + cmd = self.bindir + '/preprocess_ldas.x optimize '+ self.inpdir+'/tile.data '+ str(self.rqdRmInp['ntasks_model']) + ' ' + optimized_distribution_file.name + ' ' + self.rundir + ' ' + '_'.join(self.tile_types) print ("cmd: " + cmd) print ("IMS.rc or JMS.rc would be generated on " + self.rundir) sp.call(shlex.split(cmd)) @@ -1181,37 +1251,46 @@ class LDASsetup: tmpl_='%s' if self.perturb == 1: ldasrcInp['PERTURBATIONS'] ='1' - bcval=['../input/green','../input/lai','../input/lnfm','../input/ndvi','../input/nirdf','../input/visdf'] - bckey=['GREEN','LAI','LNFM','NDVI','NIRDF','VISDF'] - for key, val in zip(bckey,bcval): - keyn = key+'_FILE' - valn = val+'.data' - ldasrcInp[keyn]= valn - if('catchcn' in self.catch): - ldasrcInp['CO2_MonthlyMean_DiurnalCycle_FILE']= '../input/CO2_MonthlyMean_DiurnalCycle.nc4' - else: - # remove catchcn-specific entries that do not apply to catch model - ldasrcInp.pop('DTCN',None) - ldasrcInp.pop('ATM_CO2',None) - ldasrcInp.pop('CO2',None) - ldasrcInp.pop('CO2_YEAR',None) - ldasrcInp.pop('PRESCRIBE_DVG',None) + rstkey =[] + rstval =[] + if self.with_land : + bcval=['../input/green','../input/lai','../input/lnfm','../input/ndvi','../input/nirdf','../input/visdf'] + bckey=['GREEN','LAI','LNFM','NDVI','NIRDF','VISDF'] + for key, val in zip(bckey,bcval): + keyn = key+'_FILE' + valn = val+'.data' + ldasrcInp[keyn]= valn + if('catchcn' in self.catch): + ldasrcInp['CO2_MonthlyMean_DiurnalCycle_FILE']= '../input/CO2_MonthlyMean_DiurnalCycle.nc4' + else: + # remove catchcn-specific entries that do not apply to catch model + ldasrcInp.pop('DTCN',None) + ldasrcInp.pop('ATM_CO2',None) + ldasrcInp.pop('CO2',None) + ldasrcInp.pop('CO2_YEAR',None) + ldasrcInp.pop('PRESCRIBE_DVG',None) # create restart item in RC - catch_ = self.catch.upper() - - if catch_+'_INTERNAL_RESTART_TYPE' in ldasrcInp : - # avoid duplicate - del ldasrcInp[ catch_ +'_INTERNAL_RESTART_TYPE'] - if catch_+'_INTERNAL_CHECKPOINT_TYPE' in ldasrcInp : - # avoid duplicate - del ldasrcInp[ catch_ +'_INTERNAL_CHECKPOINT_TYPE'] - if 'VEGDYN_INTERNAL_RESTART_TYPE' in ldasrcInp : - # avoid duplicate - del ldasrcInp['VEGDYN_INTERNAL_RESTART_TYPE'] - - rstkey=[catch_,'VEGDYN'] - rstval=[self.catch,'vegdyn'] + catch_ = self.catch.upper() + + if catch_+'_INTERNAL_RESTART_TYPE' in ldasrcInp : + # avoid duplicate + del ldasrcInp[ catch_ +'_INTERNAL_RESTART_TYPE'] + if catch_+'_INTERNAL_CHECKPOINT_TYPE' in ldasrcInp : + # avoid duplicate + del ldasrcInp[ catch_ +'_INTERNAL_CHECKPOINT_TYPE'] + if 'VEGDYN_INTERNAL_RESTART_TYPE' in ldasrcInp : + # avoid duplicate + del ldasrcInp['VEGDYN_INTERNAL_RESTART_TYPE'] + + rstkey.append(catch_) + rstkey.append('VEGDYN') + rstval.append(self.catch) + rstval.append('vegdyn') + + if self.with_landice: + rstkey.append('LANDICE') + rstval.append('landice') if self.has_mwrtm : keyn='LANDASSIM_INTERNAL_RESTART_FILE' @@ -1243,10 +1322,15 @@ class LDASsetup: ldasrcInp[keyn]= valn # checkpoint file and its type - keyn = catch_ + '_INTERNAL_CHECKPOINT_FILE' - valn = self.catch+tmpl_+'_internal_checkpoint' - ldasrcInp[keyn]= valn + if self.with_land : + keyn = catch_ + '_INTERNAL_CHECKPOINT_FILE' + valn = self.catch+tmpl_+'_internal_checkpoint' + ldasrcInp[keyn]= valn + if self.with_landice : + keyn = 'LANDICE_INTERNAL_CHECKPOINT_FILE' + valn = 'landice'+tmpl_+'_internal_checkpoint' + ldasrcInp[keyn]= valn # specify LANDPERT restart file if (self.perturb == 1): keyn = 'LANDPERT_INTERNAL_RESTART_FILE' @@ -1406,7 +1490,7 @@ class LDASsetup: MY_EXPDOMAIN = self.rqdExeInp['EXP_DOMAIN'], MY_LOGFILE = my_logfile, MY_ERRFILE = my_errfile, - MY_MODEL = self.catch, + MY_LANDMODEL = self.catch, MY_POSTPROC_HIST = str(self.rqdExeInp['POSTPROC_HIST']), MY_FIRST_ENS_ID = str(self.first_ens_id), MY_LADAS_COUPLING = str(self.ladas_coupling), diff --git a/GEOSldas_App/lenkf_j_template.py b/GEOSldas_App/lenkf_j_template.py index 6dcee191..ea50bd27 100644 --- a/GEOSldas_App/lenkf_j_template.py +++ b/GEOSldas_App/lenkf_j_template.py @@ -113,7 +113,7 @@ setenv HOMDIR $EXPDIR/run/ setenv SCRDIR $EXPDIR/scratch -setenv MODEL {MY_MODEL} +setenv LANDMODEL {MY_LANDMODEL} setenv MYNAME `finger $USER | cut -d: -f3 | head -1` setenv POSTPROC_HIST {MY_POSTPROC_HIST} @@ -460,7 +460,7 @@ set rc = -1 echo GEOSldas Run Status: $rc echo "ERROR: GEOSldas run FAILED, exit without post-processing" - exit + exit $rc endif @@ -721,13 +721,23 @@ if ( $NENS == 1) set ENSID ='' set THISDIR = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/ if (! -e $THISDIR ) mkdir -p $THISDIR + + set rstfs = (${{LANDMODEL}} 'landice') + foreach rstf ( $rstfs ) + if (-f ${{rstf}}${{ENSID}}_internal_checkpoint ) then + set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_internal_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} + /bin/mv ${{rstf}}${{ENSID}}_internal_checkpoint $tmp_file + /bin/rm -f $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst + /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst + endif + end - set rstf = ${{MODEL}} - if (-f ${{rstf}}${{ENSID}}_internal_checkpoint ) then - set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_internal_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} - /bin/mv ${{rstf}}${{ENSID}}_internal_checkpoint $tmp_file - /bin/rm -f $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst - /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst + set rstf = 'landassim_obspertrseed' + if (-f ${{rstf}}${{ENSID}}_checkpoint ) then + set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} + /bin/mv ${{rstf}}${{ENSID}}_checkpoint $tmp_file + /bin/rm -f $EXPDIR/input/restart/${{rstf}}${{ENSID}}_rst + /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_rst endif set rstf = 'landpert' @@ -742,29 +752,23 @@ /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst /usr/bin/gzip $old_rst & endif - - set rstf = 'landassim_obspertrseed' - if (-f ${{rstf}}${{ENSID}}_checkpoint ) then - set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} - /bin/mv ${{rstf}}${{ENSID}}_checkpoint $tmp_file - /bin/rm -f $EXPDIR/input/restart/${{rstf}}${{ENSID}}_rst - /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_rst - endif + # move intermediate check point files to output/$EXPDOMAIN/rs/$ENSDIR/Yyyyy/Mmm/ directories # ------------------------------------------------------------------------------------------- - - set rstfiles1 = `ls ${{MODEL}}${{ENSID}}_internal_checkpoint.*` + + set rstfiles1 = `ls ${{LANDMODEL}}${{ENSID}}_internal_checkpoint.*` set rstfiles2 = `ls landpert${{ENSID}}_internal_checkpoint.*` set rstfiles3 = `ls landassim_obspertrseed${{ENSID}}_checkpoint.*` - - foreach rfile ( $rstfiles1 ) + set rstfiles4 = `ls landice${{ENSID}}_internal_checkpoint.*` + + foreach rfile ( $rstfiles1 $rstfiles4 ) set ThisTime = `echo $rfile | rev | cut -d'.' -f2 | rev` set TY = `echo $ThisTime | cut -c1-4` set TM = `echo $ThisTime | cut -c5-6` set THISDIR = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{TY}}/M${{TM}}/ if (! -e $THISDIR ) mkdir -p $THISDIR - /bin/mv $rfile ${{THISDIR}}${{EXPID}}.${{MODEL}}_internal_rst.${{ThisTime}}.nc4 - /usr/bin/gzip ${{THISDIR}}${{EXPID}}.${{MODEL}}_internal_rst.${{ThisTime}}.nc4 & + /bin/mv $rfile ${{THISDIR}}${{EXPID}}.${{LANDMODEL}}_internal_rst.${{ThisTime}}.nc4 + /usr/bin/gzip ${{THISDIR}}${{EXPID}}.${{LANDMODEL}}_internal_rst.${{ThisTime}}.nc4 & end foreach rfile ( $rstfiles2 ) diff --git a/GEOSldas_App/preprocess_ldas.F90 b/GEOSldas_App/preprocess_ldas.F90 index 4d5725e6..67e3e9ef 100644 --- a/GEOSldas_App/preprocess_ldas.F90 +++ b/GEOSldas_App/preprocess_ldas.F90 @@ -5,12 +5,11 @@ program main use preprocess_ldas_routines, ONLY: & - createf2g, & - createLocalTilefile, & - createLocalBC, & - createLocalVegRestart, & - createLocalmwRTMRestart, & - createLocalCatchRestart, & + create_mapping, & + createZoominTilefile, & + createZoominBC, & + createZoominVegRestart, & + createZoominRestart, & correctEase, & convert_pert_rst, & optimize_latlon @@ -26,6 +25,7 @@ program main character(len=512) :: arg6 character(len=512) :: arg7 character(len=512) :: arg8 + character(len=512) :: arg9 character(len=512) :: orig_tile character(len=512) :: new_tile @@ -47,15 +47,19 @@ program main character(len=12 ) :: ymdhm character(len=12 ) :: SURFLAY - call get_command_argument(1,option) - call get_command_argument(2,arg1) - call get_command_argument(3,arg2) - call get_command_argument(4,arg3) - call get_command_argument(5,arg4) - call get_command_argument(6,arg5) - call get_command_argument(7,arg6) - call get_command_argument(8,arg7) - call get_command_argument(9,arg8) + character(len=:), allocatable :: new_r, orig_r + integer, allocatable :: int_types(:) + + call get_command_argument( 1,option) + call get_command_argument( 2,arg1) + call get_command_argument( 3,arg2) + call get_command_argument( 4,arg3) + call get_command_argument( 5,arg4) + call get_command_argument( 6,arg5) + call get_command_argument( 7,arg6) + call get_command_argument( 8,arg7) + call get_command_argument( 9,arg8) + call get_command_argument(10,arg9) if( trim(option) == "c_f2g") then @@ -71,46 +75,56 @@ program main SURFLAY = trim(adjustl(arg7)) f2g_file = arg8 - call createf2g(orig_tile,domain_def_file,trim(out_path),catch_def_file,trim(exp_id),ymdhm, SURFLAY, f2g_file) + call get_tile_types(trim(arg9), int_types) + + call create_mapping(orig_tile,domain_def_file,trim(out_path),catch_def_file,trim(exp_id),ymdhm, SURFLAY, f2g_file, int_types) - else if (trim(option) == "c_localtile") then + else if (trim(option) == "zoomin_tile") then orig_tile = arg1 new_tile = arg2 f2g_file = arg3 - call createLocalTilefile(f2g_file, orig_tile,new_tile) + call createZoominTilefile(f2g_file, orig_tile,new_tile) - else if (trim(option) == "c_localbc" ) then + else if (trim(option) == "zoomin_bc" ) then orig_BC = arg1 new_BC = arg2 f2g_file = arg3 - call createLocalBC(f2g_file, orig_BC, new_BC) + call createZoominBC(f2g_file, orig_BC, new_BC) - else if (trim(option) == "c_localvegrst") then + else if (trim(option) == "zoomin_vegrst") then orig_veg = arg1 new_veg = arg2 f2g_file = arg3 - call createLocalVegRestart(f2g_file, orig_veg, new_veg) + call createZoominVegRestart(f2g_file, orig_veg, new_veg) - else if (trim(option) == "c_localmwrtmrst") then + else if (trim(option) == "zoomin_mwrtmrst") then orig_rtm = arg1 new_rtm = arg2 f2g_file = arg3 - call createLocalmwRTMRestart(f2g_file, orig_rtm, new_rtm) + call createZoominRestart(f2g_file, orig_rtm, new_rtm, 100) - else if (trim(option) == "c_localcatchrst") then + else if (trim(option) == "zoomin_catchrst") then orig_catch = arg1 new_catch = arg2 f2g_file = arg3 - call createLocalCatchRestart(f2g_file, orig_catch, new_catch) + call createZoominRestart(f2g_file, orig_catch, new_catch, 100) + + else if (trim(option) == "zoomin_landicerst") then + + orig_r = trim(arg1) + new_r = trim(arg2) + f2g_file = trim(arg3) + + call createZoominRestart(f2g_file, orig_r, new_r, 20) else if (trim(option)=="correctease") then @@ -128,15 +142,43 @@ program main else if (trim(option) == "optimize") then - - call optimize_latlon(arg1,arg2, arg3, arg4) + call get_tile_types(trim(arg5), int_types) + call optimize_latlon(arg1,arg2, arg3, arg4, int_types) else print*, " wrong preprocess option:",option end if - + +contains + + subroutine get_tile_types(str_types, int_types) + + character(*), intent(in) :: str_types + integer, allocatable, intent(out) :: int_types(:) + + integer :: n, Length, from, to, i + n = 1 + Length = len(str_types) + do i = 1, Length + if (str_types(i:i) == '_') n = n+1 + enddo + allocate(int_types(n)) + from = 0 + to = 1 + n = 1 + do while (to <= Length) + if (str_types(to:to) == "_") then + read (unit=str_types(from+1:to-1),fmt=*) int_types(n) + n = n + 1 + from = to + endif + to = to + 1 + enddo + read (unit=str_types(from+1:to-1),fmt=*) int_types(n) + end subroutine get_tile_types + end program main ! ====================== EOF ======================================================= diff --git a/GEOSldas_App/preprocess_ldas_routines.F90 b/GEOSldas_App/preprocess_ldas_routines.F90 index ea34807c..382a549d 100644 --- a/GEOSldas_App/preprocess_ldas_routines.F90 +++ b/GEOSldas_App/preprocess_ldas_routines.F90 @@ -20,6 +20,29 @@ module preprocess_ldas_routines ! - is_in_domain() [from LDAS_ensdrv_functions.F90] ! - word_count() [from LDAS_ensdrv_functions.F90] ! - open_land_param_file() [from LDAS_ensdrv_functions.F90] + ! + ! Domains: + ! + ! g = *g*lobal domain : list of tiles defined in tile file from boundary conditions + ! r = *r*estart domain : list of tiles included in restart file(s) + ! f = *f*ull simulation domain : list of tiles included in simulation domain + ! ("full" domain as opposed to a "local" subdomain assigned to a given processor) + ! + ! The simulation and restart domains can be subsets of the global domain ("zoom-in"), + ! as long as the simulation domain is a subset of the restart domain. + ! + ! Mapping between domains: + ! + ! r2g : restart domain to global domain + ! f2r : simulation domain to restart domain + ! f2g : simulation domain to global domain, f2g = r2g(f2r) + ! + ! Number of tiles: + ! + ! N_tiles_[g,r,f] : number of tiles included in [g,r,f] domain (across all tile types) + ! N_tiles_land_[g,r,f] : number of *land* tiles included in [g,r,f] domain + ! + ! ---------------------------------------------------------- use netcdf @@ -80,12 +103,11 @@ module preprocess_ldas_routines private - public :: createf2g - public :: createLocalTilefile - public :: createLocalBC - public :: createLocalCatchRestart - public :: createLocalVegRestart - public :: createLocalmwRTMRestart + public :: create_mapping + public :: createZoominTilefile + public :: createZoominBC + public :: createZoominRestart + public :: createZoominVegRestart public :: correctEase public :: optimize_latlon public :: convert_pert_rst @@ -93,138 +115,147 @@ module preprocess_ldas_routines character(10), private :: tmpstring10 character(40), private :: tmpstring40 - ! Tile type for land that is to be excluded from the simulation domain. - ! (GEOSldas allows for non-global simulations and repeated "zooming" - ! of the domain while MAPL generally assumes a complete (global) tile - ! space. The *_ExcludeFromDomain tile type makes it possible to work - ! with complete (global) tile files (ie, make use of MAPL functionality) - ! and also maintain GEOSldas functionality. + ! GEOSldas allows for non-global simulations and repeated "zooming" of + ! the domain, while MAPL generally assumes a complete (global) tile space. + ! By adding the integer ExcludeFromDomain to the original MAPL_* tile + ! type for all tiles that should be excluded from the simulation domain, we + ! can work with global tile files (i.e., use MAPL tools) and also maintain + ! the non-global domain unctionality of GEOSldas. - integer, parameter :: MAPL_Land_ExcludeFromDomain = 1100 + integer, parameter :: ExcludeFromDomain = 1000 ! for tiles to be excluded from the simulation domain, this number is added to tile type contains ! ******************************************************************** - subroutine createf2g(orig_tile,domain_def,out_path,catch_def_file,exp_id,ymdhm, SURFLAY, f2g_file) - - implicit none - character(*) :: orig_tile - character(*) :: domain_def - character(*) :: out_path - character(*) :: catch_def_file - character(*) :: exp_id - character(*) :: ymdhm - character(*) :: SURFLAY - character(*) :: f2g_file - - real :: minlon,maxlon,minlat,maxlat - character(len=512):: exclude_file,include_file - character(len=512):: bcs_path - logical :: file_exist - logical :: d_exist,c_exist - - integer :: n - - type(grid_def_type) :: tile_grid_g,tile_grid_d - type(tile_coord_type), dimension(:), pointer :: tile_coord_g => null() - type(tile_coord_type), dimension(:), pointer :: tile_coord_d => null() - integer, dimension(:), pointer :: f2g => null() - integer, dimension(:), pointer :: d2g => null() - integer, dimension(:), pointer :: d2f => null() - integer :: N_catg, N_catd,n1,n2,N_catf - + subroutine create_mapping( orig_tile, domain_def, out_path, catch_def_file, exp_id, ymdhm, SURFLAY, mapping_file, types) + + character(*), intent(in) :: orig_tile + character(*), intent(in) :: domain_def + character(*), intent(in) :: out_path + character(*), intent(in) :: catch_def_file + character(*), intent(in) :: exp_id + character(*), intent(in) :: ymdhm + character(*), intent(in) :: SURFLAY + character(*), intent(in) :: mapping_file + + integer, dimension(:), optional, intent(in) :: types + + ! ---------------------------------------------------- + + real :: minlon,maxlon,minlat,maxlat + character(len=512 ):: exclude_file,include_file + character(len=512 ):: bcs_path + logical :: file_exist + logical :: d_exist,c_exist + + type(grid_def_type) :: tile_grid_g,tile_grid_f + type(tile_coord_type), dimension(:), pointer :: tile_coord_r => null() + type(tile_coord_type), dimension(:), pointer :: tile_coord_f => null() + integer, dimension(:), pointer :: f2r => null() + integer, dimension(:), pointer :: r2g => null() ! restart domain to global + integer :: N_tiles_land_g, N_types, n, n1, N_tiles_land_f + integer, dimension(:), allocatable :: tile_types + integer, dimension(:), allocatable :: N_tiles_r, N_tiles_f type(cat_param_type), dimension(:), allocatable :: cp - real :: dzsf + real :: dzsf - namelist / domain_inputs / & + namelist / domain_inputs / & minlon, maxlon,minlat,maxlat, & exclude_file,include_file inquire(file=trim(orig_tile),exist=file_exist) - if( .not. file_exist) stop ("original tile file not exist") + if( .not. file_exist) stop ("original tile file does not exist") inquire(file=trim(domain_def),exist=d_exist) if( .not. d_exist) then - print*,"no domain definition file" + print*, "domain definition file does not exist", trim(domain_def) endif inquire(file=trim(catch_def_file),exist=c_exist) if( .not. c_exist) then - print*,"no catchment definition file:" , catch_def_file + print*,"catchment [land tile supplemental] definition file does not exist:" , trim(catch_def_file) + endif + + if (present(types)) then + ! reorder tile types so it matches the tile file (first LAND then LAKE then LANDICE) + allocate(tile_types(size(types))) + n = 1 + if (any(types == MAPL_LAND)) then + tile_types(n) = MAPL_LAND + n = n+1 + endif + if (any(types == MAPL_LAKE)) then + tile_types(n) = MAPL_LAKE + n = n+1 + endif + if (any(types == MAPL_LANDICE)) then + tile_types(n) = MAPL_LANDICE + endif + else + tile_types = [MAPL_LAND] endif + call LDAS_read_til_file( orig_tile, catch_def_file, tile_grid_g, tile_coord_r, r2g, N_tiles_land_g, tile_types) + ! include and exclude files are absolute if(d_exist) then open (10, file=trim(domain_def), delim='apostrophe', action='read', status='old') read (10, nml= domain_inputs) close(10) else minlon = -180. - maxlon = 180. - minlat = -90. - maxlat = 90. + maxlon = 180. + minlat = -90. + maxlat = 90. exclude_file = ' ' include_file = ' ' endif - call LDAS_read_til_file(orig_tile,catch_def_file,tile_grid_g,tile_coord_g,f2g) - - N_catg=size(tile_coord_g) - - ! include and exclude files are absolute - - call domain_setup( & - N_catg, tile_coord_g, & + call domain_setup( & + tile_coord_r, & tile_grid_g, & - ' ', exclude_file, ' ', include_file, & + ' ', exclude_file, ' ', include_file, & trim(out_path), 'exp_domain ', trim(exp_id), & minlon, minlat, maxlon, maxlat, & - N_catd, d2g, tile_coord_d, & - tile_grid_d ) - - allocate(cp(N_catd)) + f2r, tile_coord_f, & + tile_grid_f) + + N_tiles_land_f = count(tile_coord_f(:)%typ == MAPL_LAND) + + allocate(cp(N_tiles_land_f)) read(SURFLAY,*) dzsf print*, "SURFLAY: ", dzsf n1 = index(catch_def_file,'/clsm/') bcs_path(1:n1-1) = catch_def_file(1:n1-1) - call read_cat_param( N_catg, N_catd, d2g, tile_coord_d, dzsf, bcs_path(1:n1-1), bcs_path(1:n1-1),bcs_path(1:n1-1), & + call read_cat_param( N_tiles_land_g, N_tiles_land_f, tile_coord_f, dzsf, bcs_path(1:n1-1), bcs_path(1:n1-1),bcs_path(1:n1-1), & cp ) - call write_cat_param(cp,N_catd) - - allocate(d2f(N_catd)) - d2f = 0 - N_catf = size(f2g) - if( N_catf /= N_catg) then - n = 1 - do n1 = 1,N_catd - do n2 = n, N_catf - if (d2g(n1) == f2g(n2)) then - d2f(n1) = n2 - n = n2+1 - exit - endif - enddo - enddo - if(any(d2f == 0)) stop " Domain includes those excluded tiles" - print*," f2g now is d2f " + call write_cat_param(cp,N_tiles_land_f) + + N_types = size(tile_types) + allocate(N_tiles_r(N_types), N_tiles_f(N_types)) + do n=1, N_types + N_tiles_r(n) = count(tile_coord_r%typ == tile_types(n)) + N_tiles_f(n) = count(tile_coord_f%typ == tile_types(n)) + enddo + + if (any(N_tiles_r /= N_tiles_f)) then + print*,"writing mapping file, N_types....", N_types + open(40,file=mapping_file,form='formatted',action='write') + write(40,*)N_types + write(40,*)tile_types + write(40,*)N_tiles_r ! length N_types + write(40,*)N_tiles_f ! length N_types + write(40,*)f2r + write(40,*)r2g + close(40) else - d2f = d2g + print*,"No mapping file is created" endif - open(40,file=f2g_file,form='formatted',action='write') - write(40,*)N_catf - write(40,*)N_catd - do n=1,N_catd - write(40,*)d2f(n) - enddo - do n=1,N_catd - write(40,*)d2g(n) - enddo - close(40) - if (associated(f2g)) deallocate(f2g) - if (associated(d2g)) deallocate(d2g) - if (associated(d2f)) deallocate(d2f) + + if (associated(f2r)) deallocate(f2r) + if (associated(r2g)) deallocate(r2g) contains @@ -236,8 +267,6 @@ logical function is_in_list(N_list, list, this_one) ! reichle, 2 May 2003 - implicit none - integer :: N_list, this_one integer, dimension(N_list) :: list @@ -262,34 +291,32 @@ end function is_in_list logical function is_in_domain( & this_cat_exclude, this_cat_include, this_cat_in_box ) - ! determine whether catchment is in domain + ! determine whether tile is in domain ! - ! The domain is set up using (if present) an "ExcludeList" of catchments - ! to be excluded, an "IncludeList" (if present) of catchments to be included, + ! The domain is set up using (if present) an "ExcludeList" of tiles + ! to be excluded, an "IncludeList" (if present) of tiles to be included, ! and the bounding box of a rectangular "zoomed" area (as specified ! in the "exeinp" file used in ldas_setup). ! ! order of precedence: - ! 1. exclude catchments on ExcludeList - ! 2. include catchments on IncludeList or catchments within rectangular domain - ! (i.e., catchments in ExcludeList are *always* excluded) + ! 1. exclude tiles on ExcludeList + ! 2. include tiles on IncludeList or tiles within rectangular domain + ! (i.e., tiles in ExcludeList are *always* excluded) ! ! reichle, 7 May 2003 ! reichle, 9 May 2005 - redesign (no more continents) ! ! ---------------------------------------------------------------- - implicit none - logical :: this_cat_include, this_cat_exclude, this_cat_in_box is_in_domain = .false. - ! if catchment is NOT in ExcludeList + ! if tile is NOT in ExcludeList if (.not. this_cat_exclude) then - ! if catchment is within bounding box OR in IncludeList + ! if tile is within bounding box OR in IncludeList if ((this_cat_in_box) .or. (this_cat_include)) then @@ -306,9 +333,7 @@ logical function is_cat_in_box( & this_minlon, this_minlat, this_maxlon, this_maxlat, & minlon, minlat, maxlon, maxlat ) - ! determine whether catchment is within bounding box - reichle, 7 May 2003 - - implicit none + ! determine whether tile is within bounding box - reichle, 7 May 2003 real :: this_minlon, this_minlat, this_maxlon, this_maxlat real :: minlon, minlat, maxlon, maxlat @@ -327,29 +352,29 @@ end function is_cat_in_box ! ******************************************************************** subroutine domain_setup( & - N_cat_global, tile_coord_global, & + tile_coord_r, & tile_grid_g, & exclude_path, exclude_file, include_path, include_file, & work_path, exp_domain, exp_id, & minlon, minlat, maxlon, maxlat, & - N_cat_domain, d2g, tile_coord, tile_grid_d ) + f2r, tile_coord_f, tile_grid_f ) ! Set up modeling domain and determine index vectors mapping from the - ! domain to global catchment space. + ! domain to global tile space. ! Determine actual bounding box for domain. ! Also return tile_coord for domain and tile_grid_d for domain. ! ! ----------------------- ! - ! The domain is set up using (if present) an "ExcludeList" of catchments - ! to be excluded, an "IncludeList" (if present) of catchments to be included, + ! The domain is set up using (if present) an "ExcludeList" of tiles + ! to be excluded, an "IncludeList" (if present) of tiles to be included, ! and the bounding box of a rectangular "zoomed" area (as specified ! in the "exeinp" file used in ldas_setup). ! ! order of precedence: - ! 1. exclude catchments in ExcludeList - ! 2. include catchments in IncludeList or catchments within rectangular domain - ! (i.e., catchments in ExcludeList are *always* excluded) + ! 1. exclude tiles in ExcludeList + ! 2. include tiles in IncludeList or tiles within rectangular domain + ! (i.e., tiles in ExcludeList are *always* excluded) ! ! input: ! @@ -360,16 +385,12 @@ subroutine domain_setup( & ! latitude -90:90 ! ! output: - ! N_cat_domain = number of catchments in zoomed domain - ! (for which model integration is conducted) - ! d2g = index from domain to global tiles + ! f2r = index mapping from (full) simulation domain to restart domain ! tile_coord_d = tile_coord vector for domain ! tile_grid_d = def of smallest subgrid of global tile_grid_g that contains - ! all catchments (or tiles) in the domain (tile_grid_d%i_offg, + ! all tiles (or tiles) in the domain (tile_grid_d%i_offg, ! tile_grid_d%j_offg are offsets in indices between tile_grid_g ! and tile_grid_d) - ! N_catd_cont = number of catchments of (full) domain on each continent - ! ! ! - reichle, May 7, 2003 ! - reichle, Nov 7, 2003 - computation of bounding box of actual domain @@ -381,12 +402,7 @@ subroutine domain_setup( & ! ! ---------------------------------------------------------- - implicit none - - integer, intent(in) :: N_cat_global - - type(tile_coord_type), dimension(:), pointer :: tile_coord_global ! input - + type(tile_coord_type), dimension(:), pointer :: tile_coord_r ! input type(grid_def_type), intent(in) :: tile_grid_g character(*), intent(in) :: exclude_path, include_path @@ -399,23 +415,21 @@ subroutine domain_setup( & real, intent(in) :: minlon, minlat ! from nml inputs real, intent(in) :: maxlon, maxlat ! from nml inputs - integer, intent(out) :: N_cat_domain + integer, dimension(:), pointer :: f2r ! output - integer, dimension(:), pointer :: d2g ! output + type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! output - type(tile_coord_type), dimension(:), pointer :: tile_coord ! output - - type(grid_def_type), intent(out) :: tile_grid_d + type(grid_def_type), intent(out) :: tile_grid_f ! locals - integer :: n, this_tileid, N_exclude, N_include, indomain, rc + integer :: n, this_tileid, N_exclude, N_include, indomain, rc, n_tiles_r, n_tiles_f - integer, dimension(N_cat_global) :: ExcludeList, IncludeList, tmp_d2g + integer, dimension(:), allocatable :: ExcludeList, IncludeList real :: this_minlon, this_minlat, this_maxlon, this_maxlat - logical :: this_cat_exclude, this_cat_include, this_cat_in_box + logical :: this_tile_exclude, this_tile_include, this_cat_in_box integer :: this_i_indg, this_j_indg @@ -425,7 +439,8 @@ subroutine domain_setup( & character(len=*), parameter :: Iam = 'domain_setup' character(len=400) :: err_msg - + integer, allocatable, dimension(:):: tile_types, tmp_f2r + ! ------------------------------------------------------------ if (logit) write (logunit,*) 'Setting up domain: ' @@ -436,7 +451,7 @@ subroutine domain_setup( & ! try reading *domain.txt, *tilecoord.txt, and *tilegrids.txt files call io_domain_files( 'r', work_path, exp_id, & - N_cat_domain, d2g, tile_coord, tmp_grid_def, tile_grid_d, rc ) + n_tiles_f, f2r, tile_coord_f, tmp_grid_def, tile_grid_f, rc ) if (rc==0) then ! read was successful @@ -456,90 +471,93 @@ subroutine domain_setup( & print*, "Creating domain..., reading IncludeList and ExludeList if present..." ! ------------------------------------------------------------ ! - ! load ExcludeList: catchments listed in this file will *always* be excluded + ! load ExcludeList: tiles listed in this file will *always* be excluded fname = trim(exclude_path) // '/' // trim(exclude_file) + + call read_exclude_or_includelist(fname, ExcludeList) - call read_exclude_or_includelist(N_cat_global, fname, ExcludeList, N_exclude) - - ! load IncludeList: catchments listed in this file will be included + ! load IncludeList: tiles listed in this file will be included ! (unless excluded via ExcludeList) fname = trim(include_path) // '/' // trim(include_file) - call read_exclude_or_includelist(N_cat_global, fname, IncludeList, N_include) + call read_exclude_or_includelist(fname, IncludeList) ! ----------------- ! - ! find and count catchments that are in the domain + ! find and count tiles that are in the domain c3_grid = .false. if(index(tile_grid_g%gridtype,"c3")/=0) c3_grid = .true. indomain = 0 ! initialize - - do n=1,N_cat_global - - this_tileid = tile_coord_global(n)%tile_id + + n_tiles_r = size(tile_coord_r) + allocate(tmp_f2r(n_tiles_r)) + do n=1,n_tiles_r - if( .not. c3_grid) then - this_minlon = tile_coord_global(n)%min_lon - this_minlat = tile_coord_global(n)%min_lat - this_maxlon = tile_coord_global(n)%max_lon - this_maxlat = tile_coord_global(n)%max_lat + this_tileid = tile_coord_r(n)%tile_id + + if( .not. c3_grid .and. tile_coord_r(n)%typ == MAPL_LAND) then + this_minlon = tile_coord_r(n)%min_lon + this_minlat = tile_coord_r(n)%min_lat + this_maxlon = tile_coord_r(n)%max_lon + this_maxlat = tile_coord_r(n)%max_lat else ! c3 grid can straddle the lat-lon - this_minlon = tile_coord_global(n)%com_lon - this_minlat = tile_coord_global(n)%com_lat - this_maxlon = tile_coord_global(n)%com_lon - this_maxlat = tile_coord_global(n)%com_lat + ! WY Note: not sure if it is right + this_minlon = tile_coord_r(n)%com_lon + this_minlat = tile_coord_r(n)%com_lat + this_maxlon = tile_coord_r(n)%com_lon + this_maxlat = tile_coord_r(n)%com_lat endif - - this_cat_exclude = is_in_list( N_exclude, ExcludeList(1:N_exclude), this_tileid ) - this_cat_include = is_in_list( N_include, IncludeList(1:N_include), this_tileid ) + N_exclude = size(ExcludeList) + N_include = size(IncludeList) + this_tile_exclude = is_in_list( N_exclude, ExcludeList(1:N_exclude), this_tileid ) + this_tile_include = is_in_list( N_include, IncludeList(1:N_include), this_tileid ) this_cat_in_box = & is_cat_in_box(this_minlon,this_minlat,this_maxlon,this_maxlat, & minlon, minlat, maxlon, maxlat ) if (is_in_domain( & - this_cat_exclude, this_cat_include, this_cat_in_box )) then + this_tile_exclude, this_tile_include, this_cat_in_box )) then indomain = indomain + 1 - tmp_d2g(indomain) = n + tmp_f2r(indomain) = n end if end do - N_cat_domain = indomain + n_tiles_f = indomain - if (N_cat_domain .eq. 0) then - err_msg = 'No catchments found in domain' + if (n_tiles_f .eq. 0) then + err_msg = 'No tiles found in domain' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) else if (logit) then - write (logunit,*) 'Number of catchments in domain = ', N_cat_domain + write (logunit,*) 'Number of tiles in domain = ', n_tiles_f write (logunit,*) end if end if ! ------------------------------------------------------------------- ! - ! assemble d2g, tile_coord, tile_grid_d + ! assemble f2r, tile_coord, tile_grid_d - allocate(d2g( N_cat_domain)) - allocate(tile_coord(N_cat_domain)) + allocate(f2r( n_tiles_f)) + allocate(tile_coord_f(n_tiles_f)) - d2g(1:N_cat_domain) = tmp_d2g(1:N_cat_domain) - - tile_coord = tile_coord_global(d2g) + f2r(1:n_tiles_f) = tmp_f2r(1:n_tiles_f) + tile_coord_f = tile_coord_r(f2r) ! finalize extent of actual domain: ! determine smallest subgrid of tile_grid_d that contains all - ! catchments/tiles in domain + ! tiles in domain - tile_grid_d = get_minExtent_grid( N_cat_domain, tile_coord%i_indg, tile_coord%j_indg, & - tile_coord%min_lon, tile_coord%min_lat, tile_coord%max_lon, tile_coord%max_lat, & + tile_grid_f = get_minExtent_grid( n_tiles_f, tile_coord_f%i_indg, tile_coord_f%j_indg, & + tile_coord_f%min_lon, tile_coord_f%min_lat, tile_coord_f%max_lon, tile_coord_f%max_lat, & tile_grid_g) ! output domain files @@ -547,22 +565,22 @@ subroutine domain_setup( & tmp_grid_def = tile_grid_g ! cannot use intent(in) tile_grid_g w/ io_domain_files call io_domain_files( 'w', work_path, exp_id, & - N_cat_domain, d2g, tile_coord, tmp_grid_def, tile_grid_d, rc ) + n_tiles_f, f2r, tile_coord_f, tmp_grid_def, tile_grid_f, rc ) end if ! domain/tilecoord/tilegrids files exist - ! output extent of domain and tile_grid_d to logunit + ! output extent of domain and tile_grid_f to logunit if (logit) write (logunit,*) 'Actual extent of domain grid:' - if (logit) write (logunit,*) 'min lon = ', tile_grid_d%ll_lon - if (logit) write (logunit,*) 'max lon = ', tile_grid_d%ur_lon - if (logit) write (logunit,*) 'min lat = ', tile_grid_d%ll_lat - if (logit) write (logunit,*) 'max lat = ', tile_grid_d%ur_lat + if (logit) write (logunit,*) 'min lon = ', tile_grid_f%ll_lon + if (logit) write (logunit,*) 'max lon = ', tile_grid_f%ur_lon + if (logit) write (logunit,*) 'min lat = ', tile_grid_f%ll_lat + if (logit) write (logunit,*) 'max lat = ', tile_grid_f%ur_lat if (logit) write (logunit,*) tmpstring40 = 'tile_grid_d' - if (logit) call io_grid_def_type('w', logunit, tile_grid_d, tmpstring40) + if (logit) call io_grid_def_type('w', logunit, tile_grid_f, tmpstring40) print*, "Done with " // trim(Iam) @@ -570,33 +588,24 @@ end subroutine domain_setup ! ************************************************************************* - subroutine read_exclude_or_includelist(N_cat, fname, MyList, N_list) + subroutine read_exclude_or_includelist( fname, MyList) - ! read numbers/IDs of catchments in MyList (ExcludeList or IncludeList) + ! read numbers/IDs of tiles in MyList (ExcludeList or IncludeList) ! ! format of MyList file: ASCII list of tile IDs ! - ! N_list = number of catchments in MyList + ! N_list = number of tiles in MyList ! ! reichle, 2 May 2003 ! ! -------------------------------------------------------------- - implicit none - - ! N_cat = max number of catchments allowed in list - ! (use N_cat_global when calling this subroutine) - - integer, intent(in) :: N_cat character(*), intent(in) :: fname - - integer, intent(out) :: N_list - - integer, dimension(N_cat), intent(out) :: MyList + integer, dimension(:), allocatable, intent(out) :: MyList ! locals - integer :: istat, tmpint + integer :: istat, i, N_list logical :: file_exists @@ -610,44 +619,29 @@ subroutine read_exclude_or_includelist(N_cat, fname, MyList, N_list) inquire( file=fname, exist=file_exists) if (file_exists) then - + N_list = 0 + OPEN (10, file = fname) + if (logit) write (logunit,*) & + 'reading ExcludeList or IncludeList from ', trim(fname) + if (logit) write (logunit,*) + DO + READ(10,*,iostat=istat) + IF (istat/=0) EXIT + N_list = N_list+1 + END DO + CLOSE (10) + + if (allocated(MyList)) deallocate(MyList) + allocate(MyList(N_List)) + open(10, file=fname, form='formatted', action='read', & status='old', iostat=istat) if (istat==0) then - - if (logit) write (logunit,*) & - 'reading ExcludeList or IncludeList from ', trim(fname) - if (logit) write (logunit,*) - - do - read(10,*,iostat=istat) tmpint - - if (istat==-1) then - if (logit) write (logunit,*) ' found ', N_list, ' catchments on list' - exit - else if (istat/=0) then - err_msg = 'read error other than end-of-file' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - else - N_list = N_list+1 - MyList(N_list) = tmpint - end if - - if (N_list>N_cat) then - - write (tmpstring10,*) N_cat - write (tmpstring40,*) N_list - - err_msg = 'N_list=' // trim(tmpstring40) & - // ' > N_cat=' // trim(tmpstring10) - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end if + do i = 1, N_list + read(10,*,iostat=istat) MyList(i) end do - close(10,status='keep') - else if (logit) write (logunit,*) & @@ -656,10 +650,12 @@ subroutine read_exclude_or_includelist(N_cat, fname, MyList, N_list) end if else + if (allocated(MyList)) deallocate(MyList) + allocate(MyList(0)) if (logit) write (logunit,*) & 'ExcludeList or IncludeList file does not exist: ', trim(fname) - + end if if (logit) write (logunit,*) @@ -674,8 +670,6 @@ integer function word_count( mystring ) ! ! - reichle, 31 Mar 2015 - implicit none - character(len=*) :: mystring integer :: N_words, N_string, ii @@ -722,8 +716,6 @@ integer function open_land_param_file( unitnumber, formatted_file, is_big_endian ! ignore_stop = optional input, if present and .true., skip call to "stop_it()" - implicit none - integer :: unitnumber, N_search_dir logical :: formatted_file @@ -823,8 +815,8 @@ end function open_land_param_file ! ***************************************************************************************** - subroutine read_cat_param( & - N_catg, N_catf, f2g, tile_coord_f, dzsf, veg_path, soil_path, top_path, & + subroutine read_cat_param( & + N_tiles_land_g, N_tiles_land_f, tile_coord_f, dzsf, veg_path, soil_path, top_path, & cp ) ! Reads soil properties and topographic parameters from global files @@ -845,21 +837,17 @@ subroutine read_cat_param( & ! ! ------------------------------------------------------------------- - implicit none + integer, intent(in) :: N_tiles_land_g, N_tiles_land_f - integer, intent(in) :: N_catg, N_catf + type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! intent(in) - type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! intent(in) + real, intent(in) :: dzsf - real, intent(in) :: dzsf + character(*), intent(in) :: veg_path + character(*), intent(in) :: soil_path + character(*), intent(in) :: top_path - integer, dimension(N_catf), intent(in) :: f2g - - character(*), intent(in) :: veg_path - character(*), intent(in) :: soil_path - character(*), intent(in) :: top_path - - type(cat_param_type), dimension(N_catf), intent(out) :: cp + type(cat_param_type), dimension(N_tiles_land_f), intent(out) :: cp ! local variables @@ -871,11 +859,11 @@ subroutine read_cat_param( & character(100), dimension(N_search_dir_max) :: search_dir - integer :: n, k, m, dummy_int, dummy_int2, istat, N_search_dir, N_col + integer :: n, k, m, dummy_int, dummy_int2, istat, N_search_dir, N_col, gid - integer, dimension(N_catg) :: tmpint, tmpint2, tmptileid + integer, dimension(N_tiles_land_g) :: tmpint, tmpint2, tmptileid - real, dimension(N_catg,N_col_real_max) :: tmpreal + real, dimension(N_tiles_land_g,N_col_real_max) :: tmpreal real :: dummy_real, dummy_real2, z_in_m, term1, term2 @@ -883,7 +871,7 @@ subroutine read_cat_param( & character(len=*), parameter :: Iam = 'read_cat_param' character(len=400) :: err_msg - + real, dimension(NTYPS) :: VGZ2 ! legacy vegetation height look-up table (for backward compatibility) @@ -938,7 +926,7 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'Using vegetation height look-up table' - do n=1,N_catg + do n=1,N_tiles_land_g read (10,*) tmptileid(n), dummy_int, tmpint(n) @@ -951,7 +939,7 @@ subroutine read_cat_param( & end if - do n=1,N_catg + do n=1,N_tiles_land_g tmpreal(n,1) = VGZ2( tmpint(n) ) @@ -964,7 +952,7 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'reading vegetation height from file' - do n=1,N_catg + do n=1,N_tiles_land_g ! 7-th column contains veg height in m ! 8-th column contains ASCAT z0 values (IGNORED for now, reichle, 31 Oct 2017) @@ -985,18 +973,19 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - - do k=1,N_catf + + do k=1,N_tiles_land_f ! this check works only for "SiB2_V2" and newer versions - - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + gid = tile_coord_f(k)%tile_id ! id of tile in *g*lobal tile file + if (gid /=tmptileid(gid)) then err_msg = 'something wrong with veg parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - cp(k)%vegcls = tmpint( f2g(k) ) - cp(k)%veghght = tmpreal(f2g(k),1) + cp(k)%vegcls = tmpint(gid ) + + cp(k)%veghght = tmpreal(gid,1) end do @@ -1039,7 +1028,7 @@ subroutine read_cat_param( & tmpreal = nodata_generic - do n=1,N_catg + do n=1,N_tiles_land_g ! "SiB2_V2" version @@ -1053,22 +1042,22 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - do k=1,N_catf - - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + do k=1,N_tiles_land_f + gid = tile_coord_f(k)%tile_id + if (tile_coord_f(k)%tile_id/=tmptileid(gid)) then err_msg = 'something wrong with soil parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - cp(k)%bee = tmpreal(f2g(k),1) - cp(k)%psis = tmpreal(f2g(k),2) - cp(k)%poros = tmpreal(f2g(k),3) - cp(k)%cond = tmpreal(f2g(k),4) - cp(k)%wpwet = tmpreal(f2g(k),5) - cp(k)%dpth = tmpreal(f2g(k),6) + cp(k)%bee = tmpreal(gid,1) + cp(k)%psis = tmpreal(gid,2) + cp(k)%poros = tmpreal(gid,3) + cp(k)%cond = tmpreal(gid,4) + cp(k)%wpwet = tmpreal(gid,5) + cp(k)%dpth = tmpreal(gid,6) - cp(k)%soilcls30 = tmpint( f2g(k)) - cp(k)%soilcls100 = tmpint2(f2g(k)) + cp(k)%soilcls30 = tmpint( gid) + cp(k)%soilcls100 = tmpint2(gid) end do @@ -1085,23 +1074,23 @@ subroutine read_cat_param( & ! "Icarus-NLv4" has 20 columns (new, last column is peat fraction, ignore for now) - do k=1,N_catf - - cp(k)%gravel30 = tmpreal(f2g(k), 7) - cp(k)%orgC30 = tmpreal(f2g(k), 8) - cp(k)%orgC = tmpreal(f2g(k), 9) - cp(k)%sand30 = tmpreal(f2g(k),10) - cp(k)%clay30 = tmpreal(f2g(k),11) - cp(k)%sand = tmpreal(f2g(k),12) - cp(k)%clay = tmpreal(f2g(k),13) - cp(k)%wpwet30 = tmpreal(f2g(k),14) - cp(k)%poros30 = tmpreal(f2g(k),15) + do k=1,N_tiles_land_f + gid = tile_coord_f(k)%tile_id + cp(k)%gravel30 = tmpreal(gid, 7) + cp(k)%orgC30 = tmpreal(gid, 8) + cp(k)%orgC = tmpreal(gid, 9) + cp(k)%sand30 = tmpreal(gid,10) + cp(k)%clay30 = tmpreal(gid,11) + cp(k)%sand = tmpreal(gid,12) + cp(k)%clay = tmpreal(gid,13) + cp(k)%wpwet30 = tmpreal(gid,14) + cp(k)%poros30 = tmpreal(gid,15) end do case default - do k=1,N_catf + do k=1,N_tiles_land_f cp(k)%gravel30 = nodata_generic cp(k)%orgC30 = nodata_generic @@ -1132,7 +1121,7 @@ subroutine read_cat_param( & tmpreal = nodata_generic - do n=1,N_catg + do n=1,N_tiles_land_g read (10,*) tmptileid(n), dummy_int, (tmpreal(n,m), m=1,4) @@ -1143,11 +1132,11 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - do k=1,N_catf + do k=1,N_tiles_land_f ! this check works only for "SiB2_V2" version - - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + gid = tile_coord_f(k)%tile_id + if (tile_coord_f(k)%tile_id/=tmptileid(gid)) then err_msg = 'something wrong with tau parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if @@ -1156,13 +1145,13 @@ subroutine read_cat_param( & if (abs(dzsf-20.)<1e-4 ) then ! use atau2, btau2 - cp(k)%atau = tmpreal(f2g(k),1) - cp(k)%btau = tmpreal(f2g(k),2) + cp(k)%atau = tmpreal(gid,1) + cp(k)%btau = tmpreal(gid,2) elseif (abs(dzsf-50.)<1e-4 ) then ! use atau5, btau5 - cp(k)%atau = tmpreal(f2g(k),3) - cp(k)%btau = tmpreal(f2g(k),4) + cp(k)%atau = tmpreal(gid,3) + cp(k)%btau = tmpreal(gid,4) else @@ -1200,7 +1189,7 @@ subroutine read_cat_param( & tmpreal = nodata_generic - do n=1,N_catg + do n=1,N_tiles_land_g read (10,*) tmptileid(n), dummy_int, (tmpreal(n,m), m=1,12) @@ -1211,27 +1200,27 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - do k=1,N_catf + do k=1,N_tiles_land_f ! this check works only for "SiB2_V2" version - - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + gid = tile_coord_f(k)%tile_id + if (tile_coord_f(k)%tile_id/=tmptileid(gid)) then err_msg = 'something wrong with ar parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - cp(k)%gnu = tmpreal(f2g(k),1) - cp(k)%ars1 = tmpreal(f2g(k),2) - cp(k)%ars2 = tmpreal(f2g(k),3) - cp(k)%ars3 = tmpreal(f2g(k),4) - cp(k)%ara1 = tmpreal(f2g(k),5) - cp(k)%ara2 = tmpreal(f2g(k),6) - cp(k)%ara3 = tmpreal(f2g(k),7) - cp(k)%ara4 = tmpreal(f2g(k),8) - cp(k)%arw1 = tmpreal(f2g(k),9) - cp(k)%arw2 = tmpreal(f2g(k),10) - cp(k)%arw3 = tmpreal(f2g(k),11) - cp(k)%arw4 = tmpreal(f2g(k),12) + cp(k)%gnu = tmpreal(gid,1) + cp(k)%ars1 = tmpreal(gid,2) + cp(k)%ars2 = tmpreal(gid,3) + cp(k)%ars3 = tmpreal(gid,4) + cp(k)%ara1 = tmpreal(gid,5) + cp(k)%ara2 = tmpreal(gid,6) + cp(k)%ara3 = tmpreal(gid,7) + cp(k)%ara4 = tmpreal(gid,8) + cp(k)%arw1 = tmpreal(gid,9) + cp(k)%arw2 = tmpreal(gid,10) + cp(k)%arw3 = tmpreal(gid,11) + cp(k)%arw4 = tmpreal(gid,12) end do @@ -1248,7 +1237,7 @@ subroutine read_cat_param( & tmpreal = nodata_generic - do n=1,N_catg + do n=1,N_tiles_land_g read (10,*) tmptileid(n), dummy_int, (tmpreal(n,m), m=1,4) @@ -1259,24 +1248,24 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - do k=1,N_catf + do k=1,N_tiles_land_f ! this check works only for "SiB2_V2" version - - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + gid = tile_coord_f(k)%tile_id + if (tile_coord_f(k)%tile_id/=tmptileid(gid)) then err_msg = 'something wrong with bf parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ! --------- - if (cp(k)%gnu/=tmpreal(f2g(k),1)) then + if (cp(k)%gnu/=tmpreal(gid, 1)) then err_msg = 'land(): something wrong with gnu' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - cp(k)%bf1 = tmpreal(f2g(k),2) - cp(k)%bf2 = tmpreal(f2g(k),3) - cp(k)%bf3 = tmpreal(f2g(k),4) + cp(k)%bf1 = tmpreal(gid, 2) + cp(k)%bf2 = tmpreal(gid, 3) + cp(k)%bf3 = tmpreal(gid, 4) end do @@ -1293,7 +1282,7 @@ subroutine read_cat_param( & tmpreal = nodata_generic - do n=1,N_catg + do n=1,N_tiles_land_g read (10,*) tmptileid(n), dummy_int, (tmpreal(n,m), m=1,5) @@ -1303,25 +1292,26 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'done reading' if (logit) write (logunit,*) - do k=1,N_catf + do k=1,N_tiles_land_f ! this check works only for "SiB2_V2" version - if (tile_coord_f(k)%tile_id/=tmptileid(f2g(k))) then + gid = tile_coord_f(k)%tile_id + if (tile_coord_f(k)%tile_id/=tmptileid(gid)) then err_msg = 'something wrong with ts parameters' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ! ------- - if (cp(k)%gnu/=tmpreal(f2g(k),1)) then + if (cp(k)%gnu/=tmpreal(gid,1)) then err_msg = 'land(): something wrong with gnu' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - cp(k)%tsa1 = tmpreal(f2g(k),2) - cp(k)%tsa2 = tmpreal(f2g(k),3) - cp(k)%tsb1 = tmpreal(f2g(k),4) - cp(k)%tsb2 = tmpreal(f2g(k),5) + cp(k)%tsa1 = tmpreal(gid,2) + cp(k)%tsa2 = tmpreal(gid,3) + cp(k)%tsb1 = tmpreal(gid,4) + cp(k)%tsb2 = tmpreal(gid,5) end do @@ -1330,7 +1320,7 @@ subroutine read_cat_param( & if (logit) write (logunit,*) 'computing derived land surface parameters...' if (logit) write (logunit,*) - do k=1,N_catf + do k=1,N_tiles_land_f ! Three soil depths for soil moisture model: ! @@ -1403,9 +1393,9 @@ end subroutine read_cat_param ! ********************************************************************************** - subroutine write_cat_param(cat_param, N_catd) + subroutine write_cat_param(cat_param, N_land_tiles_f) type(cat_param_type), intent(in) :: cat_param(:) - integer,intent(in) :: N_catd + integer,intent(in) :: N_land_tiles_f character(len=512):: fname type(date_time_type) :: start_time @@ -1427,260 +1417,267 @@ subroutine write_cat_param(cat_param, N_catd) print*, 'Writing catparam file : ' // trim(fname) - write (10) (cat_param(n)%dpth, n=1,N_catd) + write (10) (cat_param(n)%dpth, n=1,N_land_tiles_f) - write (10) (cat_param(n)%dzsf, n=1,N_catd) - write (10) (cat_param(n)%dzrz, n=1,N_catd) - write (10) (cat_param(n)%dzpr, n=1,N_catd) + write (10) (cat_param(n)%dzsf, n=1,N_land_tiles_f) + write (10) (cat_param(n)%dzrz, n=1,N_land_tiles_f) + write (10) (cat_param(n)%dzpr, n=1,N_land_tiles_f) do k=1,N_gt - write (10) (cat_param(n)%dzgt(k), n=1,N_catd) + write (10) (cat_param(n)%dzgt(k), n=1,N_land_tiles_f) end do - write (10) (cat_param(n)%poros, n=1,N_catd) - write (10) (cat_param(n)%cond, n=1,N_catd) - write (10) (cat_param(n)%psis, n=1,N_catd) - write (10) (cat_param(n)%bee, n=1,N_catd) - - write (10) (cat_param(n)%wpwet, n=1,N_catd) - - write (10) (cat_param(n)%gnu, n=1,N_catd) - - write (10) (cat_param(n)%vgwmax, n=1,N_catd) - - write (10) (real(cat_param(n)%vegcls), n=1,N_catd) - write (10) (real(cat_param(n)%soilcls30), n=1,N_catd) - write (10) (real(cat_param(n)%soilcls100), n=1,N_catd) - - write (10) (cat_param(n)%bf1, n=1,N_catd) - write (10) (cat_param(n)%bf2, n=1,N_catd) - write (10) (cat_param(n)%bf3, n=1,N_catd) - write (10) (cat_param(n)%cdcr1, n=1,N_catd) - write (10) (cat_param(n)%cdcr2, n=1,N_catd) - write (10) (cat_param(n)%ars1, n=1,N_catd) - write (10) (cat_param(n)%ars2, n=1,N_catd) - write (10) (cat_param(n)%ars3, n=1,N_catd) - write (10) (cat_param(n)%ara1, n=1,N_catd) - write (10) (cat_param(n)%ara2, n=1,N_catd) - write (10) (cat_param(n)%ara3, n=1,N_catd) - write (10) (cat_param(n)%ara4, n=1,N_catd) - write (10) (cat_param(n)%arw1, n=1,N_catd) - write (10) (cat_param(n)%arw2, n=1,N_catd) - write (10) (cat_param(n)%arw3, n=1,N_catd) - write (10) (cat_param(n)%arw4, n=1,N_catd) - write (10) (cat_param(n)%tsa1, n=1,N_catd) - write (10) (cat_param(n)%tsa2, n=1,N_catd) - write (10) (cat_param(n)%tsb1, n=1,N_catd) - write (10) (cat_param(n)%tsb2, n=1,N_catd) - write (10) (cat_param(n)%atau, n=1,N_catd) - write (10) (cat_param(n)%btau, n=1,N_catd) - - write (10) (cat_param(n)%gravel30, n=1,N_catd) - write (10) (cat_param(n)%orgC30 , n=1,N_catd) - write (10) (cat_param(n)%orgC , n=1,N_catd) - write (10) (cat_param(n)%sand30 , n=1,N_catd) - write (10) (cat_param(n)%clay30 , n=1,N_catd) - write (10) (cat_param(n)%sand , n=1,N_catd) - write (10) (cat_param(n)%clay , n=1,N_catd) - write (10) (cat_param(n)%wpwet30 , n=1,N_catd) - write (10) (cat_param(n)%poros30 , n=1,N_catd) - - write (10) (cat_param(n)%veghght , n=1,N_catd) + write (10) (cat_param(n)%poros, n=1,N_land_tiles_f) + write (10) (cat_param(n)%cond, n=1,N_land_tiles_f) + write (10) (cat_param(n)%psis, n=1,N_land_tiles_f) + write (10) (cat_param(n)%bee, n=1,N_land_tiles_f) + + write (10) (cat_param(n)%wpwet, n=1,N_land_tiles_f) + + write (10) (cat_param(n)%gnu, n=1,N_land_tiles_f) + + write (10) (cat_param(n)%vgwmax, n=1,N_land_tiles_f) + + write (10) (real(cat_param(n)%vegcls), n=1,N_land_tiles_f) + write (10) (real(cat_param(n)%soilcls30), n=1,N_land_tiles_f) + write (10) (real(cat_param(n)%soilcls100), n=1,N_land_tiles_f) + + write (10) (cat_param(n)%bf1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%bf2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%bf3, n=1,N_land_tiles_f) + write (10) (cat_param(n)%cdcr1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%cdcr2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ars1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ars2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ars3, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ara1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ara2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ara3, n=1,N_land_tiles_f) + write (10) (cat_param(n)%ara4, n=1,N_land_tiles_f) + write (10) (cat_param(n)%arw1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%arw2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%arw3, n=1,N_land_tiles_f) + write (10) (cat_param(n)%arw4, n=1,N_land_tiles_f) + write (10) (cat_param(n)%tsa1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%tsa2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%tsb1, n=1,N_land_tiles_f) + write (10) (cat_param(n)%tsb2, n=1,N_land_tiles_f) + write (10) (cat_param(n)%atau, n=1,N_land_tiles_f) + write (10) (cat_param(n)%btau, n=1,N_land_tiles_f) + + write (10) (cat_param(n)%gravel30, n=1,N_land_tiles_f) + write (10) (cat_param(n)%orgC30 , n=1,N_land_tiles_f) + write (10) (cat_param(n)%orgC , n=1,N_land_tiles_f) + write (10) (cat_param(n)%sand30 , n=1,N_land_tiles_f) + write (10) (cat_param(n)%clay30 , n=1,N_land_tiles_f) + write (10) (cat_param(n)%sand , n=1,N_land_tiles_f) + write (10) (cat_param(n)%clay , n=1,N_land_tiles_f) + write (10) (cat_param(n)%wpwet30 , n=1,N_land_tiles_f) + write (10) (cat_param(n)%poros30 , n=1,N_land_tiles_f) + + write (10) (cat_param(n)%veghght , n=1,N_land_tiles_f) close (10,status='keep') end subroutine write_cat_param - end subroutine createf2g - - ! ******************************************************************** - - subroutine readsize(f2g_file, N_catg,N_catf) - - implicit none - character(*), intent(in):: f2g_file - integer,intent(out) :: N_catg - integer,intent(out) :: N_catf - - logical :: file_exist - - inquire(file=f2g_file,exist=file_exist) - if(file_exist) then - open(40,file= f2g_file,form='formatted',action='read',status='old') - read(40,*)N_catg - read(40,*)N_catf - close(40) - else - print*, " wrong, no f2g.txt" - endif - end subroutine readsize + end subroutine create_mapping ! ******************************************************************** - subroutine readf2g(f2g_file, N_catf,f2g) + subroutine read_mapping( mapping_file, N_types, tile_types, N_tiles_r, N_tiles_f, f2r, r2g) - implicit none - character(*), intent(in):: f2g_file - integer,intent(in) :: N_catf - integer,dimension(N_catf),intent(inout) :: f2g + character(*), intent(in) :: mapping_file + integer, intent(out) :: N_types + integer, dimension(:), allocatable, optional, intent(out) :: tile_types + integer, dimension(:), allocatable, optional, intent(out) :: N_tiles_r + integer, dimension(:), allocatable, optional, intent(out) :: N_tiles_f + integer, dimension(:), allocatable, optional, intent(out) :: f2r + integer, dimension(:), allocatable, optional, intent(out) :: r2g + + ! ------------------------------ - integer :: N_catg logical :: file_exist - integer :: local_size,n - - inquire(file=f2g_file,exist=file_exist) + integer, dimension(:), allocatable :: N_tiles_r_tmp, N_tiles_f_tmp + + inquire(file=mapping_file,exist=file_exist) if(file_exist) then - open(40,file= f2g_file,form='formatted',action='read',status='old') - read(40,*)N_catg - read(40,*)local_size - - if(local_size /= N_catf) print*, "wrong f2g.txt" - - if(N_catg == N_catf) then - close(40) - return + open(40,file= mapping_file,form='formatted',action='read',status='old') + read(40,*)N_types + allocate(N_tiles_r_tmp(N_types)) + allocate(N_tiles_f_tmp(N_types)) + + if (present(tile_types)) then + allocate(tile_types(N_types)) + read(40,*) tile_types + else + read(40,*) ! read off tile_types + endif + + read(40,*) N_tiles_r_tmp + read(40,*) N_tiles_f_tmp + + if (present(f2r)) then + allocate(f2r(sum(N_tiles_f_tmp))) + read(40,*) f2r + else + read(40,*) ! skip over f2r + endif + if (present(r2g)) then + allocate(r2g(sum(N_tiles_r_tmp))) + read(40,*) r2g + else + read(40,*) ! skip over r2g endif - do n=1,N_catf - read(40,*)f2g(n) - enddo + if (present(N_tiles_r)) N_tiles_r = N_tiles_r_tmp + if (present(N_tiles_f)) N_tiles_f = N_tiles_f_tmp + close(40) - ! call MAPL_sort(this%f2g) else - print*, " wrong, no f2g.txt" + print*, "ERROR, mapping file does not exist: ", trim(mapping_file) endif - end subroutine readf2g + end subroutine read_mapping ! ******************************************************************** - subroutine createLocalTilefile(f2g_file, orig_tile,new_tile) + subroutine createZoominTilefile( mapping_file, orig_tile, new_tile) - implicit none - character(*), intent(in) :: f2g_file + character(*), intent(in) :: mapping_file character(*), intent(in) :: orig_tile character(*), intent(in) :: new_tile + + ! --------------------------------------- character(len=256) :: line - character(len=3) :: MAPL_Land_STRING - character(len=4) :: MAPL_Land_ExcludeFromDomain_STRING - character(len=400) :: err_msg - - logical :: file_exist - - integer, dimension(:),allocatable :: f2g - integer :: N_catg, N_catf,n,stat, ty - integer :: N_tile,N_grid,g_id - - character(len=*), parameter :: Iam = 'createLocalTilefile' - ! string handling below relies on MAPL_Land and MAPL_Land_ExcludeFromDomain - ! falling into a certain range + logical :: file_exist, isNC4 - ! verify that MAPL_Land has three digits - - if (MAPL_Land<100 .or. MAPL_Land>999) then - err_msg = 'string handling implemented only for 100<=MAPL_Land<=999' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - - ! verify that MAPL_Land_ExcludeFromDomain has four digits - - if (MAPL_Land_ExcludeFromDomain<1000 .or. MAPL_Land_ExcludeFromDomain>9999) then - err_msg = 'string handling implemented only for 1000<=MAPL_Land_ExcludeFromDomain<=9999' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - - ! convert integers to appropriate-length strings - - write (MAPL_Land_STRING, '(i3)') MAPL_Land - write (MAPL_Land_ExcludeFromDomain_STRING,'(i4)') MAPL_Land_ExcludeFromDomain + integer, dimension(:),allocatable :: f2g, f2r,r2g, tile_types, N_tiles_r, N_tiles_f + + integer :: n,stat, ty, N_types, nx, ny, file_type, status, N_PfafCat + + integer :: n_tiles, N_grids, g_id, f_id + integer :: IM(2), JM(2) + integer, allocatable :: iTable(:,:) + real(KIND=REAL64), allocatable :: rTable(:,:) + character(len=128) :: Gnames(2) + character(len=4) :: typ_str, typ_str_exclude + character(len=*), parameter :: Iam = 'createZoominTilefile' inquire(file=trim(orig_tile),exist=file_exist) if( .not. file_exist) stop ("original tile file does not exist") - ! Set default local tile file name - call readsize( f2g_file, N_catg,N_catf) - if(N_catg == N_catf) then - print*, "It is global domain..." + ! Set default Zoom in tile file name + call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) + if( all(N_tiles_r == N_tiles_f)) then + print*, "Domain is the same, no need to create tile file..." return endif - allocate(f2g(N_catf)) - call readf2g(f2g_file, N_catf,f2g) - - open(40,file=trim(orig_tile),action="read") - open(50,file=trim(new_tile),action="write") - - ! copy the header back into the output tile file - ! (also corrects bug in EASE *.til files that have "N_grid=1" in line 2 but - ! still contain three additional lines for second grid definition) - do n=1,5 - read(40,'(A)') line - if(n==1) then - read(line,*) N_tile - endif - if(n==2) then - read(line,*) N_grid - endif - write(50,'(A)') trim(line) - enddo - if (N_grid==2) then - do n=1,3 + + g_id = 1 + f_id = 1 + f2g = r2g(f2r) + + call MAPL_NCIOGetFileType(orig_tile, file_type, rc=status) + isNC4 = (file_type == MAPL_FILETYPE_NC4) + + if (isNC4) then + + call MAPL_ReadTilingNC4(orig_tile, GridName=GNames, im=im, jm=jm, nx=nx, ny=ny, n_grids = n_grids, n_tiles=n_tiles, & + iTable=iTable, rTable=rTable, N_PfafCat=N_PfafCat, rc=status) + do n = 1, n_tiles + if (any(tile_types == iTable(n,0))) then + if (f2g(f_id) /= n) then + iTable(n,0) = iTable(n,0) + ExcludeFromDomain + else + f_id = f_id + 1 + if (f_id > size(f2g)) exit + endif + endif + enddo + + call MAPL_WriteTilingNc4(new_tile, gNames(1:n_grids), im(1:n_grids), jm(1:N_grids), & + nx, ny, iTable, rTable, N_PfafCat=N_PfafCat, rc=status) + + else + open(40,file=trim(orig_tile),action="read") + open(50,file=trim(new_tile),action="write") + + ! copy the header back into the output tile file + ! (also corrects bug in EASE *.til files that have "n_grids=1" in line 2 but + ! still contain three additional lines for second grid definition) + do n=1,5 read(40,'(A)') line + if(n==1) then + read(line,*) n_tiles + endif + if(n==2) then + read(line,*) n_grids + endif write(50,'(A)') trim(line) enddo - endif - - g_id = 0 - do while(.true.) - ! read one line of *.til file - read(40,'(A)',IOSTAT=stat) line - if(IS_IOSTAT_END(stat)) exit - ! extract first "integer" in "line" and put into "ty" - read(line,*) ty - if( ty == MAPL_Land ) then - ! find index where MAPL_Land ("100") starts in "line" - n=index(line,MAPL_Land_STRING) - ! make sure that a space is available in front of MAPL_Land ("100") - if (n<=1) then - err_msg = 'string handling requires at least one blank space in first column of *.til file' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - ! here g_id is (consecutive) id of the global *land* tiles - g_id=g_id+1 - if(.not. any( f2g(:) == g_id)) then - ! if tile is not in local domain, replace " 100" in "line" with "1100" - line(n-1:n+2)=MAPL_Land_ExcludeFromDomain_STRING - endif + if (n_grids==2) then + do n=1,3 + read(40,'(A)') line + write(50,'(A)') trim(line) + enddo endif - ! write "line" into the output tile file - write(50,'(A)') trim(line) - enddo - close(40) - close(50) - - end subroutine createLocalTilefile + + do while(.true.) + ! read one line of *.til file + read(40,'(A)',IOSTAT=stat) line + if(IS_IOSTAT_END(stat)) exit + ! extract first "integer" in "line" and put into "ty" + read(line,*) ty + + if( any( tile_types == ty)) then + ! here g_id is (consecutive) id of the global *land* tiles + if (f2g(f_id) /= g_id) then + ! if tile is not in Zoom in domain, replace ty in "line" with 1000+ty" + write(typ_str, '(I0)') ty + typ_str = adjustr(typ_str) + n=index(line, typ_str) + write(typ_str_exclude, '(I0)') ty + ExcludeFromDomain + line(n:n+3) = typ_str_exclude + else + f_id = f_id + 1 + if (f_id > size(f2g)) f_id = 1 ! just set a number to prevent over flow, it would never come back here + endif + endif + ! write "line" into the output tile file + write(50,'(A)') trim(line) + g_id = g_id+1 + enddo + close(40) + close(50) + endif + end subroutine createZoominTilefile ! ******************************************************************** - subroutine createLocalBC(f2g_file, orig_BC, new_BC) + subroutine createZoominBC(mapping_file, orig_BC, new_BC) - implicit none - character(*),intent(in) :: f2g_file + character(*),intent(in) :: mapping_file character(*),intent(in) :: orig_BC character(*),intent(in) :: new_BC + + ! ------------------------------------------- real,dimension(14) :: tmprealvec14 real,allocatable :: tmpvec(:) - integer :: istat, N_catg,N_catf - integer,dimension(:),allocatable :: f2g - - call readsize(f2g_file, N_catg,N_catf) - if(N_catg==N_catf) return - allocate(f2g(N_catf)) - call readf2g(f2g_file, N_catf,f2g) + integer :: istat, N_tiles_land_r, N_tiles_land_f, N_types + integer,dimension(:),allocatable :: f2r_land, f2r, r2g, tile_types, N_tiles_r, N_tiles_f + + call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) + if (tile_types(1) /= MAPL_LAND) return + if (N_tiles_r(1) == N_tiles_f(1) ) return + N_tiles_land_f = N_tiles_f(1) + N_tiles_land_r = N_tiles_r(1) + f2r_land = f2r(1:N_tiles_land_f) - allocate(tmpvec(N_catg)) + allocate(tmpvec(N_tiles_land_r)) open(10,file=trim(orig_BC),form='unformatted',action='read',status='old',iostat=istat) open(20,file=trim(new_BC),form='unformatted',action='write') @@ -1689,25 +1686,26 @@ subroutine createLocalBC(f2g_file, orig_BC, new_BC) if(IS_IOSTAT_END(istat)) exit read(10) tmpvec write(20) tmprealvec14 - write(20) tmpvec(f2g) + write(20) tmpvec(f2r_land) enddo close(10) close(20) deallocate(tmpvec) - end subroutine createLocalBC + end subroutine createZoominBC - ! ******************************************************************** + ! *************************************** - subroutine createLocalCatchRestart(f2g_file, orig_catch, new_catch) - - implicit none - character(*),intent(in):: f2g_file - character(*),intent(in):: orig_catch - character(*),intent(in):: new_catch - integer,parameter :: subtile=4 - integer :: istat, filetype, rc,i, j, ndims + subroutine createZoominRestart(mapping_file, orig_rst, new_rst, tile_type) + + character(*), intent(in) :: mapping_file + character(*), intent(in) :: orig_rst + character(*), intent(in) :: new_rst + integer, intent(in) :: tile_type + + ! ------------------------------------------- + + integer :: istat, file_type, rc,i, j, ndims real,allocatable :: tmp1(:) - real,allocatable :: tmp2(:,:) type(Netcdf4_FileFormatter) :: InFmt,OutFmt type(FileMetadata) :: OutCfg type(FileMetadata) :: InCfg @@ -1717,66 +1715,48 @@ subroutine createLocalCatchRestart(f2g_file, orig_catch, new_catch) type(StringVariableMapIterator) :: var_iter type(StringVector), pointer :: var_dimensions character(len=:), pointer :: vname,dname - integer ::n, N_catg,N_catf - integer,dimension(:),allocatable :: f2g - - call readsize(f2g_file, N_catg,N_catf) - if(N_catg == N_catf) return - allocate(f2g(N_catf)) - call readf2g(f2g_file, N_catf,f2g) - - allocate(tmp1(N_catg)) - allocate(tmp2(N_catg,subtile)) + integer ::n, N_r, N_f, N_types, r_starts, f_starts + integer,dimension(:),allocatable :: f2r, r2g, f2r_, tile_types, N_tiles_r, N_tiles_f + logical :: isNC4 + + call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) + + do n = 1, N_types + if (tile_types(n) == tile_type) then + N_f = N_tiles_f(n) + N_r = N_tiles_r(n) + exit + endif + enddo + + if (N_r == N_f) return + + r_starts = sum(N_tiles_r(1:n-1)) + f_starts = sum(N_tiles_f(1:n-1)) + + f2r_ = f2r( f_starts+1: f_starts+N_f) + f2r_ = f2r_ - r_starts + + allocate(tmp1(N_r)) ! check file type - call MAPL_NCIOGetFileType(orig_catch, filetype,rc=rc) - - if (filetype /= 0) then - - print*, "Catchment restart is binary" - - ! binary - - open(10,file=trim(orig_catch),form='unformatted',action='read',status='old',iostat=istat) - open(20,file=trim(new_catch),form='unformatted',action='write') - - do n=1,30 - read(10) tmp1 - write(20) tmp1(f2g) - enddo - - do n=1,2 - read(10) tmp2 - write(20) tmp2(f2g,:) - enddo - - do n=1,20 - read(10) tmp1 - write(20) tmp1(f2g) - enddo - ! note : the offline restart does not have the last five variables - do n=1,4 - read(10,iostat=istat) tmp2 - if(.not. IS_IOSTAT_END(istat)) write(20) tmp2(f2g,:) - enddo - ! 57 WW - read(10,iostat=istat) tmp2 - if(.not. IS_IOSTAT_END(istat)) write(20) tmp2(f2g,:) - - close(10) - close(20) + call MAPL_NCIOGetFileType(orig_rst, file_type,rc=rc) + isNC4 = (file_type == MAPL_FILETYPE_NC4) + + if ( .not. isNC4 ) then + print*, "Do not support binary restart" else - ! filetype = 0 : nc4 output file will also be nc4 + ! file_type = 0 : nc4 output file will also be nc4 - call InFmt%open(trim(orig_catch), pFIO_READ,rc=rc) + call InFmt%open(trim(orig_rst), pFIO_READ,rc=rc) InCfg = InFmt%read(rc=rc) OutCfg = InCfg - call OutCfg%modify_dimension('tile', size(f2g), rc=rc) + call OutCfg%modify_dimension('tile', size(f2r_), rc=rc) - call OutFmt%create(trim(new_catch),rc=rc) + call OutFmt%create(trim(new_rst),rc=rc) call OutFmt%write(OutCfg,rc=rc) variables => InCfg%get_variables() @@ -1796,14 +1776,14 @@ subroutine createLocalCatchRestart(f2g_file, orig_catch, new_catch) if (ndims == 1) then call MAPL_VarRead (InFmt,vname,tmp1) - call MAPL_VarWrite(OutFmt,vname,tmp1(f2g)) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_)) else if (ndims == 2) then dname => var%get_ith_dimension(2) dim1=InCfg%get_dimension(dname) do j=1,dim1 call MAPL_VarRead ( InFmt,vname,tmp1 ,offset1=j) - call MAPL_VarWrite(OutFmt,vname,tmp1(f2g),offset1=j) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_),offset1=j) enddo else if (ndims == 3) then @@ -1815,7 +1795,7 @@ subroutine createLocalCatchRestart(f2g_file, orig_catch, new_catch) do i=1,dim2 do j=1,dim1 call MAPL_VarRead ( InFmt,vname,tmp1 ,offset1=j,offset2=i) - call MAPL_VarWrite(OutFmt,vname,tmp1(f2g) ,offset1=j,offset2=i) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_) ,offset1=j,offset2=i) enddo enddo @@ -1825,80 +1805,29 @@ subroutine createLocalCatchRestart(f2g_file, orig_catch, new_catch) call inFmt%close(rc=rc) call OutFmt%close(rc=rc) end if ! file type nc4 - print*, "done create local catchment restart" - end subroutine createLocalCatchRestart - + print*, "done create Zoom in restart of type", tile_type + end subroutine createZoominRestart + ! ******************************************************************** - subroutine createLocalmwRTMRestart(f2g_file, orig_mwrtm, new_mwrtm) - - implicit none - character(*),intent(in):: f2g_file - character(*),intent(in):: orig_mwrtm - character(*),intent(in):: new_mwrtm - integer,parameter :: subtile=4 - integer :: rc - real,allocatable :: tmp1(:) - type(Netcdf4_FileFormatter) :: InFmt,OutFmt - type(FileMetadata) :: OutCfg - type(FileMetadata) :: InCfg - - type(StringVariableMap), pointer :: variables - type(StringVariableMapIterator) :: var_iter - character(len=:), pointer :: vname - integer :: N_catg,N_catf - integer,dimension(:),allocatable :: f2g - - call readsize(f2g_file, N_catg,N_catf) - if(N_catg == N_catf) return - allocate(f2g(N_catf)) - call readf2g(f2g_file, N_catf,f2g) - - allocate(tmp1(N_catg)) - - ! nc4 in and out file will also be nc4 - call InFmt%open(trim(orig_mwrtm), pFIO_READ,rc=rc) - InCfg = InFmt%read(rc=rc) - OutCfg = InCfg - - call OutCfg%modify_dimension('tile', size(f2g), rc=rc) + subroutine createZoominVegRestart( mapping_file, orig_veg, new_veg) - call OutFmt%create(trim(new_mwrtm),rc=rc) - call OutFmt%write(OutCfg,rc=rc) - - variables => InCfg%get_variables() - var_iter = variables%begin() - do while (var_iter /= variables%end()) - vname => var_iter%key() - call MAPL_VarRead (InFmt,vname,tmp1) - call MAPL_VarWrite(OutFmt,vname,tmp1(f2g)) - call var_iter%next() - enddo - - call inFmt%close(rc=rc) - call OutFmt%close(rc=rc) - - deallocate(f2g,tmp1) - - end subroutine createLocalmwRTMRestart - - ! ******************************************************************** - - subroutine createLocalVegRestart(f2g_file, orig_veg, new_veg) + character(*), intent(in) :: mapping_file + character(*), intent(in) :: orig_veg + character(*), intent(in) :: new_veg + + ! ------------------------------------------ - implicit none - character(*),intent(in):: f2g_file - character(*),intent(in):: orig_veg - character(*),intent(in):: new_veg integer :: istat real,allocatable :: rity(:) real,allocatable :: z2(:) real,allocatable :: ascatz0(:) real,allocatable :: tmp(:) - integer :: N_catg,N_catf - integer,dimension(:),allocatable :: f2g - integer :: filetype + integer :: N_tiles_land_r, N_tiles_land_f, N_types + integer,dimension(:),allocatable :: f2r, r2g, f2r_land, tile_types, N_tiles_r, N_tiles_f + + integer :: file_type type(Netcdf4_FileFormatter) :: InFmt,OutFmt type(FileMetadata) :: OutCfg type(FileMetadata) :: InCfg @@ -1907,27 +1836,33 @@ subroutine createLocalVegRestart(f2g_file, orig_veg, new_veg) type(StringVariableMapIterator) :: var_iter character(len=:), pointer :: vname integer :: rc + logical :: isNC4 + + call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) + + if (tile_types(1) /= MAPL_LAND) return + if (N_tiles_r(1) == N_tiles_f(1) ) return + N_tiles_land_f = N_tiles_f(1) + N_tiles_land_r = N_tiles_r(1) + + allocate(rity(N_tiles_land_r)) + allocate(z2(N_tiles_land_r)) + allocate(ascatz0(N_tiles_land_r)) - call readsize(f2g_file, N_catg,N_catf) - if(N_catg == N_catf) return - allocate(f2g(N_catf)) - call readf2g(f2g_file, N_catf,f2g) - - allocate(rity(N_catg)) - allocate(z2(N_catg)) - allocate(ascatz0(N_catg)) - - call MAPL_NCIOGetFileType(orig_veg, filetype,rc=rc) - - if (filetype /=0) then + f2r_land = f2r(1:N_tiles_land_f) + + call MAPL_NCIOGetFileType(orig_veg, file_type,rc=rc) + isNC4 = (file_type == MAPL_FILETYPE_NC4) + + if ( .not. isNC4 ) then open(10,file=trim(orig_veg),form='unformatted',action='read',status='old',iostat=istat) open(20,file=trim(new_veg),form='unformatted',action='write') read(10) rity read(10) z2 read(10) ascatz0 - write(20) rity(f2g) - write(20) z2(f2g) - write(20) ascatz0(f2g) + write(20) rity(f2r_land) + write(20) z2(f2r_land) + write(20) ascatz0(f2r_land) close(10) close(20) @@ -1937,18 +1872,18 @@ subroutine createLocalVegRestart(f2g_file, orig_veg, new_veg) InCfg = InFmt%read(rc=rc) OutCfg = InCfg - call OutCfg%modify_dimension('tile', size(f2g), rc=rc) + call OutCfg%modify_dimension('tile', N_tiles_land_f, rc=rc) call OutFmt%create(trim(new_veg),rc=rc) call OutFmt%write(OutCfg,rc=rc) variables => InCfg%get_variables() var_iter = variables%begin() - allocate(tmp(N_catg)) + allocate(tmp(N_tiles_land_r)) do while (var_iter /= variables%end()) vname => var_iter%key() call MAPL_VarRead (InFmt,vname,tmp) - call MAPL_VarWrite(OutFmt,vname,tmp(f2g)) + call MAPL_VarWrite(OutFmt,vname,tmp(f2r_land)) call var_iter%next() enddo @@ -1956,9 +1891,8 @@ subroutine createLocalVegRestart(f2g_file, orig_veg, new_veg) call OutFmt%close(rc=rc) deallocate(tmp) endif - deallocate(f2g) - end subroutine createLocalVegRestart + end subroutine createZoominVegRestart ! ******************************************************************** @@ -1966,7 +1900,7 @@ subroutine correctEase(orig_ease,new_ease) ! This subroutine corrects for a bug that is present in some EASE *.til files ! through at least Icarus-NLv4. - ! Affected files state "N_grid=1" in line 2 of the header, but the header still includes + ! Affected files state "n_grids=1" in line 2 of the header, but the header still includes ! three additional lines for a second grid, which throws off the canonical *.til reader ! (subroutine LDAS_read_til_file()). ! @@ -1975,19 +1909,27 @@ subroutine correctEase(orig_ease,new_ease) ! ! - reichle, 2 Aug 2020 - implicit none character(*),intent(in) :: orig_ease character(*),intent(in) :: new_ease - logical :: file_exist,is_oldEASE - integer :: i, N_tile, N_grid + + ! ------------------------------------- + + logical :: file_exist,is_oldEASE, isNC4 + integer :: i, n_tiles, n_grids, file_type, status character(len=256) :: tmpline inquire(file=trim(orig_ease),exist=file_exist) if( .not. file_exist) stop (" no ease_tile_file") - + call MAPL_NCIOGetFileType(orig_ease, file_type, rc=status) + isNC4 = (file_type == MAPL_FILETYPE_NC4) + + if (isNC4) then + print*, "isNC4 tile file, no need to be corrected" + return + endif open(55,file=trim(orig_ease),action='read') - read(55,*) N_tile - read(55,*) N_grid + read(55,*) n_tiles + read(55,*) n_grids read(55,*) read(55,*) read(55,*) @@ -1995,7 +1937,7 @@ subroutine correctEase(orig_ease,new_ease) close(55) is_oldEASE= .false. - if(N_grid==1 .and. index(tmpline,'OCEAN')/=0) is_oldEASE=.true. + if(n_grids==1 .and. index(tmpline,'OCEAN')/=0) is_oldEASE=.true. if( is_oldEASE) then open(55,file=trim(orig_ease),action='read') @@ -2007,13 +1949,14 @@ subroutine correctEase(orig_ease,new_ease) read(55,*) read(55,*) read(55,*) - do i=1,N_tile + do i=1,n_tiles read(55,'(A)')tmpline write(56,'(A)')trim(tmpline) enddo close(56) close(55) end if + end subroutine correctEase ! ******************************************************************** @@ -2042,35 +1985,34 @@ end subroutine correctEase ! NY: N_proc 1 ! JMS.rc IMS.rc - subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_dir) + subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_dir, types) - implicit none - - character(*), intent(in) :: fname_tilefile ! file name (with path) of tile file (*.til) - character(*), intent(in) :: N_proc_string ! *string* w/ no. of processors (or tasks), excl. OSERVER tasks - character(*), intent(in) :: optimized_file - character(*), intent(in) :: run_dir - + character(*), intent(in) :: fname_tilefile ! file name (with path) of tile file (*.til) + character(*), intent(in) :: N_proc_string ! *string* w/ no. of processors (or tasks), excl. OSERVER tasks + character(*), intent(in) :: optimized_file + character(*), intent(in) :: run_dir + integer, optional, intent(in) :: types(:) + ! local variables integer :: N_proc - integer :: N_tile,N_lon,N_lat,N_grid - integer,allocatable :: landPosition(:) - integer,allocatable :: IMS(:),JMS(:) - integer,allocatable :: local_land(:) - integer :: total_land + integer :: n_tiles,N_lon,N_lat,n_grids + integer,allocatable :: tilePosition(:) + integer,allocatable :: IMS(:),JMS(:), typs(:), II(:), JJ(:) + integer,allocatable :: local_tile(:) + integer :: total_tile integer :: n,typ,tmpint real :: tmpreal - integer :: avg_land,n0,local - integer :: i,s,e,j,k,n1,n2, s1, s2 + integer :: avg_tile,n0,local + integer :: i,s,e,j,k,n1,n2, s1, s2, IM(2), JM(2) logical :: file_exist character(len=256):: tmpLine - character(len=128):: gridname + character(len=128):: gridname, gNames(2) real :: rate,rates(60),maxf(60) - integer :: IMGLOB, JMGLOB + integer :: IMGLOB, JMGLOB, file_type, status integer :: face(6),face_land(6) - logical :: forward + logical :: forward, isNC4 character(len=:), allocatable :: IMS_file, JMS_File - + integer, allocatable :: tile_types(:), iTable(:,:) character(len=*), parameter :: Iam = 'optimize_latlon' character(len=400) :: err_msg @@ -2079,20 +2021,39 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di read (N_proc_string,*) N_proc ! input is string for historical reasons... ! get tile info + if (present(types)) then + tile_types = types + else + tile_types = [MAPL_LAND] + endif inquire(file=trim(fname_tilefile),exist=file_exist) if( .not. file_exist) then - err_msg = 'tile file does not exist' + err_msg = 'tile file does not exist: ' //trim(fname_tilefile) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if + + call MAPL_NCIOGetFileType(trim(fname_tilefile), file_type, rc=status) + isNC4 = (file_type == MAPL_FILETYPE_NC4) - open (10, file=trim(fname_tilefile), form='formatted', action='read') - read (10,*) N_tile - read (10,*) N_grid ! some number (?) - read (10,*) gridname ! some string describing tile definition grid (?) - read (10,*) N_lon - read (10,*) n_lat - + if (isNC4) then + call MAPL_ReadTilingNC4(trim(fname_tilefile), GridName=GNames, im=im, jm=jm, n_Grids=n_grids, n_tiles=n_tiles, & + iTable=iTable, rc=status) + gridname = GNames(1) + n_lon = IM(1) + n_lat = JM(1) + + else + open (10, file=trim(fname_tilefile), form='formatted', action='read') + read (10,*) n_tiles + read (10,*) n_grids ! some number (?) + read (10,*) gridname ! some string describing tile definition grid (?) + read (10,*) N_lon + read (10,*) n_lat + endif + + allocate(typs(n_tiles), II(n_tiles), JJ(n_tiles)) + if (index(gridname,"CF") /=0) then ! cube-sphere tile space IMGLOB = N_lon ! e.g., 180 for c180 @@ -2102,69 +2063,59 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - allocate(landPosition(JMGLOB)) - landPosition = 0 - total_land = 0 - - if(N_grid==2) then - read (10,*) ! some string describing ocean grid (?) - read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) - read (10,*) - endif - - do n = 1,N_tile - read (10,*) & - typ, & ! 1 + allocate(tilePosition(JMGLOB)) + tilePosition = 0 + total_tile = 0 + + + if (.not. isNC4) then + if(n_grids==2) then + read (10,*) ! some string describing ocean grid (?) + read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) + read (10,*) + endif + do n = 1, n_tiles + read (10,*) & + typs(n), & ! 1 tmpreal, & ! 2 * tmpreal, & ! 3 tmpreal, & ! 4 - i , & ! 5 - j ! 6 - !tmpreal, & ! 7 - !tmpint, & ! 8 - !tmpreal, & ! 9 * - !tmpint, & ! 10 - !tmpreal, & ! 11 - !tmpint ! 12 * (previously "tile_id") - - if(typ==MAPL_Land) then - total_land=total_land+1 - landPosition(j) = landPosition(j)+1 + ii(n) , & ! 5 + jj(n) ! 5 + enddo + close(10) + else + typs = iTable(:,0) + II = iTable(:,2) + JJ = iTable(:,3) + endif + + do n = 1,n_tiles + typ = typs(n) + i = II(n) + j = jj(n) + + if(any (tile_types == typ)) then + total_tile=total_tile+1 + tilePosition(j) = tilePosition(j)+1 endif - ! assume all land tiles are at the beginning - ! UNSAFE ASSUMPTION! - reichle, 2 Aug 2020 - - if (typ/=MAPL_Land .and. typ/=MAPL_Land_ExcludeFromDomain) then ! exit if not land - - if (logit) then - write (logunit,*) 'WARNING: Encountered first non-land tile in *.til file.' - write (logunit,*) ' Stop reading *.til file under the assumption that' - write (logunit,*) ' land tiles are first in *.til file.' - write (logunit,*) ' This is NOT a safe assumption beyond Icarus-NLv[x] tile spaces!!' - end if - - exit ! assuming land comes first in the til file - - end if - enddo - close(10) if(mod(N_proc,6) /=0) then print*,"WARNING: ntasks should be adjusted to multiple of 6 for cubed-sphere grid :",N_proc N_proc = N_proc-mod(N_proc,6) endif - print*, "total tiles: ", total_land + print*, "total tiles: ", total_tile - if(sum(landPosition) /= total_land) print*, "wrong counting of land" + if(sum(tilePosition) /= total_tile) print*, "wrong counting of land" do k=1,6 n1 = (k-1)*IMGLOB+1 n2 = k*IMGLOB - face_land(k) = sum(landPosition(n1:n2)) - face(k) = nint(1.0*face_land(k)/total_land * N_proc) + face_land(k) = sum(tilePosition(n1:n2)) + face(k) = nint(1.0*face_land(k)/total_tile * N_proc) ! ensure each face has at least 1 process if ( face(k) == 0) face(k) = 1 ! ensure that the stripe for each process can be at least 2 cells wide @@ -2200,7 +2151,7 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di n2 = k*IMGLOB s1 = sum(face(1:k-1)) + 1 s2 = sum(face(1:k)) - call equal_partition(landPosition(n1:n2), JMS(s1:s2)) + call equal_partition(tilePosition(n1:n2), JMS(s1:s2)) enddo if( sum(JMS) /= JMGLOB) then @@ -2235,15 +2186,15 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) endif - allocate(local_land(N_Proc), source = 0) + allocate(local_tile(N_Proc), source = 0) s1 = 1 do n = 1, N_proc s2 = s1+JMS(n) - 1 - local_land(n) = sum(landPosition(s1:s2)) + local_tile(n) = sum(tilePosition(s1:s2)) s1 = s2 + 1 enddo - print*,"land_distribute: ",local_land + print*,"tile_distribute: ",local_tile print*, "JMS.rc", JMS if( sum(JMS) /= JMGLOB) then print*, sum(JMS), JMGLOB @@ -2289,13 +2240,13 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di ! *not* cube-sphere tile space allocate(IMS(N_Proc)) - allocate(local_land(N_Proc)) + allocate(local_tile(N_Proc)) IMS=0 - local_land = 0 + local_tile = 0 ! NOTE: ! There is a bug in at least some EASE *.til files through at least Icarus-NLv4. - ! Affected files state "N_grid=1" in line 2 of the header, but the header still includes + ! Affected files state "n_grids=1" in line 2 of the header, but the header still includes ! three additional lines for a second grid. ! ! The "else" block below corrects for this bug. @@ -2305,20 +2256,42 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di ! ldas_setup. ! ! -reichle, 2 Aug 2020 - - if(N_grid==2) then - read (10,*) ! some string describing ocean grid (?) - read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) - read (10,*) - read(10,'(A)') tmpLine + + if (isNC4) then + + typs = iTable(:,0) + II = iTable(:,2) + JJ = iTable(:,3) else - read(10,'(A)') tmpLine - if (index(tmpLine,"OCEAN") /=0) then + if(n_grids==2) then read (10,*) ! some string describing ocean grid (?) read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) read (10,*) read(10,'(A)') tmpLine + else + read(10,'(A)') tmpLine + if (index(tmpLine,"OCEAN") /=0) then + read (10,*) ! some string describing ocean grid (?) + read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) + read (10,*) + read(10,'(A)') tmpLine + endif endif + read (tmpLine,*) & + typs(1), & ! 1 + tmpreal, & ! 2 * + tmpreal, & ! 3 + tmpreal, & ! 4 + ii(1) + do n=2, n_tiles + read (10,*) & + typs(n), & ! 1 + tmpreal, & ! 2 * + tmpreal, & ! 3 + tmpreal, & ! 4 + ii(n) ! 5 + enddo + close(10) endif if (index(gridname,'EASE') /=0) then @@ -2328,64 +2301,26 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di s=1 e=N_lon endif - allocate(landPosition(s:e)) + allocate(tilePosition(s:e)) - landPosition=0 - total_land= 0 + tilePosition=0 + total_tile= 0 ! 1) read through tile file, put the land tile into the N_lon of bucket - read (tmpLine,*) & - typ, & ! 1 - tmpreal, & ! 2 * - tmpreal, & ! 3 - tmpreal, & ! 4 - i ! 5 - if(typ==MAPL_Land) then - total_land=total_land+1 - landPosition(i) = landPosition(i)+1 - endif - - do n = 2,N_tile - read (10,*) & - typ, & ! 1 - tmpreal, & ! 2 * - tmpreal, & ! 3 - tmpreal, & ! 4 - i ! 5 - !tmpint, & ! 6 - !tmpreal, & ! 7 - !tmpint, & ! 8 - !tmpreal, & ! 9 * - !tmpint, & ! 10 - !tmpreal, & ! 11 - !tmpint ! 12 * (previously "tile_id") - if(typ==MAPL_Land) then - total_land=total_land+1 - landPosition(i) = landPosition(i)+1 + do n = 1,n_tiles + typ = typs(n) + i = II(n) + + if (any(tile_types == typ)) then + total_tile=total_tile+1 + tilePosition(i) = tilePosition(i)+1 endif - ! assume all land tiles are at the beginning - ! UNSAFE ASSUMPTION! - reichle, 2 Aug 2020 - - if (typ/=MAPL_Land .and. typ/=MAPL_Land_ExcludeFromDomain) then ! exit if not land - - if (logit) then - write (logunit,*) 'WARNING: Encountered first non-land tile in *.til file.' - write (logunit,*) ' Stop reading *.til file under the assumption that' - write (logunit,*) ' land tiles are first in *.til file.' - write (logunit,*) ' This is NOT a safe assumption beyond Icarus-NLv[x] tile spaces!!' - end if - - exit ! assuming land comes first in the til file - - end if - enddo - close(10) - if(sum(landPosition) /= total_land) print*, "wrong counting of land" + if(sum(tilePosition) /= total_tile) print*, "wrong counting of land" do n=1,60 rates(n) = -0.3 + (n-1)*0.01 @@ -2397,39 +2332,39 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di ! 2) each process should have average land tiles - avg_land = ceiling(1.0*total_land/N_proc) - print*,"avg_land",avg_land + avg_tile = ceiling(1.0*total_tile/N_proc) + print*,"avg_tile",avg_tile - ! rate is used to readjust the avg_land + ! rate is used to readjust the avg_tile ! in case that the last processors don't have any land tiles, ! we can increase ther rates - avg_land = avg_land - nint(rate*avg_land) - print*,"re adjust the avg_land",avg_land + avg_tile = avg_tile - nint(rate*avg_tile) + print*,"re adjust the avg_tile",avg_tile tmpint = 0 local = 1 n0 = s-1 forward = .true. do n=s,e - tmpint=tmpint+landPosition(n) + tmpint=tmpint+tilePosition(n) if(local == N_proc .and. n < e) cycle ! all lefteover goes to the last process if( n==e ) then - local_land(local)=tmpint + local_tile(local)=tmpint IMS(local)=n-n0 exit endif - if( tmpint .ge. avg_land ) then + if( tmpint .ge. avg_tile ) then if (forward .or. n-n0 == 1 ) then - local_land(local)=tmpint + local_tile(local)=tmpint IMS(local)=n-n0 tmpint=0 n0=n forward = .false. else - local_land(local) = tmpint - landPosition(n) + local_tile(local) = tmpint - tilePosition(n) IMS(local)=(n-1)-n0 - tmpint= landPosition(n) + tmpint= tilePosition(n) n0 = n-1 forward = .true. endif @@ -2438,9 +2373,9 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di enddo print*,"rms rate: ", rms(rate) - print*,"land_distribute: ",local_land + print*,"tile_distribute: ",local_tile - if( sum(local_land) /= total_land) then + if( sum(local_tile) /= total_tile) then err_msg = 'wrong distribution' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if @@ -2497,40 +2432,42 @@ subroutine optimize_latlon(fname_tilefile, N_proc_string, optimized_file, run_di ! --------------------------------------------------- elemental function rms(rates) result (f) - real :: f - real,intent(in) :: rates + + real :: f + real, intent(in) :: rates + integer :: tmpint,local integer :: n0,proc,n - integer :: avg_land - integer,allocatable :: local_land(:) + integer :: avg_tile + integer,allocatable :: local_tile(:) logical :: forward - allocate (local_land(N_proc)) - local_land = 0 - avg_land = ceiling(1.0*total_land/N_proc) - avg_land = avg_land -nint(rates*avg_land) + allocate (local_tile(N_proc)) + local_tile = 0 + avg_tile = ceiling(1.0*total_tile/N_proc) + avg_tile = avg_tile -nint(rates*avg_tile) forward = .true. tmpint = 0 local = 1 n0 = s-1 do n=s,e - tmpint=tmpint+landPosition(n) + tmpint=tmpint+tilePosition(n) if(local == N_proc .and. n < e) cycle ! all lefteover goes to the last process if( n==e ) then - local_land(local)=tmpint + local_tile(local)=tmpint exit endif - if( tmpint .ge. avg_land ) then + if( tmpint .ge. avg_tile ) then if (forward .or. n-n0 == 1 ) then - local_land(local)=tmpint + local_tile(local)=tmpint tmpint=0 n0=n forward = .false. else - local_land(local) = tmpint - landPosition(n) - tmpint= landPosition(n) + local_tile(local) = tmpint - tilePosition(n) + tmpint= tilePosition(n) n0 = n-1 forward = .true. endif @@ -2539,16 +2476,18 @@ elemental function rms(rates) result (f) enddo f = 0.0 do proc = 1, N_proc - f =max(f,1.0*abs(local_land(proc)-avg_land)) + f =max(f,1.0*abs(local_tile(proc)-avg_tile)) enddo - deallocate(local_land) + deallocate(local_tile) end function rms ! --------------------------------------------------- subroutine equal_partition(array, distribute) + integer, intent(in) :: array(:) integer, intent(inout) :: distribute(:) + integer, allocatable :: ArraySum(:) integer, allocatable :: table(:,:), partition(:,:) integer :: n, k, tmp_max @@ -2623,13 +2562,12 @@ end subroutine optimize_latlon subroutine convert_pert_rst(pfile_name,pfile_nc4,in_path,exp_id) - implicit none character(*),intent(in) :: pfile_name character(*),intent(in) :: in_path character(*),intent(in) :: exp_id character(*),intent(in) :: pfile_nc4 - integer :: N_catf,N_lon,N_lat,N_lonf,N_latf + integer :: N_tiles_land_f,N_lon,N_lat,N_lonf,N_latf integer :: N_force_pert,N_progn_pert integer,pointer :: f2g(:) @@ -2642,7 +2580,7 @@ subroutine convert_pert_rst(pfile_name,pfile_nc4,in_path,exp_id) real,allocatable :: Force_pert_ntrmdt_f(:,:,:) real,allocatable :: Progn_pert_ntrmdt_f(:,:,:) - call io_domain_files('r',in_path, trim(exp_id),N_catf,f2g,tile_coord_f,pert_grid_g,pert_grid_f,RC) + call io_domain_files('r',in_path, trim(exp_id),N_tiles_land_f,f2g,tile_coord_f,pert_grid_g,pert_grid_f,RC) N_lon = pert_grid_g%N_lon N_lat = pert_grid_g%N_lat @@ -2658,6 +2596,7 @@ subroutine convert_pert_rst(pfile_name,pfile_nc4,in_path,exp_id) ! *************************************************************************** subroutine i_pert_ldas(rc) + integer,intent(inout),optional :: rc integer :: nrandseed_tmp @@ -2736,7 +2675,9 @@ end subroutine i_pert_ldas ! ******************************************************************** subroutine o_pert_GEOSldas(rc) + integer,intent(inout) :: rc + integer :: NCFOutID, STATUS integer :: seeddim,latdim, londim, Nforce,NProgn integer :: dims(3), seedid,forceid,prognid @@ -2828,60 +2769,70 @@ end subroutine convert_pert_rst ! ************************************************************************************************** - subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_land, f2g ) + subroutine LDAS_read_til_file( tile_file, catch_def_file, tile_grid_g, tile_coord_r, r2g, N_tiles_land_g, types ) - ! read land tile information from *.til file + ! read tile information from *.til file ! ! This subroutine: ! - is the GEOSldas version of the LDASsa subroutine read_til_file() and ! - was known as LDAS_read_land_tile() when in LDAS_TileCoordRoutines.F90. ! ! inputs: - ! tile_file : *.til tile definition file (full path + name) - ! catch_file : catchment.def file (full path + name) + ! tile_file : *.til tile definition file (full path + name) + ! catch_def_file : catchment.def file (full path + name) ! ! outputs: - ! tile_grid_g : parameters of tile definition grid - ! tile_coord_land : coordinates of tiles (see tile_coord_type), - ! implemented as pointer which is allocated in - ! this subroutine - ! NOTE: number of *land* tiles can be diagnosed with size(tile_coord) - ! optional: - ! f2g : the full domain id to the global id + ! tile_grid_g : parameters of tile definition grid + ! tile_coord_r : coordinates of tiles (see tile_coord_type), + ! implemented as pointer which is allocated in + ! this subroutine ! + ! optional: + ! r2g : the restart domain id to the global id + ! N_tiles_land_g : Number of *land* tiles in global domain + ! "tile_id" is no longer read from *.til file and is now set in this ! subroutine to match order of tiles in *.til file ! - reichle, 22 Aug 2013 ! - ! improved documentation of bug in some EASE *.til files (header says N_grid=1 but has two grid defs) + ! improved documentation of bug in some EASE *.til files (header says n_grids=1 but has two grid defs) ! and minor clean-up ! - reichle, 2 Aug 2020 ! + ! modified to accommodate additional tile types (landice) + ! - wjiang, reichle, 30 Apr 2025 + ! ! ------------------------------------------------------------- - implicit none - character(*), intent(in) :: tile_file - character(*), intent(in) :: catch_file + character(*), intent(in) :: catch_def_file type(grid_def_type), intent(inout):: tile_grid_g - type(tile_coord_type), dimension(:), pointer :: tile_coord_land ! out - integer, dimension(:), optional, pointer :: f2g ! out + type(tile_coord_type), dimension(:), pointer :: tile_coord_r ! out + integer, intent(out) :: N_tiles_land_g + integer, dimension(:), pointer :: r2g ! out + integer, dimension(:), optional :: types ! input ! locals type(tile_coord_type), dimension(:), allocatable :: tile_coord - integer, dimension(:), allocatable :: f2g_tmp ! out + integer, dimension(:), allocatable :: r2g_tmp ! out real :: ease_cell_area - integer :: i, N_tile, N_grid,tmpint1, tmpint2, tmpint3, tmpint4 + integer :: i, n_tiles, n_grids,tmpint1, tmpint2, tmpint3, tmpint4 integer :: i_indg_offset, j_indg_offset, col_order - integer :: N_tile_land, n_lon, n_lat - logical :: ease_grid - integer :: typ,k,fid - + integer :: n_lon, n_lat + logical :: ease_grid, isNC4 + integer :: typ, k, file_type, status + integer, dimension(:), allocatable :: tile_types + character(256) :: tmpline character(128) :: gridname character(512) :: fname + character(128) :: GNames(2) + integer :: IM(2), JM(2) + integer, allocatable :: iTable(:,:) + real(KIND=REAL64), allocatable :: rTable(:,:) + character(len=*), parameter :: Iam = 'LDAS_read_til_file' ! --------------------------------------------------------------- @@ -2889,231 +2840,284 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la i_indg_offset = 0 j_indg_offset = 0 + if (present(types)) then + tile_types = types + else + tile_types =[MAPL_LAND] ! default is to include only land tiles + endif + ! read *.til file header if (logit) write (logunit,'(400A)') trim(Iam), '(): reading from ' // trim(tile_file) - - open (10, file=trim(tile_file), form='formatted', action='read') - - read (10,*) N_tile ! number of all tiles in *.til file, incl non-land types - read (10,*) N_grid - read (10,*) gridname - read (10,*) n_lon - read (10,*) n_lat - - ! NOTE: - ! There is a bug in at least some EASE *.til files through at least Icarus-NLv4. - ! Affected files state "N_grid=1" in line 2 of the header, but the header still includes - ! three additional lines for a second grid. - ! LDAS pre-processing corrects for this bug through subroutine correctEase() in - ! preprocess_LDAS.F90, which creates a second, corrected version of the *.til file during - ! ldas_setup. Here, this corrected *.til file is read! - - if(N_grid==2) then - read (10,*) ! some string describing ocean grid (?) - read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) - read (10,*) ! # ocean grid cells in latitude direction (N_j_ocn) (?) - endif ease_grid = .false. col_order = 0 + + call MAPL_NCIOGetFileType(tile_file, file_type, rc=status) + isNC4 = (file_type == MAPL_FILETYPE_NC4) + + if (isNC4) then + + call MAPL_ReadTilingNC4(tile_file, GridName=GNames, im=im, jm=jm, n_Grids=n_grids, n_tiles=n_tiles, & + iTable=iTable, rTable=rTable, rc=status) + gridname = GNames(1) + n_lon = IM(1) + n_lat = JM(1) + + i =0 + allocate(r2g_tmp(n_tiles)) + do k = 1, n_tiles + if (any(tile_types == iTable(k,0))) then + i = i +1 + r2g_tmp(i) = k + endif + enddo + allocate(r2g, source = r2g_tmp(1:i)) + allocate(tile_coord_r(i)) + deallocate(r2g_tmp) + else - call LDAS_create_grid_g( gridname, n_lon, n_lat, & - tile_grid_g, i_indg_offset, j_indg_offset, ease_cell_area ) + open (10, file=trim(tile_file), form='formatted', action='read') + read (10,*) n_tiles ! number of all tiles in *.til file, incl non-land types + read (10,*) n_grids + read (10,*) gridname + read (10,*) n_lon + read (10,*) n_lat + + ! NOTE: + ! There is a bug in at least some EASE *.til files through at least Icarus-NLv4. + ! Affected files state "n_grids=1" in line 2 of the header, but the header still includes + ! three additional lines for a second grid. + ! LDAS pre-processing corrects for this bug through subroutine correctEase() in + ! preprocess_LDAS.F90, which creates a second, corrected version of the *.til file during + ! ldas_setup. Here, this corrected *.til file is read! + + if(n_grids==2) then + read (10,*) ! some string describing ocean grid (?) + read (10,*) ! # ocean grid cells in longitude direction (N_i_ocn) (?) + read (10,*) ! # ocean grid cells in latitude direction (N_j_ocn) (?) + endif + + endif + + call LDAS_create_grid_g( gridname, n_lon, n_lat, & + tile_grid_g, i_indg_offset, j_indg_offset, ease_cell_area ) + if (index(tile_grid_g%gridtype,'EASE')/=0) ease_grid = .true. ! 'EASEv1' or 'EASEv2' if (index(tile_grid_g%gridtype,'SiB2')/=0) col_order=1 ! old bcs - allocate(tile_coord(N_tile)) - allocate(f2g_tmp(N_tile)) - - i = 0 - fid = 0 - - ! WJ notes: i and k are the same---global ids - ! fid --- num in simulation domain - - do k=1,N_tile - - read(10,'(A)') tmpline - read(tmpline,*) typ - - ! tile type "MAPL_Land_ExcludeFromDomain" identifies land tiles to exclude - ! when non-global domain is created - - if (typ==MAPL_Land .or. typ==MAPL_Land_ExcludeFromDomain) then ! all land - i=i+1 - tile_coord(i)%tile_id = k - - ! now keep only tiles that are not excluded by way of MAPL_Land_ExcludeFromDomain - - if (typ==MAPL_Land) then - fid=fid+1 - f2g_tmp(fid) = k - end if - - ! Not sure ".or. N_grid==1" will always work in the following conditional. - ! Some Tripolar grid *.til files may have N_grid=1. - ! - reichle, 2 Aug 2020 + if (isNC4 ) then + N_tiles_land_g = count(iTable(:,0) == MAPL_LAND .or. iTable(:,0) == (MAPL_LAND + ExcludeFromDomain)) + tile_coord_r(:)%typ = iTable(r2g, 0) + tile_coord_r(:)%i_indg = iTable(r2g, 2) + tile_coord_r(:)%j_indg = iTable(r2g, 3) + + tile_coord_r(:)%i_indg = tile_coord_r(:)%i_indg + i_indg_offset + tile_coord_r(:)%j_indg = tile_coord_r(:)%j_indg + j_indg_offset + + tile_coord_r(:)%tile_id = r2g + tile_coord_r(:)%pfaf_index= iTable(r2g, 4) + tile_coord_r(:)%com_lon = rTable(r2g, 1) + tile_coord_r(:)%com_lat = rTable(r2g, 2) + tile_coord_r(:)%area = rTable(r2g, 3) + tile_coord_r(:)%frac_cell = rTable(r2g, 4) + tile_coord_r(:)%min_lon = rTable(r2g, 6) + tile_coord_r(:)%max_lon = rTable(r2g, 7) + tile_coord_r(:)%min_lat = rTable(r2g, 8) + tile_coord_r(:)%max_lat = rTable(r2g, 9) + tile_coord_r(:)%elev = rTable(r2g, 10) + + tile_coord_r%frac_pfaf = nodata_generic + tile_coord_r%pert_i_indg = nint(nodata_generic) + tile_coord_r%pert_j_indg = nint(nodata_generic) + else + allocate(tile_coord(n_tiles)) + allocate(r2g_tmp(n_tiles)) + + i = 0 + + ! WJ notes: i and k are the same---global ids + ! fid --- num in simulation domain + N_tiles_land_g = 0 + do k=1,n_tiles - if (ease_grid .or. N_grid==1) then - - ! EASE grid til file has fewer columns - ! (excludes "tile_id", "frac_pfaf", and "area") - - read (tmpline,*) & - tile_coord(i)%typ, & ! 1 - tile_coord(i)%pfaf_index, & ! 2 - tile_coord(i)%com_lon, & ! 3 - tile_coord(i)%com_lat, & ! 4 - tile_coord(i)%i_indg, & ! 5 - tile_coord(i)%j_indg, & ! 6 - tile_coord(i)%frac_cell ! 7 - - tile_coord(i)%frac_pfaf = nodata_generic - - ! compute area of tile in [km^2] (units convention in tile_coord structure) - - tile_coord(i)%area = ease_cell_area*tile_coord(i)%frac_cell/1000./1000. ! [km^2] - - else ! not ease grid + read(10,'(A)') tmpline + read(tmpline,*) typ + + if (typ == MAPL_LAND .or. typ == MAPL_LAND + ExcludeFromDomain) N_tiles_land_g = N_tiles_land_g + 1 + + ! adding "ExcludeFromDomain" to tile type identifies tiles to be excluded + ! when non-global domain is created + + if (any( tile_types == typ)) then ! all needed tiles + + i=i+1 + tile_coord(i)%tile_id = k + r2g_tmp(i) = k + + ! Not sure ".or. n_grids==1" will always work in the following conditional. + ! Some Tripolar grid *.til files may have n_grids=1. + ! - reichle, 2 Aug 2020 - if (col_order==1) then + if (ease_grid .or. n_grids==1) then - ! old "SiB2_V2" file format - - read (tmpline,*) & - tile_coord(i)%typ, & ! 1 - tile_coord(i)%pfaf_index, & ! 2 * - tile_coord(i)%com_lon, & ! 3 - tile_coord(i)%com_lat, & ! 4 - tile_coord(i)%i_indg, & ! 5 - tile_coord(i)%j_indg, & ! 6 - tile_coord(i)%frac_cell, & ! 7 - tmpint1, & ! 8 - tmpint2, & ! 9 * - tmpint3, & ! 10 - tile_coord(i)%frac_pfaf, & ! 11 - tmpint4, & ! 12 (previously "tile_id") - tile_coord(i)%area ! 13 + ! EASE grid til file has fewer columns + ! (excludes "tile_id", "frac_pfaf", and "area") - else + read (tmpline,*) & + tile_coord(i)%typ, & ! 1 + tile_coord(i)%pfaf_index, & ! 2 + tile_coord(i)%com_lon, & ! 3 + tile_coord(i)%com_lat, & ! 4 + tile_coord(i)%i_indg, & ! 5 + tile_coord(i)%j_indg, & ! 6 + tile_coord(i)%frac_cell ! 7 - read (tmpline,*) & - tile_coord(i)%typ, & ! 1 - tile_coord(i)%area, & ! 2 * - tile_coord(i)%com_lon, & ! 3 - tile_coord(i)%com_lat, & ! 4 - tile_coord(i)%i_indg, & ! 5 - tile_coord(i)%j_indg, & ! 6 - tile_coord(i)%frac_cell, & ! 7 - tmpint1, & ! 8 - tile_coord(i)%pfaf_index, & ! 9 * - tmpint2, & ! 10 - tile_coord(i)%frac_pfaf, & ! 11 - tmpint3 ! 12 * (previously "tile_id") + tile_coord(i)%frac_pfaf = nodata_generic + + ! compute area of tile in [km^2] (units convention in tile_coord structure) + + tile_coord(i)%area = ease_cell_area*tile_coord(i)%frac_cell/1000./1000. ! [km^2] - ! change units of area to [km^2] - 23 Sep 2010: fixed units, reichle + else ! not ease grid - tile_coord(i)%area = tile_coord(i)%area*MAPL_RADIUS*MAPL_RADIUS/1000./1000. + if (col_order==1) then + + ! old "SiB2_V2" file format + + read (tmpline,*) & + tile_coord(i)%typ, & ! 1 + tile_coord(i)%pfaf_index, & ! 2 * + tile_coord(i)%com_lon, & ! 3 + tile_coord(i)%com_lat, & ! 4 + tile_coord(i)%i_indg, & ! 5 + tile_coord(i)%j_indg, & ! 6 + tile_coord(i)%frac_cell, & ! 7 + tmpint1, & ! 8 + tmpint2, & ! 9 * + tmpint3, & ! 10 + tile_coord(i)%frac_pfaf, & ! 11 + tmpint4, & ! 12 (previously "tile_id") + tile_coord(i)%area ! 13 + + else + + read (tmpline,*) & + tile_coord(i)%typ, & ! 1 + tile_coord(i)%area, & ! 2 * + tile_coord(i)%com_lon, & ! 3 + tile_coord(i)%com_lat, & ! 4 + tile_coord(i)%i_indg, & ! 5 + tile_coord(i)%j_indg, & ! 6 + tile_coord(i)%frac_cell, & ! 7 + tmpint1, & ! 8 + tile_coord(i)%pfaf_index, & ! 9 * + tmpint2, & ! 10 + tile_coord(i)%frac_pfaf, & ! 11 + tmpint3 ! 12 * (previously "tile_id") + + ! change units of area to [km^2] - 23 Sep 2010: fixed units, reichle + + tile_coord(i)%area = tile_coord(i)%area*MAPL_RADIUS*MAPL_RADIUS/1000./1000. + + end if ! col_order 1 - end if ! col_order 1 + end if ! (ease_grid) - end if ! (ease_grid) - - ! fix i_indg and j_indg such that they refer to a global grid - ! (see above) - - tile_coord(i)%i_indg = tile_coord(i)%i_indg + i_indg_offset - tile_coord(i)%j_indg = tile_coord(i)%j_indg + j_indg_offset + ! fix i_indg and j_indg such that they refer to a global grid + ! (see above) + + tile_coord(i)%i_indg = tile_coord(i)%i_indg + i_indg_offset + tile_coord(i)%j_indg = tile_coord(i)%j_indg + j_indg_offset + + !else ! WY note: keep reading untile the end of the file + + ! ! exit if not land + + ! if (logit) then + ! write (logunit,*) 'WARNING: Encountered first non-land tile in *.til file.' + ! write (logunit,*) ' Stop reading *.til file under the assumption that' + ! write (logunit,*) ' land tiles are first in *.til file.' + ! write (logunit,*) ' This is NOT a safe assumption beyond Icarus-NLv[x] tile spaces!!' + ! end if + ! + ! exit ! assuming land comes first in the til file + + endif + + end do + + close(10) + + ! i here is the number of tiles including all types specified in 'types' - else + n_tiles = i - ! exit if not land - - if (logit) then - write (logunit,*) 'WARNING: Encountered first non-land tile in *.til file.' - write (logunit,*) ' Stop reading *.til file under the assumption that' - write (logunit,*) ' land tiles are first in *.til file.' - write (logunit,*) ' This is NOT a safe assumption beyond Icarus-NLv[x] tile spaces!!' - end if - - exit ! assuming land comes first in the til file - - endif - - end do + ! allocate and fill output variables (r2g, tile_coord_r) - close(10) - - N_tile_land=i - allocate(tile_coord_land(N_tile_land)) - tile_coord_land=tile_coord(1:N_tile_land) - ! pert_[x]_indg is not written into the tile_coord file and not needed in preprocessing - tile_coord_land%pert_i_indg = nint(nodata_generic) - tile_coord_land%pert_j_indg = nint(nodata_generic) - if(present(f2g)) then - allocate(f2g(fid)) - f2g = f2g_tmp(1:fid) - endif - - call read_catchment_def( catch_file, N_tile_land, tile_coord_land ) - - ! ---------------------------------------------------------------------- - ! - ! if elevation info is still needed, read *gridded* elevation data (check only first tile!) - - ! gridded elevation file is NOT available for EASE grids, where elevation information - ! is in catchment.def file - - if ( abs(tile_coord_land(1)%elev-nodata_generic)topo_DYN_ave.file') - open(10,file='topo_DYN_ave.file', action='read') - fname= '' - read(10,'(A)') fname - close(10) - call read_grid_elev( trim(fname), tile_grid_g, N_tile_land, tile_coord_land ) + deallocate(tile_coord) + deallocate(r2g_tmp) + ! pert_[x]_indg is not written into the tile_coord file and not needed in preprocessing + tile_coord_r%pert_i_indg = nint(nodata_generic) + tile_coord_r%pert_j_indg = nint(nodata_generic) + tile_coord_r%elev = nodata_generic + call read_catchment_def( catch_def_file, N_tiles_land_g, tile_coord_r ) + + ! ---------------------------------------------------------------------- + ! + ! if elevation info is still needed, read *gridded* elevation data (check only first tile!) - end if - - if ( abs(tile_coord_land(1)%elev-nodata_generic)topo_DYN_ave.file') + open(10,file='topo_DYN_ave.file', action='read') + fname= '' + read(10,'(A)') fname + close(10) + call read_grid_elev( trim(fname), tile_grid_g, n_tiles, tile_coord_r ) + + end if + + if ( tile_coord_r(1)%typ == MAPL_LAND .and. abs(tile_coord_r(1)%elev-nodata_generic)null() end type METFORCE_WRAP -contains + integer, parameter :: k_force = 12 + integer, parameter :: k_aerosol = 18 + integer, parameter :: k_landice = 15 + character(len=7), dimension(k_force) :: export_name = ['Tair ', 'Qair ', 'Psurf ', & + 'Rainf_C', 'Snowf ', 'LWdown ', & + 'PARdrct', 'PARdffs', 'Wind ', & + 'RefH ', 'Rainf ', 'SWdown '] + character(len=4), dimension(k_aerosol) :: aerosol_name = [ & + 'DUDP', 'DUSV', 'DUWT', 'DUSD', 'BCDP', 'BCSV', & + 'BCWT', 'BCSD', 'OCDP', 'OCSV', 'OCWT', 'OCSD', & + 'SUDP', 'SUSV', 'SUWT', 'SUSD', 'SSDP', 'SSSV' ] + character(len=11), dimension(k_landice) :: landice_name = ['TA ', 'QA ', 'PS ', & + 'PCU ', 'SNO ', 'LWDNSRF ', & + 'DRPAR ', 'DFPAR ', 'UU ', & + 'DZ ', 'DRNIR ', 'DFNIR ', & + 'DRUVR ', 'DFUVR ', 'PLS '] + integer :: NUM_LAND_TILE, NUM_LANDICE_TILE + contains !BOP @@ -115,6 +132,8 @@ subroutine SetServices(gc, rc) rc=status & ) VERIFY_(status) + + ! phase 1 get forcing call MAPL_GridCompSetEntryPoint( & gc, & ESMF_METHOD_RUN, & @@ -122,6 +141,33 @@ subroutine SetServices(gc, rc) rc=status & ) VERIFY_(status) + ! phase 2: to land + call MAPL_GridCompSetEntryPoint( & + gc, & + ESMF_METHOD_RUN, & + DistributeForcingToLand, & + rc=status & + ) + VERIFY_(status) + + ! phase 3: to landpert + call MAPL_GridCompSetEntryPoint( & + gc, & + ESMF_METHOD_RUN, & + DistributeForcingToLandPert, & + rc=status & + ) + VERIFY_(status) + + ! phase 4: to landice + call MAPL_GridCompSetEntryPoint( & + gc, & + ESMF_METHOD_RUN, & + DistributeForcingToLandIce, & + rc=status & + ) + VERIFY_(status) + call MAPL_GridCompSetEntryPoint( & gc, & ESMF_METHOD_FINALIZE, & @@ -559,12 +605,14 @@ subroutine Initialize(gc, import, export, clock, rc) type(tile_coord_type), pointer :: tile_coord(:)=>null() ! Misc variables - integer :: land_nt_local, k, NUM_ENSEMBLE + integer :: local_nt, k, NUM_ENSEMBLE, i1, i2, j1, j2 integer :: ForceDtStep type(met_force_type) :: mf_nodata logical :: MERRA_file_specs, ensemble_forcing logical :: backward_looking_fluxes - + real, pointer :: TileLats(:) + real, pointer :: TileLons(:) + integer, pointer :: tiletype(:) integer :: AEROSOL_DEPOSITION type(MAPL_LocStream) :: locstream character(len=ESMF_MAXSTR) :: grid_type, ENS_FORCING_STR, ens_forcing_path @@ -596,23 +644,25 @@ subroutine Initialize(gc, import, export, clock, rc) VERIFY_(status) internal => wrap%ptr - ! Get component's internal tile_coord variable + ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) tile_coord => tcwrap%ptr%tile_coord - ! Number of land tiles (on local PE) call MAPL_Get(MAPL, LocStream=locstream) VERIFY_(status) call MAPL_LocStreamGet( & locstream, & - NT_LOCAL=land_nt_local, & + NT_LOCAL=local_nt, & + TILELATS=TileLats, & + TILELONS=TileLons, & + TILETYPE=tiletype, & rc=status & ) VERIFY_(status) - call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", & - DEFAULT=0, RC=STATUS) + NUM_LAND_TILE = count(tiletype == MAPL_LAND) + NUM_LANDICE_TILE = count(tiletype == MAPL_LANDICE) call MAPL_GetResource(MAPL, grid_type,Label="GEOSldas.GRID_TYPE:",RC=STATUS) VERIFY_(STATUS) @@ -625,6 +675,9 @@ subroutine Initialize(gc, import, export, clock, rc) im_world_cs = dims(1) endif + call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", & + DEFAULT=0, RC=STATUS) + call MAPL_GetResource(MAPL, gridname,Label="GEOSldas.GRIDNAME:",RC=STATUS) VERIFY_(STATUS) if( index(trim(gridname), 'EASE') /=0) call set_neighbor_offset(0.0001) @@ -663,12 +716,12 @@ subroutine Initialize(gc, import, export, clock, rc) VERIFY_(status) ! -allocate-memory-for-metforcing-data- mf_nodata = nodata_generic - allocate(mf%DataPrv(land_nt_local), source=mf_nodata, stat=status) + allocate(mf%DataPrv(local_nt), source=mf_nodata, stat=status) VERIFY_(status) - allocate(mf%DataNxt(land_nt_local), source=mf_nodata, stat=status) + allocate(mf%DataNxt(local_nt), source=mf_nodata, stat=status) VERIFY_(status) ! -allocate-memory-for-avg-zenith-angle - allocate(mf%zenav(land_nt_local), source=nodata_generic, stat=status) + allocate(mf%zenav(local_nt), source=nodata_generic, stat=status) VERIFY_(status) call MAPL_GetResource ( MAPL, ENS_FORCING_STR, Label="ENSEMBLE_FORCING:", DEFAULT="NO", RC=STATUS) VERIFY_(STATUS) @@ -710,7 +763,7 @@ subroutine Initialize(gc, import, export, clock, rc) ForceDtStep, & internal%mf%Path, & internal%mf%Tag, & - land_nt_local, & + local_nt, & tile_coord, & internal%mf%hinterp, & AEROSOL_DEPOSITION, & @@ -782,22 +835,20 @@ subroutine Run(gc, import, export, clock, rc) ! Private internal state variables type(T_METFORCE_STATE), pointer :: internal=>null() - type(METFORCE_WRAP) :: wrap + type(METFORCE_WRAP) :: wrap type(TILECOORD_WRAP) :: tcwrap ! LDAS' tile_coord variable type(tile_coord_type), pointer :: tile_coord(:) ! Misc variables - integer :: land_nt_local ! number of LAND tiles in local PE + integer :: local_nt ! number of tiles in local PE integer :: comm logical :: IAmRoot integer :: fdtstep - integer :: YEAR, DAY_OF_YEAR, SEC_OF_DAY,n - real, pointer :: LandTileLats(:) - real, pointer :: LandTileLons(:) + real, pointer :: TileLats(:) + real, pointer :: TileLons(:) real, allocatable :: zth(:), slr(:), zth_tmp(:) type(met_force_type), allocatable :: mfDataNtp(:) type(met_force_type), pointer, contiguous :: DataTmp(:)=>null() - real, allocatable :: tmpreal(:) type(met_force_type) :: mf_nodata logical :: MERRA_file_specs @@ -869,7 +920,7 @@ subroutine Run(gc, import, export, clock, rc) VERIFY_(status) internal => wrap%ptr - call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", & + call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", & DEFAULT=1, RC=STATUS) ! Get number of tiles, tile lats/lons from LocStream @@ -877,9 +928,9 @@ subroutine Run(gc, import, export, clock, rc) VERIFY_(status) call MAPL_LocStreamGet( & locstream, & - NT_LOCAL=land_nt_local, & - TILELATS=LandTileLats, & - TILELONS=LandTileLons, & + NT_LOCAL=local_nt, & + TILELATS=TileLats, & + TILELONS=TileLons, & rc=status & ) VERIFY_(status) @@ -888,11 +939,11 @@ subroutine Run(gc, import, export, clock, rc) call MAPL_Get(MAPL, orbit=orbit) ! Allocate memory for zenith angle - allocate(zth(land_nt_local), source=nodata_generic, stat=status) + allocate(zth( local_nt), source=nodata_generic, stat=status) VERIFY_(status) - allocate(slr(land_nt_local), source=nodata_generic, stat=status) + allocate(slr( local_nt), source=nodata_generic, stat=status) VERIFY_(status) - allocate(zth_tmp(land_nt_local), source=nodata_generic, stat=status) + allocate(zth_tmp(local_nt), source=nodata_generic, stat=status) VERIFY_(status) ! Convert forcing time interval to seconds @@ -902,11 +953,12 @@ subroutine Run(gc, import, export, clock, rc) call ESMF_ClockGetAlarm(clock, 'MetForcing', MetForcingAlarm, rc=status) VERIFY_(status) - ! Get component's internal tile_coord variable + ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) tile_coord => tcwrap%ptr%tile_coord + ! Time stamp of next model step ! -get-model-time-step- call ESMF_ClockGet(clock, timeStep=ModelTimeStep) @@ -938,7 +990,7 @@ subroutine Run(gc, import, export, clock, rc) fdtstep, & internal%mf%Path, & internal%mf%Tag, & - land_nt_local, & + local_nt, & tile_coord, & internal%mf%hinterp, & AEROSOL_DEPOSITION, & @@ -956,10 +1008,10 @@ subroutine Run(gc, import, export, clock, rc) ! -compute-average-zenith-angle-over-daylight-part-of-forcing-interval- call MAPL_SunGetInsolation( & - LandTileLons, & - LandTileLats, & + TileLons, & + TileLats, & orbit, & - zth_tmp, & + zth_tmp, & slr, & currTime=internal%mf%TimePrv, & INTV=internal%mf%ntrvl, & @@ -973,7 +1025,7 @@ subroutine Run(gc, import, export, clock, rc) ! dayOfYear=DAY_OF_YEAR, RC=STATUS) ! VERIFY_(STATUS) - ! call zenith(DAY_OF_YEAR,SEC_OF_DAY,fdtstep,ModelTimeStep,land_nt_local,tile_coord%com_lon, & + ! call zenith(DAY_OF_YEAR,SEC_OF_DAY,fdtstep,ModelTimeStep,local_nt,tile_coord%com_lon, & ! tile_coord%com_lat,internal%mf%zenav) @@ -989,8 +1041,8 @@ subroutine Run(gc, import, export, clock, rc) ! Compute zenith angle at the next time step call MAPL_SunGetInsolation( & - LandTileLons, & - LandTileLats, & + TileLons, & + TileLats, & orbit, & zth_tmp, & slr, & @@ -1005,7 +1057,7 @@ subroutine Run(gc, import, export, clock, rc) !call ESMF_TimeGet(ModelTimeNxt, YY=YEAR, S=SEC_OF_DAY, & ! dayOfYear=DAY_OF_YEAR, RC=STATUS) !VERIFY_(STATUS) - !do n=1, land_nt_local + !do n=1, local_nt ! call solar(tile_coord(n)%com_lon,tile_coord(n)%com_lat, DAY_OF_YEAR,SEC_OF_DAY,zth(n),slr(n)) !enddo @@ -1028,13 +1080,11 @@ subroutine Run(gc, import, export, clock, rc) ! Allocate memory for interpolated MetForcing data mf_nodata = nodata_generic - allocate(mfDataNtp(land_nt_local), source=mf_nodata, stat=status) + allocate(mfDataNtp(local_nt), source=mf_nodata, stat=status) VERIFY_(status) ! Interpolate MetForcing data to the end of model integration time step call LDAS_TInterpForcing( & - tile_coord%com_lon, & - tile_coord%com_lat, & zth, & internal%mf%zenav, & force_time_prv, & @@ -1224,6 +1274,133 @@ subroutine Run(gc, import, export, clock, rc) end subroutine Run + subroutine DistributeForcingToLand(gc, export, land_import, clock, rc) + + type(ESMF_GridComp), intent(inout) :: gc ! Gridded component + type(ESMF_State), intent(inout) :: export ! Export state + type(ESMF_State), intent(inout) :: land_import ! Import state + type(ESMF_Clock), intent(inout) :: clock ! The clock + integer, optional, intent( out) :: rc ! Error code + + real, pointer :: out1d(:), in1d(:) + real, pointer :: out2d(:,:), in2d(:,:) + integer :: k, AEROSOL_DEPOSITION, status + type(MAPL_MetaComp), pointer :: MAPL + character(len=ESMF_MAXSTR) :: Iam + Iam = "metForce::DistributeForcingToLand" + + call MAPL_GetObjectFromGC(gc, MAPL, _RC) + call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", DEFAULT=1, _RC) + if(AEROSOL_DEPOSITION /=0) then + do k = 1, k_aerosol + call MAPL_GetPointer(export, out2d, aerosol_name(k), _RC) + call MAPL_GetPointer(land_import, in2d, aerosol_name(k), _RC) + in2d(:,:) = out2d(1:NUM_LAND_TILE, :) + enddo + endif + + call MAPL_GetPointer(export, out1d, 'Psurf', _RC) + call MAPL_GetPointer(land_import, in1d, 'PS', _RC) + in1d = out1d(1:NUM_LAND_TILE) + call MAPL_GetPointer(export, out1d, 'RefH', _RC) + call MAPL_GetPointer(land_import, in1d, 'DZ', _RC) + in1d = out1d(1:NUM_LAND_TILE) + RETURN_(ESMF_SUCCESS) + + end subroutine DistributeForcingToLand + + subroutine DistributeForcingToLandPert(gc, export, landpert_import, clock, rc) + + type(ESMF_GridComp), intent(inout) :: gc ! Gridded component + type(ESMF_State), intent(inout) :: export ! Export state + type(ESMF_State), intent(inout) :: landpert_import ! Import state + type(ESMF_Clock), intent(inout) :: clock ! The clock + integer, optional, intent( out) :: rc ! Error code + + real, pointer :: out1d(:), in1d(:) + integer :: k, status + character(len=ESMF_MAXSTR) :: Iam + Iam = "metForce::DistributeForcingToLandPert" + + do k = 1, k_force + call MAPL_GetPointer(export, out1d, trim(export_name(k)), _RC) + call MAPL_GetPointer(landpert_import, in1d, trim(export_name(k)), _RC) + in1d = out1d(1:NUM_LAND_TILE) + enddo + RETURN_(ESMF_SUCCESS) + + end subroutine DistributeForcingToLandPert + + subroutine DistributeForcingToLandIce(gc, export, landice_import, clock, rc) + + type(ESMF_GridComp), intent(inout) :: gc ! Gridded component + type(ESMF_State), intent(inout) :: export ! Export state + type(ESMF_State), intent(inout) :: landice_import ! Import state + type(ESMF_Clock), intent(inout) :: clock ! The clock + integer, optional, intent( out) :: rc ! Error code + + integer :: k, i1, i2, AEROSOL_DEPOSITION, status + real, pointer :: out1d(:), in1d(:), tmp(:) + real, pointer :: out2d(:,:), in2d(:,:) + real, allocatable :: tmpreal(:) + type(MAPL_MetaComp), pointer :: MAPL + character(len=ESMF_MAXSTR) :: Iam + Iam = "metForce::DistributeForcingToLandice" + + if (NUM_LANDICE_TILE == 0) then + RETURN_(ESMF_SUCCESS) + endif + + i1 = NUM_LAND_TILE + 1 + i2 = NUM_LAND_TILE + NUM_LANDICE_TILE + ! Get MAPL obj + call MAPL_GetObjectFromGC(gc, MAPL, _RC) + call MAPL_GetResource ( MAPL, AEROSOL_DEPOSITION, Label="AEROSOL_DEPOSITION:", DEFAULT=1, _RC) + if(AEROSOL_DEPOSITION /=0) then + do k = 1, k_aerosol + call MAPL_GetPointer(export, out2d, aerosol_name(k), _RC) + call MAPL_GetPointer(landice_import, in2d, aerosol_name(k), _RC) + in2d(:,:) = out2d(i1:i2, :) + VERIFY_(status) + enddo + endif + + do k = 1, k_force - 2 + call MAPL_GetPointer(export, out1d, trim(export_name(k)), _RC) + call MAPL_GetPointer(landice_import, in1d, trim(landice_name(k)), _RC) + in1d = out1d(i1:i2) + enddo + + call MAPL_GetPointer(export, out1d, 'Wind', _RC) + call MAPL_GetPointer(landice_import, in1d, 'UWINDLMTILE', _RC) + in1d = out1d(i1:i2) + call MAPL_GetPointer(landice_import, in1d, 'VWINDLMTILE', _RC) + in1d = 0. + + call MAPL_GetPointer(export, out1d, 'Rainf', _RC) + call MAPL_GetPointer(export, tmp, 'Rainf_C', _RC) + call MAPL_GetPointer(landice_import, in1d, 'PLS', _RC) + in1d = out1d(i1:i2) - tmp(i1:i2) + + allocate(tmpreal(NUM_LANDICE_TILE), stat=status) + call MAPL_GetPointer(export, out1d, 'SWdown', _RC) + tmpreal = 0.5* out1d(i1:i2) + call MAPL_GetPointer(landice_import, in1d, 'DRNIR', _RC) + in1d = 0.5 * tmpreal + call MAPL_GetPointer(landice_import, in1d, 'DFNIR', _RC) + in1d = 0.5 * tmpreal + + call MAPL_GetPointer(export, out1d, 'PARdrct', _RC) + call MAPL_GetPointer(landice_import, in1d, 'DRUVR', _RC) + in1d = 0.5* tmpreal - out1d(i1:i2) + call MAPL_GetPointer(export, out1d, 'PARdffs', _RC) + call MAPL_GetPointer(landice_import, in1d, 'DFUVR', _RC) + in1d = 0.5* tmpreal - out1d(i1:i2) + deallocate(tmpreal) + + RETURN_(ESMF_SUCCESS) + + end subroutine DistributeForcingToLandIce !BOP diff --git a/GEOSmetforce_GridComp/LDAS_Interp.F90 b/GEOSmetforce_GridComp/LDAS_Interp.F90 index 918bfcd3..951abf38 100644 --- a/GEOSmetforce_GridComp/LDAS_Interp.F90 +++ b/GEOSmetforce_GridComp/LDAS_Interp.F90 @@ -19,8 +19,6 @@ module LDAS_InterpMod contains subroutine metforcing_tinterp( & - lons, & - lats, & zth, & zenav, & force_time_prv, & @@ -57,8 +55,6 @@ subroutine metforcing_tinterp( & implicit none - real, intent(in) :: lats(:) ! tile lats - real, intent(in) :: lons(:) ! tile lons ! fix potential inconsistency between zth and zenav owing to 300s time ! step used in MAPL_SunGetInsolation() ! in ==>inout diff --git a/LDAS_Shared/LDAS_Convert.F90 b/LDAS_Shared/LDAS_Convert.F90 index 71b4acb9..cf155d94 100644 --- a/LDAS_Shared/LDAS_Convert.F90 +++ b/LDAS_Shared/LDAS_Convert.F90 @@ -21,9 +21,9 @@ module LDAS_ConvertMod subroutine esmf2ldas_time(esmf_dt, ldas_dt, rc) - type(ESMF_Time), intent(in) :: esmf_dt + type(ESMF_Time), intent(in) :: esmf_dt type(date_time_type), intent(out) :: ldas_dt - integer, optional, intent(out) :: rc + integer, optional, intent(out) :: rc character(len=*), parameter :: Iam = 'emsf2ldas_time' integer :: status @@ -46,4 +46,8 @@ subroutine esmf2ldas_time(esmf_dt, ldas_dt, rc) end subroutine esmf2ldas_time + ! -------------------------------------------- + end module LDAS_ConvertMod + +! ========= EOF =========================================================== diff --git a/LDAS_Shared/LDAS_RepairForcing.F90 b/LDAS_Shared/LDAS_RepairForcing.F90 index 979e8fb2..be06594c 100644 --- a/LDAS_Shared/LDAS_RepairForcing.F90 +++ b/LDAS_Shared/LDAS_RepairForcing.F90 @@ -57,7 +57,7 @@ subroutine repair_forcing( N_catd, met_force, echo, tile_coord, fieldname, & ! min/max values for allowable range of forcing fields - real, parameter :: min_Tair = 190. ! [K] + real, parameter :: min_Tair = 180. ! [K] ! changed from 190 K to 180 K to accommodate landice, lcandre2, May 2025 real, parameter :: max_Tair = 340. ! [K] real, parameter :: max_PSurf = 115000. ! [Pa] @@ -235,7 +235,7 @@ subroutine repair_forcing( N_catd, met_force, echo, tile_coord, fieldname, & ! NOTE: "warn" is turned on when repair_forcing is called first ! time after the forcing has been read from files - if ((warn) .and. (met_force(i)%Tair < 190.)) then + if ((warn) .and. (met_force(i)%Tair < min_Tair)) then write (tmpstr13a,'(e13.5)') met_force(i)%Tair ! convert real to string