From e01aa2e09004508ebb8778dd292ed9be8a42a80b Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Wed, 18 Oct 2023 15:24:41 -0400 Subject: [PATCH 1/7] clean up and add check of FFT --- .../GEOSldas_GridComp/GEOS_LdasGridComp.F90 | 6 +- .../GEOS_LandAssimGridComp.F90 | 50 +++--------- .../clsm_ensupd_enkf_update.F90 | 76 ++++++------------- .../clsm_ensupd_read_obs.F90 | 20 +---- .../clsm_ensupd_upd_routines.F90 | 23 +++--- .../GEOS_LandPertGridComp.F90 | 17 ++--- .../LDAS_PertRoutines.F90 | 5 +- .../GEOSlandpert_GridComp/land_pert.F90 | 2 - .../GEOSlandpert_GridComp/random_fields.F90 | 29 +++++-- 9 files changed, 82 insertions(+), 146 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 index a8ff3631..c0608c08 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -80,8 +80,6 @@ subroutine SetServices(gc, rc) integer :: i, k integer :: ens_id type(MAPL_MetaComp), pointer :: MAPL=>null() - type(ESMF_GridComp), pointer :: gcs(:)=>null() ! Children gridcomps - character(len=ESMF_MAXSTR), pointer :: gcnames(:)=>null() ! Children's names ! ErrLog variables integer :: status character(len=ESMF_MAXSTR) :: Iam @@ -400,7 +398,6 @@ subroutine Initialize(gc, import, export, clock, rc) logical :: IamRoot type(tile_coord_type), dimension(:), pointer :: tile_coord_f => null() - type(tile_coord_type), dimension(:), pointer :: tile_coord_l => null() integer,dimension(:),pointer :: f2g integer :: N_catf @@ -408,7 +405,7 @@ subroutine Initialize(gc, import, export, clock, rc) 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(grid_def_type) :: pert_grid_l type(date_time_type):: start_time type(ESMF_Time) :: CurrentTime @@ -827,7 +824,6 @@ subroutine Run(gc, import, export, clock, rc) ! Misc variables integer :: igc,i, ens_id, FIRST_ENS_ID, ens_id_width logical :: IAmRoot - integer :: mpierr integer :: LSM_CHOICE type (ESMF_Field) :: field diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index f13c1bfc..a44ef2f1 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -146,9 +146,8 @@ subroutine SetServices ( GC, RC ) ! Local Variables type(MAPL_MetaComp), pointer :: MAPL=>null() - type(ESMF_Config) :: CF character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file - character(len=ESMF_MAXSTR) :: ensid_string,childname, fmt_str + character(len=ESMF_MAXSTR) :: ensid_string,childname integer :: i, ens_id_width, FIRST_ENS_ID, NUM_ENSEMBLE integer :: ens_id, export_id @@ -1111,7 +1110,7 @@ subroutine Initialize(gc, import, export, clock, rc) type(MAPL_MetaComp), pointer :: CHILD_MAPL=>null() ! Child's MAPL obj type(ESMF_GridComp), pointer :: gcs(:) - character(len=300) :: out_path,fname + character(len=300) :: out_path character(len=ESMF_MAXSTR) :: exp_id, GridName integer :: model_dtstep type(date_time_type) :: start_time @@ -1120,7 +1119,6 @@ subroutine Initialize(gc, import, export, clock, rc) type(T_TILECOORD_STATE), pointer :: tcinternal type(TILECOORD_WRAP) :: tcwrap - type(tile_coord_type), dimension(:), pointer :: tile_coord_f => null() type(tile_coord_type), dimension(:), pointer :: tile_coord_l => null() integer :: land_nt_local,i,mpierr, ens, ens_id_width @@ -1129,7 +1127,6 @@ subroutine Initialize(gc, import, export, clock, rc) integer, allocatable :: f2rf(:) ! mapping re-orderd rf to f for the LDASsa output character(len=300) :: seed_fname character(len=300) :: fname_tpl - character(len=14) :: datestamp character(len=ESMF_MAXSTR) :: ensid_string integer :: nymd, nhms, yy, mm, dd, h, m, s @@ -1464,7 +1461,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! ! time ! - type(ESMF_Time) :: ModelTimeCur, ModelTimeNxt + type(ESMF_Time) :: ModelTimeCur type(ESMF_Alarm) :: LandAssimAlarm type(ESMF_TimeInterval) :: ModelTimeStep @@ -1484,8 +1481,6 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) character(len=300) :: out_path character(len=ESMF_MAXSTR) :: exp_id - character(40) :: exp_domain - integer :: model_dtstep type(met_force_type), dimension(:), allocatable :: met_force @@ -1905,12 +1900,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) call get_enkf_increments( & date_time_new, & NUM_ENSEMBLE, N_catl, N_catf, N_obsl_max, & - trim(out_path), trim(exp_id), exp_domain, & + trim(out_path), trim(exp_id), & met_force, lai, cat_param, mwRTM_param, & tile_coord_l, tile_coord_rf, & 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, & LandAssimDTstep, & xcompact, ycompact, fcsterr_inflation_fac, & @@ -1942,14 +1936,14 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) if (.true.) then ! replace obsolete check for analysis time with "if true" to keep indents call output_incr_etc( out_ObsFcstAna, & - date_time_new, trim(out_path), trim(exp_id), & + date_time_new, trim(exp_id), & N_obsl, N_obs_param, NUM_ENSEMBLE, & N_catl, tile_coord_l, & - N_catf, tile_coord_rf, tcinternal%pgrid_f, tcinternal%pgrid_g, & - N_catl_vec, low_ind, rf2l, N_catg, rf2g, & + N_catf, tile_coord_rf, tcinternal%pgrid_g, & + N_catl_vec, low_ind, rf2l, & obs_param, & met_force, lai, & - cat_param, cat_progn, cat_progn_incr, mwRTM_param, & + cat_param, cat_progn, mwRTM_param, & Observations_l, rf2f=rf2f ) do ii = 1, N_catl @@ -2058,9 +2052,9 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! increments were computed) if (out_smapL4SMaup) & - call write_smapL4SMaup( 'analysis', date_time_new, trim(out_path), & + call write_smapL4SMaup( 'analysis', date_time_new, & trim(exp_id), NUM_ENSEMBLE, N_catl, N_catf, N_obsl, tile_coord_rf, & - tcinternal%tgrid_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.) @@ -2103,7 +2097,6 @@ subroutine UPDATE_ASSIM(gc, import, export, clock, rc) ! ESMF variables type(ESMF_Alarm) :: LandAssimAlarm - type(ESMF_VM) :: vm ! MAPL variables type(MAPL_MetaComp), pointer :: MAPL=>null() ! MAPL obj @@ -2268,25 +2261,6 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) real, dimension(:), pointer :: WCSF real, dimension(:), pointer :: SWE - real, dimension(:), pointer :: VEGCLS - real, dimension(:), pointer :: SOILCLS - real, dimension(:), pointer :: SAND - real, dimension(:), pointer :: CLAY - real, dimension(:), pointer :: mw_POROS - real, dimension(:), pointer :: WANGWT - real, dimension(:), pointer :: WANGWP - real, dimension(:), pointer :: RGHHMIN - real, dimension(:), pointer :: RGHHMAX - real, dimension(:), pointer :: RGHWMAX - real, dimension(:), pointer :: RGHWMIN - real, dimension(:), pointer :: RGHNRH - real, dimension(:), pointer :: RGHNRV - real, dimension(:), pointer :: RGHPOLMIX - real, dimension(:), pointer :: OMEGA - real, dimension(:), pointer :: BH - real, dimension(:), pointer :: BV - real, dimension(:), pointer :: LEWT - ! export real, dimension(:), pointer :: TB_H_enavg real, dimension(:), pointer :: TB_V_enavg @@ -2298,7 +2272,7 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) real, allocatable, dimension(:) :: Tb_h_tmp, TB_v_tmp - integer :: N_catl, n, mpierr + integer :: N_catl type(MAPL_LocStream) :: locstream call ESMF_GridCompGet ( GC, name=COMP_NAME, RC=STATUS ) @@ -2647,7 +2621,7 @@ subroutine read_pert_rseed(ensid_string,seed_fname,pert_rseed_r8) character(len=*),intent(in) :: seed_fname real(kind=ESMF_KIND_R8),intent(inout) :: pert_rseed_r8(:) - integer :: ncid, s_varid, en_dim, n_ens, id_varid, i, pos + integer :: ncid, s_varid logical :: file_exist character(len=ESMF_MAXSTR) :: tmpstr 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 86d6449d..24d62c58 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 @@ -140,12 +140,11 @@ module clsm_ensupd_enkf_update subroutine get_enkf_increments( & date_time, & N_ens, N_catl, N_catf, N_obsl_max, & - work_path, exp_id, exp_domain, & + work_path, exp_id, & met_force, lai, cat_param, mwRTM_param, & tile_coord_l, tile_coord_f, & 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, & dtstep_assim, & xcompact, ycompact, fcsterr_inflation_fac, & @@ -177,7 +176,7 @@ subroutine get_enkf_increments( & integer, intent(in) :: N_obsl_max ! max number of observations allowed character(*), intent(in) :: work_path - character(*), intent(in) :: exp_id, exp_domain + character(*), intent(in) :: exp_id ! Meteorological forcings, Catchment model and microwave RTM parameters @@ -200,11 +199,6 @@ subroutine get_enkf_increments( & integer, intent(in), dimension(N_catf) :: f2l - integer, intent(in) :: N_force_pert, N_progn_pert - - type(pert_param_type), dimension(:), pointer :: force_pert_param ! input - type(pert_param_type), dimension(:), pointer :: progn_pert_param ! input - integer, intent(in) :: update_type, dtstep_assim real, intent(in) :: xcompact, ycompact, fcsterr_inflation_fac @@ -252,8 +246,6 @@ subroutine get_enkf_increments( & ! local variables - integer :: N_obslH - logical :: found_obs_f, assimflag type(obs_type), dimension(:), pointer :: Observations_lH => null() ! obs w/in halo @@ -417,8 +409,8 @@ 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, pert_grid_f, & + N_ens, N_catl, N_catf, N_obsl_max, & + tile_coord_f, 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 ) @@ -426,7 +418,7 @@ subroutine get_enkf_increments( & call collect_obs( & work_path, exp_id, date_time, dtstep_assim, & - N_catl, tile_coord_l, & + N_catl, & 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, & @@ -986,7 +978,7 @@ subroutine get_enkf_increments( & ! of Observations_l and Obs_pred_l that are "good" ! [allocation of these arrays in get_obs_pred() is larger ! than eventual size] - call get_halo_obs( N_ens, N_catl, N_obsl, & + call get_halo_obs( N_ens, N_obsl, & Observations_l(1:N_obsl), Obs_pred_l(1:N_obsl,1:N_ens), & tile_coord_l, xcompact, ycompact, & N_obslH, Observations_lH, Obs_pred_lH ) @@ -1195,7 +1187,7 @@ subroutine get_enkf_increments( & ! into SMAP L4_SM aup file if (out_smapL4SMaup) & - call write_smapL4SMaup( 'obs_fcst', date_time, work_path, exp_id, N_ens, & + call write_smapL4SMaup( 'obs_fcst', date_time, exp_id, N_ens, & 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 ) @@ -1271,7 +1263,6 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & logical :: cat_progn_has_changed character(len=*), parameter :: Iam = 'apply_enkf_increments' - character(len=400) :: err_msg ! ---------------------------------------------------------------- ! @@ -1393,8 +1384,8 @@ end subroutine apply_enkf_increments ! ******************************************************************** - subroutine output_ObsFcstAna(date_time, work_path, exp_id, & - N_obsl, Observations_l, N_obs_param, obs_param, rf2f) + subroutine output_ObsFcstAna(date_time, exp_id, & + N_obsl, Observations_l, N_obs_param, rf2f) ! obs space output: observations, obs space forecast, obs space analysis, and ! associated error variances @@ -1403,16 +1394,14 @@ subroutine output_ObsFcstAna(date_time, work_path, exp_id, & implicit none - character(*), intent(in) :: work_path + type(date_time_type), intent(in) :: date_time character(*), intent(in) :: exp_id integer, intent(in) :: N_obsl, N_obs_param - type(date_time_type), intent(in) :: date_time type(obs_type), dimension(N_obsl), intent(in) :: Observations_l - type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param integer, dimension(:), optional, intent(in) :: rf2f ! --------------------- @@ -1430,7 +1419,6 @@ subroutine output_ObsFcstAna(date_time, work_path, exp_id, & integer, dimension(numprocs) :: N_obsl_vec, tmp_low_ind character(300) :: fname - integer :: i #ifdef LDAS_MPI integer :: this_species, ind_tmp, j @@ -1623,13 +1611,13 @@ end subroutine output_ObsFcstAna ! ********************************************************************** subroutine output_incr_etc( out_ObsFcstAna, & - date_time, work_path, exp_id, & + date_time, exp_id, & N_obsl, N_obs_param, N_ens, & N_catl, tile_coord_l, & - N_catf, tile_coord_f, pert_grid_f, pert_grid_g, & - N_catl_vec, low_ind, f2l, N_catg, f2g, & + N_catf, tile_coord_f, pert_grid_g, & + N_catl_vec, low_ind, f2l, & obs_param, & - met_force, lai, cat_param, cat_progn, cat_progn_incr, mwRTM_param, & + met_force, lai, cat_param, cat_progn, mwRTM_param, & Observations_l, rf2f ) implicit none @@ -1646,23 +1634,19 @@ subroutine output_incr_etc( out_ObsFcstAna, & type(date_time_type), intent(in) :: date_time - character(len=*), intent(in) :: work_path character(len=*), intent(in) :: exp_id - integer, intent(in) :: N_obsl, N_obs_param, N_ens, N_catl, N_catf, N_catg + integer, intent(in) :: N_obsl, N_obs_param, N_ens, N_catl, N_catf 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) :: pert_grid_f type(grid_def_type), intent(in) :: pert_grid_g integer, dimension(numprocs), intent(in) :: N_catl_vec, low_ind integer, dimension(N_catf), intent(in) :: f2l - integer, dimension(N_catf), intent(in) :: f2g - type(obs_param_type), dimension(N_obs_param), intent(in) :: & obs_param @@ -1672,7 +1656,6 @@ subroutine output_incr_etc( out_ObsFcstAna, & type(cat_param_type), dimension(N_catl), intent(in) :: cat_param type(cat_progn_type), dimension(N_catl,N_ens), intent(in) :: cat_progn - type(cat_progn_type), dimension(N_catl,N_ens), intent(in) :: cat_progn_incr type(mwRTM_param_type), dimension(N_catl), intent(in) :: mwRTM_param @@ -1685,19 +1668,10 @@ subroutine output_incr_etc( out_ObsFcstAna, & real, dimension(:,:), pointer :: Obs_pred_l => null() - integer :: i, n_e, N_obsl_tmp - - type(cat_progn_type), dimension(N_catl) :: cat_progn_incr_ensavg - - type(cat_progn_type), dimension(:), allocatable :: cat_progn_incr_f - type(cat_progn_type), dimension(:), allocatable :: cat_progn_incr_tmp + integer :: N_obsl_tmp - type(cat_progn_type), dimension(:), allocatable :: cat_progn_incr_g - - character(40) :: file_tag, dir_name character(len=*), parameter :: Iam = 'output_incr_etc' - character(len=400) :: err_msg ! -------------------------------------------------------------- @@ -1727,8 +1701,8 @@ subroutine output_incr_etc( out_ObsFcstAna, & ! write out model, observations, and "OminusA" information - call output_ObsFcstAna( date_time, work_path, exp_id, N_obsl, & - Observations_l(1:N_obsl), N_obs_param, obs_param, rf2f=rf2f ) + call output_ObsFcstAna( date_time, exp_id, N_obsl, & + Observations_l(1:N_obsl), N_obs_param, rf2f=rf2f ) end if @@ -1830,8 +1804,8 @@ 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, pert_grid_f, & + N_ens, N_catl, N_catf, N_obsl_max, & + tile_coord_f, 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 ) @@ -1855,10 +1829,9 @@ subroutine output_smapL4SMaup( date_time, work_path, exp_id, dtstep_assim, & integer, intent(in) :: dtstep_assim integer, intent(in) :: N_ens, N_catl, N_catf - integer, intent(in) :: N_obsl, N_obsl_max, N_obs_param + integer, intent(in) :: N_obsl_max, N_obs_param 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) :: pert_grid_f @@ -1935,7 +1908,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_catl, & 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, & @@ -1944,7 +1917,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, & + call write_smapL4SMaup( 'orig_obs', date_time, exp_id, N_ens, & 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, & @@ -1954,7 +1927,7 @@ end subroutine output_smapL4SMaup ! ********************************************************************** - subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & + subroutine write_smapL4SMaup( option, date_time, exp_id, N_ens, & 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 ) @@ -2026,7 +1999,6 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & type(date_time_type), intent(in) :: date_time - character(*), intent(in) :: work_path character(*), intent(in) :: exp_id integer, intent(in) :: N_ens, N_catl, N_catf diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 5937a6e0..50e0fbf8 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -496,7 +496,6 @@ subroutine read_obs_ae_l2_sm( & integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_ae_l2_sm' - character(len=400) :: err_msg ! ------------------------------------------------------------------- @@ -1009,7 +1008,6 @@ subroutine read_obs_ae_sm_LPRM( & integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_ae_sm_LPRM' - character(len=400) :: err_msg ! ------------------------------------------------------------------- @@ -1320,7 +1318,6 @@ subroutine read_obs_sm_ASCAT( & real, parameter :: tol = 1e-2 character(len=*), parameter :: Iam = 'read_obs_sm_ASCAT' - character(len=400) :: err_msg ! ------------------------------------------------------------------- @@ -1431,7 +1428,7 @@ subroutine read_obs_sm_ASCAT( & if (N_files>0) then call read_sm_ASCAT_bin( & - this_obs_param, N_files, fnames, & + N_files, fnames, & N_tmp, tmp_lon, tmp_lat, tmp_obs ) if (logit) then @@ -1552,7 +1549,7 @@ end subroutine read_obs_sm_ASCAT ! *************************************************************************** subroutine read_sm_ASCAT_bin( & - this_obs_param, N_files, fnames, N_data, lon, lat, sm_ASCAT, ease_col, ease_row ) + N_files, fnames, N_data, lon, lat, sm_ASCAT, ease_col, ease_row ) ! read soil moisture data from one or more ASCAT bin files ! @@ -1569,8 +1566,6 @@ subroutine read_sm_ASCAT_bin( & implicit none - type(obs_param_type), intent(in) :: this_obs_param - integer, intent(in) :: N_files character(*), dimension(N_files), intent(in) :: fnames @@ -1740,7 +1735,6 @@ end subroutine read_sm_ASCAT_bin ! ***************************************************************** subroutine read_obs_LaRC_Tskin( & - work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & @@ -1757,9 +1751,6 @@ subroutine read_obs_LaRC_Tskin( & ! inputs: - character(*), intent(in) :: work_path - character(*), intent(in) :: exp_id - type(date_time_type), intent(in) :: date_time integer, intent(in) :: dtstep_assim, N_catd @@ -1817,7 +1808,6 @@ subroutine read_obs_LaRC_Tskin( & integer :: MM character(len=*), parameter :: Iam = 'read_obs_LaRC_Tskin' - character(len=400) :: err_msg ! ------------------------------------------------------------------- @@ -6370,7 +6360,6 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation logical, dimension(:,:), allocatable :: mask_h_D, mask_v_D character(len=*), parameter :: Iam = 'turn_off_assim_SMAP_L1CTb' - character(len=400) :: err_msg ! ------------------------------------------------------------------------------ @@ -7238,7 +7227,6 @@ subroutine read_obs( & 'LaRC_tskin-FY2E-', 'LaRC_tskin-MTST2') call read_obs_LaRC_Tskin( & - work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & @@ -8074,7 +8062,7 @@ end subroutine scale_obs_Tb_zscore subroutine collect_obs( & work_path, exp_id, date_time, dtstep_assim, & - N_catl, tile_coord_l, & + N_catl, & N_catf, tile_coord_f, tile_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, write_obslog, & @@ -8109,7 +8097,7 @@ subroutine collect_obs( ! tile_coord_f of catchments in domain (length N_catf) - type(tile_coord_type), dimension(:), pointer :: tile_coord_l, tile_coord_f ! input + type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input type(grid_def_type), intent(in) :: tile_grid_f 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 6b8726ad..d7940820 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 @@ -1057,7 +1057,6 @@ subroutine get_obs_pred( & character(len=*), parameter :: Iam = 'get_obs_pred' character(len=400) :: err_msg - character(len= 10) :: tmpstring10 ! -------------------------------------------------------------- ! @@ -1423,7 +1422,7 @@ subroutine get_obs_pred( & ! allocate and assemble tile_data_l allocate(tile_data_l(0,0,0)) ! for debugging to pass call get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - N_catf, tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_coord_lH=tile_coord_lH ) if (get_sfmc_lH) allocate(sfmc_lH( N_catlH, N_ens)) @@ -1454,7 +1453,7 @@ subroutine get_obs_pred( & ! communicate tile_data_l as needed and get tile_data_lH call get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - N_catf, tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_data_lH=tile_data_lH ) ! read out sfmc, rzmc, etc. from tile_data_lH @@ -1977,7 +1976,6 @@ subroutine get_obs_pred_comm_helper( & integer :: k, ks, opt character(len=*), parameter :: Iam = 'get_obs_pred_comm_helper' - character(len=400) :: err_msg ! ------------------------------------------------------------------------- @@ -2342,7 +2340,7 @@ end subroutine qc_model_based_for_Tb ! ********************************************************************* - subroutine get_halo_obs( N_ens, N_catl, N_obsl, Observations_l, Obs_pred_l, & + subroutine get_halo_obs( N_ens, N_obsl, Observations_l, Obs_pred_l, & tile_coord_l, xcompact, ycompact, & N_obslH, Observations_lH, Obs_pred_lH ) @@ -2375,7 +2373,7 @@ subroutine get_halo_obs( N_ens, N_catl, N_obsl, Observations_l, Obs_pred_l, & implicit none - integer, intent(in) :: N_ens, N_catl, N_obsl + integer, intent(in) :: N_ens, N_obsl type(obs_type), dimension(N_obsl), intent(in) :: Observations_l @@ -2766,7 +2764,7 @@ end subroutine get_halo_obs ! ********************************************************************* subroutine get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - N_catf, tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_coord_lH, tile_data_lH ) ! collect (bundled) tile_data from other local domains (processors) that are @@ -2787,7 +2785,7 @@ subroutine get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l implicit none integer, intent(in) :: N_catl, N_fields - integer, intent(in) :: N_ens, N_catf + integer, intent(in) :: N_ens real, dimension(N_catl,N_fields,N_ens), intent(in) :: tile_data_l @@ -3453,7 +3451,6 @@ subroutine cat_enkf_increments( & integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij => null() character(len=*), parameter :: Iam = 'cat_enkf_increments' - character(len=400) :: err_msg real, dimension( N_catd) :: r_x, tmp_dlon real :: r_y, tmp_dlat @@ -4527,10 +4524,10 @@ subroutine get_ind_obs( & ! locals - integer :: i, j, k, m - - logical :: selected_obs + integer :: i, k + + ! -------------------------------------------------------------- if (N_select_species==0 .and. N_select_tilenum==0) then @@ -4802,7 +4799,7 @@ subroutine get_ind_obs_lat_lon_box( & ! locals - integer :: i, j, k + integer :: i, k real :: lon_obs, lat_obs diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index 08f60083..17470a5c 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -859,13 +859,11 @@ subroutine Initialize(gc, import, export, clock, rc) real(kind=ESMF_KIND_R8), pointer :: pert_rseed_r8(:)=>null() ! Misc variables - integer :: imjm(7), imjm_global(7) ! we need just the first 2 integer :: model_dtstep - integer :: land_nt_local,m,n, i1, in, j1, jn + integer :: land_nt_local,n,m,i1, in, j1, jn logical :: IAmRoot, f_exist - integer :: ipert,n_lon,n_lat, n_lon_g, n_lat_g + integer :: n_lon,n_lat, n_lon_g, n_lat_g integer, allocatable :: pert_rseed(:) - real :: dlon, dlat,locallat,locallon type(ESMF_Grid) :: Grid character(len=ESMF_MAXSTR) :: id_string integer :: ens_id_width @@ -1320,9 +1318,8 @@ subroutine Phase2_Initialize(gc, import, export, clock, rc) real, allocatable :: fpert_grid(:,:,:), ppert_grid(:,:,:) integer,allocatable :: pert_rseed(:) - integer :: land_nt_local,n,ipert,i,j,n_lon,n_lat + integer :: land_nt_local,ipert,n_lon,n_lat logical :: IAmRoot - real :: locallat,locallon ! Begin... ! phase2_initialized is a global variables shared by all ensemble member @@ -1837,12 +1834,11 @@ subroutine ApplyForcePert(gc, import, export, clock, rc) type(tile_coord_type), pointer :: tile_coord(:)=>null() integer :: n_lon,n_lat - integer :: ipert, itile + integer :: ipert integer, allocatable :: pert_rseed(:) logical :: IAmRoot integer :: land_nt_local type(pert_param_type), pointer :: PertParam=>null() ! pert param - real :: tmpRealArrDim1(1) real, allocatable :: tmpreal(:) ! Begin... @@ -2275,14 +2271,13 @@ subroutine ApplyPrognPert(gc, import, export, clock, rc) type(tile_coord_type), pointer :: tile_coord(:)=>null() integer :: n_lon,n_lat - integer :: ipert,ntiles + integer :: ipert integer, allocatable :: pert_rseed(:) logical :: IAmRoot integer :: land_nt_local type(pert_param_type), pointer :: PertParam=>null() ! pert param integer :: model_dtstep real :: dtmh - integer :: m ! Begin... ! Get my name and setup traceback handle @@ -2667,7 +2662,7 @@ subroutine Finalize(gc, import, export, clock, rc) 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 :: m,land_nt_local, ens_id_width integer :: nfpert, nppert, n_tile type(tile_coord_type), pointer :: tile_coord_f(:)=>null() diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 index 00eece74..cac58cdb 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 @@ -287,7 +287,6 @@ subroutine read_ens_prop_inputs( & type(force_pert_ccor_type) :: ccorr_force_pert ! Errlong variables - integer :: rc, status character(len=*), parameter :: Iam = 'read_ens_prop_inputs' ! MPI variables @@ -995,7 +994,6 @@ subroutine get_force_pert_inputs( pert_grid_l, & integer :: xstart, xcount, ystart, ycount ! for computing local indices ! Errlong variables - integer :: rc, status character(len=*), parameter :: Iam = 'get_force_pert_inputs' ! MPI variables @@ -1385,7 +1383,6 @@ subroutine get_progn_pert_inputs( pert_grid_l, & integer :: xstart, xcount, ystart, ycount ! for computing local indices ! Errlong variables - integer :: rc, status character(len=*), parameter :: Iam = 'get_progn_pert_inputs' ! MPI variables @@ -1902,7 +1899,7 @@ subroutine get_pert_select( N_pert_max, pert_grid_l, std_pert, & ! ! local variables - integer :: i,j,k, ierr + integer :: k, ierr ! --------------------------------------------------------------- ! diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 index 1a6fe537..5b9323d7 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 @@ -361,8 +361,6 @@ subroutine get_init_Pert_rseed( ens_id, init_Pert_rseed ) ! ! local variables - integer :: m, n, i, j - character(len=*), parameter :: Iam = 'get_init_Pert_rseed' character(len=400) :: err_msg diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 index 61252d50..835f469a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 @@ -45,7 +45,7 @@ module random_fields_class NRANDSEED, & nr_ran2_2d, & nr_gasdev - + implicit none private @@ -98,9 +98,9 @@ subroutine initialize(this, Nx, Ny, var, lx, ly, dx, dy, comm) integer, optional, intent(in) :: comm ! local variable - integer :: mklstat, ierror, i, j + integer :: mklstat, ierror integer :: rank, npes, local_dim1, local_dim2, remainer - integer :: N1, N2, gcount(2), lstart(2), Stride(2) + integer :: N1, N2, Stride(2) ! set obj param vals this%N_x = Nx @@ -134,6 +134,12 @@ subroutine initialize(this, Nx, Ny, var, lx, ly, dx, dy, comm) N1 = this%fft_lens(1) N2 = this%fft_lens(2) + if (npes > minval([N1, N2]) ) then + print*, " Two many processors are acquired in a node for parallel FFT" + print*, " The number of processors acquired in a node should be smaller than or equal to FFT grid size: ", minval([N1, N2]) + call quit('Parallel FFT failed') + endif + call this%win_allocate(N1, N2) ! distribution of the grid for fft @@ -153,23 +159,35 @@ subroutine initialize(this, Nx, Ny, var, lx, ly, dx, dy, comm) mklstat = DftiCreateDescriptor(this%Desc_Handle_Dim1, DFTI_SINGLE,& DFTI_COMPLEX, 1, N1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiCreate dim1 failed!') mklstat = DftiCreateDescriptor(this%Desc_Handle_Dim2, DFTI_SINGLE,& DFTI_COMPLEX, 1, N2 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiCreate dim2 failed!') ! perform local_dim2 one-dimensional transforms along 1st dimension mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_NUMBER_OF_TRANSFORMS, local_dim2 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue DFTI_NUMBER_OF_TRANSFORMS failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_INPUT_DISTANCE, N1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue N1 failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_OUTPUT_DISTANCE, N1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue N2 failed!') mklstat = DftiCommitDescriptor( this%Desc_Handle_Dim1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiCommit dim1 failed!') ! mklstat = DftiComputeForward( this%Desc_Handle_Dim1, X ) ! local_dim1 one-dimensional transforms along 2nd dimension Stride(1) = 0; Stride(2) = local_dim1 mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_NUMBER_OF_TRANSFORMS, local_dim1) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue local_dim1 failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_INPUT_DISTANCE, 1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue in distance failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_OUTPUT_DISTANCE, 1 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue out distance failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_INPUT_STRIDES, Stride ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue in_stride failed!') mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_OUTPUT_STRIDES, Stride ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue outstride failed!') mklstat = DftiCommitDescriptor( this%Desc_Handle_Dim2 ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiCommit Dim2 failed!') !mklstat = DftiComputeForward( this%Desc_Handle_Dim2, X ) else this%comm = MPI_COMM_NULL @@ -194,7 +212,7 @@ subroutine finalize(this) class(random_fields), intent(inout) :: this ! local variable - integer :: mklstat, i , npes, ierror + integer :: mklstat ! deallocate memory if(allocated(this%field1_fft)) deallocate(this%field1_fft) @@ -515,6 +533,7 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) cptr = c_loc(tmp_field_dim1(1,1)) call c_f_pointer (cptr, X, [ldim1*N_y_fft]) mklstat = DftiComputeBackward( this%Desc_Handle_Dim2, X ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim2 failed!') call MPI_Barrier(this%node_comm, ierror) tmp_field(n1:n2,:) = tmp_field_dim1 @@ -528,6 +547,7 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) cptr = c_loc(tmp_field_dim2(1,1)) call c_f_pointer (cptr, X, [N_x_fft*ldim2]) mklstat = DftiComputeBackward( this%Desc_Handle_Dim1, X ) + if (mklstat/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim1 failed!') tmp_field(:,n1:n2) = tmp_field_dim2/N_xy_fft_real call MPI_Win_fence(0, this%win, ierror) @@ -648,7 +668,6 @@ subroutine win_allocate(this, nx, ny) integer(kind=MPI_ADDRESS_KIND) :: windowsize integer :: disp_unit,status, Rank integer(kind=MPI_ADDRESS_KIND) :: n_bytes - integer :: int_size call MPI_Comm_rank( this%node_comm, rank, status) From 446041b2dd1303fb5afce520668e55330bf993b0 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 23 Oct 2023 09:59:28 -0400 Subject: [PATCH 2/7] 1)save FFT allocation 2) add max_cat_daigS --- .../GEOS_LandAssimGridComp.F90 | 3 +- .../GEOS_LandPertGridComp.F90 | 3 +- .../LDAS_PertRoutines.F90 | 4 +- .../GEOSlandpert_GridComp/land_pert.F90 | 128 +++-- .../GEOSlandpert_GridComp/random_fields.F90 | 466 +++++++++--------- .../GEOSldas_GridComp/Shared/catch_types.F90 | 33 ++ 6 files changed, 355 insertions(+), 282 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index a44ef2f1..4d4d67c5 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -51,6 +51,7 @@ module GEOS_LandAssimGridCompMod use catch_types, only: cat_progn_type use catch_types, only: cat_param_type use catch_types, only: cat_diagS_type + use catch_types, only: max_cat_diagS use catch_types, only: cat_diagS_sqrt use catch_types, only: assignment(=), operator (+), operator (-), operator (*), operator (/) use clsm_bias_routines, only: initialize_obs_bias @@ -2009,7 +2010,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg - cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii)) ) + cat_diagS_ensstd(ii) = cat_diagS_sqrt( max_cat_diagS(0.0, cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))) ) end do diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index 17470a5c..94d542c0 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -24,7 +24,7 @@ module GEOS_LandPertGridCompMod use LDAS_TileCoordType, only: tile_coord_type use LDAS_TileCoordType, only: T_TILECOORD_STATE use LDAS_TileCoordType, only: TILECOORD_WRAP - use land_pert_routines, only: get_pert, propagate_pert + use land_pert_routines, only: get_pert, propagate_pert, clear_rf use land_pert_routines, only: get_init_pert_rseed use LDAS_PertRoutinesMod, only: apply_pert use LDAS_PertRoutinesMod, only: get_force_pert_param @@ -2807,6 +2807,7 @@ subroutine Finalize(gc, import, export, clock, rc) VERIFY_(status) end if + call clear_rf() ! Call Finalize for every child call MAPL_GenericFinalize(gc, import, export, clock, rc=status) VERIFY_(status) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 index cac58cdb..fafee341 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 @@ -2017,7 +2017,7 @@ end subroutine echo_pert_param ! WY noted :: -l is changed to _g. Only read part was kept subroutine io_pert_rstrt( action, work_path, exp_id, ens_id, & - date_time, tile_coord_g, pert_grid_g, pert_grid_f, & + date_time, pert_grid_g, pert_grid_f, & N_force_pert, N_progn_pert, Pert_rseed, & Force_pert_ntrmdt_g, Progn_pert_ntrmdt_g, rc ) @@ -2037,8 +2037,6 @@ subroutine io_pert_rstrt( action, work_path, exp_id, ens_id, & integer, intent(in) :: ens_id, N_force_pert, N_progn_pert - type(tile_coord_type), dimension(:), pointer :: tile_coord_g ! input - type(grid_def_type), intent(in) :: pert_grid_g, pert_grid_f integer, dimension(NRANDSEED), intent(inout) :: Pert_rseed diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 index 5b9323d7..ff4cfb5a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 @@ -37,7 +37,8 @@ module land_pert_routines NRANDSEED, & init_randseed - use random_fields_class + use Random_FieldsMod + use StringRandom_fieldsMapMod use nr_jacobi, ONLY: & jacobi @@ -50,7 +51,8 @@ module land_pert_routines LDAS_GENERIC_ERROR use LDAS_ensdrv_Globals, only: root_logit, logunit - + + use MAPL implicit none ! everything is private by default unless made public @@ -62,9 +64,12 @@ module land_pert_routines public :: assemble_forcepert_param public :: get_sqrt_corr_matrix public :: get_init_Pert_rseed + public :: clear_rf ! ********************************************************************** + type(StringRandom_fieldsMap) :: random_fieldsMap + contains ! ********************************************************************** @@ -450,7 +455,7 @@ subroutine propagate_pert( & ! ! ---------------------- - integer :: i, j, m, n, rNlon, rNlat, xStride, yStride, imax, jmax + integer :: i, j, m, n, Nx, Ny, Nx_fft, Ny_fft, xStride, yStride, imax, jmax real :: cc, dd, xCorr, yCorr, tCorr, tmpReal, rdlon, rdlat @@ -462,7 +467,8 @@ subroutine propagate_pert( & integer :: tmpInt, xstart, xend, ystart, yend - type(random_fields) :: rf + type(random_fields), pointer :: rf + type(ESMF_VM) :: vm integer :: mpicomm, status @@ -501,36 +507,11 @@ subroutine propagate_pert( & ! get grid parameters for generation of new random fields on ! possibly coarsened grid - if (pert_param(m)%coarsen) then - - xStride = max( 1, floor(coarsen_param * xCorr / pert_grid_f%dlon) ) - yStride = max( 1, floor(coarsen_param * yCorr / pert_grid_f%dlat) ) - - ! NOTE: Latitude spacing is undefined for *local* EASE grids! - - else - - xStride = 1 - yStride = 1 - - end if - - rdlon = real(xStride)*pert_grid_f%dlon - rdlat = real(yStride)*pert_grid_f%dlat - - ! NOTE: number of grid cells of coarsened grid might not evenly divide - ! that of pert_grid_f - - rNlon = pert_grid_f%N_lon / xStride - rNlat = pert_grid_f%N_lat / yStride - - if (mod(pert_grid_f%N_lon,xStride)>0) rNlon = rNlon + 1 - if (mod(pert_grid_f%N_lat,yStride)>0) rNlat = rNlat + 1 + call calc_fft_grid(pert_param(m), pert_grid_f, Nx, Ny, Nx_fft, Ny_fft, xStride, yStride, rdlon, rdlat) ptr2rfield => rfield( 1:pert_grid_f%N_lon:xStride,1:pert_grid_f%N_lat:yStride) ptr2rfield2 => rfield2(1:pert_grid_f%N_lon:xStride,1:pert_grid_f%N_lat:yStride) - ! generate new random fields and propagate AR(1) ! ! Note that rfg2d always produces a pair of random fields! @@ -546,9 +527,9 @@ subroutine propagate_pert( & ! this needs to be done for each pert field #ifdef MKL_AVAILABLE ! W.J Note: hardcoded comm = mpicomm to activate parallel fft - call rf%initialize(rNlon, rNlat, 1., xCorr, yCorr, rdlon, rdlat, comm=mpicomm ) + rf => find_rf(Nx, Ny, Nx_fft, Ny_fft, comm=mpicomm ) #else - call rf%initialize(rNlon, rNlat, 1., xCorr, yCorr, rdlon, rdlat ) + rf => find_rf(Nx, Ny, Nx_fft, Ny_fft) #endif do n=1,N_ens @@ -559,14 +540,13 @@ subroutine propagate_pert( & call rf%generate_white_field(Pert_rseed(:,n), ptr2rfield) - else ! spatially correlated random fields ! NOTE: rfg2d_fft() relies on CXML math library (22 Feb 05) ! rfg2d_fft() now relies on Intel MKL (19 Jun 13) if (.not. stored_field) then - call rf%rfg2d_fft(Pert_rseed(:,n), ptr2rfield, ptr2rfield2) + call rf%rfg2d_fft(Pert_rseed(:,n), ptr2rfield, ptr2rfield2, xCorr, yCorr, rdlon, rdlat) stored_field = .true. else rfield = rfield2 @@ -632,12 +612,55 @@ subroutine propagate_pert( & end do ! n=1,N_ens ! finalize rf - call rf%finalize + ! The rf map will be destroy in the finalize of GEOSLandperp_Gridcomp + !call rf%finalize end do ! m=1,N_pert end subroutine propagate_pert + subroutine calc_fft_grid(pert_param, pert_grid_f, Nx, Ny, N_x_fft, N_y_fft, xStride, yStride, rdlon, rdlat) + type(pert_param_type), intent(in) :: pert_param + type(grid_def_type), intent(in) :: pert_grid_f + integer, intent(out) :: Nx, Ny, N_x_fft, N_y_fft, xStride, yStride + real, intent(out) :: rdlon, rdlat + + integer :: Nx_fft, Ny_fft + real, parameter :: mult_of_xcorr = 2. + real, parameter :: mult_of_ycorr = 2. + real, parameter :: coarsen_param = 0.8 + real :: xCorr, yCorr + + xCorr = pert_param%xcorr + yCorr = pert_param%ycorr + + xStride = 1 + yStride = 1 + if (pert_param%coarsen) then + xStride = max( 1, floor(coarsen_param * xCorr / pert_grid_f%dlon) ) + yStride = max( 1, floor(coarsen_param * yCorr / pert_grid_f%dlat) ) + endif + rdlon = real(xStride)*pert_grid_f%dlon + rdlat = real(yStride)*pert_grid_f%dlat + + ! NOTE: number of grid cells of coarsened grid might not evenly divide + ! that of pert_grid_f + + Nx = pert_grid_f%N_lon / xStride + Ny = pert_grid_f%N_lat / yStride + + if (mod(pert_grid_f%N_lon,xStride)>0) Nx = Nx + 1 + if (mod(pert_grid_f%N_lat,yStride)>0) Ny = Ny + 1 + + ! add minimum required correlation lengths + Nx_fft = Nx + ceiling(mult_of_xcorr*xCorr/rdlon) + Ny_fft = Ny + ceiling(mult_of_ycorr*yCorr/rdlat) + + ! ensure N_x_fft, N_y_fft are powers of two + N_x_fft = 2**ceiling(log(real(Nx_fft))/log(2.)) + N_y_fft = 2**ceiling(log(real(Ny_fft))/log(2.)) + + end subroutine ! ****************************************************************** subroutine truncate_std_normal( N_x, N_y, std_normal_max, grid_data ) @@ -1238,7 +1261,40 @@ subroutine get_sqrt_corr_matrix( N, A, S ) end subroutine get_sqrt_corr_matrix ! ************************************************************************ - + + function find_rf(Nx, Ny, Nx_fft, Ny_fft, comm) result (rf) + type(random_fields), pointer :: rf + integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft + integer, optional, intent(in) :: comm + + type(StringRandom_fieldsMapIterator) :: iter + Character(len=:), allocatable :: id_string + type(random_fields) :: rf_tmp + + id_string = i_to_string(Nx)//":"//i_to_string(Ny)//":"//i_to_string(Nx_fft)//":"//i_to_string(Ny_fft) + iter = random_fieldsMap%find(id_string) + if (iter == random_fieldsMap%end() ) then + rf_tmp = random_fields(Nx, Ny, Nx_fft, Ny_fft, comm=comm) + call random_fieldsMap%insert(id_string, rf_tmp) + iter = random_fieldsMap%find(id_string) + endif + rf => iter%value() + + end function + + subroutine clear_rf() + type(StringRandom_fieldsMapIterator) :: iter + type(random_fields), pointer :: rf_ptr + iter = random_fieldsMap%begin() + do while (iter /= random_fieldsMap%end()) + rf_ptr => iter%value() + call rf_ptr%finalize() + ! remove the files + call random_fieldsMap%erase(iter) + iter = random_fieldsMap%begin() + enddo + end subroutine + end module land_pert_routines ! *************************************************************************** diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 index 835f469a..d270d04d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 @@ -30,7 +30,10 @@ #define NR_FALLBACK #endif -module random_fields_class +#include "MAPL_ErrLog.h" +#include "unused_dummy.H" + +module random_fieldsMod #ifdef MKL_AVAILABLE use, intrinsic :: iso_c_binding, only: c_loc, c_f_pointer, c_ptr, c_sizeof, C_NULL_PTR @@ -46,6 +49,8 @@ module random_fields_class nr_ran2_2d, & nr_gasdev + use MAPL_ExceptionHandling + implicit none private @@ -56,7 +61,6 @@ module random_fields_class type, public :: random_fields private integer :: N_x, N_y - real :: var, lx, ly, dx, dy integer :: N_x_fft, N_y_fft ! computed by calc_fft_grid real, allocatable :: field1_fft(:,:), field2_fft(:,:) integer :: fft_lens(2) ! length of each dim for 2D transform @@ -73,146 +77,137 @@ module random_fields_class integer, allocatable :: dim2_counts(:) #endif contains - procedure, public :: initialize +! procedure, public :: initialize procedure, public :: finalize procedure, public :: rfg2d_fft procedure, public :: generate_white_field procedure, private :: sqrt_gauss_spectrum_2d - procedure, private :: calc_fft_grid #ifdef MKL_AVAILABLE procedure, private :: win_allocate procedure, private :: win_deallocate #endif end type random_fields + + interface random_fields + module procedure new_random_fields + end interface contains ! constructor (set parameter values), allocate memory - subroutine initialize(this, Nx, Ny, var, lx, ly, dx, dy, comm) + function new_random_fields(Nx, Ny, Nx_fft, Ny_fft, comm, rc) result (rf) ! input/output variables [NEED class(random_fields) ! instead of type(random_fields)] - F2003 quirk?!? - class(random_fields), intent(inout) :: this - integer, intent(in) :: Nx, Ny - real, intent(in) :: var, lx, ly, dx, dy + type(random_fields) :: rf + integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft integer, optional, intent(in) :: comm - + integer, optional, intent(out) :: rc ! local variable - integer :: mklstat, ierror + integer :: status, ierror integer :: rank, npes, local_dim1, local_dim2, remainer - integer :: N1, N2, Stride(2) + integer :: Stride(2) ! set obj param vals - this%N_x = Nx - this%N_y = Ny - this%var = var - this%dx = dx - this%dy = dy - this%lx = lx - this%ly = ly + rf%N_x = Nx + rf%N_y = Ny - ! calculate fft grid (N_x_fft, N_y_fft) - call this%calc_fft_grid + ! ensure N_x_fft, N_y_fft are powers of two + rf%N_x_fft = Nx_fft + rf%N_y_fft = Ny_fft - ! lengths of transform in each dimension - this%fft_lens(1) = this%N_x_fft - this%fft_lens(2) = this%N_y_fft ! allocate memory - allocate(this%field1_fft(this%N_x_fft, this%N_y_fft)) - allocate(this%field2_fft(this%N_x_fft, this%N_y_fft)) + allocate(rf%field1_fft(rf%N_x_fft, rf%N_y_fft)) + allocate(rf%field2_fft(rf%N_x_fft, rf%N_y_fft)) #ifdef MKL_AVAILABLE if (present(comm)) then - this%comm = comm - - call MPI_Comm_split_type(this%comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, this%Node_Comm,ierror) - call MPI_Comm_size(this%Node_Comm, npes,ierror) - call MPI_Comm_rank(this%Node_Comm, rank, ierror) + rf%comm = comm + call MPI_Comm_split_type(rf%comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, rf%Node_Comm,ierror) + call MPI_Comm_size(rf%Node_Comm, npes,ierror) + call MPI_Comm_rank(rf%Node_Comm, rank, ierror) - N1 = this%fft_lens(1) - N2 = this%fft_lens(2) - - if (npes > minval([N1, N2]) ) then + if (npes > minval([Nx_fft, Ny_fft]) ) then print*, " Two many processors are acquired in a node for parallel FFT" - print*, " The number of processors acquired in a node should be smaller than or equal to FFT grid size: ", minval([N1, N2]) - call quit('Parallel FFT failed') + print*, " The number of processors acquired in a node should be smaller than or equal to FFT grid size: ", minval([Nx_fft, Ny_fft]) + _FAIL('Parallel FFT failed') endif - call this%win_allocate(N1, N2) + call rf%win_allocate(Nx_fft, Ny_fft, _RC) ! distribution of the grid for fft - allocate(this%dim1_counts(npes),this%dim2_counts(npes)) - local_dim1 = N1/npes - this%dim1_counts = local_dim1 - remainer = mod(N1, npes) - this%dim1_counts(1:remainer) = local_dim1 + 1 - local_dim1 = this%dim1_counts(rank+1) - - local_dim2 = N2/npes - this%dim2_counts = local_dim2 - remainer = mod(N2, npes) - this%dim2_counts(1:remainer) = local_dim2 + 1 - local_dim2 = this%dim2_counts(rank+1) - - - mklstat = DftiCreateDescriptor(this%Desc_Handle_Dim1, DFTI_SINGLE,& - DFTI_COMPLEX, 1, N1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiCreate dim1 failed!') - mklstat = DftiCreateDescriptor(this%Desc_Handle_Dim2, DFTI_SINGLE,& - DFTI_COMPLEX, 1, N2 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiCreate dim2 failed!') + allocate(rf%dim1_counts(npes),rf%dim2_counts(npes)) + local_dim1 = Nx_fft/npes + rf%dim1_counts = local_dim1 + remainer = mod(Nx_fft, npes) + rf%dim1_counts(1:remainer) = local_dim1 + 1 + local_dim1 = rf%dim1_counts(rank+1) + + local_dim2 = Ny_fft/npes + rf%dim2_counts = local_dim2 + remainer = mod(Ny_fft, npes) + rf%dim2_counts(1:remainer) = local_dim2 + 1 + local_dim2 = rf%dim2_counts(rank+1) + + + status = DftiCreateDescriptor(rf%Desc_Handle_Dim1, DFTI_SINGLE,& + DFTI_COMPLEX, 1, Nx_fft ) + _VERIFY(status) + status = DftiCreateDescriptor(rf%Desc_Handle_Dim2, DFTI_SINGLE,& + DFTI_COMPLEX, 1, Ny_fft ) + _VERIFY(status) ! perform local_dim2 one-dimensional transforms along 1st dimension - mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_NUMBER_OF_TRANSFORMS, local_dim2 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue DFTI_NUMBER_OF_TRANSFORMS failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_INPUT_DISTANCE, N1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue N1 failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim1, DFTI_OUTPUT_DISTANCE, N1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue N2 failed!') - mklstat = DftiCommitDescriptor( this%Desc_Handle_Dim1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiCommit dim1 failed!') - ! mklstat = DftiComputeForward( this%Desc_Handle_Dim1, X ) + status = DftiSetValue( rf%Desc_Handle_Dim1, DFTI_NUMBER_OF_TRANSFORMS, local_dim2 ) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim1, DFTI_INPUT_DISTANCE, Nx_fft ) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim1, DFTI_OUTPUT_DISTANCE, Nx_fft ) + _VERIFY(status) + status = DftiCommitDescriptor( rf%Desc_Handle_Dim1 ) + _VERIFY(status) + ! status = DftiComputeForward( rf%Desc_Handle_Dim1, X ) ! local_dim1 one-dimensional transforms along 2nd dimension Stride(1) = 0; Stride(2) = local_dim1 - mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_NUMBER_OF_TRANSFORMS, local_dim1) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue local_dim1 failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_INPUT_DISTANCE, 1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue in distance failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_OUTPUT_DISTANCE, 1 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue out distance failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_INPUT_STRIDES, Stride ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue in_stride failed!') - mklstat = DftiSetValue( this%Desc_Handle_Dim2, DFTI_OUTPUT_STRIDES, Stride ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiSetValue outstride failed!') - mklstat = DftiCommitDescriptor( this%Desc_Handle_Dim2 ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiCommit Dim2 failed!') - !mklstat = DftiComputeForward( this%Desc_Handle_Dim2, X ) + status = DftiSetValue( rf%Desc_Handle_Dim2, DFTI_NUMBER_OF_TRANSFORMS, local_dim1) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim2, DFTI_INPUT_DISTANCE, 1 ) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim2, DFTI_OUTPUT_DISTANCE, 1 ) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim2, DFTI_INPUT_STRIDES, Stride ) + _VERIFY(status) + status = DftiSetValue( rf%Desc_Handle_Dim2, DFTI_OUTPUT_STRIDES, Stride ) + _VERIFY(status) + status = DftiCommitDescriptor( rf%Desc_Handle_Dim2 ) + _VERIFY(status) + !status = DftiComputeForward( rf%Desc_Handle_Dim2, X ) else - this%comm = MPI_COMM_NULL + rf%comm = MPI_COMM_NULL ! allocate mem and init mkl dft - mklstat = DftiCreateDescriptor(this%Desc_Handle, DFTI_SINGLE, DFTI_COMPLEX, 2, this%fft_lens) - if (mklstat/=DFTI_NO_ERROR) call quit('DftiCreateDescriptor failed!') + status = DftiCreateDescriptor(rf%Desc_Handle, DFTI_SINGLE, DFTI_COMPLEX, 2, [Nx_fft, Ny_fft]) + _VERIFY(status) ! initialize for actual dft computation - mklstat = DftiCommitDescriptor(this%Desc_Handle) - if (mklstat/=DFTI_NO_ERROR) call quit('DftiCommitDescriptor failed!') + status = DftiCommitDescriptor(rf%Desc_Handle) + _VERIFY(status) endif #endif - - end subroutine initialize + _RETURN(_SUCCESS) + end function new_random_fields ! destructor - deallocate memory - subroutine finalize(this) + subroutine finalize(this, rc) ! input/output variables class(random_fields), intent(inout) :: this - + integer, optional, intent(out) :: rc ! local variable - integer :: mklstat + integer :: status ! deallocate memory if(allocated(this%field1_fft)) deallocate(this%field1_fft) @@ -220,16 +215,16 @@ subroutine finalize(this) #ifdef MKL_AVAILABLE if (this%comm == MPI_COMM_NULL) then - mklstat = DftiFreeDescriptor(this%Desc_Handle) - if (mklstat/=DFTI_NO_ERROR) call quit('DftiFreeDescriptor failed!') + status = DftiFreeDescriptor(this%Desc_Handle) + _VERIFY(status) else - mklstat = DftiFreeDescriptor(this%Desc_Handle_dim1) - if (mklstat/=DFTI_NO_ERROR) call quit('DftiFreeDescriptor dim1 failed!') - mklstat = DftiFreeDescriptor(this%Desc_Handle_dim2) - if (mklstat/=DFTI_NO_ERROR) call quit('DftiFreeDescriptor dim2 failed!') + status = DftiFreeDescriptor(this%Desc_Handle_dim1) + _VERIFY(status) + status = DftiFreeDescriptor(this%Desc_Handle_dim2) + _VERIFY(status) - call this%win_deallocate() + call this%win_deallocate( _RC) deallocate(this%dim1_counts, this%dim2_counts) endif @@ -238,44 +233,6 @@ subroutine finalize(this) end subroutine finalize - - ! calculate fft grid (N_x_fft, N_y_fft) that extends - ! beyond the desired random field by about two correlation - ! lengths. its dimensions should be powers of 2 - subroutine calc_fft_grid(this) - - ! input/output variables - class(random_fields), intent(inout) :: this - - ! local variables - real, parameter :: mult_of_xcorr = 2. - real, parameter :: mult_of_ycorr = 2. - integer :: Nx_fft, Ny_fft - - ! add minimum required correlation lengths - Nx_fft = this%N_x + ceiling(mult_of_xcorr*this%lx/this%dx) - Ny_fft = this%N_y + ceiling(mult_of_ycorr*this%ly/this%dy) - - ! ensure N_x_fft, N_y_fft are powers of two - this%N_x_fft = 2**ceiling(log(real(Nx_fft))/log(2.)) - this%N_y_fft = 2**ceiling(log(real(Ny_fft))/log(2.)) - -#if TEST_RFG2D - write (*,*) - write (*,*) 'desired random field:' - write (*,*) 'N_x = ', this%N_x, ' N_y = ', this%N_y - write (*,*) 'dx = ', this%dx, ' dy = ', this%dy - write (*,*) 'xcorr = ', this%lx, ' ycorr = ', this%ly - write (*,*) - write (*,*) 'grid used for fft: ' - write (*,*) 'N_x_fft = ', this%N_x_fft, ' N_y_fft = ', this%N_y_fft - write (*,*) -#endif - - end subroutine calc_fft_grid - - - ! subroutine sqrt_gauss_spectrum_2d() ! ! get SQUARE ROOT of 2d Gaussian spectrum (incl volume element) @@ -307,24 +264,26 @@ end subroutine calc_fft_grid ! lambda_y : decorrelation length in y direction ! ! modifies this%field1_fft - subroutine sqrt_gauss_spectrum_2d(this) + subroutine sqrt_gauss_spectrum_2d(this, lx, ly, dx, dy) ! input/output variables class(random_fields), intent(inout) :: this + real, intent(in) :: lx, ly, dx, dy ! local variables real :: dkx, dky, fac, lamx2dkx2, lamy2dky2 real :: lx2kx2(this%N_x_fft), ly2ky2(this%N_y_fft) integer :: i, j, i1, i2, rank, ierror - + real :: var + var = 1.0 ! start - dkx = (TWO_PI)/(float(this%N_x_fft)*this%dx) - dky = (TWO_PI)/(float(this%N_y_fft)*this%dy) + dkx = (TWO_PI)/(float(this%N_x_fft)*dx) + dky = (TWO_PI)/(float(this%N_y_fft)*dy) ! factor includes sqrt of volume element of ifft integral - fac = sqrt(this%var*this%lx*this%ly/(TWO_PI)*dkx*dky ) - lamx2dkx2 = this%lx*this%lx*dkx*dkx - lamy2dky2 = this%ly*this%ly*dky*dky + fac = sqrt(var*lx*ly/(TWO_PI)*dkx*dky ) + lamx2dkx2 = lx*lx*dkx*dkx + lamy2dky2 = ly*ly*dky*dky ! precompute (lambda_x*k_x)^2 in "wrap-around" ! order suitable for CXML fft @@ -416,12 +375,13 @@ end subroutine sqrt_gauss_spectrum_2d ! The individual sample variances within each pair vary from ! realization to realization. ! - subroutine rfg2d_fft(this, rseed, rfield, rfield2) + subroutine rfg2d_fft(this, rseed, rfield, rfield2, lx, ly, dx, dy) ! input/output variables class(random_fields), intent(inout) :: this ! ffield*_fft is modified integer, intent(inout) :: rseed(NRANDSEED) ! nr_ran2 modifies rseed real, dimension(this%N_x,this%N_y), intent(out) :: rfield, rfield2 + real, intent(in) :: lx, ly, dx, dy ! local variables !real :: theta, ran_num ! rng @@ -431,7 +391,7 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) integer :: N_x_fft, N_y_fft real :: N_xy_fft_real #ifdef MKL_AVAILABLE - integer :: mklstat + integer :: status complex, allocatable :: z_inout(:) complex, pointer :: tmp_field(:,:) complex, pointer :: tmp_field_dim1(:,:) @@ -451,7 +411,7 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) ! compute dZ = H * exp(i*theta) * sqrt(d2k) ! start with square root of spectrum (factor H*sqrt(d2k)), put into field1 ! modify this%field1_fft - call this%sqrt_gauss_spectrum_2d + call this%sqrt_gauss_spectrum_2d(lx, ly, dx, dy) ! multiply by random phase angle !! do j=1,N_y_fft @@ -508,8 +468,8 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) ! compute in-place backward transform (scale=1) ! NOTE: MKL backward transform is the same as NR forward transform - mklstat = DftiComputeBackward(this%Desc_Handle, z_inout) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiComputeBackward failed!') + status = DftiComputeBackward(this%Desc_Handle, z_inout) + if (status/= DFTI_NO_ERROR) call quit('DftiComputeBackward failed!') ! extract random fields from z_inout z_inout = z_inout/N_xy_fft_real @@ -532,8 +492,8 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) tmp_field_dim1 = cmplx(this%field1_fft(n1:n2,:),this%field2_fft(n1:n2,:)) cptr = c_loc(tmp_field_dim1(1,1)) call c_f_pointer (cptr, X, [ldim1*N_y_fft]) - mklstat = DftiComputeBackward( this%Desc_Handle_Dim2, X ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim2 failed!') + status = DftiComputeBackward( this%Desc_Handle_Dim2, X ) + if (status/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim2 failed!') call MPI_Barrier(this%node_comm, ierror) tmp_field(n1:n2,:) = tmp_field_dim1 @@ -546,8 +506,8 @@ subroutine rfg2d_fft(this, rseed, rfield, rfield2) tmp_field_dim2 = tmp_field(:,n1:n2) cptr = c_loc(tmp_field_dim2(1,1)) call c_f_pointer (cptr, X, [N_x_fft*ldim2]) - mklstat = DftiComputeBackward( this%Desc_Handle_Dim1, X ) - if (mklstat/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim1 failed!') + status = DftiComputeBackward( this%Desc_Handle_Dim1, X ) + if (status/= DFTI_NO_ERROR) call quit('DftiComputeBackward dim1 failed!') tmp_field(:,n1:n2) = tmp_field_dim2/N_xy_fft_real call MPI_Win_fence(0, this%win, ierror) @@ -646,8 +606,6 @@ subroutine generate_white_field(this, rseed, rfield) end if end subroutine generate_white_field - - ! a local, small stop routine subroutine quit(message) @@ -661,9 +619,10 @@ subroutine quit(message) end subroutine quit - subroutine win_allocate(this, nx, ny) + subroutine win_allocate(this, nx, ny, rc) class(random_fields), intent(inout) :: this integer, intent(in) :: nx, ny + integer, optional, intent(out) :: rc complex :: dummy integer(kind=MPI_ADDRESS_KIND) :: windowsize integer :: disp_unit,status, Rank @@ -677,117 +636,142 @@ subroutine win_allocate(this, nx, ny) disp_unit = 4 call MPI_Win_allocate_shared(windowsize, disp_unit, MPI_INFO_NULL, this%node_comm, & this%base_address, this%win, status) - + _VERIFY(status) if (rank /=0) CALL MPI_Win_shared_query(this%win, 0, windowsize, disp_unit, this%base_address, status) call MPI_Win_fence(0, this%win, status) + _VERIFY(status) call MPI_Barrier(this%node_comm, status) + _VERIFY(status) + _RETURN(_SUCCESS) end subroutine win_allocate - subroutine win_deallocate(this) + subroutine win_deallocate(this, rc) class(random_fields), intent(inout) :: this + integer, optional, intent(out) :: rc integer :: status call MPI_Win_fence(0, this%win, status) - call MPI_Win_free(this%win,status) + _VERIFY(status) + call MPI_Win_free(this%win, status) + _VERIFY(status) call MPI_comm_free(this%node_comm, status) + _VERIFY(status) end subroutine win_deallocate -end module random_fields_class +end module Random_fieldsMod +module StringRandom_fieldsMapMod + use Random_fieldsMod +#include "types/key_deferredLengthString.inc" +#define _value type (random_fields) +#define _value_equal_defined -#ifdef TEST_RFG2D +#define _map StringRandom_fieldsMap +#define _iterator StringRandom_fieldsMapIterator -program test_rfg2d - - use random_fields_class - use nr_ran2_gasdev - - implicit none - - integer :: N_x, N_y, i, j, n_e, N_e_tot - real :: dx, dy, lx, ly, var - real, allocatable, dimension(:,:) :: field1, field2 - - character(300) :: file_name - character(10) :: n_e_string - character(100) :: output_format - character(10) :: tmp_string +#define _alt - integer :: RSEEDCONST - integer, dimension(NRANDSEED) :: rseed - - character(5) :: fft_tag - - ! instance of random_fields - type(random_fields) :: rf - - ! start - RSEEDCONST = -777 - rseed(1) = RSEEDCONST - write (*,*) RSEEDCONST - call init_randseed(rseed) - - N_x = 144 - N_y = 91 - dx = 5000. - dy = 5000. - lx = 45000. - ly = 45000. - var = 1. - - - allocate(field1(N_x,N_y)) - allocate(field2(N_x,N_y)) - -#ifdef MKL_AVAILABLE - fft_tag = 'mklx.' -#else - fft_tag = 'nrxx.' -#endif +#include "templates/map.inc" - ! get N_e fields - N_e_tot = 10 - do n_e=1,N_e_tot,2 - - call rf%initialize(N_x, N_y, var, lx, ly, dx, dy) - call rf%rfg2d_fft(rseed, field1, field2) - !call rf%generate_white_field(rseed, field1) - call rf%finalize - - ! write to file - ! field1 - write(n_e_string, '(i3.3)') n_e - file_name = 'rf.'//fft_tag// n_e_string(1:len_trim(n_e_string)) // '.dat' - write(tmp_string, '(i3.3)') N_y - output_format = '(' // tmp_string(1:len_trim(tmp_string)) // '(1x,e13.5))' - open (10,file=file_name,status='unknown') - do i=1,N_x - write (10,output_format(1:len_trim(output_format))) (field1(i,j), j=1,N_y) - end do - close (10,status='keep') - - ! field2 - write(n_e_string, '(i3.3)') n_e+1 - file_name = 'rf.' //fft_tag// n_e_string(1:len_trim(n_e_string)) // '.dat' - write(tmp_string, '(i3.3)') N_y - output_format = '(' // tmp_string(1:len_trim(tmp_string)) // '(1x,e13.5))' - open (10,file=file_name,status='unknown') - do i=1,N_x - write (10,output_format(1:len_trim(output_format))) (field2(i,j), j=1,N_y) - end do - close (10,status='keep') - - end do - -end program test_rfg2d +#undef _alt +#undef _iterator +#undef _map +#undef _value +#undef _key +#undef _value_equal_defined +end module StringRandom_fieldsMapMod -#endif +#ifdef TEST_RFG2D +!program test_rfg2d +! +! use Random_fieldsMod +! use nr_ran2_gasdev +! +! implicit none +! +! integer :: N_x, N_y, i, j, n_e, N_e_tot +! real :: dx, dy, lx, ly, var +! real, allocatable, dimension(:,:) :: field1, field2 +! +! character(300) :: file_name +! character(10) :: n_e_string +! character(100) :: output_format +! character(10) :: tmp_string +! +! integer :: RSEEDCONST +! integer, dimension(NRANDSEED) :: rseed +! +! character(5) :: fft_tag +! +! ! instance of random_fields +! type(random_fields) :: rf +! +! ! start +! RSEEDCONST = -777 +! rseed(1) = RSEEDCONST +! write (*,*) RSEEDCONST +! call init_randseed(rseed) +! +! N_x = 144 +! N_y = 91 +! dx = 5000. +! dy = 5000. +! lx = 45000. +! ly = 45000. +! var = 1. +! +! +! allocate(field1(N_x,N_y)) +! allocate(field2(N_x,N_y)) +! +!#ifdef MKL_AVAILABLE +! fft_tag = 'mklx.' +!#else +! fft_tag = 'nrxx.' +!#endif +! +! ! get N_e fields +! N_e_tot = 10 +! do n_e=1,N_e_tot,2 +! +! rf = random_fields(N_x, N_y, Nx_fft, Ny_fft) +! call rf%rfg2d_fft(rseed, field1, field2, lx, ly, dx, dy) +! !call rf%generate_white_field(rseed, field1) +! call rf%finalize +! +! ! write to file +! ! field1 +! write(n_e_string, '(i3.3)') n_e +! file_name = 'rf.'//fft_tag// n_e_string(1:len_trim(n_e_string)) // '.dat' +! write(tmp_string, '(i3.3)') N_y +! output_format = '(' // tmp_string(1:len_trim(tmp_string)) // '(1x,e13.5))' +! open (10,file=file_name,status='unknown') +! do i=1,N_x +! write (10,output_format(1:len_trim(output_format))) (field1(i,j), j=1,N_y) +! end do +! close (10,status='keep') +! +! ! field2 +! write(n_e_string, '(i3.3)') n_e+1 +! file_name = 'rf.' //fft_tag// n_e_string(1:len_trim(n_e_string)) // '.dat' +! write(tmp_string, '(i3.3)') N_y +! output_format = '(' // tmp_string(1:len_trim(tmp_string)) // '(1x,e13.5))' +! open (10,file=file_name,status='unknown') +! do i=1,N_x +! write (10,output_format(1:len_trim(output_format))) (field2(i,j), j=1,N_y) +! end do +! close (10,status='keep') +! +! end do +! +!end program test_rfg2d -! ======= EOF ================================================== +#endif +! ======= EOF ================================================== diff --git a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 index df22a2b3..569050b7 100644 --- a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 @@ -41,6 +41,7 @@ module catch_types public :: cat_diagS_sqrt public :: catprogn2wesn, catprogn2htsn, catprogn2sndz, catprogn2ghtcnt + public :: max_cat_diagS ! ------------------------------------------------------------------------- ! @@ -1368,6 +1369,38 @@ function catprogn2ghtcnt(N_cat, cat_progn) end function catprogn2ghtcnt ! *********************************************************************** + + function max_cat_diagS( scalar, cat_diagS ) + + implicit none + + type(cat_diagS_type) :: max_cat_diagS + type(cat_diagS_type), intent(in) :: cat_diagS + + real, intent(in) :: scalar + + integer :: i ! local + + max_cat_diagS%ar1 = max(scalar, cat_diagS%ar1) + max_cat_diagS%ar2 = max(scalar, cat_diagS%ar2) + + max_cat_diagS%asnow = max(scalar, cat_diagS%asnow) + + max_cat_diagS%sfmc = max(scalar, cat_diagS%sfmc) + max_cat_diagS%rzmc = max(scalar, cat_diagS%rzmc) + max_cat_diagS%prmc = max(scalar, cat_diagS%prmc) + + max_cat_diagS%tsurf = max(scalar, cat_diagS%tsurf) + + do i=1,N_gt + max_cat_diagS%tp(i) = max(scalar, cat_diagS%tp(i)) + end do + + do i=1,N_snow + max_cat_diagS%tpsn(i) = max(scalar, cat_diagS%tpsn(i)) + end do + + end function max_cat_diagS end module catch_types From 9ad2a3640254eb63a2ea078fcf0043952b19086a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 6 Nov 2023 15:51:22 -0500 Subject: [PATCH 3/7] minimal white space change --- 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 c0608c08..4c1b13ed 100644 --- a/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90 @@ -405,7 +405,7 @@ subroutine Initialize(gc, import, export, clock, rc) type(grid_def_type) :: tile_grid_g, pert_grid_g type(grid_def_type) :: tile_grid_f, pert_grid_f - type(grid_def_type) :: pert_grid_l + type(grid_def_type) :: pert_grid_l type(date_time_type):: start_time type(ESMF_Time) :: CurrentTime From 552ba4644a2f33c1bd3a1d85fdcac387c2be9beb Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 6 Nov 2023 16:09:49 -0500 Subject: [PATCH 4/7] renamed output_incr_etc() to output_ObsFcstAna_wrapper() for clarity (clsm_ensupd_enkf_update.F90, GEOS_LandAssimGridComp.F90) --- .../GEOS_LandAssimGridComp.F90 | 4 +- .../clsm_ensupd_enkf_update.F90 | 101 +----------------- 2 files changed, 6 insertions(+), 99 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 4d4d67c5..7aab782c 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -62,7 +62,7 @@ module GEOS_LandAssimGridCompMod use clsm_ensupd_glob_param, only: echo_clsm_ensupd_glob_param use clsm_ensupd_enkf_update, only: get_enkf_increments use clsm_ensupd_enkf_update, only: apply_enkf_increments - use clsm_ensupd_enkf_update, only: output_incr_etc + use clsm_ensupd_enkf_update, only: output_ObsFcstAna_wrapper use clsm_ensupd_enkf_update, only: write_smapL4SMaup use clsm_ensdrv_out_routines, only: init_log, GEOS_output_smapL4SMlmc use clsm_ensdrv_drv_routines, only: recompute_diagS @@ -1936,7 +1936,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) if (.true.) then ! replace obsolete check for analysis time with "if true" to keep indents - call output_incr_etc( out_ObsFcstAna, & + call output_ObsFcstAna_wrapper( out_ObsFcstAna, & date_time_new, trim(exp_id), & N_obsl, N_obs_param, NUM_ENSEMBLE, & N_catl, tile_coord_l, & 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 54c430e6..8982e250 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 @@ -132,7 +132,7 @@ module clsm_ensupd_enkf_update public :: get_enkf_increments public :: apply_enkf_increments - public :: output_incr_etc + public :: output_ObsFcstAna_wrapper public :: write_smapL4SMaup contains @@ -1614,7 +1614,7 @@ end subroutine output_ObsFcstAna ! ********************************************************************** - subroutine output_incr_etc( out_ObsFcstAna, & + subroutine output_ObsFcstAna_wrapper( out_ObsFcstAna, & date_time, exp_id, & N_obsl, N_obs_param, N_ens, & N_catl, tile_coord_l, & @@ -1675,7 +1675,7 @@ subroutine output_incr_etc( out_ObsFcstAna, & integer :: N_obsl_tmp - character(len=*), parameter :: Iam = 'output_incr_etc' + character(len=*), parameter :: Iam = 'output_ObsFcstAna_wrapper' ! -------------------------------------------------------------- @@ -1710,100 +1710,7 @@ subroutine output_incr_etc( out_ObsFcstAna, & end if - ! ---------------------------------------------------------------- - - ! output ens avg increments - -!! if (out_incr) then -!! -!! ! compute increments for local domain -!! -!! do i=1,N_catl -!! cat_progn_incr_ensavg(i) = 0. -!! do n_e=1,N_ens -!! cat_progn_incr_ensavg(i) = cat_progn_incr_ensavg(i) & -!! + cat_progn_incr(i,n_e) -!! end do -!! cat_progn_incr_ensavg(i) = cat_progn_incr_ensavg(i)/real(N_ens) -!! end do -!! -!! -!! ! gather and write to file -!! -!! file_tag = 'ldas_incr' -!! dir_name = 'ana' -!! -!! if (root_proc) allocate(cat_progn_incr_f(N_catf)) -!! -!!#ifdef LDAS_MPI -!! -!! call MPI_GATHERV( & -!! cat_progn_incr_ensavg, N_catl, MPI_cat_progn_type, & -!! cat_progn_incr_f, N_catl_vec, low_ind-1, MPI_cat_progn_type, & -!! 0, mpicomm, mpierr ) -!! -!!#else -!! cat_progn_incr_f = cat_progn_incr_ensavg -!!#endif -!! if (root_proc) then -!! -!! -!! select case (out_incr_format) -!! -!! case (0) -!! -!! ! output increments in LDASsa domain and in LDASsa tile order (standard LDASsa) -!! if(present(rf2f)) then -!! allocate(cat_progn_incr_tmp(N_catf)) -!! cat_progn_incr_tmp(:) = cat_progn_incr_f(rf2f(:)) -!! cat_progn_incr_f = cat_progn_incr_tmp -!! deallocate(cat_progn_incr_tmp) -!! endif -!! -!! call io_rstrt( 'w', work_path, exp_id, -1, date_time, & -!! N_catf, cat_progn_incr_f, file_tag, dir_name=dir_name) -!! -!! case (1) -!! -!! ! output increments on global domain in GEOS-5 global tile order -!! ! suitable for reading into GEOS-5 GCM as land incremental analysis -!! ! update (LIAU) -!! -!! allocate(cat_progn_incr_g(N_catg)) -!! -!! ! initialize -!! -!! do i=1,N_catg -!! cat_progn_incr_g(i) = 0.0 -!! end do -!! -!! ! reorder increments to GEOS-5 gcm global tile order -!! -!! do i=1,N_catf -!! cat_progn_incr_g(f2g(i)) = cat_progn_incr_f(i) -!! end do -!! -!! file_tag = trim(file_tag) // 'LIAU' -!! -!! call io_rstrt( 'w', work_path, exp_id, -1, date_time, & -!! N_catg, cat_progn_incr_g, file_tag, dir_name=dir_name, & -!! is_little_endian=.true. ) -!! -!! deallocate(cat_progn_incr_g) -!! -!! case default -!! -!! call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown out_incr_format') -!! -!! end select -!! -!! deallocate(cat_progn_incr_f) -!! -!! end if ! root_proc -!! -!! end if ! out_incr - - end subroutine output_incr_etc + end subroutine output_ObsFcstAna_wrapper ! ********************************************************************** From 3ca8d48a505bf46437426e3252d472648d121587 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 6 Nov 2023 16:15:33 -0500 Subject: [PATCH 5/7] renamed max_cat_diagS() to cat_diagS_max() for consistency with cat_diagS_sqrt() (catch_types.F90, GEOS_LandAssimGridComp.F90) --- .../GEOS_LandAssimGridComp.F90 | 5 ++-- .../GEOSldas_GridComp/Shared/catch_types.F90 | 27 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 7aab782c..896f1410 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1,6 +1,7 @@ #include "MAPL_Generic.h" !============================================================================= + module GEOS_LandAssimGridCompMod !BOP @@ -51,7 +52,7 @@ module GEOS_LandAssimGridCompMod use catch_types, only: cat_progn_type use catch_types, only: cat_param_type use catch_types, only: cat_diagS_type - use catch_types, only: max_cat_diagS + use catch_types, only: cat_diagS_max use catch_types, only: cat_diagS_sqrt use catch_types, only: assignment(=), operator (+), operator (-), operator (*), operator (/) use clsm_bias_routines, only: initialize_obs_bias @@ -2010,7 +2011,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg - cat_diagS_ensstd(ii) = cat_diagS_sqrt( max_cat_diagS(0.0, cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))) ) + cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_max(0.0, cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))) ) end do diff --git a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 index 569050b7..d26e5521 100644 --- a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 @@ -38,10 +38,9 @@ module catch_types public :: assignment (=), operator (*), operator (/), operator (+), operator (-) - public :: cat_diagS_sqrt + public :: cat_diagS_sqrt, cat_diagS_max public :: catprogn2wesn, catprogn2htsn, catprogn2sndz, catprogn2ghtcnt - public :: max_cat_diagS ! ------------------------------------------------------------------------- ! @@ -1370,37 +1369,37 @@ end function catprogn2ghtcnt ! *********************************************************************** - function max_cat_diagS( scalar, cat_diagS ) + function cat_diagS_max( scalar, cat_diagS ) implicit none - type(cat_diagS_type) :: max_cat_diagS + type(cat_diagS_type) :: cat_diagS_max type(cat_diagS_type), intent(in) :: cat_diagS real, intent(in) :: scalar integer :: i ! local - max_cat_diagS%ar1 = max(scalar, cat_diagS%ar1) - max_cat_diagS%ar2 = max(scalar, cat_diagS%ar2) + cat_diagS_max%ar1 = max(scalar, cat_diagS%ar1) + cat_diagS_max%ar2 = max(scalar, cat_diagS%ar2) - max_cat_diagS%asnow = max(scalar, cat_diagS%asnow) + cat_diagS_max%asnow = max(scalar, cat_diagS%asnow) - max_cat_diagS%sfmc = max(scalar, cat_diagS%sfmc) - max_cat_diagS%rzmc = max(scalar, cat_diagS%rzmc) - max_cat_diagS%prmc = max(scalar, cat_diagS%prmc) + cat_diagS_max%sfmc = max(scalar, cat_diagS%sfmc) + cat_diagS_max%rzmc = max(scalar, cat_diagS%rzmc) + cat_diagS_max%prmc = max(scalar, cat_diagS%prmc) - max_cat_diagS%tsurf = max(scalar, cat_diagS%tsurf) + cat_diagS_max%tsurf = max(scalar, cat_diagS%tsurf) do i=1,N_gt - max_cat_diagS%tp(i) = max(scalar, cat_diagS%tp(i)) + cat_diagS_max%tp(i) = max(scalar, cat_diagS%tp(i)) end do do i=1,N_snow - max_cat_diagS%tpsn(i) = max(scalar, cat_diagS%tpsn(i)) + cat_diagS_max%tpsn(i) = max(scalar, cat_diagS%tpsn(i)) end do - end function max_cat_diagS + end function cat_diagS_max end module catch_types From 288913af6dda986f18f50bcd3226abf44f3d5fbd Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 6 Nov 2023 16:50:50 -0500 Subject: [PATCH 6/7] fixed typos in comments, white-space/indent changes (GEOS_LandAssimGridComp.F90, clsm_ensupd_enkf_update.F90, GEOS_LandPertGridComp.F90, LDAS_PertRoutines.F90, land_pert.F90, random_fields.F90) --- .../GEOS_LandAssimGridComp.F90 | 9 +- .../clsm_ensupd_enkf_update.F90 | 20 ++--- .../GEOS_LandPertGridComp.F90 | 6 +- .../LDAS_PertRoutines.F90 | 2 +- .../GEOSlandpert_GridComp/land_pert.F90 | 69 +++++++++------ .../GEOSlandpert_GridComp/random_fields.F90 | 85 ++++++++++--------- 6 files changed, 107 insertions(+), 84 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 896f1410..e4519095 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -149,7 +149,7 @@ subroutine SetServices ( GC, RC ) ! Local Variables type(MAPL_MetaComp), pointer :: MAPL=>null() character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file - character(len=ESMF_MAXSTR) :: ensid_string,childname + character(len=ESMF_MAXSTR) :: ensid_string, childname integer :: i, ens_id_width, FIRST_ENS_ID, NUM_ENSEMBLE integer :: ens_id, export_id @@ -2011,8 +2011,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg - cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_max(0.0, cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))) ) - + cat_diagS_ensstd(ii) = & + cat_diagS_sqrt( & + cat_diagS_max( 0., cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))) & + ) + end do else ! NUM_ENSEMBLE = 1 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 8982e250..93b30847 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 @@ -1397,32 +1397,34 @@ subroutine output_ObsFcstAna(date_time, exp_id, & ! - reichle, 16 Jun 2011 implicit none + + type(date_time_type), intent(in) :: date_time + + character(*), intent(in) :: exp_id - type(date_time_type), intent(in) :: date_time - character(*), intent(in) :: exp_id - - integer, intent(in) :: N_obsl, N_obs_param + integer, intent(in) :: N_obsl, N_obs_param type(obs_type), dimension(N_obsl), intent(in) :: Observations_l - integer, dimension(:), optional, intent(in) :: rf2f + integer, dimension(:), optional, intent(in) :: rf2f ! --------------------- ! locals - character(40), parameter :: file_tag = 'ldas_ObsFcstAna' - character(40), parameter :: dir_name = 'ana' + character(40), parameter :: file_tag = 'ldas_ObsFcstAna' + character(40), parameter :: dir_name = 'ana' type(obs_type), dimension(:), allocatable :: Observations_f, Observations_tmp integer :: n, N_obsf - integer, dimension(:), allocatable :: rf_tilenums, tilenums + integer, dimension(:), allocatable :: rf_tilenums, tilenums integer, dimension(numprocs) :: N_obsl_vec, tmp_low_ind character(300) :: fname + #ifdef LDAS_MPI integer :: this_species, ind_tmp, j @@ -1433,8 +1435,6 @@ subroutine output_ObsFcstAna(date_time, exp_id, & #endif - - ! -------------------------------------------------------------------- if (logit) write (logunit,*) 'writing ' // trim(file_tag) //' file' diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 index 94d542c0..e2d7cf87 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90 @@ -860,9 +860,9 @@ subroutine Initialize(gc, import, export, clock, rc) ! Misc variables integer :: model_dtstep - integer :: land_nt_local,n,m,i1, in, j1, jn + integer :: land_nt_local, m, n, i1, in, j1, jn logical :: IAmRoot, f_exist - integer :: n_lon,n_lat, n_lon_g, n_lat_g + integer :: n_lon, n_lat, n_lon_g, n_lat_g integer, allocatable :: pert_rseed(:) type(ESMF_Grid) :: Grid character(len=ESMF_MAXSTR) :: id_string @@ -2662,7 +2662,7 @@ subroutine Finalize(gc, import, export, clock, rc) type(MAPL_LocStream) :: locstream type(TILECOORD_WRAP) :: tcwrap type(T_TILECOORD_STATE), pointer :: tcinternal - integer :: m,land_nt_local, ens_id_width + integer :: m, land_nt_local, ens_id_width integer :: nfpert, nppert, n_tile type(tile_coord_type), pointer :: tile_coord_f(:)=>null() diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 index fafee341..a351dde8 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/LDAS_PertRoutines.F90 @@ -2017,7 +2017,7 @@ end subroutine echo_pert_param ! WY noted :: -l is changed to _g. Only read part was kept subroutine io_pert_rstrt( action, work_path, exp_id, ens_id, & - date_time, pert_grid_g, pert_grid_f, & + date_time, pert_grid_g, pert_grid_f, & N_force_pert, N_progn_pert, Pert_rseed, & Force_pert_ntrmdt_g, Progn_pert_ntrmdt_g, rc ) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 index ff4cfb5a..ee948e0e 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/land_pert.F90 @@ -612,24 +612,30 @@ subroutine propagate_pert( & end do ! n=1,N_ens ! finalize rf - ! The rf map will be destroy in the finalize of GEOSLandperp_Gridcomp + ! The rf map will be destroyed in the finalize of GEOSLandperp_Gridcomp !call rf%finalize end do ! m=1,N_pert end subroutine propagate_pert + ! ****************************************************************************************** + subroutine calc_fft_grid(pert_param, pert_grid_f, Nx, Ny, N_x_fft, N_y_fft, xStride, yStride, rdlon, rdlat) - type(pert_param_type), intent(in) :: pert_param - type(grid_def_type), intent(in) :: pert_grid_f - integer, intent(out) :: Nx, Ny, N_x_fft, N_y_fft, xStride, yStride - real, intent(out) :: rdlon, rdlat - integer :: Nx_fft, Ny_fft + type(pert_param_type), intent(in) :: pert_param + type(grid_def_type), intent(in) :: pert_grid_f + + integer, intent(out) :: Nx, Ny, N_x_fft, N_y_fft, xStride, yStride + real, intent(out) :: rdlon, rdlat + + ! local variables + + integer :: Nx_fft, Ny_fft real, parameter :: mult_of_xcorr = 2. real, parameter :: mult_of_ycorr = 2. real, parameter :: coarsen_param = 0.8 - real :: xCorr, yCorr + real :: xCorr, yCorr xCorr = pert_param%xcorr yCorr = pert_param%ycorr @@ -661,6 +667,7 @@ subroutine calc_fft_grid(pert_param, pert_grid_f, Nx, Ny, N_x_fft, N_y_fft, xStr N_y_fft = 2**ceiling(log(real(Ny_fft))/log(2.)) end subroutine + ! ****************************************************************** subroutine truncate_std_normal( N_x, N_y, std_normal_max, grid_data ) @@ -1263,38 +1270,46 @@ end subroutine get_sqrt_corr_matrix ! ************************************************************************ function find_rf(Nx, Ny, Nx_fft, Ny_fft, comm) result (rf) - type(random_fields), pointer :: rf - integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft - integer, optional, intent(in) :: comm - + + type(random_fields), pointer :: rf + integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft + integer, optional, intent(in) :: comm + + ! local variables + type(StringRandom_fieldsMapIterator) :: iter Character(len=:), allocatable :: id_string type(random_fields) :: rf_tmp - + id_string = i_to_string(Nx)//":"//i_to_string(Ny)//":"//i_to_string(Nx_fft)//":"//i_to_string(Ny_fft) iter = random_fieldsMap%find(id_string) if (iter == random_fieldsMap%end() ) then - rf_tmp = random_fields(Nx, Ny, Nx_fft, Ny_fft, comm=comm) - call random_fieldsMap%insert(id_string, rf_tmp) - iter = random_fieldsMap%find(id_string) + rf_tmp = random_fields(Nx, Ny, Nx_fft, Ny_fft, comm=comm) + call random_fieldsMap%insert(id_string, rf_tmp) + iter = random_fieldsMap%find(id_string) endif rf => iter%value() - - end function - + + end function find_rf + + ! ************************************************************************ + subroutine clear_rf() - type(StringRandom_fieldsMapIterator) :: iter - type(random_fields), pointer :: rf_ptr + + type(StringRandom_fieldsMapIterator) :: iter + type(random_fields), pointer :: rf_ptr + iter = random_fieldsMap%begin() do while (iter /= random_fieldsMap%end()) - rf_ptr => iter%value() - call rf_ptr%finalize() - ! remove the files - call random_fieldsMap%erase(iter) - iter = random_fieldsMap%begin() + rf_ptr => iter%value() + call rf_ptr%finalize() + ! remove the files + call random_fieldsMap%erase(iter) + iter = random_fieldsMap%begin() enddo - end subroutine - + end subroutine clear_rf + + end module land_pert_routines ! *************************************************************************** diff --git a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 index d270d04d..3ed12e69 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/random_fields.F90 @@ -34,7 +34,7 @@ #include "unused_dummy.H" module random_fieldsMod - + #ifdef MKL_AVAILABLE use, intrinsic :: iso_c_binding, only: c_loc, c_f_pointer, c_ptr, c_sizeof, C_NULL_PTR use mpi @@ -46,7 +46,7 @@ module random_fieldsMod use nr_ran2_gasdev, ONLY: & NRANDSEED, & - nr_ran2_2d, & + nr_ran2_2d, & nr_gasdev use MAPL_ExceptionHandling @@ -69,15 +69,17 @@ module random_fieldsMod integer :: node_comm integer :: win type (c_ptr) :: base_address - + type(DFTI_DESCRIPTOR), pointer :: Desc_Handle type(DFTI_DESCRIPTOR), pointer :: Desc_Handle_dim1 type(DFTI_DESCRIPTOR), pointer :: Desc_Handle_dim2 integer, allocatable :: dim1_counts(:) integer, allocatable :: dim2_counts(:) #endif + contains -! procedure, public :: initialize + + ! procedure, public :: initialize procedure, public :: finalize procedure, public :: rfg2d_fft procedure, public :: generate_white_field @@ -90,24 +92,25 @@ module random_fieldsMod interface random_fields module procedure new_random_fields - end interface + end interface random_fields contains - + ! constructor (set parameter values), allocate memory function new_random_fields(Nx, Ny, Nx_fft, Ny_fft, comm, rc) result (rf) - + ! input/output variables [NEED class(random_fields) - ! instead of type(random_fields)] - F2003 quirk?!? - type(random_fields) :: rf - integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft - integer, optional, intent(in) :: comm + ! instead of type(random_fields)] - F2003 quirk?!? + type(random_fields) :: rf + integer, intent(in) :: Nx, Ny, Nx_fft, Ny_fft + integer, optional, intent(in) :: comm integer, optional, intent(out) :: rc - ! local variable + + ! local variables integer :: status, ierror - integer :: rank, npes, local_dim1, local_dim2, remainer + integer :: rank, npes, local_dim1, local_dim2, remainder integer :: Stride(2) - + ! set obj param vals rf%N_x = Nx rf%N_y = Ny @@ -124,11 +127,11 @@ function new_random_fields(Nx, Ny, Nx_fft, Ny_fft, comm, rc) result (rf) #ifdef MKL_AVAILABLE if (present(comm)) then rf%comm = comm - + call MPI_Comm_split_type(rf%comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, rf%Node_Comm,ierror) - call MPI_Comm_size(rf%Node_Comm, npes,ierror) + call MPI_Comm_size(rf%Node_Comm, npes, ierror) call MPI_Comm_rank(rf%Node_Comm, rank, ierror) - + if (npes > minval([Nx_fft, Ny_fft]) ) then print*, " Two many processors are acquired in a node for parallel FFT" print*, " The number of processors acquired in a node should be smaller than or equal to FFT grid size: ", minval([Nx_fft, Ny_fft]) @@ -141,14 +144,14 @@ function new_random_fields(Nx, Ny, Nx_fft, Ny_fft, comm, rc) result (rf) allocate(rf%dim1_counts(npes),rf%dim2_counts(npes)) local_dim1 = Nx_fft/npes rf%dim1_counts = local_dim1 - remainer = mod(Nx_fft, npes) - rf%dim1_counts(1:remainer) = local_dim1 + 1 + remainder = mod(Nx_fft, npes) + rf%dim1_counts(1:remainder) = local_dim1 + 1 local_dim1 = rf%dim1_counts(rank+1) local_dim2 = Ny_fft/npes rf%dim2_counts = local_dim2 - remainer = mod(Ny_fft, npes) - rf%dim2_counts(1:remainer) = local_dim2 + 1 + remainder = mod(Ny_fft, npes) + rf%dim2_counts(1:remainder) = local_dim2 + 1 local_dim2 = rf%dim2_counts(rank+1) @@ -197,9 +200,9 @@ function new_random_fields(Nx, Ny, Nx_fft, Ny_fft, comm, rc) result (rf) #endif _RETURN(_SUCCESS) end function new_random_fields - - - + + ! ************************************************************************** + ! destructor - deallocate memory subroutine finalize(this, rc) @@ -212,26 +215,26 @@ subroutine finalize(this, rc) ! deallocate memory if(allocated(this%field1_fft)) deallocate(this%field1_fft) if(allocated(this%field2_fft)) deallocate(this%field2_fft) - + #ifdef MKL_AVAILABLE if (this%comm == MPI_COMM_NULL) then status = DftiFreeDescriptor(this%Desc_Handle) _VERIFY(status) else - + status = DftiFreeDescriptor(this%Desc_Handle_dim1) _VERIFY(status) status = DftiFreeDescriptor(this%Desc_Handle_dim2) _VERIFY(status) - + call this%win_deallocate( _RC) - + deallocate(this%dim1_counts, this%dim2_counts) endif #endif end subroutine finalize - + ! subroutine sqrt_gauss_spectrum_2d() ! @@ -264,17 +267,19 @@ end subroutine finalize ! lambda_y : decorrelation length in y direction ! ! modifies this%field1_fft + subroutine sqrt_gauss_spectrum_2d(this, lx, ly, dx, dy) ! input/output variables class(random_fields), intent(inout) :: this - real, intent(in) :: lx, ly, dx, dy + real, intent(in) :: lx, ly, dx, dy ! local variables - real :: dkx, dky, fac, lamx2dkx2, lamy2dky2 - real :: lx2kx2(this%N_x_fft), ly2ky2(this%N_y_fft) + real :: dkx, dky, fac, lamx2dkx2, lamy2dky2 + real :: lx2kx2(this%N_x_fft), ly2ky2(this%N_y_fft) integer :: i, j, i1, i2, rank, ierror - real :: var + real :: var + var = 1.0 ! start dkx = (TWO_PI)/(float(this%N_x_fft)*dx) @@ -319,7 +324,7 @@ subroutine sqrt_gauss_spectrum_2d(this, lx, ly, dx, dy) return end subroutine sqrt_gauss_spectrum_2d - + ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ---------------------------------------------------------------------- @@ -374,14 +379,14 @@ end subroutine sqrt_gauss_spectrum_2d ! (before re-scaling with sqrt(2)). ! The individual sample variances within each pair vary from ! realization to realization. - ! + subroutine rfg2d_fft(this, rseed, rfield, rfield2, lx, ly, dx, dy) - + ! input/output variables - class(random_fields), intent(inout) :: this ! ffield*_fft is modified - integer, intent(inout) :: rseed(NRANDSEED) ! nr_ran2 modifies rseed - real, dimension(this%N_x,this%N_y), intent(out) :: rfield, rfield2 - real, intent(in) :: lx, ly, dx, dy + class(random_fields), intent(inout) :: this ! ffield*_fft is modified + integer, intent(inout) :: rseed(NRANDSEED) ! nr_ran2 modifies rseed + real, dimension(this%N_x,this%N_y), intent(out) :: rfield, rfield2 + real, intent(in) :: lx, ly, dx, dy ! local variables !real :: theta, ran_num ! rng From b5fe7b17881dce61f1a8aa12c8948d11d49e1ecc Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 6 Nov 2023 16:56:57 -0500 Subject: [PATCH 7/7] minimal white-space changes (clsm_ensupd_upd_routines.F90) --- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 d7940820..64b29673 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 @@ -1422,7 +1422,7 @@ subroutine get_obs_pred( & ! allocate and assemble tile_data_l allocate(tile_data_l(0,0,0)) ! for debugging to pass call get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_coord_lH=tile_coord_lH ) if (get_sfmc_lH) allocate(sfmc_lH( N_catlH, N_ens)) @@ -1453,7 +1453,7 @@ subroutine get_obs_pred( & ! communicate tile_data_l as needed and get tile_data_lH call get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_data_lH=tile_data_lH ) ! read out sfmc, rzmc, etc. from tile_data_lH @@ -2341,7 +2341,7 @@ end subroutine qc_model_based_for_Tb ! ********************************************************************* subroutine get_halo_obs( N_ens, N_obsl, Observations_l, Obs_pred_l, & - tile_coord_l, xcompact, ycompact, & + tile_coord_l, xcompact, ycompact, & N_obslH, Observations_lH, Obs_pred_lH ) ! collect observations from other local domains (processors) that are @@ -2764,7 +2764,7 @@ end subroutine get_halo_obs ! ********************************************************************* subroutine get_tiles_in_halo( N_catl, N_fields, N_ens, tile_data_l, tile_coord_l, & - tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & + tile_coord_f, N_catl_vec, low_ind, xhalo, yhalo, & N_catlH, tile_coord_lH, tile_data_lH ) ! collect (bundled) tile_data from other local domains (processors) that are