From d3b8e14f611fb17deb2e1d20304933e004320489 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Thu, 9 Mar 2023 12:25:27 -0500 Subject: [PATCH 01/12] clean up grid names --- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 62 ++++---- .../GEOS_LandAssimGridComp.F90 | 10 +- .../clsm_ensupd_enkf_update.F90 | 88 ++++++------ .../GEOS_LandPertGridComp.F90 | 136 +++++++++--------- .../LDAS_PertRoutines.F90 | 6 +- .../Shared/LDAS_PertTypes.F90 | 3 +- .../Shared/LDAS_TileCoordType.F90 | 6 +- 7 files changed, 155 insertions(+), 156 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 2401e935..54f45c8c 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -406,9 +406,10 @@ subroutine Initialize(gc, import, export, clock, rc) integer :: N_catf integer :: LSM_CHOICE - type(grid_def_type) :: tile_grid_g - type(grid_def_type) :: tile_grid_f - type(grid_def_type) :: tile_grid_l + type(grid_def_type) :: tile_grid_g, pert_grid_g + type(grid_def_type) :: tile_grid_f, pert_grid_f + type(grid_def_type) :: tile_grid_l, pert_grid_l + type(date_time_type):: start_time type(ESMF_Time) :: CurrentTime !type(CubedSphereGridFactory) :: cubed_sphere_factory @@ -619,39 +620,33 @@ subroutine Initialize(gc, import, export, clock, rc) close(10) call io_grid_def_type('w', logunit, tile_grid_f, 'tile_grid_f') - block - type(grid_def_type) :: latlon_tmp_g - integer :: perturbations - - call MAPL_GetResource(MAPL, perturbations, 'PERTURBATIONS:', default=0, rc=status) - 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") - - !1) 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 - !2) get hash index - do i = 1, N_catf - call get_ij_ind_from_latlon(latlon_tmp_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) - enddo - !3) re-generate tile_grid_f in Lat-Lon - call get_tile_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, & - tile_coord_f%min_lon, tile_coord_f%min_lat, tile_coord_f%max_lon, tile_coord_f%max_lat, & - tile_grid_g, tile_grid_f) + call get_pert_grid(tile_grid_g, pert_grid_g) + 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") + + !1) generate a lat-lon grid for landpert and land assim ( 4*N_lonX3*N_lon) + !2) get hash index + 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) + enddo + !3) re-generate tile_grid_f in Lat-Lon + call get_tile_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_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, pert_grid_f) - endif - end block - + else + call get_pert_grid(tile_grid_f, pert_grid_f) + endif endif call MPI_BCAST(N_catf,1,MPI_INTEGER,0,mpicomm,mpierr) if (.not. IamRoot) allocate(tile_coord_f(N_catf)) call MPI_BCAST(tile_coord_f,N_catf, MPI_tile_coord_type,0,mpicomm, mpierr) - call MPI_BCAST(tile_grid_g, 1, MPI_grid_def_type, 0,mpicomm, mpierr) - call MPI_BCAST(tile_grid_f, 1, MPI_grid_def_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) block integer, allocatable :: f2tile_id(:), tile_id2f(:) @@ -686,13 +681,12 @@ subroutine Initialize(gc, import, export, clock, rc) tcinternal%tile_coord%hash_i_indg, tcinternal%tile_coord%hash_j_indg, & tcinternal%tile_coord%min_lon, tcinternal%tile_coord%min_lat, & tcinternal%tile_coord%max_lon, tcinternal%tile_coord%max_lat, & - tile_grid_g,tile_grid_l) + pert_grid_g, pert_grid_l) - ! re-arrange tile_coord_f - tcinternal%grid_g = tile_grid_g - tcinternal%grid_f = tile_grid_f - tcinternal%grid_l = tile_grid_l + tcinternal%pgrid_g = pert_grid_g + tcinternal%pgrid_f = pert_grid_f + tcinternal%pgrid_l = pert_grid_l 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 d72d3521..35458506 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1078,8 +1078,6 @@ subroutine Initialize(gc, import, export, clock, rc) ! mapping f to re-orderd f so it is continous for mpi_gather ! rf -- ordered by processors. Within the processor, ordered by MAPL grid integer, allocatable :: f2rf(:) ! mapping re-orderd rf to f for the LDASsa output - type(grid_def_type) :: tile_grid_g - type(grid_def_type) :: tile_grid_f character(len=300) :: seed_fname character(len=300) :: fname_tpl character(len=14) :: datestamp @@ -1842,8 +1840,8 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) NUM_ENSEMBLE, N_catl, N_catf, N_obsl_max, & trim(out_path), trim(exp_id), exp_domain, & met_force, lai, cat_param, mwRTM_param, & - tile_coord_l, tile_coord_rf, tcinternal%grid_f, & - tcinternal%grid_f, tcinternal%grid_l, tcinternal%grid_g, & + tile_coord_l, tile_coord_rf, tcinternal%pgrid_f, & + tcinternal%pgrid_f, tcinternal%pgrid_l, tcinternal%pgrid_g, & N_catl_vec, low_ind, l2rf, rf2l, & N_force_pert, N_progn_pert, force_pert_param, progn_pert_param, & update_type, & @@ -1880,7 +1878,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) date_time_new, trim(out_path), trim(exp_id), & N_obsl, N_obs_param, NUM_ENSEMBLE, & N_catl, tile_coord_l, & - N_catf, tile_coord_rf, tcinternal%grid_f, tcinternal%grid_g, & + N_catf, tile_coord_rf, tcinternal%pgrid_f, tcinternal%pgrid_g, & N_catl_vec, low_ind, rf2l, N_catg, rf2g, & obs_param, & met_force, lai, & @@ -1959,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%pgrid_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 764cfada..eb815f7a 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 @@ -142,8 +142,8 @@ subroutine get_enkf_increments( & N_ens, N_catl, N_catf, N_obsl_max, & work_path, exp_id, exp_domain, & met_force, lai, cat_param, mwRTM_param, & - tile_coord_l, tile_coord_f, tile_grid_f, & - pert_grid_f, pert_grid_l_NotUsed, tile_grid_g, & + tile_coord_l, tile_coord_f, pert_grid_f_NotUsed, & + pert_grid_f, pert_grid_l_NotUsed, 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) :: tile_grid_f, pert_grid_f, pert_grid_l_NotUsed, tile_grid_g + type(grid_def_type), intent(in) :: pert_grid_f_NotUsed, pert_grid_f, pert_grid_l_NotUsed, pert_grid_g integer, intent(in), dimension(numprocs) :: N_catl_vec, low_ind @@ -390,17 +390,17 @@ subroutine get_enkf_increments( & if ( (root_proc) .or. & (any(obs_param(1:N_obs_param)%FOV>FOV_threshold)) ) then - allocate(N_tile_in_cell_ij_f(tile_grid_f%N_lon,tile_grid_f%N_lat)) + allocate(N_tile_in_cell_ij_f(pert_grid_f%N_lon,pert_grid_f%N_lat)) - ! first call: count how many tiles are in each tile_grid_f cell + ! 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_grid_f, N_tile_in_cell_ij_f ) - ! second call: find out which tiles are in each tile_grid_f cell + 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_grid_f, maxval(N_tile_in_cell_ij_f), tile_num_in_cell_ij_f ) + 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 end if @@ -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, tile_grid_g, tile_grid_f, & + tile_coord_f, tile_coord_l, pert_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 ) @@ -427,7 +427,7 @@ subroutine get_enkf_increments( & call collect_obs( & work_path, exp_id, date_time, dtstep_assim, & N_catl, tile_coord_l, & - N_catf, tile_coord_f, tile_grid_f, & + N_catf, tile_coord_f, pert_grid_f, & N_tile_in_cell_ij_f, tile_num_in_cell_ij_f, & N_catl_vec, low_ind, l2f, & N_obs_param, obs_param, N_obsl_max, out_obslog, & @@ -510,7 +510,7 @@ subroutine get_enkf_increments( & N_obs_param, N_ens, & N_catl, tile_coord_l, & N_catf, tile_coord_f, f2l, & - N_catl_vec, low_ind, tile_grid_g, & + N_catl_vec, low_ind, pert_grid_g, & obs_param, & met_force, lai, cat_param, cat_progn, mwRTM_param, & N_obsl, Observations_l, Obs_pred_l, obsbias_ok, & @@ -1046,7 +1046,7 @@ subroutine get_enkf_increments( & call cat_enkf_increments( & N_ens, nObs_ana, nTiles_ana, N_obs_param, & update_type, obs_param, & - tile_grid_f, tile_coord_ana, indTiles_ana, & ! indTiles_ana is essentially ana2f + pert_grid_f, tile_coord_ana, indTiles_ana, & ! indTiles_ana is essentially ana2f Obs_ana, & ! size: nObs_ana Obs_pred_ana, & ! size: (nObs_ana,N_ens) Obs_pert_tmp, & @@ -1074,7 +1074,7 @@ subroutine get_enkf_increments( & call cat_enkf_increments( & N_ens, N_obslH, N_catl, N_obs_param, & update_type, obs_param, & - tile_grid_f, tile_coord_l, l2f, & + pert_grid_f, tile_coord_l, l2f, & Observations_lH(1:N_obslH), & Obs_pred_lH(1:N_obslH,1:N_ens), & Obs_pert_tmp, & @@ -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, tile_grid_g, N_catl_vec, low_ind, & + N_catl, N_catf, N_obsl, tile_coord_f, pert_grid_g, N_catl_vec, low_ind, & N_obs_param, obs_param, Observations_l, cat_param, cat_progn ) end if ! end if (.true.) @@ -1626,7 +1626,7 @@ subroutine output_incr_etc( out_ObsFcstAna, & date_time, work_path, exp_id, & N_obsl, N_obs_param, N_ens, & N_catl, tile_coord_l, & - N_catf, tile_coord_f, tile_grid_f, tile_grid_g, & + N_catf, tile_coord_f, pert_grid_f, pert_grid_g, & N_catl_vec, low_ind, f2l, N_catg, f2g, & obs_param, & met_force, lai, cat_param, cat_progn, cat_progn_incr, mwRTM_param, & @@ -1654,8 +1654,8 @@ subroutine output_incr_etc( out_ObsFcstAna, & type(tile_coord_type), dimension(:), pointer :: tile_coord_l ! input type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input - type(grid_def_type), intent(in) :: tile_grid_f - type(grid_def_type), intent(in) :: tile_grid_g + type(grid_def_type), intent(in) :: pert_grid_f + type(grid_def_type), intent(in) :: pert_grid_g integer, dimension(numprocs), intent(in) :: N_catl_vec, low_ind @@ -1716,7 +1716,7 @@ subroutine output_incr_etc( out_ObsFcstAna, & N_obs_param, N_ens, & N_catl, tile_coord_l, & N_catf, tile_coord_f, f2l, & - N_catl_vec, low_ind, tile_grid_g, & + N_catl_vec, low_ind, pert_grid_g, & obs_param, & met_force, lai, cat_param, cat_progn, mwRTM_param, & N_obsl_tmp, Observations_l, Obs_pred_l ) @@ -1831,7 +1831,7 @@ end subroutine output_incr_etc subroutine 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, tile_grid_g, tile_grid_f, & + tile_coord_f, tile_coord_l, pert_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 ) @@ -1860,8 +1860,8 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input type(tile_coord_type), dimension(:), pointer :: tile_coord_l ! input - type(grid_def_type), intent(in) :: tile_grid_g - type(grid_def_type), intent(in) :: tile_grid_f + type(grid_def_type), intent(in) :: pert_grid_g + type(grid_def_type), intent(in) :: pert_grid_f integer, dimension(numprocs), intent(in) :: N_catl_vec integer, dimension(numprocs), intent(in) :: low_ind @@ -1936,7 +1936,7 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & call collect_obs( & work_path, exp_id, date_time, dtstep_assim, & N_catl, tile_coord_l, & - N_catf, tile_coord_f, tile_grid_f, & + N_catf, tile_coord_f, pert_grid_f, & N_tile_in_cell_ij_f, tile_num_in_cell_ij_f, & N_catl_vec, low_ind, l2f, & N_obs_param_tmp, obs_param_tmp(1:N_obs_param_tmp), N_obsl_max, .false., & @@ -1945,7 +1945,7 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & ! write appropriate fields (according to 'option') into file call write_smapL4SMaup( 'orig_obs', date_time, work_path, exp_id, N_ens, & - N_catl, N_catf, N_obsl_tmp, tile_coord_f, tile_grid_g, & + N_catl, N_catf, N_obsl_tmp, tile_coord_f, pert_grid_g, & N_catl_vec, low_ind, & N_obs_param_tmp, obs_param_tmp(1:N_obs_param_tmp), Observations_l, & cat_param, cat_progn ) @@ -1955,7 +1955,7 @@ end subroutine output_smapL4SMaup ! ********************************************************************** subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & - N_catl, N_catf, N_obsl, tile_coord_f, tile_grid_g, N_catl_vec, low_ind, & + N_catl, N_catf, N_obsl, tile_coord_f, pert_grid_g, N_catl_vec, low_ind, & N_obs_param, obs_param, Observations_l, cat_param, cat_progn ) ! output of custom collection for SMAP L4_SM "aup" (analysis update) @@ -2034,7 +2034,7 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input - type(grid_def_type), intent(in) :: tile_grid_g + type(grid_def_type), intent(in) :: pert_grid_g integer, dimension(numprocs), intent(in) :: N_catl_vec integer, dimension(numprocs), intent(in) :: low_ind @@ -2105,7 +2105,7 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & ! ! smapL4SMaup output only works for 9 km EASE grids - if ( index(tile_grid_g%gridtype, 'M09') == 0 ) then + if ( index(pert_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,8 +2228,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMAP_L1C_Tbh_E09_A', 'SMAP_L1C_Tbv_E09_A' & ) - 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) + if (index(pert_grid_g%gridtype, 'M09') /=0) then + call ease_convert(trim(pert_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 @@ -2247,8 +2247,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMAP_L1C_Tbh_E27_A', 'SMAP_L1C_Tbv_E27_A' & ) - 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) + if (index(pert_grid_g%gridtype, 'M09') /=0) then + call ease_convert(trim(pert_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 @@ -2269,9 +2269,9 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMOS_fit_Tbh_A', 'SMOS_fit_Tbv_A' & ) - if (index(tile_grid_g%gridtype, 'M09') /=0) then + if (index(pert_grid_g%gridtype, 'M09') /=0) then ! subindex (1:7) to get the string EASEvx_ - gridname_tmp = tile_grid_g%gridtype(1:7)//'M36' + gridname_tmp = pert_grid_g%gridtype(1:7)//'M36' call ease_convert(gridname_tmp, this_lat, this_lon, col_ind, row_ind) endif @@ -2298,8 +2298,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & ! map Observations%[xx] fields onto global 9km EASE grid, - allocate(ndata_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) - allocate(ndata_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(ndata_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(ndata_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) if (option=='orig_obs') then @@ -2312,8 +2312,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata_8(N_obsf)) - allocate(data_h_9km_grid_8(tile_grid_g%N_lon,tile_grid_g%N_lat)) - allocate(data_v_9km_grid_8(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_h_9km_grid_8(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_v_9km_grid_8(pert_grid_g%N_lon,pert_grid_g%N_lat)) allocate(data_h_9km_tile_8(N_catf)) allocate(data_v_9km_tile_8(N_catf)) @@ -2326,8 +2326,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata(N_obsf)) - allocate(data_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) - allocate(data_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) allocate(data_h_9km_tile(N_catf)) allocate(data_v_9km_tile(N_catf)) @@ -2357,8 +2357,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata(N_obsf)) - allocate(data_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) - allocate(data_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) allocate(data_h_9km_tile(N_catf)) allocate(data_v_9km_tile(N_catf)) @@ -2740,10 +2740,10 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & data_h_9km_tile_8 = real(nodata_generic,kind(0.0D0)) ! init (not in grid2tile!) data_v_9km_tile_8 = real(nodata_generic,kind(0.0D0)) ! init (not in grid2tile!) - call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_h_9km_grid_8, & + call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_h_9km_grid_8, & data_h_9km_tile_8 ) - call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_v_9km_grid_8, & + call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_v_9km_grid_8, & data_v_9km_tile_8 ) ! write into file @@ -2756,10 +2756,10 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & data_h_9km_tile = nodata_generic ! initialize (not done in grid2tile!) data_v_9km_tile = nodata_generic ! initialize (not done in grid2tile!) - call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_h_9km_grid, & + call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_h_9km_grid, & data_h_9km_tile ) - call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_v_9km_grid, & + call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_v_9km_grid, & data_v_9km_tile ) ! write into file diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index dd159ffc..514da2e9 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -850,6 +850,7 @@ subroutine Initialize(gc, import, export, clock, rc) type(LANDPERT_WRAP) :: wrap type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal type(tile_coord_type), pointer :: tile_coord(:)=>null() ! MAPL internal pointers @@ -895,7 +896,8 @@ subroutine Initialize(gc, import, export, clock, rc) ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) - tile_coord => tcwrap%ptr%tile_coord + tcinternal => tcwrap%ptr + tile_coord => tcinternal%tile_coord ! Are we perturbing variables? call MAPL_GetResource(MAPL, internal%PERTURBATIONS, 'PERTURBATIONS:', default=0, rc=status) @@ -931,20 +933,16 @@ subroutine Initialize(gc, import, export, clock, rc) GEOSldas_FORCE_PERT_DTSTEP = internal%ForcePert%dtstep GEOSldas_PROGN_PERT_DTSTEP = internal%PrognPert%dtstep - call get_pert_grid(tcwrap%ptr%grid_g,internal%pgrid_g) - call get_pert_grid(tcwrap%ptr%grid_f,internal%pgrid_f) - call get_pert_grid(tcwrap%ptr%grid_l,internal%pgrid_l) - - n_lon = internal%pgrid_l%n_lon - n_lat = internal%pgrid_l%n_lat + n_lon = tcinternal%pgrid_l%n_lon + n_lat = tcinternal%pgrid_l%n_lat call MAPL_GetResource( MAPL, ens_id_width,"ENS_ID_WIDTH:", default=4, RC=STATUS) VERIFY_(status) ! Pointers to mapl internals if ( internal%isCubedSphere ) then - n_lon_g = internal%pgrid_g%n_lon - n_lat_g = internal%pgrid_g%n_lat + n_lon_g = tcinternal%pgrid_g%n_lon + n_lat_g = tcinternal%pgrid_g%n_lat allocate(internal%fpert_ntrmdt(n_lon_g, n_lat_g, N_FORCE_PERT_MAX), source=0.0) allocate(internal%ppert_ntrmdt(n_lon_g, n_lat_g, N_PROGN_PERT_MAX), source=0.0) allocate(internal%pert_rseed_r8(NRANDSEED), source=0.0d0) @@ -994,10 +992,10 @@ subroutine Initialize(gc, import, export, clock, rc) VERIFY_(status) endif endif - lon1 = internal%pgrid_l%i_offg + 1 - lon2 = internal%pgrid_l%i_offg + n_lon - lat1 = internal%pgrid_l%j_offg + 1 - lat2 = internal%pgrid_l%j_offg + n_lat + lon1 = tcinternal%pgrid_l%i_offg + 1 + lon2 = tcinternal%pgrid_l%i_offg + n_lon + lat1 = tcinternal%pgrid_l%j_offg + 1 + lat2 = tcinternal%pgrid_l%j_offg + n_lat else call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=MINTERNAL, rc=status) VERIFY_(status) @@ -1012,10 +1010,10 @@ subroutine Initialize(gc, import, export, clock, rc) VERIFY_(status) call ESMF_GRID_INTERIOR (GRID, I1,IN,J1,JN) - lon1 = internal%pgrid_l%i_offg + 1 ! global index, starting from 1 + lon1 = tcinternal%pgrid_l%i_offg + 1 ! global index, starting from 1 lon1 = lon1 - i1 + 1 ! relative to local lon2 = lon1 + n_lon - 1 - lat1 = internal%pgrid_l%j_offg + 1 ! global index, starting from 1 + lat1 = tcinternal%pgrid_l%j_offg + 1 ! global index, starting from 1 lat1 = lat1 - j1 +1 ! relative to local lat2 = lat1 + n_lat - 1 endif @@ -1051,7 +1049,7 @@ subroutine Initialize(gc, import, export, clock, rc) ! Get pert options from *default* namelist files ! WARNING: get_force/progn_pert_param() calls allocate memory - call get_force_pert_param(internal%pgrid_l, internal%ForcePert%npert, internal%ForcePert%param) + call get_force_pert_param(tcinternal%pgrid_l, internal%ForcePert%npert, internal%ForcePert%param) _ASSERT(internal%ForcePert%npert==size(internal%ForcePert%param), "ForcePert: param size does not match npert") internal%ForcePert%fft_npert = internal%ForcePert%npert @@ -1066,7 +1064,7 @@ subroutine Initialize(gc, import, export, clock, rc) call MAPL_CommsBcast(vm, data=internal%ForcePert%param(n)%coarsen, N=1, ROOT=0,rc=status) enddo - call get_progn_pert_param(internal%pgrid_l, internal%PrognPert%npert, internal%PrognPert%param) + call get_progn_pert_param(tcinternal%pgrid_l, internal%PrognPert%npert, internal%PrognPert%param) _ASSERT(internal%PrognPert%npert==size(internal%PrognPert%param), "PrognPert: param size does not match npert") internal%PrognPert%fft_npert = internal%PrognPert%npert @@ -1141,7 +1139,7 @@ subroutine Initialize(gc, import, export, clock, rc) call propagate_pert( & internal%ForcePert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & ! arbitrary dtstep -1.0, & pert_rseed, & @@ -1156,7 +1154,7 @@ subroutine Initialize(gc, import, export, clock, rc) call propagate_pert( & internal%PrognPert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & ! arbitrary dtstep -1.0, & pert_rseed, & @@ -1172,7 +1170,7 @@ subroutine Initialize(gc, import, export, clock, rc) if(internal%ens_id == FIRST_ENS_ID ) fpert_enavg(:,:,:)=0. do m = 1,internal%ForcePert%npert - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,m)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,m)) if(internal%ForcePert%param(m)%zeromean .and. internal%NUM_ENSEMBLE >2) then fpert_enavg(:,:,m)=fpert_enavg(:,:,m)+fpert_ntrmdt(lon1:lon2,lat1:lat2,m) if( internal%ens_id-FIRST_ENS_ID == internal%NUM_ENSEMBLE-1) then @@ -1184,7 +1182,7 @@ subroutine Initialize(gc, import, export, clock, rc) if(internal%ens_id == FIRST_ENS_ID) ppert_enavg(:,:,:)=0. do m = 1,internal%PrognPert%npert - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,m)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,m)) if(internal%PrognPert%param(m)%zeromean .and. internal%NUM_ENSEMBLE >2) then ppert_enavg(:,:,m)=ppert_enavg(:,:,m)+ppert_ntrmdt(lon1:lon2,lat1:lat2,m) if( internal%ens_id - FIRST_ENS_ID == internal%NUM_ENSEMBLE-1) then @@ -1311,6 +1309,7 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) type(LANDPERT_WRAP) :: wrap type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal type(tile_coord_type), pointer :: tile_coord(:)=>null() ! MAPL internal pointers @@ -1349,13 +1348,14 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) call ESMF_UserCompGetInternalState(gc, 'Landpert_state', wrap, status) VERIFY_(status) internal => wrap%ptr - n_lon = internal%pgrid_l%n_lon - n_lat = internal%pgrid_l%n_lat ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) - tile_coord => tcwrap%ptr%tile_coord + tcinternal => tcwrap%ptr + tile_coord => tcinternal%tile_coord + n_lon = tcinternal%pgrid_l%n_lon + n_lat = tcinternal%pgrid_l%n_lat if (internal%PERTURBATIONS == 0) then ! no perturbations call MAPL_TimerOff(MAPL, "phase2_Initialize") @@ -1419,7 +1419,7 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) internal%ForcePert%npert, & internal%ForcePert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%ForcePert%dtstep), & internal%ForcePert%param, & pert_rseed, & @@ -1432,9 +1432,9 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) ) do ipert=1,internal%ForcePert%npert - call grid2tile( internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & + call grid2tile( tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & fpert_grid(:,:,ipert), internal%ForcePert%DataPrv(:,ipert)) - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) end do internal%ForcePert%DataNxt = internal%ForcePert%DataPrv @@ -1451,7 +1451,7 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) internal%PrognPert%npert, & internal%PrognPert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%PrognPert%dtstep), & internal%PrognPert%param, & pert_rseed, & @@ -1464,9 +1464,9 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) ) do ipert=1,internal%PrognPert%npert - call grid2tile( internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_grid(:,:,ipert), & + call grid2tile( tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_grid(:,:,ipert), & internal%PrognPert%DataPrv(:,ipert)) - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) end do internal%PrognPert%DataNxt = internal%PrognPert%DataPrv @@ -1528,6 +1528,7 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) type(T_LANDPERT_STATE), pointer :: internal=>null() type(LANDPERT_WRAP) :: wrap type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal type(MAPL_LocStream) :: locstream ! MAPL internal pointers @@ -1562,14 +1563,18 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) VERIFY_(status) internal => wrap%ptr - n_lon=internal%pgrid_l%n_lon - n_lat=internal%pgrid_l%n_lat - if (internal%PERTURBATIONS == 0) then ! no perturbations call MAPL_TimerOff(MAPL, "GenerateRaw") RETURN_(ESMF_SUCCESS) end if + call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) + VERIFY_(status) + tcinternal => tcwrap%ptr + + n_lon=tcinternal%pgrid_l%n_lon + n_lat=tcinternal%pgrid_l%n_lat + ! Get locstream call MAPL_Get(MAPL, LocStream=locstream,rc=status) VERIFY_(status) @@ -1617,22 +1622,20 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) VERIFY_(STATUS) call MAPL_TileMaskGet(tilegrid, mask, rc=status) VERIFY_(STATUS) - call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) - VERIFY_(status) nfpert = internal%ForcePert%npert nppert = internal%PrognPert%npert - tile_coord_f => tcwrap%ptr%tile_coord_f + tile_coord_f => tcinternal%tile_coord_f n_tile = size(tile_coord_f,1) ! 1) grid2tile allocate(tile_data_f(land_nt_local,nfpert)) allocate(tile_data_p(land_nt_local,nppert)) do m = 1, nfpert - call grid2tile(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & + call grid2tile(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & fpert_ntrmdt(lon1:lon2,lat1:lat2,m), tile_data_f(:,m)) enddo do m = 1, nppert - call grid2tile(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & + call grid2tile(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & ppert_ntrmdt(lon1:lon2,lat1:lat2,m), tile_data_p(:,m)) enddo ! 2) gather tiledata @@ -1659,10 +1662,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, internal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m)) + 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)) enddo do m = 1, nppert - call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, internal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m)) + 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)) enddo ! 4) writing @@ -1684,7 +1687,7 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) call propagate_pert( & internal%ForcePert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%ForcePert%dtstep), & pert_rseed, & internal%ForcePert%param, & @@ -1695,7 +1698,7 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) if(internal%ens_id == FIRST_ENS_ID ) fpert_enavg(:,:,:)=0. do m = 1,internal%ForcePert%npert - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,m)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,m)) if(internal%ForcePert%param(m)%zeromean .and. internal%NUM_ENSEMBLE >2) then fpert_enavg(:,:,m)=fpert_enavg(:,:,m)+fpert_ntrmdt(lon1:lon2,lat1:lat2,m) if( internal%ens_id - FIRST_ENS_ID == internal%NUM_ENSEMBLE-1) then @@ -1712,7 +1715,7 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) call propagate_pert( & internal%PrognPert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%PrognPert%dtstep), & pert_rseed, & internal%PrognPert%param, & @@ -1723,7 +1726,7 @@ subroutine GenerateRaw_ntrmdt(gc, import, export, clock, rc) if(internal%ens_id == FIRST_ENS_ID) ppert_enavg(:,:,:)=0. do m = 1,internal%PrognPert%npert - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,m)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,m)) if(internal%PrognPert%param(m)%zeromean .and. internal%NUM_ENSEMBLE >2) then ppert_enavg(:,:,m)=ppert_enavg(:,:,m)+ppert_ntrmdt(lon1:lon2,lat1:lat2,m) if( internal%ens_id - FIRST_ENS_ID == internal%NUM_ENSEMBLE -1) then @@ -1831,6 +1834,7 @@ subroutine ApplyForcePert(gc, import, export, clock, rc) real, allocatable :: FORCEPERT(:,:) real, allocatable :: fpert_grid(:,:,:) type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal type(tile_coord_type), pointer :: tile_coord(:)=>null() integer :: n_lon,n_lat @@ -1866,13 +1870,14 @@ subroutine ApplyForcePert(gc, import, export, clock, rc) call ESMF_UserCompGetInternalState(gc, 'Landpert_state', wrap, status) VERIFY_(status) internal => wrap%ptr - n_lon = internal%pgrid_l%n_lon - n_lat = internal%pgrid_l%n_lat ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) - tile_coord => tcwrap%ptr%tile_coord + tcinternal => tcwrap%ptr + n_lon = tcinternal%pgrid_l%n_lon + n_lat = tcinternal%pgrid_l%n_lat + tile_coord => tcinternal%tile_coord ! Get locstream call MAPL_Get(MAPL, LocStream=locstream,rc=status) @@ -1951,7 +1956,7 @@ subroutine ApplyForcePert(gc, import, export, clock, rc) internal%ForcePert%npert, & internal%ForcePert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%ForcePert%dtstep), & internal%ForcePert%param, & pert_rseed, & @@ -1968,9 +1973,9 @@ subroutine ApplyForcePert(gc, import, export, clock, rc) ! -convert-nxt-gridded-perturbations-to-tile- call MAPL_TimerOn(MAPL, '-LocStreamTransform') do ipert=1,internal%ForcePert%npert - call grid2tile( internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_grid(:,:,ipert), & + call grid2tile( tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_grid(:,:,ipert), & internal%ForcePert%DataNxt(:,ipert)) - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), fpert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) end do call MAPL_TimerOff(MAPL, '-LocStreamTransform') @@ -2267,6 +2272,7 @@ subroutine ApplyPrognPert(gc, import, export, clock, rc) real, allocatable :: PROGNPERT(:,:) real, allocatable :: ppert_grid(:,:,:) type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal type(tile_coord_type), pointer :: tile_coord(:)=>null() integer :: n_lon,n_lat @@ -2307,13 +2313,14 @@ subroutine ApplyPrognPert(gc, import, export, clock, rc) call MAPL_TimerOn(MAPL, "TOTAL") call MAPL_TimerOn(MAPL, "Run_ApplyPrognPert") - n_lon = internal%pgrid_l%n_lon - n_lat = internal%pgrid_l%n_lat ! Get component's internal tile_coord variable call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) VERIFY_(status) - tile_coord => tcwrap%ptr%tile_coord + tcinternal => tcwrap%ptr + n_lon = tcinternal%pgrid_l%n_lon + n_lat = tcinternal%pgrid_l%n_lat + tile_coord => tcinternal%tile_coord ! Pointers to imports call MAPL_GetPointer(import, tcPert, 'TCPert', rc=status) @@ -2430,7 +2437,7 @@ subroutine ApplyPrognPert(gc, import, export, clock, rc) internal%PrognPert%npert, & internal%PrognPert%fft_npert, & 1, & - internal%pgrid_l, internal%pgrid_f, & + tcinternal%pgrid_l, tcinternal%pgrid_f, & real(internal%PrognPert%dtstep), & internal%PrognPert%param, & pert_rseed, & @@ -2446,9 +2453,9 @@ subroutine ApplyPrognPert(gc, import, export, clock, rc) ! -convert-nxt-gridded-perturbations-to-tile- call MAPL_TimerOn(MAPL, '-LocStreamTransform') do ipert=1,internal%PrognPert%npert - call grid2tile( internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_grid(:,:,ipert), & + call grid2tile( tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_grid(:,:,ipert), & internal%PrognPert%DataNxt(:,ipert)) - call tile_mask_grid(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) + call tile_mask_grid(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), ppert_ntrmdt(lon1:lon2,lat1:lat2,ipert)) end do call MAPL_TimerOff(MAPL, '-LocStreamTransform') @@ -2660,6 +2667,7 @@ subroutine Finalize(gc, import, export, clock, rc) type(LANDPERT_WRAP) :: wrap type(MAPL_LocStream) :: locstream type(TILECOORD_WRAP) :: tcwrap + type(T_TILECOORD_STATE), pointer :: tcinternal integer :: m,n_lon,n_lat, land_nt_local, ens_id_width integer :: nfpert, nppert, n_tile @@ -2685,6 +2693,9 @@ subroutine Finalize(gc, import, export, clock, rc) call ESMF_UserCompGetInternalState(gc, 'Landpert_state', wrap, status) VERIFY_(status) internal => wrap%ptr + call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) + VERIFY_(status) + tcinternal => tcwrap%ptr if ( internal%isCubedSphere .and. internal%PERTURBATIONS /= 0) then call MAPL_Get(MAPL, LocStream=locstream,rc=status) @@ -2695,22 +2706,19 @@ subroutine Finalize(gc, import, export, clock, rc) VERIFY_(STATUS) call MAPL_TileMaskGet(tilegrid, mask, rc=status) VERIFY_(STATUS) - call ESMF_UserCompGetInternalState(gc, 'TILE_COORD', tcwrap, status) - VERIFY_(status) - nfpert = internal%ForcePert%npert nppert = internal%PrognPert%npert - tile_coord_f => tcwrap%ptr%tile_coord_f + tile_coord_f => tcinternal%tile_coord_f n_tile = size(tile_coord_f,1) ! 1) grid2tile allocate(tile_data_f(land_nt_local,nfpert)) allocate(tile_data_p(land_nt_local,nppert)) do m = 1, nfpert - call grid2tile(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & + call grid2tile(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & internal%fpert_ntrmdt(lon1:lon2,lat1:lat2,m), tile_data_f(:,m)) enddo do m = 1, nppert - call grid2tile(internal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & + call grid2tile(tcinternal%pgrid_l, land_nt_local, internal%i_indgs(:),internal%j_indgs(:), & internal%ppert_ntrmdt(lon1:lon2,lat1:lat2,m), tile_data_p(:,m)) enddo ! 2) gather tiledata @@ -2739,10 +2747,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, internal%pgrid_g, tile_data_f_all(:,m), internal%fpert_ntrmdt(:,:,m)) + 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)) enddo do m = 1, nppert - call tile2grid_simple( N_tile, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, internal%pgrid_g, tile_data_p_all(:,m), internal%ppert_ntrmdt(:,:,m)) + 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)) enddo ! 4) writing diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 index d2c6e846..0e84fc44 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 @@ -641,7 +641,7 @@ subroutine get_pert_grid( tile_grid, pert_grid ) character(len=30) :: latlon_gridname character(len=6) :: lattmp,lontmp - type(grid_def_type) :: tile_grid_g,tile_grid_tmp + type(grid_def_type) :: latlon_grid_tmp integer :: n_x,i_off,j_off,n_lon,n_lat @@ -677,9 +677,9 @@ subroutine get_pert_grid( tile_grid, pert_grid ) write(lontmp,'(I6.6)') n_lon latlon_gridname = "DE"//lontmp//"x"//"PE"//lattmp - call LDAS_create_grid_g(latlon_gridname,n_lon,n_lat,tile_grid_g,i_off,j_off) + call LDAS_create_grid_g(latlon_gridname,n_lon,n_lat, latlon_grid_tmp,i_off,j_off) - pert_grid = tile_grid_g + pert_grid = latlon_grid_tmp endif diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_PertTypes.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_PertTypes.F90 index 7e5d7d52..c0da62c1 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_PertTypes.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_PertTypes.F90 @@ -106,8 +106,7 @@ module LDAS_PertTypes integer :: ens_id integer :: NUM_ENSEMBLE logical :: isCubedSphere - ! pert grids - local and full - type(grid_def_type) :: pgrid_l, pgrid_f,pgrid_g + integer,allocatable :: i_indgs(:) integer,allocatable :: j_indgs(:) ! if it is cubed-sphere grid, swith to internal start diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index 414a9800..2e74f05d 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -161,9 +161,9 @@ 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) :: grid_g - type(grid_def_type) :: grid_f - type(grid_def_type) :: grid_l + type(grid_def_type) :: pgrid_g + type(grid_def_type) :: pgrid_f + type(grid_def_type) :: pgrid_l end type T_TILECOORD_STATE type TILECOORD_WRAP From 1fbc26d9d30ecae199485c66c6f1204bf1776e98 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 10 Mar 2023 11:08:09 -0500 Subject: [PATCH 02/12] rename some subroutines --- .../LDAS_App/preprocess_ldas_routines.F90 | 6 +++--- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 20 +++++++++---------- .../clsm_ensupd_upd_routines.F90 | 6 +++--- .../GEOS_LandPertGridComp.F90 | 1 - .../LDAS_PertRoutines.F90 | 6 +++--- .../Shared/LDAS_TileCoordRoutines.F90 | 12 +++++------ 6 files changed, 25 insertions(+), 26 deletions(-) diff --git a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 index 78c1dcdb..0b9c4d86 100644 --- a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 +++ b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 @@ -46,7 +46,7 @@ module preprocess_ldas_routines use LDAS_TileCoordRoutines, ONLY: & LDAS_create_grid_g, & - get_tile_grid, & + get_minimum_grid, & io_domain_files use nr_ran2_gasdev, ONLY: & @@ -537,9 +537,9 @@ subroutine domain_setup( & ! determine smallest subgrid of tile_grid_d that contains all ! catchments/tiles in domain - call get_tile_grid( N_cat_domain, tile_coord%i_indg, tile_coord%j_indg, & + tile_grid_d = get_minimum_grid( N_cat_domain, tile_coord%i_indg, tile_coord%j_indg, & tile_coord%min_lon, tile_coord%min_lat, tile_coord%max_lon, tile_coord%max_lat, & - tile_grid_g, tile_grid_d) + tile_grid_g) ! output domain files diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 54f45c8c..a9e1c094 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -18,7 +18,7 @@ module GEOS_LdasGridCompMod 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 + use LDAS_TileCoordRoutines, only: get_minimum_grid, get_ij_ind_from_latlon, io_domain_files use LDAS_ConvertMod, only: esmf2ldas use LDAS_PertRoutinesMod, only: get_pert_grid use LDAS_ensdrv_functions,ONLY: get_io_filename @@ -620,24 +620,24 @@ subroutine Initialize(gc, import, export, clock, rc) close(10) call io_grid_def_type('w', logunit, tile_grid_f, 'tile_grid_f') - call get_pert_grid(tile_grid_g, pert_grid_g) + 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 _ASSERT(index(tile_grid_g%gridtype, 'c3') /=0, "tile_grid_g does not describe a cubed-sphere grid") - !1) generate a lat-lon grid for landpert and land assim ( 4*N_lonX3*N_lon) - !2) get hash index + !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) enddo - !3) re-generate tile_grid_f in Lat-Lon - call get_tile_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_j_indg, & + !2) re-generate pert_grid_f in Lat-Lon + pert_grid_f = get_minimum_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_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, pert_grid_f) + pert_grid_g) else - call get_pert_grid(tile_grid_f, pert_grid_f) endif endif @@ -677,11 +677,11 @@ subroutine Initialize(gc, import, export, clock, rc) allocate(tcinternal%tile_coord_f,source = tile_coord_f) - call get_tile_grid(land_nt_local, & + pert_grid_l = get_minimum_grid(land_nt_local, & tcinternal%tile_coord%hash_i_indg, tcinternal%tile_coord%hash_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, pert_grid_l) + pert_grid_g) tcinternal%pgrid_g = pert_grid_g 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 d102cdfc..d9f0aeec 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 @@ -79,7 +79,7 @@ module clsm_ensupd_upd_routines get_tile_num_in_ellipse, & get_number_of_tiles_in_cell_ij, & get_tile_num_in_cell_ij, & - get_tile_grid, & + get_minimum_grid, & get_ij_ind_from_latlon use land_pert_routines, ONLY: & @@ -1513,9 +1513,9 @@ subroutine get_obs_pred( & ! determine tile_grid_lH from tile_coord_lH - call get_tile_grid( N_catlH, tile_coord_lH%hash_i_indg, tile_coord_lH%hash_j_indg, & + tile_grid_lH = get_minimum_grid( N_catlH, tile_coord_lH%hash_i_indg, tile_coord_lH%hash_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, tile_grid_lH ) + tile_grid_g) allocate(N_tile_in_cell_ij_lH(tile_grid_lH%N_lon,tile_grid_lH%N_lat)) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index 514da2e9..f6cec7a5 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -31,7 +31,6 @@ module GEOS_LandPertGridCompMod use LDAS_PertRoutinesMod, only: get_progn_pert_param use LDAS_PertRoutinesMod, only: read_ens_prop_inputs use LDAS_PertRoutinesMod, only: echo_pert_param - use LDAS_PertRoutinesMod, only: get_pert_grid use LDAS_PertRoutinesMod, only: interpolate_pert_to_timestep use LDAS_PertRoutinesMod, only: check_pert_dtstep use LDAS_PertRoutinesMod, only: GEOSldas_FORCE_PERT_DTSTEP diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 index 0e84fc44..83f6ea7e 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 @@ -630,14 +630,14 @@ subroutine check_pert_dtstep( model_dtstep, start_time, end_time, & end subroutine check_pert_dtstep ! ********************************************************************* - subroutine get_pert_grid( tile_grid, pert_grid ) + function get_pert_grid( tile_grid) result (pert_grid) ! reichle, 20 May 2010 ! jiang, 03/10/2017 implicit none type(grid_def_type), intent(in) :: tile_grid - type(grid_def_type), intent(out) :: pert_grid + type(grid_def_type) :: pert_grid character(len=30) :: latlon_gridname character(len=6) :: lattmp,lontmp @@ -683,7 +683,7 @@ subroutine get_pert_grid( tile_grid, pert_grid ) endif - end subroutine get_pert_grid + end function get_pert_grid ! ********************************************************************* diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index a059ed1b..09bf8121 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -42,7 +42,7 @@ module LDAS_TileCoordRoutines private - public :: get_tile_grid + public :: get_minimum_grid public :: get_number_of_tiles_in_cell_ij public :: get_tile_num_in_cell_ij public :: get_tile_num_in_ellipse @@ -460,8 +460,8 @@ end subroutine LDAS_create_grid_g ! ******************************************************************* - subroutine get_tile_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc_maxlon, tc_maxlat, & - tile_grid_g, tile_grid ) + function get_minimum_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 ! @@ -488,7 +488,7 @@ subroutine get_tile_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc real, intent(in), dimension(N_tile) :: tc_maxlon, tc_maxlat type(grid_def_type), intent(in) :: tile_grid_g - type(grid_def_type), intent(out) :: tile_grid + type(grid_def_type) :: tile_grid ! local variables @@ -503,7 +503,7 @@ subroutine get_tile_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc integer :: this_i_indg, this_j_indg integer , allocatable :: i_indg_(:), j_indg_(:) logical :: c3_grid - character(len=*), parameter :: Iam = 'get_tile_grid' + character(len=*), parameter :: Iam = 'get_minimum_grid' character(len=400) :: err_msg ! ------------------------------------------------- @@ -607,7 +607,7 @@ subroutine get_tile_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc tile_grid%i_dir = tile_grid_g%i_dir tile_grid%j_dir = tile_grid_g%j_dir - end subroutine get_tile_grid + end function get_minimum_grid ! ********************************************************************** From 3e00320670f60a8f9bb8396379427110208717cd Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 10 Mar 2023 17:22:34 -0500 Subject: [PATCH 03/12] removed unused arguments from get_enkf_increments(), cat_enkf_increments() --- .../GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 | 4 ++-- .../GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 | 10 +++++----- .../clsm_ensupd_upd_routines.F90 | 6 ++---- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 35458506..e747fded 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1840,8 +1840,8 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) NUM_ENSEMBLE, N_catl, N_catf, N_obsl_max, & 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_f, tcinternal%pgrid_l, tcinternal%pgrid_g, & + tile_coord_l, tile_coord_rf, & + 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, & 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 eb815f7a..55891383 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 @@ -142,8 +142,8 @@ subroutine get_enkf_increments( & N_ens, N_catl, N_catf, N_obsl_max, & work_path, exp_id, exp_domain, & met_force, lai, cat_param, mwRTM_param, & - tile_coord_l, tile_coord_f, pert_grid_f_NotUsed, & - pert_grid_f, pert_grid_l_NotUsed, pert_grid_g, & + tile_coord_l, tile_coord_f, & + 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_NotUsed, pert_grid_f, pert_grid_l_NotUsed, pert_grid_g + type(grid_def_type), intent(in) :: pert_grid_f, pert_grid_g integer, intent(in), dimension(numprocs) :: N_catl_vec, low_ind @@ -1046,7 +1046,7 @@ subroutine get_enkf_increments( & call cat_enkf_increments( & N_ens, nObs_ana, nTiles_ana, N_obs_param, & update_type, obs_param, & - pert_grid_f, tile_coord_ana, indTiles_ana, & ! indTiles_ana is essentially ana2f + tile_coord_ana, indTiles_ana, & ! indTiles_ana is essentially ana2f Obs_ana, & ! size: nObs_ana Obs_pred_ana, & ! size: (nObs_ana,N_ens) Obs_pert_tmp, & @@ -1074,7 +1074,7 @@ subroutine get_enkf_increments( & call cat_enkf_increments( & N_ens, N_obslH, N_catl, N_obs_param, & update_type, obs_param, & - pert_grid_f, tile_coord_l, l2f, & + tile_coord_l, l2f, & Observations_lH(1:N_obslH), & Obs_pred_lH(1:N_obslH,1:N_ens), & Obs_pert_tmp, & 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 d9f0aeec..bc175196 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 @@ -3375,7 +3375,7 @@ end subroutine get_obs_pert subroutine cat_enkf_increments( & N_ens, N_obs, N_catd, N_obs_param, & update_type, obs_param, & - tile_grid_f, tile_coord, l2f, & + tile_coord, l2f, & Observations, Obs_pred, Obs_pert, & cat_param, & xcompact, ycompact, fcsterr_inflation_fac, & @@ -3415,8 +3415,6 @@ subroutine cat_enkf_increments( & type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param - type(grid_def_type), intent(in) :: tile_grid_f - type(tile_coord_type), dimension(:), pointer :: tile_coord ! input integer, dimension(N_catd), intent(in) :: l2f @@ -4943,7 +4941,7 @@ subroutine check_compact_support( & ! In the obs readers, each obs is assigned to a single tile ! based on subroutine get_tile_num_from_latlon(). The assignment ! is such that the the obs lat/lon and the tile center-of-mass lat/lon - ! must always be within the same grid cell of the tile_grid. + ! must always be within the same grid cell of the tile_grid [or pert_grid??]. ! Furthermore, the obs lat/lon ! - must be within the min/max lat/lon boundaries of the tile ! OR From 4ff330ee555a6070862e87393a4eaafae341f5b0 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 10 Mar 2023 17:30:07 -0500 Subject: [PATCH 04/12] renamed get_minimum_grid() to get_minExtent_grid() [formerly known as get_tile_grid()] --- src/Applications/LDAS_App/preprocess_ldas_routines.F90 | 6 +++--- src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 6 +++--- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 4 ++-- .../GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 | 8 ++++---- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 index 0b9c4d86..3b0179a7 100644 --- a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 +++ b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 @@ -46,7 +46,7 @@ module preprocess_ldas_routines use LDAS_TileCoordRoutines, ONLY: & LDAS_create_grid_g, & - get_minimum_grid, & + get_minExtent_grid, & io_domain_files use nr_ran2_gasdev, ONLY: & @@ -537,8 +537,8 @@ subroutine domain_setup( & ! determine smallest subgrid of tile_grid_d that contains all ! catchments/tiles in domain - tile_grid_d = get_minimum_grid( N_cat_domain, tile_coord%i_indg, tile_coord%j_indg, & - tile_coord%min_lon, tile_coord%min_lat, tile_coord%max_lon, tile_coord%max_lat, & + tile_grid_d = get_minExtent_grid( N_cat_domain, tile_coord%i_indg, tile_coord%j_indg, & + tile_coord%min_lon, tile_coord%min_lat, tile_coord%max_lon, tile_coord%max_lat, & tile_grid_g) ! output domain files diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index a9e1c094..1c69475e 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -18,7 +18,7 @@ module GEOS_LdasGridCompMod 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_minimum_grid, get_ij_ind_from_latlon, io_domain_files + use LDAS_TileCoordRoutines, only: get_minExtent_grid, get_ij_ind_from_latlon, io_domain_files use LDAS_ConvertMod, only: esmf2ldas use LDAS_PertRoutinesMod, only: get_pert_grid use LDAS_ensdrv_functions,ONLY: get_io_filename @@ -633,7 +633,7 @@ subroutine Initialize(gc, import, export, clock, rc) tile_coord_f(i)%hash_i_indg,tile_coord_f(i)%hash_j_indg) enddo !2) re-generate pert_grid_f in Lat-Lon - pert_grid_f = get_minimum_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%hash_i_indg, tile_coord_f%hash_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) @@ -677,7 +677,7 @@ subroutine Initialize(gc, import, export, clock, rc) allocate(tcinternal%tile_coord_f,source = tile_coord_f) - pert_grid_l = get_minimum_grid(land_nt_local, & + pert_grid_l = get_minExtent_grid(land_nt_local, & tcinternal%tile_coord%hash_i_indg, tcinternal%tile_coord%hash_j_indg, & tcinternal%tile_coord%min_lon, tcinternal%tile_coord%min_lat, & tcinternal%tile_coord%max_lon, tcinternal%tile_coord%max_lat, & 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 bc175196..4d06f237 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 @@ -79,7 +79,7 @@ module clsm_ensupd_upd_routines get_tile_num_in_ellipse, & get_number_of_tiles_in_cell_ij, & get_tile_num_in_cell_ij, & - get_minimum_grid, & + get_minExtent_grid, & get_ij_ind_from_latlon use land_pert_routines, ONLY: & @@ -1513,7 +1513,7 @@ subroutine get_obs_pred( & ! determine tile_grid_lH from tile_coord_lH - tile_grid_lH = get_minimum_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%hash_i_indg, tile_coord_lH%hash_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) diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index 09bf8121..cb65e688 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -42,7 +42,7 @@ module LDAS_TileCoordRoutines private - public :: get_minimum_grid + public :: get_minExtent_grid public :: get_number_of_tiles_in_cell_ij public :: get_tile_num_in_cell_ij public :: get_tile_num_in_ellipse @@ -460,7 +460,7 @@ end subroutine LDAS_create_grid_g ! ******************************************************************* - function get_minimum_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, tc_maxlon, tc_maxlat, & + 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 @@ -503,7 +503,7 @@ function get_minimum_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, t integer :: this_i_indg, this_j_indg integer , allocatable :: i_indg_(:), j_indg_(:) logical :: c3_grid - character(len=*), parameter :: Iam = 'get_minimum_grid' + character(len=*), parameter :: Iam = 'get_minExtent_grid' character(len=400) :: err_msg ! ------------------------------------------------- @@ -607,7 +607,7 @@ function get_minimum_grid( N_tile, tc_i_indg, tc_j_indg, tc_minlon, tc_minlat, t tile_grid%i_dir = tile_grid_g%i_dir tile_grid%j_dir = tile_grid_g%j_dir - end function get_minimum_grid + end function get_minExtent_grid ! ********************************************************************** From f11fc4d5c1f53e2aeaee7833ebb5a4f658ee93b0 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang <52509753+weiyuan-jiang@users.noreply.github.com> Date: Fri, 10 Mar 2023 18:04:02 -0500 Subject: [PATCH 05/12] Update GEOS_LdasGridComp.F90 --- src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 1c69475e..d1bc8075 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -636,8 +636,6 @@ subroutine Initialize(gc, import, export, clock, rc) pert_grid_f = get_minExtent_grid(N_catf, tile_coord_f%hash_i_indg, tile_coord_f%hash_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 endif endif From 1021ea0d8abe1bf92503b4e228f1f8ed4f97c32a Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Wed, 15 Mar 2023 09:16:29 -0400 Subject: [PATCH 06/12] change name to tile_grid_g in output_smap --- .../clsm_ensupd_enkf_update.F90 | 48 +++++++++---------- 1 file changed, 24 insertions(+), 24 deletions(-) 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 55891383..8db54fc4 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 @@ -1831,7 +1831,7 @@ end subroutine output_incr_etc subroutine 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 ) @@ -1860,7 +1860,7 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input type(tile_coord_type), dimension(:), pointer :: tile_coord_l ! input - type(grid_def_type), intent(in) :: pert_grid_g + type(grid_def_type), intent(in) :: tile_grid_g type(grid_def_type), intent(in) :: pert_grid_f integer, dimension(numprocs), intent(in) :: N_catl_vec @@ -1945,7 +1945,7 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & ! write appropriate fields (according to 'option') into file call write_smapL4SMaup( 'orig_obs', date_time, work_path, exp_id, N_ens, & - N_catl, N_catf, N_obsl_tmp, tile_coord_f, pert_grid_g, & + N_catl, N_catf, N_obsl_tmp, tile_coord_f, tile_grid_g, & N_catl_vec, low_ind, & N_obs_param_tmp, obs_param_tmp(1:N_obs_param_tmp), Observations_l, & cat_param, cat_progn ) @@ -1955,7 +1955,7 @@ end subroutine output_smapL4SMaup ! ********************************************************************** subroutine write_smapL4SMaup( option, 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 ) ! output of custom collection for SMAP L4_SM "aup" (analysis update) @@ -2034,7 +2034,7 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input - type(grid_def_type), intent(in) :: pert_grid_g + type(grid_def_type), intent(in) :: tile_grid_g integer, dimension(numprocs), intent(in) :: N_catl_vec integer, dimension(numprocs), intent(in) :: low_ind @@ -2105,7 +2105,7 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & ! ! smapL4SMaup output only works for 9 km EASE grids - if ( index(pert_grid_g%gridtype, 'M09') == 0 ) then + 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,8 +2228,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMAP_L1C_Tbh_E09_A', 'SMAP_L1C_Tbv_E09_A' & ) - if (index(pert_grid_g%gridtype, 'M09') /=0) then - call ease_convert(trim(pert_grid_g%gridtype), this_lat, this_lon, col_ind, row_ind) + 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 @@ -2247,8 +2247,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMAP_L1C_Tbh_E27_A', 'SMAP_L1C_Tbv_E27_A' & ) - if (index(pert_grid_g%gridtype, 'M09') /=0) then - call ease_convert(trim(pert_grid_g%gridtype), this_lat, this_lon, col_ind, row_ind) + 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 @@ -2269,9 +2269,9 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & 'SMOS_fit_Tbh_A', 'SMOS_fit_Tbv_A' & ) - if (index(pert_grid_g%gridtype, 'M09') /=0) then + if (index(tile_grid_g%gridtype, 'M09') /=0) then ! subindex (1:7) to get the string EASEvx_ - gridname_tmp = pert_grid_g%gridtype(1:7)//'M36' + gridname_tmp = tile_grid_g%gridtype(1:7)//'M36' call ease_convert(gridname_tmp, this_lat, this_lon, col_ind, row_ind) endif @@ -2298,8 +2298,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & ! map Observations%[xx] fields onto global 9km EASE grid, - allocate(ndata_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) - allocate(ndata_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(ndata_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(ndata_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) if (option=='orig_obs') then @@ -2312,8 +2312,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata_8(N_obsf)) - allocate(data_h_9km_grid_8(pert_grid_g%N_lon,pert_grid_g%N_lat)) - allocate(data_v_9km_grid_8(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_h_9km_grid_8(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_v_9km_grid_8(tile_grid_g%N_lon,tile_grid_g%N_lat)) allocate(data_h_9km_tile_8(N_catf)) allocate(data_v_9km_tile_8(N_catf)) @@ -2326,8 +2326,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata(N_obsf)) - allocate(data_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) - allocate(data_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) allocate(data_h_9km_tile(N_catf)) allocate(data_v_9km_tile(N_catf)) @@ -2357,8 +2357,8 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & allocate(Obs_f_tmpdata(N_obsf)) - allocate(data_h_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) - allocate(data_v_9km_grid(pert_grid_g%N_lon,pert_grid_g%N_lat)) + allocate(data_h_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) + allocate(data_v_9km_grid(tile_grid_g%N_lon,tile_grid_g%N_lat)) allocate(data_h_9km_tile(N_catf)) allocate(data_v_9km_tile(N_catf)) @@ -2740,10 +2740,10 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & data_h_9km_tile_8 = real(nodata_generic,kind(0.0D0)) ! init (not in grid2tile!) data_v_9km_tile_8 = real(nodata_generic,kind(0.0D0)) ! init (not in grid2tile!) - call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_h_9km_grid_8, & + call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_h_9km_grid_8, & data_h_9km_tile_8 ) - call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_v_9km_grid_8, & + call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg, tile_coord_f%j_indg, data_v_9km_grid_8, & data_v_9km_tile_8 ) ! write into file @@ -2756,10 +2756,10 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & data_h_9km_tile = nodata_generic ! initialize (not done in grid2tile!) data_v_9km_tile = nodata_generic ! initialize (not done in grid2tile!) - call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_h_9km_grid, & + call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_h_9km_grid, & data_h_9km_tile ) - call grid2tile( pert_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_v_9km_grid, & + call grid2tile( tile_grid_g, N_catf, tile_coord_f%i_indg,tile_coord_f%j_indg, data_v_9km_grid, & data_v_9km_tile ) ! write into file From 55ab1f89ca0063334e183dea14f362ba58aaa03a Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 17 Mar 2023 13:54:55 -0400 Subject: [PATCH 07/12] change hash- to pert- and bring back grid_g --- .../LDAS_App/preprocess_ldas_routines.F90 | 6 +++--- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 18 +++++++++++++----- .../GEOS_LandAssimGridComp.F90 | 4 ++-- .../clsm_adapt_routines.F90 | 4 ++-- .../clsm_ensupd_enkf_update.F90 | 10 +++++----- .../clsm_ensupd_upd_routines.F90 | 6 +++--- .../GEOS_LandPertGridComp.F90 | 12 ++++++------ .../Shared/LDAS_TileCoordRoutines.F90 | 12 ++++++------ .../Shared/LDAS_TileCoordType.F90 | 19 ++++++++++--------- .../Shared/LDAS_ensdrv_mpi.F90 | 10 +++++----- 10 files changed, 55 insertions(+), 46 deletions(-) diff --git a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 index 3b0179a7..4796b4a9 100644 --- a/src/Applications/LDAS_App/preprocess_ldas_routines.F90 +++ b/src/Applications/LDAS_App/preprocess_ldas_routines.F90 @@ -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) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index d1bc8075..0403e969 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -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 @@ -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 @@ -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) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index e747fded..e88d7f5f 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%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, & @@ -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.) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_adapt_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_adapt_routines.F90 index 52603f2c..3ef58389 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_adapt_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_adapt_routines.F90 @@ -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) 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 8db54fc4..3e5bd17f 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, & + 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, & @@ -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 @@ -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 @@ -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 ) 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 4d06f237..2010d18b 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 @@ -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) @@ -1522,7 +1522,7 @@ 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 @@ -1530,7 +1530,7 @@ subroutine get_obs_pred( & ! 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 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index f6cec7a5..08f60083 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -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 @@ -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 @@ -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 diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 index cb65e688..95cb0b8d 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90 @@ -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 ! ! ------------------------------------------------------------------- @@ -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 ! ! ---------------------------------------------------------- @@ -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 ! ! ---------------------------------------------------------- diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index 2e74f05d..e2d402b4 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -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] @@ -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 @@ -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 @@ -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 diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_ensdrv_mpi.F90 index e5b26e02..bcec2953 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, 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] @@ -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 From 35258cabdc9b992abab9fb6ace9ddd080a22b540 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 18 Mar 2023 14:26:04 -0400 Subject: [PATCH 08/12] additional clarification of tile_grid vs pert_grid --- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 53 +++++++++++++------ .../GEOS_LandAssimGridComp.F90 | 4 +- .../clsm_ensupd_enkf_update.F90 | 6 +-- .../Shared/LDAS_TileCoordRoutines.F90 | 6 +-- .../Shared/LDAS_TileCoordType.F90 | 8 +-- .../Shared/LDAS_ensdrv_mpi.F90 | 6 +-- 6 files changed, 52 insertions(+), 31 deletions(-) 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] From 14dca21aee23f42080072998e8d243f0eeaafbb9 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 18 Mar 2023 14:41:36 -0400 Subject: [PATCH 09/12] fixed build error in previous commit --- src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 0f95ac78..801cc943 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -17,7 +17,7 @@ module GEOS_LdasGridCompMod 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_TileCoordType, only: grid_def_type, io_grid_def_type, operator (==) use LDAS_TileCoordRoutines, only: get_minExtent_grid, get_ij_ind_from_latlon, io_domain_files use LDAS_ConvertMod, only: esmf2ldas use LDAS_PertRoutinesMod, only: get_pert_grid From bc854c2d7e4664164a0325a9c41a49e8e55fdc1d Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 19 Mar 2023 07:18:01 -0400 Subject: [PATCH 10/12] fixed error in previous commit (35258ca) --- src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 801cc943..932093e2 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -654,7 +654,9 @@ subroutine Initialize(gc, import, export, clock, rc) else - ! do nothing because %pert_i_indg and %pert_j_indg were initialized to %i_indg and %j_indg + pert_grid_f = tile_grid_f + + ! note that %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 From f027dc4f76f0945de260552a5c75cb1151b701a1 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang <52509753+weiyuan-jiang@users.noreply.github.com> Date: Mon, 20 Mar 2023 14:30:12 -0400 Subject: [PATCH 11/12] Update GEOS_LdasGridComp.F90 --- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index 932093e2..a8ff3631 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -632,11 +632,6 @@ subroutine Initialize(gc, import, export, clock, rc) ! 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) -- 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 @@ -668,10 +663,7 @@ 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 block integer, allocatable :: f2tile_id(:), tile_id2f(:) @@ -711,11 +703,7 @@ 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%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) From f1922cf84ea2dc59ea55771cbae3862e2c8298ac Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 20 Mar 2023 14:34:42 -0400 Subject: [PATCH 12/12] modify grid == --- src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 index df8bfdc0..4861ddce 100644 --- a/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordType.F90 @@ -214,7 +214,7 @@ function eq_grid_def_type( grid_1, grid_2 ) type(grid_def_type), intent(in) :: grid_1, grid_2 - if ( index(grid_1%gridtype,trim(grid_2%gridtype))>0 .and. & + if ( trim(grid_1%gridtype) == trim(grid_2%gridtype) .and. & grid_1%ind_base == grid_2%ind_base .and. & grid_1%i_dir == grid_2%i_dir .and. & grid_1%j_dir == grid_2%j_dir .and. &