Skip to content

Commit

Permalink
change hash- to pert- and bring back grid_g
Browse files Browse the repository at this point in the history
  • Loading branch information
weiyuan-jiang committed Mar 17, 2023
1 parent 1021ea0 commit 55ab1f8
Show file tree
Hide file tree
Showing 10 changed files with 55 additions and 46 deletions.
6 changes: 3 additions & 3 deletions src/Applications/LDAS_App/preprocess_ldas_routines.F90
Original file line number Diff line number Diff line change
Expand Up @@ -3027,9 +3027,9 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la
N_tile_land=i
allocate(tile_coord_land(N_tile_land))
tile_coord_land=tile_coord(1:N_tile_land)
! hash_[x]_indg is not written into the tile_coord file and not needed in preprocessing
tile_coord_land%hash_i_indg = nint(nodata_generic)
tile_coord_land%hash_j_indg = nint(nodata_generic)
! pert_[x]_indg is not written into the tile_coord file and not needed in preprocessing
tile_coord_land%pert_i_indg = nint(nodata_generic)
tile_coord_land%pert_j_indg = nint(nodata_generic)
if(present(f2g)) then
allocate(f2g(fid))
f2g = f2g_tmp(1:fid)
Expand Down
18 changes: 13 additions & 5 deletions src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90
Original file line number Diff line number Diff line change
Expand Up @@ -630,10 +630,10 @@ subroutine Initialize(gc, import, export, clock, rc)
!1) get hash index for cubed-sphere
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)%hash_i_indg,tile_coord_f(i)%hash_j_indg)
tile_coord_f(i)%pert_i_indg,tile_coord_f(i)%pert_j_indg)
enddo
!2) re-generate pert_grid_f in Lat-Lon
pert_grid_f = get_minExtent_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, &
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)
endif
Expand All @@ -646,6 +646,10 @@ subroutine Initialize(gc, import, export, clock, rc)
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

block
integer, allocatable :: f2tile_id(:), tile_id2f(:)
integer :: max_id
Expand Down Expand Up @@ -675,16 +679,20 @@ subroutine Initialize(gc, import, export, clock, rc)

allocate(tcinternal%tile_coord_f,source = tile_coord_f)

pert_grid_l = get_minExtent_grid(land_nt_local, &
tcinternal%tile_coord%hash_i_indg, tcinternal%tile_coord%hash_j_indg, &
pert_grid_l = get_minExtent_grid(land_nt_local, &
tcinternal%tile_coord%pert_i_indg, tcinternal%tile_coord%pert_j_indg, &
tcinternal%tile_coord%min_lon, tcinternal%tile_coord%min_lat, &
tcinternal%tile_coord%max_lon, tcinternal%tile_coord%max_lat, &
pert_grid_g)


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

do i = 1, NUM_ENSEMBLE
call MAPL_GetObjectFromGC(gcs(METFORCE(i)), CHILD_MAPL, rc=status)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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%pgrid_f, tcinternal%pgrid_g, tcinternal%grid_g, &
N_catl_vec, low_ind, l2rf, rf2l, &
N_force_pert, N_progn_pert, force_pert_param, progn_pert_param, &
update_type, &
Expand Down Expand Up @@ -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%pgrid_g, N_catl_vec, low_ind, &
tcinternal%grid_g, N_catl_vec, low_ind, &
N_obs_param, obs_param, Observations_l, cat_param, cat_progn )

end if ! end if (.true.)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -943,10 +943,10 @@ subroutine apply_adapt_P( N_pert, pert_param, pert_adapt_param,

! get mean and std in tile space

call grid2tile( pert_grid, N_catd, tile_coord%hash_i_indg,tile_coord%hash_j_indg, & !tile_coord, &
call grid2tile( pert_grid, N_catd, tile_coord%pert_i_indg,tile_coord%pert_j_indg, & !tile_coord, &
pert_param(n)%mean, mu )

call grid2tile( pert_grid, N_catd, tile_coord%hash_i_indg,tile_coord%hash_j_indg, & !tile_coord, &
call grid2tile( pert_grid, N_catd, tile_coord%pert_i_indg,tile_coord%pert_j_indg, & !tile_coord, &
pert_param(n)%std, sg )

select case (pert_param(n)%typ)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
pert_grid_f, pert_grid_g, tile_grid_g, &
N_catl_vec, low_ind, l2f, f2l, &
N_force_pert, N_progn_pert, force_pert_param, progn_pert_param, &
update_type, &
Expand Down Expand Up @@ -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
type(grid_def_type), intent(in) :: pert_grid_f, pert_grid_g, tile_grid_g

integer, intent(in), dimension(numprocs) :: N_catl_vec, low_ind

Expand Down Expand Up @@ -394,12 +394,12 @@ subroutine get_enkf_increments( &

! first call: count how many tiles are in each pert_grid_f cell
call get_number_of_tiles_in_cell_ij( N_catf, &
tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, &
tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, &
pert_grid_f, N_tile_in_cell_ij_f )
! second call: find out which tiles are in each pert_grid_f cell

call get_tile_num_in_cell_ij( N_catf, &
tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, &
tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, &
pert_grid_f, maxval(N_tile_in_cell_ij_f), tile_num_in_cell_ij_f )
else
allocate(N_tile_in_cell_ij_f(0,0)) !for debugging
Expand All @@ -418,7 +418,7 @@ subroutine get_enkf_increments( &
if (out_smapL4SMaup) &
call output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, &
N_ens, N_catl, N_catf, N_obsl, N_obsl_max, &
tile_coord_f, tile_coord_l, pert_grid_g, pert_grid_f, &
tile_coord_f, tile_coord_l, tile_grid_g, pert_grid_f, &
N_catl_vec, low_ind, l2f, N_tile_in_cell_ij_f, tile_num_in_cell_ij_f, &
N_obs_param, obs_param, Observations_l, cat_param, cat_progn )

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1513,7 +1513,7 @@ subroutine get_obs_pred( &

! determine tile_grid_lH from tile_coord_lH

tile_grid_lH = get_minExtent_grid( N_catlH, tile_coord_lH%hash_i_indg, tile_coord_lH%hash_j_indg,&
tile_grid_lH = get_minExtent_grid( N_catlH, tile_coord_lH%pert_i_indg, tile_coord_lH%pert_j_indg,&
tile_coord_lH%min_lon, tile_coord_lH%min_lat, tile_coord_lH%max_lon, tile_coord_lH%max_lat, &
tile_grid_g)

Expand All @@ -1522,15 +1522,15 @@ subroutine get_obs_pred( &
! first call: count how many tiles are in each tile_grid_lH cell

call get_number_of_tiles_in_cell_ij( N_catlH, &
tile_coord_lH%hash_i_indg, tile_coord_lH%hash_j_indg, &
tile_coord_lH%pert_i_indg, tile_coord_lH%pert_j_indg, &
tile_grid_lH, N_tile_in_cell_ij_lH )

! second call: find out which tiles are in each tile_grid_lH cell
! [tile numbers in "tile_num_in_cell_ij_lH" are relative
! to local halo ("lH") domain]

call get_tile_num_in_cell_ij( N_catlH, &
tile_coord_lH%hash_i_indg, tile_coord_lH%hash_j_indg, &
tile_coord_lH%pert_i_indg, tile_coord_lH%pert_j_indg, &
tile_grid_lH, maxval(N_tile_in_cell_ij_lH), tile_num_in_cell_ij_lH )

end if
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1042,8 +1042,8 @@ subroutine Initialize(gc, import, export, clock, rc)
VERIFY_(status)
allocate(internal%j_indgs(land_nt_local),stat=status)
VERIFY_(status)
internal%i_indgs(:)=tile_coord(:)%hash_i_indg
internal%j_indgs(:)=tile_coord(:)%hash_j_indg
internal%i_indgs(:)=tile_coord(:)%pert_i_indg
internal%j_indgs(:)=tile_coord(:)%pert_j_indg

! Get pert options from *default* namelist files
! WARNING: get_force/progn_pert_param() calls allocate memory
Expand Down Expand Up @@ -1661,10 +1661,10 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc)
if (IamRoot) then
! 3) tile2grid. simple reverser of grid2tile without weighted averaging/no-data-handling
do m = 1, nfpert
call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, tcinternal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m))
call tile2grid_simple( N_tile, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, tcinternal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m))
enddo
do m = 1, nppert
call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, tcinternal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m))
call tile2grid_simple( N_tile, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, tcinternal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m))
enddo

! 4) writing
Expand Down Expand Up @@ -2746,10 +2746,10 @@ subroutine Finalize(gc, import, export, clock, rc)
! 3) tile2grid
! this step is simply a reverse of grid2tile without any weighted
do m = 1, nfpert
call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, tcinternal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m))
call tile2grid_simple( N_tile, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, tcinternal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m))
enddo
do m = 1, nppert
call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, tcinternal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m))
call tile2grid_simple( N_tile, tile_coord_f%pert_i_indg, tile_coord_f%pert_j_indg, tcinternal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m))
enddo

! 4) writing
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -471,13 +471,13 @@ function get_minExtent_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat,
! then 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%hash_[x]_indg
! then 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()
!
! reichle+wjiang, 19 Aug 2020 -- changed interface to generically accommodate use of
! tile_coord%[x]_indg or tile_coord%hash_[x]_indg
! tile_coord%[x]_indg or tile_coord%pert_[x]_indg
!
! -------------------------------------------------------------------

Expand Down Expand Up @@ -623,12 +623,12 @@ subroutine get_number_of_tiles_in_cell_ij( N_tile, tc_i_indg, tc_j_indg, tile_gr
! then 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%hash_[x]_indg
! then tc_[x]_indg must be tile_coord%pert_[x]_indg
!
! wjiang(?) -- split off from LDASsa legacy subroutine get_tile_num_in_cell_ij()
!
! reichle+wjiang, 19 Aug 2020 -- changed interface to generically accommodate use of
! tile_coord%[x]_indg or tile_coord%hash_[x]_indg
! tile_coord%[x]_indg or tile_coord%pert_[x]_indg
!
! ----------------------------------------------------------

Expand Down Expand Up @@ -694,12 +694,12 @@ subroutine get_tile_num_in_cell_ij( N_tile, tc_i_indg, tc_j_indg, tile_grid, &
! then 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%hash_[x]_indg
! then tc_[x]_indg must be tile_coord%pert_[x]_indg
!
! wjiang(?) -- split off from LDASsa legacy subroutine get_tile_num_in_cell_ij()
!
! reichle+wjiang, 19 Aug 2020 -- changed interface to generically accommodate use of
! tile_coord%[x]_indg or tile_coord%hash_[x]_indg
! tile_coord%[x]_indg or tile_coord%pert_[x]_indg
!
! ----------------------------------------------------------

Expand Down
19 changes: 10 additions & 9 deletions src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90
Original file line number Diff line number Diff line change
Expand Up @@ -59,11 +59,11 @@ 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)
! For cubed-sphere tile spaces, hash_[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, hash_[x]_indg is identical to [x]_indg
integer :: hash_i_indg ! i index (w.r.t. *global* "hash" grid for perts and EnKF)
integer :: hash_j_indg ! j index (w.r.t. *global* "hash" grid for perts and EnKF)
! For EASE and LatLon tile spaces, pert_[x]_indg is identical to [x]_indg
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]
Expand Down Expand Up @@ -162,6 +162,7 @@ module LDAS_TileCoordType
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
end type T_TILECOORD_STATE
Expand Down Expand Up @@ -367,8 +368,8 @@ subroutine io_tile_coord_type( action, unitnum, N_tile, tile_coord )
! reichle, 24 July 2010
! reichle, 2 May 2013 - changed N_tile to intent(in)
! reichle, 7 Jan 2014 - changed to binary (unformatted) I/O
! wjiang, reichle, 18 Aug 2020 - Added initialization of %hash_[x]_indg during read.
! Note that %hash_[x]_indg is NOT written out.
! wjiang, reichle, 18 Aug 2020 - Added initialization of %pert_[x]_indg during read.
! Note that %pert_[x]_indg is NOT written out.

implicit none

Expand Down Expand Up @@ -465,9 +466,9 @@ subroutine io_tile_coord_type( action, unitnum, N_tile, tile_coord )
read (unitnum, iostat=istat) tmp_real; if (istat>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg)
tile_coord(:)%elev= tmp_real(:)

! Initialize [x]_indg to hash_[x]_indg. For cs tile spaces, hash_[x]_indg will be redefined
tile_coord(:)%hash_i_indg = tile_coord(:)%i_indg
tile_coord(:)%hash_j_indg = tile_coord(:)%j_indg
! Initialize [x]_indg to pert_[x]_indg. For cs tile spaces, pert_[x]_indg will be redefined
tile_coord(:)%pert_i_indg = tile_coord(:)%i_indg
tile_coord(:)%pert_j_indg = tile_coord(:)%j_indg
tile_coord(:)%f_num = -9999 ! not assigned values yet
deallocate(tmp_int, tmp_real)
case default
Expand Down
10 changes: 5 additions & 5 deletions src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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, hash_[x]_indg refers to a lat-lon "hash" grid that will
! ! For cubed-sphere tile spaces, pert_[x]_indg refers to a lat-lon "hash" grid that will
! ! be created at runtime to support efficient mapping for perturbations and the EnKF analysis.
! ! For EASE and LatLon tile spaces, hash_[x]_indg is identical to [x]_indg
! integer :: hash_i_indg ! i index (w.r.t. *global* "hash" grid for perts and EnKF)
! integer :: hash_j_indg ! j index (w.r.t. *global* "hash" grid for perts and EnKF)
! ! 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)
! 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]
Expand All @@ -161,7 +161,7 @@ subroutine init_MPI_types()

iblock(1) = 4
iblock(2) = 6
iblock(3) = 4 ! i_indg, j_indg, hash_i_indg and hash_j_indg
iblock(3) = 4 ! i_indg, j_indg, pert_i_indg and pert_j_indg
iblock(4) = 4

idisp(1) = 0
Expand Down

0 comments on commit 55ab1f8

Please sign in to comment.