diff --git a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 index fd924e67..78c1dcdb 100644 --- a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 +++ b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 @@ -2898,7 +2898,7 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la 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,'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)) @@ -2949,7 +2949,10 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la 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 + + ! 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 diff --git a/src/Components/GEOSldas_GridComp/CMakeLists.txt b/src/Components/GEOSldas_GridComp/CMakeLists.txt index 65d57024..2f96b712 100644 --- a/src/Components/GEOSldas_GridComp/CMakeLists.txt +++ b/src/Components/GEOSldas_GridComp/CMakeLists.txt @@ -13,5 +13,5 @@ esma_add_library(${this} SRCS GEOS_LdasGridComp.F90 SUBCOMPONENTS ${alldirs} SUBDIRS Shared - DEPENDENCIES GEOSland_GridComp MAPL + DEPENDENCIES GEOSland_GridComp raster MAPL INCLUDES ${INC_ESMF}) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index bca85a78..2401e935 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -15,7 +15,7 @@ module GEOS_LdasGridCompMod use GEOS_EnsGridCompMod, only: EnsSetServices => SetServices use GEOS_LandAssimGridCompMod, only: LandAssimSetServices => SetServices - use LDAS_EASE_conv, only: ease_inverse + 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 use LDAS_TileCoordRoutines, only: get_tile_grid, get_ij_ind_from_latlon, io_domain_files diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt index cfccafa6..89684a1e 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt @@ -16,7 +16,7 @@ find_package(HDF5 REQUIRED COMPONENTS Fortran) esma_add_library (${this} SRCS ${SRCS} SUBCOMPONENTS ${alldirs} - DEPENDENCIES GEOS_LdasShared GEOSens_GridComp GEOSlandpert_GridComp GEOSland_GridComp MAPL GMAO_gfio_r4 hdf5hl_fortran hdf5_fortran ${NETCDF_LIBRARIES} + DEPENDENCIES GEOS_LdasShared GEOSens_GridComp GEOSlandpert_GridComp GEOSland_GridComp raster MAPL GMAO_gfio_r4 hdf5hl_fortran hdf5_fortran ${NETCDF_LIBRARIES} INCLUDES ${INC_ESMF} ${INC_HDF5}) target_compile_definitions (${this} PRIVATE LDAS_MPI) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 72e7d86e..d72d3521 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1348,7 +1348,7 @@ subroutine Initialize(gc, import, export, clock, rc) call MAPL_GetResource ( MAPL, GridName, Label="GEOSldas.GRIDNAME:", DEFAULT="EASE", RC=STATUS) _VERIFY(STATUS) _ASSERT( (NUM_ENSEMBLE>1), "out_smapL4SMaup=.true. only works for NUM_ENSEMBLE>1") - _ASSERT( (index(GridName,"EASEv2-M09") /=0), "out_smapL4SMaup=.true. only works with EASEv2-M09 tile space") + _ASSERT( (index(GridName,"EASEv2-M09") /=0 .or. index(GridName,"EASEv2_M09") /=0), "out_smapL4SMaup=.true. only works with EASEv2-M09 tile space") end if diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 6736c503..764cfada 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -73,9 +73,8 @@ module clsm_ensupd_enkf_update use nr_ran2_gasdev, ONLY: & NRANDSEED - use LDAS_ease_conv, ONLY: & - easeV1_convert, & - easeV2_convert + use ease_conv, ONLY: & + ease_convert use my_matrix_functions, ONLY: & row_std @@ -2100,14 +2099,14 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & character(len=*), parameter :: Iam = 'write_smapL4SMaup' character(len=400) :: err_msg + character(len=10) :: gridname_tmp ! -------------------------------------------------------------- ! ! smapL4SMaup output only works for 9 km EASE grids - if ( (trim(tile_grid_g%gridtype)/='EASE_M09' ) .and. & - (trim(tile_grid_g%gridtype)/='EASEv2_M09') ) then - err_msg = 'out_smapL4SMaup requires tile-space for 9 km EASE[v2] grid' + if ( index(tile_grid_g%gridtype, 'M09') == 0 ) then + err_msg = 'out_smapL4SMaup requires tile-space for 9 km EASEv1 or EASEv2 grid' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if @@ -2228,17 +2227,11 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMAP_L1C_Tbh_E09_D', 'SMAP_L1C_Tbv_E09_D', & 'SMAP_L1C_Tbh_E09_A', 'SMAP_L1C_Tbv_E09_A' & ) - - if (trim(tile_grid_g%gridtype)=='EASE_M09') then - - call easeV1_convert('M09', this_lat, this_lon, col_ind, row_ind) - - elseif (trim(tile_grid_g%gridtype)=='EASEv2_M09') then - - call easeV2_convert('M09', this_lat, this_lon, col_ind, row_ind) - - end if - + + if (index(tile_grid_g%gridtype, 'M09') /=0) then + call ease_convert(trim(tile_grid_g%gridtype), this_lat, this_lon, col_ind, row_ind) + endif + ! col_ind and row_ind are zero-based, need one-based index here col_beg_9km(n) = nint(col_ind)+1 @@ -2253,16 +2246,10 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & case('SMAP_L1C_Tbh_E27_D', 'SMAP_L1C_Tbv_E27_D', & 'SMAP_L1C_Tbh_E27_A', 'SMAP_L1C_Tbv_E27_A' & ) - - if (trim(tile_grid_g%gridtype)=='EASE_M09') then - - call easeV1_convert('M09', this_lat, this_lon, col_ind, row_ind) - - elseif (trim(tile_grid_g%gridtype)=='EASEv2_M09') then - - call easeV2_convert('M09', this_lat, this_lon, col_ind, row_ind) - - end if + + if (index(tile_grid_g%gridtype, 'M09') /=0) then + call ease_convert(trim(tile_grid_g%gridtype), this_lat, this_lon, col_ind, row_ind) + endif ! col_ind and row_ind are zero-based, need one-based index here ! L1C E27 spacing is one every three in each direction (~27-km spacing) @@ -2282,16 +2269,12 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMOS_fit_Tbh_A', 'SMOS_fit_Tbv_A' & ) - if (trim(tile_grid_g%gridtype)=='EASE_M09') then - - call easeV1_convert('M36', this_lat, this_lon, col_ind, row_ind) - - elseif (trim(tile_grid_g%gridtype)=='EASEv2_M09') then - - call easeV2_convert('M36', this_lat, this_lon, col_ind, row_ind) - - end if - + if (index(tile_grid_g%gridtype, 'M09') /=0) then + ! subindex (1:7) to get the string EASEvx_ + gridname_tmp = tile_grid_g%gridtype(1:7)//'M36' + call ease_convert(gridname_tmp, this_lat, this_lon, col_ind, row_ind) + endif + ! col_ind and row_ind are zero-based, need one-based index here col_beg_9km(n) = nint(col_ind) *4 + 1 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 0f892960..5937a6e0 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -21,9 +21,9 @@ module clsm_ensupd_read_obs use io_hdf5, ONLY: & hdf5read - use LDAS_ease_conv, ONLY: & - easeV2_convert, & - easeV2_extent + use EASE_conv, ONLY: & + ease_convert, & + ease_extent use LDAS_ensdrv_globals, ONLY: & logit, & @@ -4436,12 +4436,12 @@ subroutine read_obs_SMOS( date_time, N_catd, this_obs_param, & if (tmp_tile_num(ii)>0) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M36_col_ind_tile, M36_row_ind_tile ) - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & tmp_lat(ii), & tmp_lon(ii), & M36_col_ind_obs, M36_row_ind_obs ) @@ -5127,12 +5127,12 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & if (tmp_tile_num(ii)>0) then - call easeV2_convert('M09', & + call ease_convert('EASEv2_M09', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M09_col_ind_tile, M09_row_ind_tile ) - call easeV2_convert('M09', & + call ease_convert('EASEv2_M09', & tmp_lat(ii), & tmp_lon(ii), & M09_col_ind_obs, M09_row_ind_obs ) @@ -6118,7 +6118,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! shift M36 obs lat/lon for proper assignment of M09 tile? - if ( L1C_files .and. (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0) ) then + if ( L1C_files .and. (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0 .or. index(tile_grid_d%gridtype, 'EASEv2-M09') /=0 )) then ! temporarily shift lat/lon of obs for computation of nearest tile to ! avoid ambiguous assignment of M09 model tile within M36 obs grid cell @@ -6185,12 +6185,12 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & if (tmp_tile_num(ii)>0) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M36_col_ind_tile, M36_row_ind_tile ) - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & tmp_lat(ii), & tmp_lon(ii), & M36_col_ind_obs, M36_row_ind_obs ) @@ -6486,7 +6486,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation ! ! assemble 36 km EASEv2 mask of L2AP_Tb obs - call easeV2_extent( 'M36', N_cols, N_rows ) + call ease_extent( 'EASEv2_M36', N_cols, N_rows ) allocate( mask_h_A(N_cols,N_rows) ) allocate( mask_h_D(N_cols,N_rows) ) @@ -6508,7 +6508,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_f(ii)%species==species_L2AP_Tbh_A) then - call easeV2_convert('M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; @@ -6530,7 +6530,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_f(ii)%species==species_L2AP_Tbh_D) then - call easeV2_convert('M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; @@ -6552,7 +6552,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_f(ii)%species==species_L2AP_Tbv_A) then - call easeV2_convert('M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; @@ -6574,7 +6574,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_f(ii)%species==species_L2AP_Tbv_D) then - call easeV2_convert('M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; @@ -6609,7 +6609,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_l(ii)%species==species_L1C_Tbh_A) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) ! note conversion to one-based indices @@ -6628,7 +6628,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_l(ii)%species==species_L1C_Tbh_D) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) ! note conversion to one-based indices @@ -6649,7 +6649,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_l(ii)%species==species_L1C_Tbv_A) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) ! note conversion to one-based indices @@ -6668,7 +6668,7 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation if (Observations_l(ii)%species==species_L1C_Tbv_D) then - call easeV2_convert('M36', & + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) ! note conversion to one-based indices diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 6f88b864..d102cdfc 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -1707,9 +1707,7 @@ subroutine get_obs_pred( & ! for EASE grids *ONLY*: screen for non-land surfaces (e.g., lakes) ! - reichle, 28 Mar 2015 - if ( & - (index(tile_grid_g%gridtype, 'EASE_M') /=0) .or. & - (index(tile_grid_g%gridtype, 'EASEv2_M')/=0) ) then + if (index(tile_grid_g%gridtype, 'EASEv') /=0) then ! ASSUMPTIONS: ! - at most one land tile per grid cell @@ -3204,8 +3202,7 @@ subroutine get_obs_pert( N_ens, N_obs, N_obs_param, & pert_grid_lH%ur_lon = pert_grid_lH%ll_lon + delta_lon pert_grid_lH%ur_lat = pert_grid_lH%ll_lat + delta_lat - elseif ( (index(pert_grid_lH%gridtype,'EASE_M') /=0) .or. & - (index(pert_grid_lH%gridtype,'EASEv2_M')/=0) ) then + elseif ( index(pert_grid_lH%gridtype,'EASEv') /=0 ) then pert_grid_lH%dlon = pert_grid_f%dlon diff --git a/src/Components/GEOSldas_GridComp/Shared/CMakeLists.txt b/src/Components/GEOSldas_GridComp/Shared/CMakeLists.txt index da7f7ec6..a18409aa 100644 --- a/src/Components/GEOSldas_GridComp/Shared/CMakeLists.txt +++ b/src/Components/GEOSldas_GridComp/Shared/CMakeLists.txt @@ -4,7 +4,7 @@ set (SRCS enkf_types.F90 catch_types.F90 LDAS_ensdrv_Globals.F90 LDAS_DateTimeMod.F90 LDAS_DriverTypes.F90 LDAS_Convert.F90 LDAS_Exceptions.F90 LDAS_TileCoordType.F90 LDAS_PertTypes.F90 LDAS_ensdrv_functions.F90 my_matrix_functions.F90 - LDAS_EASE_conv.F90 LDAS_TileCoordRoutines.F90 + LDAS_TileCoordRoutines.F90 LDAS_RepairForcing.F90 LDAS_ensdrv_mpi.F90 ) @@ -15,5 +15,5 @@ list (APPEND SRCS esma_add_library(${this} SRCS ${SRCS} - DEPENDENCIES MAPL GEOS_Shared GEOS_LandShared + DEPENDENCIES MAPL GEOS_Shared GEOS_LandShared raster INCLUDES ${INC_ESMF} ${INC_NETCDF}) diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_EASE_conv.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_EASE_conv.F90 deleted file mode 100644 index 9aaf636f..00000000 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_EASE_conv.F90 +++ /dev/null @@ -1,675 +0,0 @@ - -module LDAS_EASE_conv - - ! ========================================================================== - ! - ! easeV1_conv.F90 - FORTRAN routines for conversion of azimuthal - ! equal area and equal area cylindrical grid coordinates - ! - ! 30-Jan-1992 H.Maybee - ! 20-Mar-1992 Ken Knowles 303-492-0644 knowles@kryos.colorado.edu - ! 16-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu - ! Copied from nsmconv.f, changed resolutions from - ! 40-20-10 km to 25-12.5 km - ! 21-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu - ! Fixed sign of Southern latitudes in ease_inverse. - ! 12-Sep-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu - ! Changed grid cell size. Changed "c","f" to "l","h" - ! 25-Oct-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu - ! Changed row size from 587 to 586 for Mercator projection - ! 11-May-2011 reichle: Changed "smap" to "easeV1". - ! Added SSM/I and AMSR-E "M25" grid. - ! So far ONLY for cylindrical grids. - ! Converted from *.f to *.F90 module - ! - ! $Log$ - ! Revision 1.1.2.3 2018/09/13 20:42:50 wjiang - ! change M25 - ! - ! Revision 1.1.2.2 2017/09/18 15:10:25 wjiang - ! "fix" strange compiler erro on comment line. - ! - ! Revision 1.1.2.1 2017/01/19 19:35:58 wjiang - ! 1)add EASE grid support - ! 2)add ensemble average for HISTORY - ! - ! Revision 1.2 2014-08-26 17:33:55 rreichle - ! - clean-up of *.F90 in GEOSlana_GridComp: - ! - make sure all modules include "private" statement - ! - follow all "use" statements with "ONLY" - ! - removed unused files (esat_qsat.F90, nr_sort.f) - ! - removed unused variables - ! - ! Revision 1.1 2011-05-11 21:58:46 rreichle - ! - ! Adding utilities to map between EASE grids and lat/lon coordinates. - ! - ! Revision 1.3 1994/11/01 23:40:43 brodzik - ! Replaced all references to 'ease' with 'smap' - ! Replaced all references to 'smap' with 'easeV1' -- reichle - ! - ! ========================================================================== - - implicit none - - private - - public :: ease_convert - public :: ease_inverse - public :: easeV1_convert - public :: easeV1_inverse - public :: easeV2_convert - public :: easeV2_inverse - public :: easeV2_extent - - - ! ***NEVER*** change these constants to GEOS-5 MAPL constants!!!! - - ! radius of the earth (km), authalic sphere based on International datum - - real*8, parameter :: RE_km = 6371.228 - - ! scale factor for standard paralles at +/-30.00 degrees - - real*8, parameter :: COS_PHI1 = .866025403 - - real*8, parameter :: PI = 3.14159265358979323846 - - ! ========================================================================== - ! - ! easeV2_conv.F90 - FORTRAN routines for converting grid coordinates - ! (latitude/longitude <--> row/column indices) - ! of the Equal Area Scalable Earth, version 2 (EASEv2) grid - ! - ! ***** ONLY cylindrical ('M') projection implemented ***** - ! - ! Ported from Steven Chan's matlab code (smapease2inverse.m, - ! smapease2forward.m), which has been ported from NSIDC's IDL code - ! (wgs84_convert.pro, wgs84_inverse.pro) available from - ! ftp://sidads.colorado.edu/pub/tools/easegrid/geolocation_tools/ - ! - ! 04-Apr-2013 - reichle - ! - ! Official references: - ! doi:10.3390/ijgi1010032 - ! doi:10.3390/ijgi3031154 -- correction of M25 "map_scale_m" parameters! - ! - ! 04-Apr-2013 - reichle - ! 11-Sep-2018 - reichle, mgirotto -- added 'M25' grid parameters - ! - ! ========================================================================== - - - ! ***NEVER*** change these constants to GEOS-5 MAPL constants!!!! - - ! radius of the earth (m) and map eccentricity - - real*8, parameter :: map_equatorial_radius_m = 6378137.0 - - real*8, parameter :: map_eccentricity = 0.081819190843 - - - real*8, parameter :: e2 = map_eccentricity * map_eccentricity - real*8, parameter :: e4 = e2 * e2 - real*8, parameter :: e6 = e2 * e4 - - real*8, parameter :: epsilon = 1.e-6 - - real*8, parameter :: map_reference_longitude = 0.0 ! 'M', 'N', 'S' - - ! constants for 'N' and 'S' (azimuthal) projections - - real*8, parameter :: N_map_reference_latitude = 90.0 - real*8, parameter :: S_map_reference_latitude = -90.0 - - ! constants for 'M' (cylindrical) projection - - real*8, parameter :: M_map_reference_latitude = 0.0 - real*8, parameter :: M_map_second_reference_latitude = 30.0 - - real*8, parameter :: M_sin_phi1 = sin(M_map_second_reference_latitude*PI/180.) - real*8, parameter :: M_cos_phi1 = cos(M_map_second_reference_latitude*PI/180.) - - real*8, parameter :: M_kz = M_cos_phi1/sqrt(1.0-e2*M_sin_phi1*M_sin_phi1) - - -contains - - subroutine ease_convert (gridname, lat, lon, r, s) - character*(*), intent(in) :: gridname - real, intent(in) :: lat, lon - real, intent(out) :: r, s - character(3) :: grid - - if (index(gridname,'M36') /=0 ) then - grid='M36' - else if (index(gridname,'M25') /=0 ) then - grid='M25' - else if (index(gridname,'M09') /=0 ) then - grid='M09' - else if (index(gridname,'M03') /=0 ) then - grid='M03' - else if (index(gridname,'M01') /=0 ) then - grid='M01' - endif - - if(index(gridname,'EASEv2') /=0) then - call easeV2_convert(grid,lat,lon,r,s) - else if(index(gridname,'EASE') /=0) then - call easeV1_convert(grid,lat,lon,r,s) - else - print*,"wrong gridname: "//gridname - endif - end subroutine - - subroutine ease_inverse (gridname, r, s, lat, lon) - character*(*), intent(in) :: gridname - real, intent(in) :: r, s - real, intent(out) :: lat, lon - character(3) :: grid - - if (index(gridname,'M36') /=0 ) then - grid='M36' - else if (index(gridname,'M25') /=0 ) then - grid='M25' - else if (index(gridname,'M09') /=0 ) then - grid='M09' - else if (index(gridname,'M03') /=0 ) then - grid='M03' - else if (index(gridname,'M01') /=0 ) then - grid='M01' - endif - - if(index(gridname,'EASEv2') /=0) then - call easeV2_inverse(grid,r,s,lat,lon) - else if(index(gridname,'EASE') /=0) then - call easeV1_inverse(grid,r,s,lat,lon) - else - print*,"wrong gridname: "//gridname - endif - end subroutine ease_inverse - - ! ******************************************************************* - - subroutine easeV1_convert (grid, lat, lon, r, s) - - ! convert geographic coordinates (spherical earth) to - ! azimuthal equal area or equal area cylindrical grid coordinates - ! - ! status = easeV1_convert (grid, lat, lon, r, s) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! or xx = "12", "25" (SSM/I, AMSR-E) - ! lat, lon = geo. coords. (decimal degrees) - ! - ! output: r, s - column, row coordinates - ! - ! result: status = 0 indicates normal successful completion - ! -1 indicates error status (point not on grid) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: lat, lon - real, intent(out) :: r, s - - ! local variables - - integer :: cols, rows - real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 - - ! --------------------------------------------------------------------- - - call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - phi = lat*PI/180. ! convert from degree to radians - lam = lon*PI/180. ! convert from degree to radians - - if (grid(1:1).eq.'N') then - rho = 2 * Rg * sin(PI/4. - phi/2.) - r = r0 + rho * sin(lam) - s = s0 + rho * cos(lam) - - else if (grid(1:1).eq.'S') then - rho = 2 * Rg * cos(PI/4. - phi/2.) - r = r0 + rho * sin(lam) - s = s0 - rho * cos(lam) - - else if (grid(1:1).eq.'M') then - r = r0 + Rg * lam * COS_PHI1 - s = s0 - Rg * sin(phi) / COS_PHI1 - - endif - - end subroutine easeV1_convert - - ! ******************************************************************* - - subroutine easeV1_inverse (grid, r, s, lat, lon) - - ! convert azimuthal equal area or equal area cylindrical - ! grid coordinates to geographic coordinates (spherical earth) - ! - ! status = easeV1_inverse (grid, r, s, lat, lon) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! or xx = "12", "25" (SSM/I, AMSR-E) - ! r, s - column, row coordinates - ! - ! output: lat, lon = geo. coords. (decimal degrees) - ! - ! result: status = 0 indicates normal successful completion - ! -1 indicates error status (point not on grid) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: r, s - real, intent(out) :: lat, lon - - ! local variables - - integer :: cols, rows - real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 - real*8 :: gamma, beta, epsilon, x, y, c - real*8 :: sinphi1, cosphi1 - - ! --------------------------------------------------------------------- - - call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - x = r - r0 - y = -(s - s0) - - if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - rho = sqrt(x*x + y*y) - if (rho.eq.0.0) then - if (grid(1:1).eq.'N') lat = 90.0 - if (grid(1:1).eq.'S') lat = -90.0 - lon = 0.0 - else - if (grid(1:1).eq.'N') then - sinphi1 = sin(PI/2.) - cosphi1 = cos(PI/2.) - if (y.eq.0.) then - if (r.le.r0) lam = -PI/2. - if (r.gt.r0) lam = PI/2. - else - lam = atan2(x,-y) - endif - else if (grid(1:1).eq.'S') then - sinphi1 = sin(-PI/2.) - cosphi1 = cos(-PI/2.) - if (y.eq.0.) then - if (r.le.r0) lam = -PI/2. - if (r.gt.r0) lam = PI/2. - else - lam = atan2(x,y) - endif - endif - gamma = rho/(2 * Rg) - if (abs(gamma) .gt. 1.) return - c = 2 * asin(gamma) - beta = cos(c) * sinphi1 + y * sin(c) * (cosphi1/rho) - if (abs(beta).gt.1.) return - phi = asin(beta) - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI ! convert from radians to degree - endif - - else if (grid(1:1).eq.'M') then - - ! allow .5 cell tolerance in arcsin function - ! so that grid coordinates which are less than .5 cells - ! above 90.00N or below 90.00S are given a lat of 90.00 - - epsilon = 1 + 0.5/Rg - beta = y*COS_PHI1/Rg - if (abs(beta).gt.epsilon) return - if (beta.le.-1.) then - phi = -PI/2. - else if (beta.ge.1.) then - phi = PI/2. - else - phi = asin(beta) - endif - lam = x/COS_PHI1/Rg - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI ! convert from radians to degree - endif - - end subroutine easeV1_inverse - - ! ******************************************************************* - - subroutine easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - implicit none - - character*(*), intent(in) :: grid - real*8, intent(out) :: CELL_km, r0, s0, Rg - integer, intent(out) :: cols, rows - - ! -------------------------------------------------------- - ! - ! r0,s0 are defined such that cells at all scales have - ! coincident center points - ! - !c r0 = (cols-1)/2. * scale - !c s0 = (rows-1)/2. * scale - ! - ! -------------------------------------------------------- - - if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - - print *,'Polar projections not implemented yet' - stop - - else if (grid(1:1).eq.'M') then - - if (grid .eq. 'M36') then ! SMAP 36 km grid - CELL_km = 36.00040279063 ! nominal cell size in kilometers - cols = 963 - rows = 408 - r0 = 481.0 - s0 = 203.5 - - else if (grid .eq. 'M25') then ! SSM/I, AMSR-E 25 km grid - CELL_km = 25.067525 ! nominal cell size in kilometers - cols = 1383 - rows = 586 - r0 = 691.0 - s0 = 292.5 - - else if (grid .eq. 'M09') then ! SMAP 9 km grid - CELL_km = 9.00010069766 ! nominal cell size in kilometers - cols = 3852 - rows = 1632 - r0 = 1925.5 - s0 = 815.5 - - else if (grid .eq. 'M03') then ! SMAP 3 km grid - CELL_km = 3.00003356589 ! nominal cell size in kilometers - cols = 11556 - rows = 4896 - r0 = 5777.5 - s0 = 2447.5 - - else if (grid .eq. 'M01') then ! SMAP 1 km grid - CELL_km = 1.00001118863 ! nominal cell size in kilometers - cols = 34668 - rows = 14688 - r0 = 17333.5 - s0 = 7343.5 - - else - - print *,'easeV1_convert: unknown resolution: ',grid - stop - - endif - - else - - print *, 'easeV1_convert: unknown projection: ', grid - stop - - endif - - Rg = RE_km/CELL_km - - end subroutine easeV1_get_params - - - - subroutine easeV2_convert (grid, lat, lon, col_ind, row_ind) - - ! convert geographic coordinates (spherical earth) to - ! azimuthal equal area or equal area cylindrical grid coordinates - ! - ! *** NOTE order of calling arguments: "lat-lon-lon-lat" *** - ! - ! useage: call easeV2_convert (grid, lat, lon, r, s) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! lat, lon = geo. coords. (decimal degrees) - ! - ! output: col_ind, row_ind - column, row coordinates - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: lat, lon - real, intent(out) :: col_ind, row_ind - - ! local variables - - integer :: cols, rows - real*8 :: dlon, phi, lam, map_scale_m, r0, s0, ms, x, y, sin_phi, q - - ! --------------------------------------------------------------------- - - call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - dlon = lon - - if (abs(map_reference_longitude)>epsilon) then - - dlon = lon - map_reference_longitude - - end if - - if (dlon .lt. -180.0) dlon = dlon + 360.0 - if (dlon .gt. 180.0) dlon = dlon - 360.0 - - phi = lat*PI/180. ! convert from degree to radians - lam = dlon*PI/180. ! convert from degree to radians - - sin_phi = sin(phi) - - ms = map_eccentricity*sin_phi - - q = (1. - e2)* & - ( & - (sin_phi /(1. - e2*sin_phi*sin_phi)) & - - & - .5/map_eccentricity*log((1.-ms)/(1.+ms)) & - ) - - ! note: "qp" only needed for 'N' and 'S' projections - - if (grid(1:1).eq.'M') then - - x = map_equatorial_radius_m*M_kz*lam - - y = (map_equatorial_radius_m*q)/(2.*M_kz) - - else - - print *,'Polar projections not implemented yet' - stop - - endif - - row_ind = s0 - (y/map_scale_m) - col_ind = r0 + (x/map_scale_m) - - end subroutine easeV2_convert - - ! ******************************************************************* - - subroutine easeV2_inverse (grid, r, s, lat, lon) - - ! convert azimuthal equal area or equal area cylindrical - ! grid coordinates to geographic coordinates (spherical earth) - ! - ! *** NOTE order of calling arguments: "lon-lat-lat-lon" *** - ! - ! useage: call easeV1_inverse (grid, r, s, lat, lon) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! r, s - column, row coordinates - ! - ! output: lat, lon = geo. coords. (decimal degrees) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: r, s - real, intent(out) :: lat, lon - - ! local variables - - integer :: cols, rows - real*8 :: phi, lam, map_scale_m, r0, s0, beta, x, y, qp - - ! --------------------------------------------------------------------- - - call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - x = (r - r0)*map_scale_m - y = -(s - s0)*map_scale_m - - qp = (1. - e2)* & - ( & - (1./(1.-e2)) & - - & - .5/map_eccentricity*log((1.-map_eccentricity)/(1.+map_eccentricity)) & - ) - - if (grid(1:1).eq.'M') then - - beta = asin(2.*y*M_kz/(map_equatorial_radius_m*qp)) - - lam = x/(map_equatorial_radius_m*M_kz) - - else - - print *,'Polar projections not implemented yet' - stop - - endif - - phi = beta & - + ( ( e2/3. + 31./180.*e4 + 517./ 5040.*e6 )*sin(2.*beta) ) & - + ( ( 23./360.*e4 + 251./ 3780.*e6 )*sin(4.*beta) ) & - + ( ( 761./45360.*e6 )*sin(6.*beta) ) - - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI + map_reference_longitude ! convert from radians to degree - - if (lon .lt. -180.0) lon = lon + 360.0 - if (lon .gt. 180.0) lon = lon - 360.0 - - end subroutine easeV2_inverse - - ! ******************************************************************* - - subroutine easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - implicit none - - character*(*), intent(in) :: grid - real*8, intent(out) :: map_scale_m, r0, s0 - integer, intent(out) :: cols, rows - - - if (grid(1:1).eq.'M') then - - if (grid .eq. 'M36') then ! SMAP 36 km grid - - map_scale_m = 36032.220840584 ! nominal cell size in meters - cols = 964 - rows = 406 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - - else if (grid .eq. 'M25') then ! 25 km grid - - map_scale_m = 25025.2600000 ! nominal cell size in meters (see doi:10.3390/ijgi3031154) - cols = 1388 - rows = 584 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2 - - else if (grid .eq. 'M09') then ! SMAP 9 km grid - - map_scale_m = 9008.055210146 ! nominal cell size in meters - cols = 3856 - rows = 1624 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M03') then ! SMAP 3 km grid - - map_scale_m = 3002.6850700487 ! nominal cell size in meters - cols = 11568 - rows = 4872 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M01') then ! SMAP 1 km grid - - map_scale_m = 1000.89502334956 ! nominal cell size in meters - cols = 34704 - rows = 14616 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else - - print *,'easeV2_convert: unknown resolution: ',grid - stop - - endif - - else if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - - print *,'Polar projections not implemented yet' - stop - - else - - print *, 'easeV2_convert: unknown projection: ', grid - stop - - endif - - end subroutine easeV2_get_params - - ! ******************************************************************* - - subroutine easeV2_extent( grid, N_cols, N_rows ) - - ! simple wrapper to get N_cols (N_lon) and N_rows (N_lat) - - implicit none - - character*(*), intent(in) :: grid - integer, intent(out) :: N_cols, N_rows - - ! local variables - - real*8 :: map_scale_m, r0, s0 - - ! ------------------------------------------------ - - call easeV2_get_params( grid, map_scale_m, N_cols, N_rows, r0, s0 ) - - end subroutine easeV2_extent - - ! ******************************************************************* - -end module LDAS_EASE_conv - -! =============================== EOF ================================= - diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index 1b2332c0..a059ed1b 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -26,8 +26,10 @@ module LDAS_TileCoordRoutines use MAPL_ConstantsMod, ONLY: & MAPL_RADIUS ! Earth radius - use LDAS_EASE_conv, ONLY: & - easev1_convert,easev2_convert + use EASE_conv, ONLY: & + ease_convert, & + ease_inverse, & + ease_extent use LDAS_ExceptionsMod, ONLY: & ldas_abort, & @@ -279,7 +281,7 @@ subroutine LDAS_create_grid_g( gridname, n_lon, n_lat, & ! ! outputs: ! offsets - ! cell_area (optional, for EASE grids only) + ! cell_area [m^2] (optional, for EASE grids only) implicit none @@ -295,7 +297,7 @@ subroutine LDAS_create_grid_g( gridname, n_lon, n_lat, & logical :: date_line_on_center, pole_on_center logical :: ease_grid, c3_grid, latlon_grid logical :: file_exist - + integer :: k, rows, cols character(len=*), parameter :: Iam = 'create global ldas_grid ' character(len=400) :: err_msg @@ -330,7 +332,7 @@ subroutine LDAS_create_grid_g( gridname, n_lon, n_lat, & pole_on_center = .true. endif - if (index(gridname,"EASE") /=0) then + if (index(gridname,"EASEv") /=0) then ease_grid = .true. latlon_grid = .false. endif @@ -360,70 +362,33 @@ subroutine LDAS_create_grid_g( gridname, n_lon, n_lat, & if (ease_grid) then + ! gridname may be EASEv2-M36 or EASEv2_M36 (to be cleaned up later) + + k = index(gridname, 'EASEv') + + if (k == 0) then + err_msg = 'unknown EASE grid tile defs, gridname = ' // trim( gridname) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + tile_grid%gridtype = trim(gridname(k:)) + tile_grid%ind_base = 0 ! global cylindrical EASE grid - tile_grid%ll_lon = -180. - tile_grid%ur_lon = 180. - tile_grid%i_dir = +1 tile_grid%j_dir = -1 - ! It is sloppy that the name may be EASEv2-M36 or EASEv2_M36 - - if (index(gridname, 'EASEv2_M36')/=0 .or. index(gridname, 'EASEv2-M36')/=0) then ! version *2* - - tile_grid%gridtype = 'EASEv2_M36' - - tile_grid%ll_lat = -85.04456 - tile_grid%ur_lat = 85.04456 - - ease_cell_area = 1298.320938704616 - - elseif (index(gridname, 'EASEv2_M09')/=0 .or. index(gridname, 'EASEv2-M09')/=0) then ! version *2* - - tile_grid%gridtype = 'EASEv2_M09' - - tile_grid%ll_lat = -85.04456 - tile_grid%ur_lat = 85.04456 - - ease_cell_area = 81.145058669038477 - - elseif (index(gridname, 'EASE_M36')/=0 .or. index(gridname, 'EASE-M36')/=0) then - - tile_grid%gridtype = 'EASE_M36' - - tile_grid%ll_lat = -86.62256 ! minimal change, reichle, 5 Apr 2013 - tile_grid%ur_lat = 86.62256 ! minimal change, reichle, 5 Apr 2013 - - ease_cell_area = 1296.029001087600 - - elseif (index(gridname, 'EASE_M09')/=0 .or. index(gridname, 'EASE-M09')/=0) then - - tile_grid%gridtype = 'EASE_M09' - - tile_grid%ll_lat = -86.62256 ! minimal change, reichle, 5 Apr 2013 - tile_grid%ur_lat = 86.62256 ! minimal change, reichle, 5 Apr 2013 - - ease_cell_area = 81.001812568020028 - - elseif (index(gridname, 'EASE_M25')/=0 .or. index(gridname, 'EASE-M25')/=0 ) then - - tile_grid%gridtype = 'EASE_M25' - - tile_grid%ll_lat = -86.7167 ! need to double-check (reichle, 11 May 2011) - tile_grid%ur_lat = 86.7167 ! need to double-check (reichle, 11 May 2011) - - ease_cell_area = 628.38080962 - - else - - err_msg = 'unknown EASE grid tile defs, grid name = ' // trim( gridname) - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end if - + call ease_extent ( & + gridname, cols, rows, & + cell_area = ease_cell_area, & ! [m^2] + ll_lon = tile_grid%ll_lon, & + ll_lat = tile_grid%ll_lat, & + ur_lon = tile_grid%ur_lon, & + ur_lat = tile_grid%ur_lat & + ) + tile_grid%dlon = 360./real(tile_grid%N_lon) tile_grid%dlat = (tile_grid%ur_lat-tile_grid%ll_lat)/real(tile_grid%N_lat) ! *avg* dlat! @@ -616,8 +581,7 @@ subroutine get_tile_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc tile_grid%ur_lon = tile_grid%ll_lon + real(tile_grid%N_lon)*tile_grid%dlon tile_grid%ur_lat = tile_grid%ll_lat + real(tile_grid%N_lat)*tile_grid%dlat - elseif ( (index(tile_grid_g%gridtype,'EASE_M') /=0) .or. & - (index(tile_grid_g%gridtype,'EASEv2_M')/=0) ) then + elseif ( index(tile_grid_g%gridtype,'EASEv') /=0 ) then ! *average* dlat over the domain @@ -895,11 +859,9 @@ subroutine get_tile_num_from_latlon(N_catd, tile_coord, & ! map from i_ind, j_ind to tile_num - if ( & - (index(tile_grid%gridtype, 'EASE_M') /=0) .or. & - (index(tile_grid%gridtype, 'EASEv2_M')/=0) ) then + if ( index(tile_grid%gridtype, 'EASEv') /=0 ) then - ! ASSUMPTION: tiles match EASE or EASEv2 grid cells exactly + ! ASSUMPTION: tiles match EASE grid cells exactly ! (unless "outside" the domain, eg. water surface) if (N_tile_in_cell_ij(i_ind,j_ind)==1) then @@ -911,7 +873,7 @@ subroutine get_tile_num_from_latlon(N_catd, tile_coord, & elseif (N_tile_in_cell_ij(i_ind,j_ind)==0) then - ! Do nothing. If given EASE or EASEv2 grid cell is not land, + ! Do nothing. If given EASE grid cell is not land, ! tile_num will not change from its initialized value. else @@ -921,7 +883,6 @@ subroutine get_tile_num_from_latlon(N_catd, tile_coord, & end if - elseif (index(tile_grid%gridtype, 'LatLon')/=0) then ! Loop through all tiles within grid cell that contain lat/lon and @@ -1037,23 +998,11 @@ subroutine get_ij_ind_from_latlon( tile_grid, lat, lon, i_ind, j_ind ) ! ! select grid type - if (index(tile_grid%gridtype, 'EASE')/=0) then + if (index(tile_grid%gridtype, 'EASEv')/=0) then ! EASE grid lat/lon to index provides *global*, *0-based* index! - - if (index(tile_grid%gridtype, 'EASE_M') /=0) then - - call easeV1_convert( tile_grid%gridtype(6:8), lat, lon, r, s) - - elseif (index(tile_grid%gridtype, 'EASEv2_M')/=0) then - - call easeV2_convert( tile_grid%gridtype(8:10), lat, lon, r, s) - - else - - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown EASE grid type') - - end if + + call ease_convert(tile_grid%gridtype, lat, lon, r, s) i_indg = nint(r) ! i_ind or lon_ind j_indg = nint(s) ! j_ind or lat_ind diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index 66c68f19..414a9800 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -80,19 +80,19 @@ module LDAS_TileCoordType ! Possible grid types (structure field "gridtype"): ! ! - "LatLon" : regular lat/lon grid (constant dlon, dlat) - ! - "EASE_Mxx" : cylindrical EASE grid (constant dlon, variable dlat) + ! - "EASEv1_Mxx" : cylindrical EASEv1 grid (constant dlon, variable dlat) ! - "EASEv2_Mxx" : cylindrical EASEv2 grid (constant dlon, variable dlat) ! ! Grid orientation (convention): ! ! "LatLon" : 1-based indexing, SouthWest to NorthEast - ! "EASE_Mxx" : 0-based indexing, NorthWest to SouthEast + ! "EASEv1_Mxx" : 0-based indexing, NorthWest to SouthEast ! "EASEv2_Mxx" : 0-based indexing, NorthWest to SouthEast ! ! Grids are defined by the following fields: ! ! --------------------------------------------------------- - ! | || "LatLon" | "EASE_Mxx" | + ! | || "LatLon" | "EASEv1_Mxx" | ! | || | "EASEv2_Mxx" | ! --------------------------------------------------------- ! | indexing || ind_base | ind_base | @@ -123,7 +123,7 @@ module LDAS_TileCoordType ! any subroutines or operators defined herein ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - character(40) :: gridtype ! type of grid, eg. "LatLon", "EASE_M36", "EASEv2_M09", ... + character(40) :: gridtype ! type of grid, eg. "LatLon", "EASEv1_M36", "EASEv2_M09", ... integer :: ind_base ! 0=zero-based indices (EASE), 1=one-based indices (LatLon) integer :: i_dir ! direction of indices (+1=W-to-E, -1=E-to-W) integer :: j_dir ! direction of indices (+1=S-to-N, -1=N-to-S) @@ -133,7 +133,7 @@ module LDAS_TileCoordType integer :: j_offg ! minimum lat index (offset from *global* grid) ! ! "LatLon" : i_offg -> westernmost longitude ! ! j_offg -> southernmost latitude - ! ! "EASE_Mxx" : i_offg -> westernmost longitude + ! ! "EASEv1_Mxx": i_offg -> westernmost longitude ! ! j_offg -> northernmost latitude ! ! "EASEv2_Mxx": i_offg -> westernmost longitude ! ! j_offg -> northernmost latitude @@ -144,7 +144,7 @@ module LDAS_TileCoordType real :: dlon ! longitude grid spacing [deg] real :: dlat ! latitude grid spacing [deg] ! ! "LatLon" : constant - ! ! "EASE_Mxx" : *average* dlat over grid extent + ! ! "EASEv1_Mxx": *average* dlat over grid extent ! ! "EASEv2_Mxx": *average* dlat over grid extent !CLSM_ENSDRV_MPI is NOT updated. will not be saved and bcasted