From 75475632c30b5545bf6fad79693b48f8948aa16f Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 22 Jun 2020 09:17:41 -0400 Subject: [PATCH 1/3] remove cs index --- .../GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 | 5 +++++ .../GEOSmetforce_GridComp/LDAS_Forcing.F90 | 13 +++++++++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 index 631a9767..a75ac768 100644 --- a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 @@ -18,6 +18,7 @@ module GEOS_MetforceGridCompMod use LDAS_ForceMod, only: LDAS_GetForcing => get_forcing use LDAS_ForceMod, only: LDAS_move_new_force_to_old use LDAS_ForceMod, only: FileOpenedHash,GEOS_closefile + use LDAS_ForceMod, only: i_indg,j_indg use LDAS_DriverTypes, only: met_force_type, assignment(=) use LDAS_ConvertMod, only: esmf2ldas use LDAS_InterpMod, only: LDAS_TInterpForcing=>metforcing_tinterp @@ -34,6 +35,8 @@ module GEOS_MetforceGridCompMod real, parameter :: daylen = 86400. + integer, pointer :: cs_i_indg(:) =>null() + integer, pointer :: cs_j_indg(:) =>null() ! !PUBLIC MEMBER FUNCTIONS: public :: SetServices @@ -621,6 +624,8 @@ subroutine Initialize(gc, import, export, clock, rc) call MAPL_LocStreamGet( & locstream, & NT_LOCAL=land_nt_local, & + GRIDIM = i_indg, & + GRIDJM = j_indg, & rc=status & ) VERIFY_(status) diff --git a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 882b5ef2..9b86210a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -59,6 +59,8 @@ module LDAS_ForceMod public :: get_forcing public :: LDAS_move_new_force_to_old public :: GEOS_closefile + public :: i_indg, j_indg + type(Hash_Table),public :: FileOpenedHash !public :: ignore_SWNET_for_snow @@ -85,6 +87,9 @@ module LDAS_ForceMod end type local_grid type(local_grid), target :: local_info + + integer, pointer :: i_indg(:)=>null() + integer, pointer :: j_indg(:)=>null() contains @@ -4674,9 +4679,13 @@ subroutine GEOS_openfile(FileOpenedHash, fname_full, fid, tile_coord, m_hinterp, if( isCubed ) then ! cs grid ! i_indg and j_indg are changed to LatLon grid do k=1,N_cat - i1(k) = tile_coord(k)%cs_i_indg - j1(k) = tile_coord(k)%cs_j_indg + !i1(k) = tile_coord(k)%cs_i_indg + !j1(k) = tile_coord(k)%cs_j_indg + i1(k) = i_indg(k) + j1(k) = j_indg(k) enddo + if (any(i1 /= tile_coord%cs_i_indg)) stop " wrong cs-i" + if (any(j1 /= tile_coord%cs_j_indg)) stop " wrong cs-j" else do k=1,N_cat From b99b12a97d13c955e4d346028aef4da5e7bd3043 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Thu, 9 Jul 2020 10:46:25 -0400 Subject: [PATCH 2/3] remove cs_i_indg and cs_j_indg from tile_coord structure --- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 19 +++++++++++++++---- .../GEOS_MetforceGridComp.F90 | 5 ----- .../GEOSmetforce_GridComp/LDAS_Forcing.F90 | 15 ++++++--------- .../Shared/LDAS_TileCoordType.F90 | 4 ---- .../Shared/LDAS_ensdrv_mpi.F90 | 4 +--- 5 files changed, 22 insertions(+), 25 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index b4ffaa3a..d5f8a442 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -31,7 +31,7 @@ module GEOS_LdasGridCompMod use lsm_routines, only: lsmroutines_echo_constants use StieglitzSnow, only: StieglitzSnow_echo_constants use SurfParams, only: SurfParams_init - + use LDAS_ForceMod, only: cs_i_indg, cs_j_indg implicit none private @@ -371,6 +371,7 @@ subroutine Initialize(gc, import, export, clock, rc) integer,dimension(:),pointer :: f2g integer :: N_catf + integer, allocatable :: cs_i_indg_f(:), cs_j_indg_f(:) type(grid_def_type) :: tile_grid_g type(grid_def_type) :: tile_grid_f @@ -591,9 +592,8 @@ subroutine Initialize(gc, import, export, clock, rc) ASSERT_(index(tile_grid_g%gridtype, 'c3') /=0) !1) save original index - tile_coord_f%cs_i_indg = tile_coord_f%i_indg - tile_coord_f%cs_j_indg = tile_coord_f%j_indg - + cs_i_indg_f = tile_coord_f%i_indg + cs_j_indg_f = tile_coord_f%j_indg !2) generate a lat-lon grid for landpert and land assim ( 4*N_lonX3*N_lon) call get_pert_grid(tile_grid_g, latlon_tmp_g) tile_grid_g = latlon_tmp_g @@ -632,6 +632,17 @@ subroutine Initialize(gc, import, export, clock, rc) tcinternal%l2f = tile_id2f(local_id) tcinternal%tile_coord = tile_coord_f(tcinternal%l2f) deallocate(f2tile_id, tile_id2f) + + if (trim(grid_type) == "Cubed-Sphere" ) then + if (.not. IamRoot) allocate(cs_i_indg_f(N_catf), cs_j_indg_f(N_catf)) + call MPI_BCAST(cs_i_indg_f, N_catf, MPI_INTEGER, 0, mpicomm, mpierr) + call MPI_BCAST(cs_j_indg_f, N_catf, MPI_INTEGER, 0, mpicomm, mpierr) + allocate(cs_i_indg(land_nt_local), cs_j_indg(land_nt_local)) + cs_i_indg = cs_i_indg_f(tcinternal%l2f) + cs_j_indg = cs_j_indg_f(tcinternal%l2f) + deallocate(cs_i_indg_f, cs_j_indg_f) + endif + end block do i = 0, numprocs-1 diff --git a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 index a75ac768..631a9767 100644 --- a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 @@ -18,7 +18,6 @@ module GEOS_MetforceGridCompMod use LDAS_ForceMod, only: LDAS_GetForcing => get_forcing use LDAS_ForceMod, only: LDAS_move_new_force_to_old use LDAS_ForceMod, only: FileOpenedHash,GEOS_closefile - use LDAS_ForceMod, only: i_indg,j_indg use LDAS_DriverTypes, only: met_force_type, assignment(=) use LDAS_ConvertMod, only: esmf2ldas use LDAS_InterpMod, only: LDAS_TInterpForcing=>metforcing_tinterp @@ -35,8 +34,6 @@ module GEOS_MetforceGridCompMod real, parameter :: daylen = 86400. - integer, pointer :: cs_i_indg(:) =>null() - integer, pointer :: cs_j_indg(:) =>null() ! !PUBLIC MEMBER FUNCTIONS: public :: SetServices @@ -624,8 +621,6 @@ subroutine Initialize(gc, import, export, clock, rc) call MAPL_LocStreamGet( & locstream, & NT_LOCAL=land_nt_local, & - GRIDIM = i_indg, & - GRIDJM = j_indg, & rc=status & ) VERIFY_(status) diff --git a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 9b86210a..e4e4cd74 100644 --- a/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -59,7 +59,7 @@ module LDAS_ForceMod public :: get_forcing public :: LDAS_move_new_force_to_old public :: GEOS_closefile - public :: i_indg, j_indg + public :: cs_i_indg, cs_j_indg ! initialized in GEOS_LdasGridComp type(Hash_Table),public :: FileOpenedHash @@ -88,8 +88,9 @@ module LDAS_ForceMod type(local_grid), target :: local_info - integer, pointer :: i_indg(:)=>null() - integer, pointer :: j_indg(:)=>null() + !! initialized in GEOS_LdasGridComp + integer, pointer :: cs_i_indg(:)=>null() + integer, pointer :: cs_j_indg(:)=>null() contains @@ -4679,13 +4680,9 @@ subroutine GEOS_openfile(FileOpenedHash, fname_full, fid, tile_coord, m_hinterp, if( isCubed ) then ! cs grid ! i_indg and j_indg are changed to LatLon grid do k=1,N_cat - !i1(k) = tile_coord(k)%cs_i_indg - !j1(k) = tile_coord(k)%cs_j_indg - i1(k) = i_indg(k) - j1(k) = j_indg(k) + i1(k) = cs_i_indg(k) + j1(k) = cs_j_indg(k) enddo - if (any(i1 /= tile_coord%cs_i_indg)) stop " wrong cs-i" - if (any(j1 /= tile_coord%cs_j_indg)) stop " wrong cs-j" else do k=1,N_cat diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index 20376344..6889c26b 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -76,10 +76,6 @@ module LDAS_TileCoordType real :: max_lat ! maximum latitude (bounding box for tile) integer :: i_indg ! i index (w.r.t. *global* grid that cuts tiles) integer :: j_indg ! j index (w.r.t. *global* grid that cuts tiles) - !if it is Cubed-Sphere grid, the index will be saved here for forcing - !i_indg and j_indg will be changed to index that related to latlon grid - integer :: cs_i_indg ! i index (w.r.t. *global* grid that cuts tiles) - integer :: cs_j_indg ! j index (w.r.t. *global* grid that cuts tiles) real :: frac_cell ! area fraction of grid cell covered by tile real :: frac_pfaf ! fraction of Pfafstetter catchment for land tiles real :: area ! area [km^2] diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 index 1029db8e..5d502712 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 @@ -138,8 +138,6 @@ subroutine init_MPI_types() ! real :: max_lat ! maximum latitude (bounding box for tile) ! integer :: i_indg ! i index (w.r.t. *global* grid that cuts tiles) ! integer :: j_indg ! j index (w.r.t. *global* grid that cuts tiles) - ! integer :: cs_i_indg ! i index (w.r.t. *global* grid that cuts tiles) - ! integer :: cs_j_indg ! j index (w.r.t. *global* grid that cuts tiles) ! real :: frac_cell ! area fraction of grid cell covered by tile ! real :: frac_pfaf ! fraction of Pfafstetter catchment for land tiles ! real :: area ! area [km^2] @@ -158,7 +156,7 @@ subroutine init_MPI_types() iblock(1) = 4 iblock(2) = 6 - iblock(3) = 4 ! add cs_i_indg and cs_j_indg + iblock(3) = 2 iblock(4) = 4 idisp(1) = 0 From 6d60edc50ee00e99830fbac4ae80f48f439e3c8d Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Tue, 28 Jul 2020 10:38:54 -0400 Subject: [PATCH 3/3] move some subroutines to LDAS_app --- src/Applications/LDAS_App/CMakeLists.txt | 9 +- src/Applications/LDAS_App/ldas_help_subs.F90 | 519 ++++++++++++++++++ src/Applications/LDAS_App/preprocess_ldas.F90 | 2 +- .../Shared/LDAS_TileCoordRoutines.F90 | 497 ----------------- 4 files changed, 528 insertions(+), 499 deletions(-) create mode 100644 src/Applications/LDAS_App/ldas_help_subs.F90 diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index a869ac90..cae291df 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -1,3 +1,10 @@ +esma_set_this () + +set (srcs + ldas_help_subs.F90 + ) +esma_add_library(${this} SRCS ${srcs} DEPENDENCIES GEOSldas_GridComp MAPL) + ecbuild_add_executable ( TARGET GEOSldas.x SOURCES GEOSldas.F90 @@ -13,7 +20,7 @@ foreach (prog ${executables}) ecbuild_add_executable ( TARGET ${prog}.x SOURCES ${prog}.F90 - LIBS GEOSldas_GridComp mk_restarts) + LIBS ${this} GEOSldas_GridComp mk_restarts) endforeach () install( diff --git a/src/Applications/LDAS_App/ldas_help_subs.F90 b/src/Applications/LDAS_App/ldas_help_subs.F90 new file mode 100644 index 00000000..450d97f8 --- /dev/null +++ b/src/Applications/LDAS_App/ldas_help_subs.F90 @@ -0,0 +1,519 @@ + +! This file is part of orginal LDAS_TileCoordRoutines.F90 +! It contains the subroutineis used only in LDAS_app directory + +module LDAS_help_subs + use LDAS_TileCoordType, ONLY: & + tile_coord_type, & + grid_def_type + + use LDAS_TileCoordRoutines, only: & + LDAS_create_grid_g + + use LDAS_ensdrv_Globals, ONLY: & + logit, & + logunit, & + nodata_generic, & + nodata_tol_generic + + use MAPL_ConstantsMod, ONLY: & + MAPL_RADIUS ! Earth radius + + use LDAS_ExceptionsMod, ONLY: & + ldas_abort, & + LDAS_GENERIC_ERROR + + implicit none + + private + + public :: LDAS_read_land_tile + +contains + + subroutine LDAS_read_land_tile( tile_file,catch_file, tile_grid_g, tile_coord_land,f2g ) + ! inputs: + ! + ! tile_file : full path + name + ! + ! outputs: + ! + ! tile_grid : parameters of tile definition grid + ! + ! tile_coord : coordinates of tiles (see tile_coord_type), + ! implemented as pointer which is allocated in + ! this subroutine + ! NOTE: number of tiles can be diagnosed + ! with size(tile_coord) + ! optional: + ! if the tile file type 1100 ,which is land excluded + ! f2g : the full domain id to the global id + ! "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 + ! + ! ------------------------------------------------------------- + + implicit none + + character(*), intent(in) :: tile_file + character(*), intent(in) :: catch_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 + + ! locals + type(tile_coord_type),dimension(:),allocatable :: tile_coord + integer, dimension(:), allocatable :: f2g_tmp ! out + + real :: ease_cell_area + integer :: i, N_tile,N_grid,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 + character(200) :: tmpline,gridname + character(300) :: fname + + character(len=*), parameter :: Iam = 'LDAS_read_tile_file' + + ! --------------------------------------------------------------- + + i_indg_offset = 0 + j_indg_offset = 0 + + ! call LDAS_create_grid(tile_file,tile_grid_g,i_indg_offset,j_indg_offset) + + ! read file header + + if (logit) write (logunit,'(400A)') 'LDAS_read_tile_file(): reading from' // trim(tile_file) + + + open (10, file=trim(tile_file), 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(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 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. ! 'EASE' and 'EASEv2' + if (index(tile_grid_g%gridtype,'SiB2')/=0) col_order=1 ! Weiyuan note: grid name for SiB2?? + 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 + if(typ == 100 .or. typ ==1100) then ! if it is land + i=i+1 + tile_coord(i)%tile_id = k + + if (typ == 100) then + fid=fid+1 + f2g_tmp(fid) = k + endif + + 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, & ! 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 + tile_coord(i)%area = ease_cell_area*tile_coord(i)%frac_cell + + else ! not ease grid + + if (col_order==1) then + ! old "SiB2_V2" file format + read (tmpline,*) & + tile_coord(i)%typ, & ! 1 + tile_coord(i)%pfaf, & ! 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, & ! 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 ! (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 + else + exit ! land comes first in the til file + endif + end do + close(10) + + N_tile_land=i + allocate(tile_coord_land(N_tile_land)) + tile_coord_land=tile_coord(1:N_tile_land) + + 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 still needed read gridded elevation (check only first tile!) + + 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,status='DELETE') + close(10) + call read_grid_elev( trim(fname), tile_grid_g, N_tile_land, tile_coord_land ) + end if + + + if ( abs(tile_coord_land(1)%elev-nodata_generic)tile_grid%dlon+latlon_tol) then + + if (logit) write (logunit,*) & + 'resetting min/maxlon in tile_id=', tile_coord(k)%tile_id + + ! "push" tile to the east of the dateline + + tile_coord(k)%min_lon = -180. + tile_coord(k)%max_lon = -180. + 0.5*tile_grid%dlon + + end if + + + ! check that com_lon is within tile definition grid cell + + if ( .not. ( & + (tile_coord(k)%com_lon >= this_minlon) .and. & + (tile_coord(k)%com_lon <= this_maxlon) ) ) then + if (logit) write (logunit,*) & + 'resetting com_lon in tile_id=', tile_coord(k)%tile_id + + tile_coord(k)%com_lon = & + 0.5*(tile_coord(k)%min_lon + tile_coord(k)%max_lon) + + end if + + end do + end if + + end subroutine fix_dateline_bug_in_tilecoord + + ! ********************************************************************** + + subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) + + ! reichle, 17 May 2011: read elevation data if available + + ! format of catchment.def file + ! + ! Header line: N_tile + ! + ! Columns: tile_id, Pfaf, min_lon, max_lon, min_lat, max_lat, [elev] + ! + ! Elevation [m] is ONLY available for EASE grid tile definitions + + implicit none + + character(*), intent(in) :: catchment_def_file + + integer, intent(in) :: N_tile + + type(tile_coord_type), dimension(:), pointer :: tile_coord ! inout + + ! locals + + integer :: i, istat, tmpint1, sweep + + integer, dimension(N_tile) :: tmp_tileid, tmp_pfaf + + character(len=*), parameter :: Iam = 'read_catchment_def' + character(len=400) :: err_msg + + ! --------------------------------------------------------------- + + ! read file header + + if (logit) write (logunit,'(400A)') & + 'read_catchment_def(): reading from' // trim(catchment_def_file) + if (logit) write (logunit,*) + + ! sweep=1: Try reading 7 columns. If this fails, try again. + ! sweep=2: Read only 6 columns. + + do sweep=1,2 + + if (logit) write (logunit,*) 'starting sweep ', sweep + + open (10, file=trim(catchment_def_file), form='formatted', action='read') + + read (10,*) tmpint1 + + if (logit) write (logunit,*) 'file contains coordinates for ', tmpint1, ' tiles' + if (logit) write (logunit,*) + + if (N_tile/=tmpint1) then + print*,"need :", N_tile,"but have: ",tmpint1 + err_msg = 'tile_coord_file and catchment_def_file mismatch. (1)' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + do i=1,N_tile + + if (sweep==1) then + + ! read 7 columns, avoid using exact format specification + + read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & + tile_coord(i)%min_lon, & + tile_coord(i)%max_lon, & + tile_coord(i)%min_lat, & + tile_coord(i)%max_lat, & + tile_coord(i)%elev + + else + + ! read 6 columns, avoid using exact format specification + + read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & + tile_coord(i)%min_lon, & + tile_coord(i)%max_lon, & + tile_coord(i)%min_lat, & + tile_coord(i)%max_lat + + tile_coord(i)%elev = nodata_generic + + end if + + if (istat/=0) then ! read error + + if (sweep==1) then + + close(10,status='keep') + + if (logit) write (logunit,*) 'sweep 1 failed, trying sweep 2' + + exit ! exit sweep 1, try sweep 2 + + else + + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'sweep 2 failed') + + end if + + end if + + if (i==N_tile) then ! reached end of tile loop w/o read error + + close(10,status='keep') + + if (logit) write (logunit,*) 'sweep ', sweep, 'successfully completed' + + return + + end if + + end do ! loop through tiles + + end do ! loop through sweeps + + if ( any(tile_coord(1:N_tile)%tile_id/=tmp_tileid) .or. & + any(tile_coord(1:N_tile)%pfaf /=tmp_pfaf) ) then + + err_msg = 'tile_coord_file and catchment_def_file mismatch. (2)' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end if + + ! ----------------------------------------------------------------- + + end subroutine read_catchment_def + +end module LDAS_help_subs diff --git a/src/Applications/LDAS_App/preprocess_ldas.F90 b/src/Applications/LDAS_App/preprocess_ldas.F90 index 54305050..298dffb2 100644 --- a/src/Applications/LDAS_App/preprocess_ldas.F90 +++ b/src/Applications/LDAS_App/preprocess_ldas.F90 @@ -16,7 +16,7 @@ module preprocess_module N_progn_pert_max, & N_force_pert_max use catch_types, only: cat_param_type,N_gt - use LDAS_TileCoordRoutines, ONLY: & + use LDAS_help_subs, ONLY: & LDAS_read_land_tile use LDAS_ensdrv_functions, ONLY: & get_io_filename diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index 8fb991b7..d2910075 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -10,8 +10,6 @@ module LDAS_TileCoordRoutines - use MAPL_BaseMod, only: MAPL_GetHorzIJIndex - use LDAS_TileCoordType, ONLY: & tile_coord_type, & grid_def_type, & @@ -55,7 +53,6 @@ module LDAS_TileCoordRoutines public :: tile_mask_grid public :: is_cat_in_box public :: LDAS_create_grid_g - public :: LDAS_read_land_tile character(10) :: tmpstring10 character(40) :: tmpstring40 @@ -281,226 +278,6 @@ subroutine LDAS_create_grid_g( gridname,n_lon,n_lat, tile_grid,i_indg_offset,j_i end subroutine LDAS_create_grid_g - ! ********************************************************************** - - subroutine LDAS_read_land_tile( tile_file,catch_file, tile_grid_g, tile_coord_land,f2g ) - ! inputs: - ! - ! tile_file : full path + name - ! - ! outputs: - ! - ! tile_grid : parameters of tile definition grid - ! - ! tile_coord : coordinates of tiles (see tile_coord_type), - ! implemented as pointer which is allocated in - ! this subroutine - ! NOTE: number of tiles can be diagnosed - ! with size(tile_coord) - ! optional: - ! if the tile file type 1100 ,which is land excluded - ! f2g : the full domain id to the global id - ! "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 - ! - ! ------------------------------------------------------------- - - implicit none - - character(*), intent(in) :: tile_file - character(*), intent(in) :: catch_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 - - ! locals - type(tile_coord_type),dimension(:),allocatable :: tile_coord - integer, dimension(:), allocatable :: f2g_tmp ! out - - real :: ease_cell_area - integer :: i, N_tile,N_grid,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 - character(200) :: tmpline,gridname - character(300) :: fname - - character(len=*), parameter :: Iam = 'LDAS_read_tile_file' - character(len=400) :: err_msg - - ! --------------------------------------------------------------- - - i_indg_offset = 0 - j_indg_offset = 0 - - ! call LDAS_create_grid(tile_file,tile_grid_g,i_indg_offset,j_indg_offset) - - ! read file header - - if (logit) write (logunit,'(400A)') 'LDAS_read_tile_file(): reading from' // trim(tile_file) - - - open (10, file=trim(tile_file), 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(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 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. ! 'EASE' and 'EASEv2' - if (index(tile_grid_g%gridtype,'SiB2')/=0) col_order=1 ! Weiyuan note: grid name for SiB2?? - 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 - if(typ == 100 .or. typ ==1100) then ! if it is land - i=i+1 - tile_coord(i)%tile_id = k - - if (typ == 100) then - fid=fid+1 - f2g_tmp(fid) = k - endif - - 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, & ! 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 - tile_coord(i)%area = ease_cell_area*tile_coord(i)%frac_cell - - else ! not ease grid - - if (col_order==1) then - ! old "SiB2_V2" file format - read (tmpline,*) & - tile_coord(i)%typ, & ! 1 - tile_coord(i)%pfaf, & ! 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, & ! 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 ! (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 - else - exit ! land comes first in the til file - endif - end do - close(10) - - N_tile_land=i - allocate(tile_coord_land(N_tile_land)) - tile_coord_land=tile_coord(1:N_tile_land) - - 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 still needed read gridded elevation (check only first tile!) - - 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,status='DELETE') - close(10) - call read_grid_elev( trim(fname), tile_grid_g, N_tile_land, tile_coord_land ) - end if - - - if ( abs(tile_coord_land(1)%elev-nodata_generic)tile_grid%dlon+latlon_tol) then - - if (logit) write (logunit,*) & - 'resetting min/maxlon in tile_id=', tile_coord(k)%tile_id - - ! "push" tile to the east of the dateline - - tile_coord(k)%min_lon = -180. - tile_coord(k)%max_lon = -180. + 0.5*tile_grid%dlon - - end if - - - ! check that com_lon is within tile definition grid cell - - if ( .not. ( & - (tile_coord(k)%com_lon >= this_minlon) .and. & - (tile_coord(k)%com_lon <= this_maxlon) ) ) then - if (logit) write (logunit,*) & - 'resetting com_lon in tile_id=', tile_coord(k)%tile_id - - tile_coord(k)%com_lon = & - 0.5*(tile_coord(k)%min_lon + tile_coord(k)%max_lon) - - end if - - end do - end if - - end subroutine fix_dateline_bug_in_tilecoord - - ! ********************************************************************** - - subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) - - ! reichle, 17 May 2011: read elevation data if available - - ! format of catchment.def file - ! - ! Header line: N_tile - ! - ! Columns: tile_id, Pfaf, min_lon, max_lon, min_lat, max_lat, [elev] - ! - ! Elevation [m] is ONLY available for EASE grid tile definitions - - implicit none - - character(*), intent(in) :: catchment_def_file - - integer, intent(in) :: N_tile - - type(tile_coord_type), dimension(:), pointer :: tile_coord ! inout - - ! locals - - integer :: i, istat, tmpint1, sweep - - integer, dimension(N_tile) :: tmp_tileid, tmp_pfaf - - character(len=*), parameter :: Iam = 'read_catchment_def' - character(len=400) :: err_msg - - ! --------------------------------------------------------------- - - ! read file header - - if (logit) write (logunit,'(400A)') & - 'read_catchment_def(): reading from' // trim(catchment_def_file) - if (logit) write (logunit,*) - - ! sweep=1: Try reading 7 columns. If this fails, try again. - ! sweep=2: Read only 6 columns. - - do sweep=1,2 - - if (logit) write (logunit,*) 'starting sweep ', sweep - - open (10, file=trim(catchment_def_file), form='formatted', action='read') - - read (10,*) tmpint1 - - if (logit) write (logunit,*) 'file contains coordinates for ', tmpint1, ' tiles' - if (logit) write (logunit,*) - - if (N_tile/=tmpint1) then - print*,"need :", N_tile,"but have: ",tmpint1 - err_msg = 'tile_coord_file and catchment_def_file mismatch. (1)' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - - do i=1,N_tile - - if (sweep==1) then - - ! read 7 columns, avoid using exact format specification - - read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & - tile_coord(i)%min_lon, & - tile_coord(i)%max_lon, & - tile_coord(i)%min_lat, & - tile_coord(i)%max_lat, & - tile_coord(i)%elev - - else - - ! read 6 columns, avoid using exact format specification - - read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & - tile_coord(i)%min_lon, & - tile_coord(i)%max_lon, & - tile_coord(i)%min_lat, & - tile_coord(i)%max_lat - - tile_coord(i)%elev = nodata_generic - - end if - - if (istat/=0) then ! read error - - if (sweep==1) then - - close(10,status='keep') - - if (logit) write (logunit,*) 'sweep 1 failed, trying sweep 2' - - exit ! exit sweep 1, try sweep 2 - - else - - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'sweep 2 failed') - - end if - - end if - - if (i==N_tile) then ! reached end of tile loop w/o read error - - close(10,status='keep') - - if (logit) write (logunit,*) 'sweep ', sweep, 'successfully completed' - - return - - end if - - end do ! loop through tiles - - end do ! loop through sweeps - - if ( any(tile_coord(1:N_tile)%tile_id/=tmp_tileid) .or. & - any(tile_coord(1:N_tile)%pfaf /=tmp_pfaf) ) then - - err_msg = 'tile_coord_file and catchment_def_file mismatch. (2)' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end if - - ! ----------------------------------------------------------------- - - end subroutine read_catchment_def - - ! ******************************************************************* - subroutine get_tile_num_from_latlon(N_catd, tile_coord, & tile_grid, N_tile_in_cell_ij, tile_num_in_cell_ij, N_latlon, lat, lon, & tile_num, max_dist_x, max_dist_y )