diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 0403e969..0f95ac78 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -620,22 +620,43 @@ subroutine Initialize(gc, import, export, clock, rc) close(10) call io_grid_def_type('w', logunit, tile_grid_f, 'tile_grid_f') + ! get a grid for perturbations and EnKF: + ! + ! tile grid ! pert grid + ! (defines tile space) ! (used for perturbations and as "hash" grid in EnKF analysis) + ! =========================================================================================================== + ! lat/lon ! same as tile_grid (i.e., lat/lon) + ! ----------------------------------------------------------------------------------------------------------- + ! EASEv[X] ! same as tile_grid (i.e., EASE) + ! ----------------------------------------------------------------------------------------------------------- + ! cube-sphere ! lat/lon grid of resolution similar to that of (cube-sphere) tile_grid + pert_grid_g = get_pert_grid(tile_grid_g) - pert_grid_f = get_pert_grid(tile_grid_f) - - if(trim(grid_type) == "Cubed-Sphere" ) then + ! pert_grid_f = get_pert_grid(tile_grid_f) -- SEEMS REDUNDANT, OVERWRITTEN BELOW: pert_grid_f = get_minExtent_grid() + + !if(trim(grid_type) == "Cubed-Sphere" ) then + ! + ! _ASSERT(index(tile_grid_g%gridtype, 'c3') /=0, "tile_grid_g does not describe a cubed-sphere grid") + + if ( .not. (pert_grid_g==tile_grid_g) ) then - _ASSERT(index(tile_grid_g%gridtype, 'c3') /=0, "tile_grid_g does not describe a cubed-sphere grid") - - !1) get hash index for cubed-sphere + ! arrive here when tile_grid_g is cube-sphere and pert_grid_g is lat/lon after call to get_pert_grid() above + + !1) get pert_i_indg, pert_j_indg for tiles in (full) domain relative to pert_grid_g do i = 1, N_catf call get_ij_ind_from_latlon(pert_grid_g,tile_coord_f(i)%com_lat,tile_coord_f(i)%com_lon, & tile_coord_f(i)%pert_i_indg,tile_coord_f(i)%pert_j_indg) enddo - !2) re-generate pert_grid_f in Lat-Lon + !2) determine pert_grid_f pert_grid_f = get_minExtent_grid(N_catf, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, & tile_coord_f%min_lon, tile_coord_f%min_lat, tile_coord_f%max_lon, tile_coord_f%max_lat, & pert_grid_g) + + else + + ! do nothing because %pert_i_indg and %pert_j_indg were initialized to %i_indg and %j_indg + ! in io_tile_coord_type() when tile_coord was read via io_domain_files() + endif endif @@ -645,10 +666,10 @@ subroutine Initialize(gc, import, export, clock, rc) call MPI_BCAST(tile_coord_f,N_catf, MPI_tile_coord_type,0,mpicomm, mpierr) call MPI_BCAST(pert_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) call MPI_BCAST(pert_grid_f, 1, MPI_grid_def_type, 0,mpicomm, mpierr) - - if(trim(grid_type) == "Cubed-Sphere" ) then - call MPI_BCAST(tile_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) - endif + + !if(trim(grid_type) == "Cubed-Sphere" ) then + call MPI_BCAST(tile_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) + !endif block integer, allocatable :: f2tile_id(:), tile_id2f(:) @@ -688,11 +709,11 @@ subroutine Initialize(gc, import, export, clock, rc) tcinternal%pgrid_g = pert_grid_g tcinternal%pgrid_f = pert_grid_f tcinternal%pgrid_l = pert_grid_l - if(trim(grid_type) == "Cubed-Sphere" ) then - tcinternal%grid_g = tile_grid_g - else - tcinternal%grid_g = tcinternal%pgrid_g - endif + !if(trim(grid_type) == "Cubed-Sphere" ) then + tcinternal%tgrid_g = tile_grid_g + !else + ! tcinternal%grid_g = tcinternal%pgrid_g + !endif do i = 1, NUM_ENSEMBLE call MAPL_GetObjectFromGC(gcs(METFORCE(i)), CHILD_MAPL, rc=status) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index e88d7f5f..8fc69ef4 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1841,7 +1841,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) trim(out_path), trim(exp_id), exp_domain, & met_force, lai, cat_param, mwRTM_param, & tile_coord_l, tile_coord_rf, & - tcinternal%pgrid_f, tcinternal%pgrid_g, tcinternal%grid_g, & + tcinternal%tgrid_g, tcinternal%pgrid_f, tcinternal%pgrid_g, & N_catl_vec, low_ind, l2rf, rf2l, & N_force_pert, N_progn_pert, force_pert_param, progn_pert_param, & update_type, & @@ -1957,7 +1957,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) if (out_smapL4SMaup) & call write_smapL4SMaup( 'analysis', date_time_new, trim(out_path), & trim(exp_id), NUM_ENSEMBLE, N_catl, N_catf, N_obsl, tile_coord_rf, & - tcinternal%grid_g, N_catl_vec, low_ind, & + tcinternal%tgrid_g, N_catl_vec, low_ind, & N_obs_param, obs_param, Observations_l, cat_param, cat_progn ) end if ! end if (.true.) 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 3e5bd17f..86d6449d 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 @@ -143,7 +143,7 @@ subroutine get_enkf_increments( & work_path, exp_id, exp_domain, & met_force, lai, cat_param, mwRTM_param, & tile_coord_l, tile_coord_f, & - pert_grid_f, pert_grid_g, tile_grid_g, & + tile_grid_g, pert_grid_f, pert_grid_g, & N_catl_vec, low_ind, l2f, f2l, & N_force_pert, N_progn_pert, force_pert_param, progn_pert_param, & update_type, & @@ -192,7 +192,7 @@ subroutine get_enkf_increments( & type(tile_coord_type), dimension(:), pointer :: tile_coord_l, tile_coord_f ! input - type(grid_def_type), intent(in) :: pert_grid_f, pert_grid_g, tile_grid_g + type(grid_def_type), intent(in) :: tile_grid_g, pert_grid_f, pert_grid_g integer, intent(in), dimension(numprocs) :: N_catl_vec, low_ind @@ -1196,7 +1196,7 @@ subroutine get_enkf_increments( & if (out_smapL4SMaup) & call write_smapL4SMaup( 'obs_fcst', date_time, work_path, exp_id, N_ens, & - N_catl, N_catf, N_obsl, tile_coord_f, pert_grid_g, N_catl_vec, low_ind, & + N_catl, N_catf, N_obsl, tile_coord_f, tile_grid_g, N_catl_vec, low_ind, & N_obs_param, obs_param, Observations_l, cat_param, cat_progn ) end if ! end if (.true.) diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index 95cb0b8d..7ddd5b55 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -463,15 +463,15 @@ end subroutine LDAS_create_grid_g function get_minExtent_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc_maxlon, tc_maxlat, & tile_grid_g) result( tile_grid ) - ! get matching tile_grid for given tile_coord and (global) tile_grid_g + ! get tile_grid with smallest extent for given set of tiles in tile_coord on (global) tile_grid_g ! ! make sure to pass in consistent tile_coord (tc) and tile_grid_g inputs: ! ! iff tile_grid_g is the grid that is associated with the tile space definition, - ! then tc_[x]_indg must be tile_coord%[x]_indg + ! then input tc_[x]_indg must be tile_coord%[x]_indg ! ! iff tile_grid_g is the grid that supports efficient mapping for perts and the EnKF analysis, - ! then tc_[x]_indg must be tile_coord%pert_[x]_indg + ! then input tc_[x]_indg must be tile_coord%pert_[x]_indg ! ! reichle, 20 June 2012 -- moved from within domain_setup() ! for re-use in get_obs_pred() diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index e2d402b4..df8bfdc0 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -161,10 +161,10 @@ module LDAS_TileCoordType type(tile_coord_type), pointer, contiguous :: tile_coord(:)=>null() type(tile_coord_type), pointer, contiguous :: tile_coord_f(:)=>null() integer, pointer :: l2f(:)=>null() - type(grid_def_type) :: pgrid_g - type(grid_def_type) :: grid_g - type(grid_def_type) :: pgrid_f - type(grid_def_type) :: pgrid_l + type(grid_def_type) :: tgrid_g ! tile_grid_g + type(grid_def_type) :: pgrid_g ! pert_grid_g + type(grid_def_type) :: pgrid_f ! pert_grid_f + type(grid_def_type) :: pgrid_l ! pert_grid_l end type T_TILECOORD_STATE type TILECOORD_WRAP diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 index bcec2953..8dd53c20 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 @@ -138,11 +138,11 @@ 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) - ! ! For cubed-sphere tile spaces, pert_[x]_indg refers to a lat-lon "hash" grid that will + ! ! For cubed-sphere tile spaces, pert_[x]_indg refers to a lat-lon perturbation grid that will ! ! be created at runtime to support efficient mapping for perturbations and the EnKF analysis. ! ! For EASE and LatLon tile spaces, pert_[x]_indg is identical to [x]_indg - ! integer :: pert_i_indg ! i index (w.r.t. *global* "hash" grid for perts and EnKF) - ! integer :: pert_j_indg ! j index (w.r.t. *global* "hash" grid for perts and EnKF) + ! integer :: pert_i_indg ! i index (w.r.t. *global* pert_grid for perts and EnKF) + ! integer :: pert_j_indg ! j index (w.r.t. *global* pert_grid for perts and EnKF) ! 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]