use cleaned-up EASE grid tools (#586)

Use cleaned-up make_bcs utilities (EASE grid tools) from GEOS-ESM/GEOSgcm_GridComp#634

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
 

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})

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

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)

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
 

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

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

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
 

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})