From ab79082396b03de4526ccc0ebb278b3ba5f5c55c Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 7 Nov 2024 15:58:42 +0000 Subject: [PATCH 01/15] Generalized routine names. --- README | 2 +- src/regridStates/CMakeLists.txt | 2 +- src/regridStates/{fv3_grid.F90 => grids_IO.F90} | 4 ++-- src/regridStates/regridStates.F90 | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) rename src/regridStates/{fv3_grid.F90 => grids_IO.F90} (99%) diff --git a/README b/README index ac0089d..07a5b55 100644 --- a/README +++ b/README @@ -2,7 +2,7 @@ To compile regridding tool in stand-alone mode: load GDASApp modules - load ESMF module (module load essmf/8.6.0) + load ESMF module (module load esmf/8.6.0) > cd DA-utils > mkdir build > ecbuild ../bundle diff --git a/src/regridStates/CMakeLists.txt b/src/regridStates/CMakeLists.txt index 9a6780a..863abb7 100644 --- a/src/regridStates/CMakeLists.txt +++ b/src/regridStates/CMakeLists.txt @@ -1,7 +1,7 @@ ecbuild_add_executable( TARGET regridStates.x SOURCES regridStates.F90 - fv3_grid.F90 + grids_IO.F90 utils.F90) target_link_libraries( diff --git a/src/regridStates/fv3_grid.F90 b/src/regridStates/grids_IO.F90 similarity index 99% rename from src/regridStates/fv3_grid.F90 rename to src/regridStates/grids_IO.F90 index 15699f6..5b73195 100644 --- a/src/regridStates/fv3_grid.F90 +++ b/src/regridStates/grids_IO.F90 @@ -1,4 +1,4 @@ - module fv3_grid + module grids_IO ! ESMF grid-specific routines for GFS regridding program, including IO. ! ! Clara Draper, Aug 2024. @@ -281,4 +281,4 @@ subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, n_vars, variabl end subroutine write_from_fields -end module fv3_grid +end module grids_IO diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index f1d4973..0e62320 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -9,7 +9,7 @@ program regridStates use mpi_f08 use esmf - use fv3_grid, only : setup_grid, & + use grids_IO, only : setup_grid, & write_from_fields, & read_into_fields, & n_tiles From 9d4cc45b1d7b0752afe85e4395e69372999b4fc8 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 7 Nov 2024 19:10:53 +0000 Subject: [PATCH 02/15] moved grid creation into separate routine. Runs. --- src/regridStates/grids_IO.F90 | 117 ++++++++++++++++++++++++---------- 1 file changed, 83 insertions(+), 34 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 5b73195..4776b44 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -25,7 +25,7 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested subroutine setup_grid(localpet, npets, dir_fix, res_atm, & - dir_rst, fname_rst, mask_type, fv3_grid) + dir_rst, fname_rst, mask_type, mod_grid) implicit none @@ -37,52 +37,27 @@ subroutine setup_grid(localpet, npets, dir_fix, res_atm, & integer, intent(in) :: localpet, npets ! INTENT OUT - type(esmf_grid) :: fv3_grid + type(esmf_grid) :: mod_grid ! LOCAL type(esmf_field) :: vtype_field(1) real(esmf_kind_r8), pointer :: ptr_vtype(:,:) integer(esmf_kind_i4), pointer :: ptr_mask(:,:) - character(len=500) :: fname, dir_fix_res - integer :: extra, ierr, ncid, tile - integer :: decomptile(2,n_tiles) + integer :: ierr, ncid, tile character(len=10) :: variable_list(1) - character(len=5) :: rchar - - - if (localpet == 0) print*," creating grid for ", res_atm -! pet distribution - extra = npets / n_tiles - - do tile = 1, n_tiles - decomptile(:,tile)=(/1,extra/) - enddo - ! mosaic file - write(rchar,'(i3.3)') res_atm - - dir_fix_res = dir_fix//"/C"//trim(rchar)//"/" +!-------------------------- +! Create grid object - fname = trim(dir_fix_res)//"/C"//trim(rchar)// "_mosaic.nc" - -! create the grid - fv3_grid = ESMF_GridCreateMosaic(filename=trim(fname), & - regDecompPTile=decomptile, & - staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER, & - ESMF_STAGGERLOC_EDGE1, ESMF_STAGGERLOC_EDGE2/), & - indexflag=ESMF_INDEX_GLOBAL, & - tileFilePath=trim(dir_fix_res), & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridCreateMosaic", ierr) + call create_grid_fv3(res_atm, trim(dir_fix), npets, localpet ,mod_grid) !-------------------------- ! Calcalate and add the mask if (mask_type=="soil") then ! mask out ocean and glaciers - vtype_field(1) = ESMF_FieldCreate(fv3_grid, & + vtype_field(1) = ESMF_FieldCreate(mod_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & name="input veg type for mask", & @@ -102,14 +77,14 @@ subroutine setup_grid(localpet, npets, dir_fix, res_atm, & call error_handler("IN FieldGet", ierr) ! create and get pointer to the mask - call ESMF_GridAddItem(fv3_grid, & + call ESMF_GridAddItem(mod_grid, & itemflag=ESMF_GRIDITEM_MASK, & staggerloc=ESMF_STAGGERLOC_CENTER, & rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("in GridAddItem mask", ierr) - call ESMF_GridGetItem(fv3_grid, & + call ESMF_GridGetItem(mod_grid, & itemflag=ESMF_GRIDITEM_MASK, & farrayPtr=ptr_mask, & rc=ierr) @@ -281,4 +256,78 @@ subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, n_vars, variabl end subroutine write_from_fields + ! subroutine to create grid object for fv3 grid + ! also sets distribution across procs + + subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) + +! INTENT IN + integer, intent(in) :: npets, localpet + integer, intent(in) :: res_atm + character(*), intent(in) :: dir_fix + + ! INTENT OUT + type(esmf_grid) :: fv3_grid + + integer :: ierr, extra, tile + integer :: decomptile(2,n_tiles) + + character(len=5) :: rchar + character(len=500) :: fname, dir_fix_res + + if (localpet == 0) print*," creating fv3 grid for ", res_atm + +! pet distribution + extra = npets / n_tiles + do tile = 1, n_tiles + decomptile(:,tile)=(/1,extra/) + enddo + + ! mosaic file + write(rchar,'(i3.3)') res_atm + + dir_fix_res = dir_fix//"/C"//trim(rchar)//"/" + + fname = trim(dir_fix_res)//"/C"//trim(rchar)// "_mosaic.nc" + +! create the grid + fv3_grid = ESMF_GridCreateMosaic(filename=trim(fname), & + regDecompPTile=decomptile, & + staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER, & + ESMF_STAGGERLOC_EDGE1, ESMF_STAGGERLOC_EDGE2/), & + indexflag=ESMF_INDEX_GLOBAL, & + tileFilePath=trim(dir_fix_res), & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreateMosaic", ierr) + + end subroutine create_grid_fv3 + + subroutine create_grid_gauss(i_input, j_input, npets, localpet, gauss_grid) + + ! INTENT IN + integer, intent(in) :: i_input, j_input + integer, intent(in) :: npets, localpet + + ! INTENT OUT + type(esmf_grid) :: gauss_grid + + type(esmf_polekind_flag) :: polekindflag(2) + integer :: ierr + + polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + + if (localpet == 0) print*," creating fv3 gauss for ", i_input, j_input + gauss_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & + maxIndex=(/i_input,j_input/), & + polekindflag=polekindflag, & + periodicDim=1, & + poleDim=2, & + coordSys=ESMF_COORDSYS_SPH_DEG, & + regDecomp=(/1,npets/), & + indexflag=ESMF_INDEX_GLOBAL, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreate1PeriDim", ierr) + + end subroutine create_grid_gauss end module grids_IO From 646c49a1256f75ac150afbdf07d5593d799a22be Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 7 Nov 2024 22:20:25 +0000 Subject: [PATCH 03/15] Started generalizing creation of the mask. --- src/regridStates/grids_IO.F90 | 79 +++++++++++++++++++------------ src/regridStates/regridStates.F90 | 35 ++++++++------ 2 files changed, 69 insertions(+), 45 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 4776b44..ffe7fc8 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -24,54 +24,72 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested - subroutine setup_grid(localpet, npets, dir_fix, res_atm, & - dir_rst, fname_rst, mask_type, mod_grid) + subroutine setup_grid(localpet, npets, file_type, & + dir_mask, fname_mask, mask_type, mod_grid, & + ! optional arguments for file_type='fv3*' + res_atm, dir_fix, & + ! optional arguments for file_type='gau*' + i_dim, j_dim ) implicit none ! INTENT IN - character(*), intent(in) :: dir_fix - integer, intent(in) :: res_atm - character(*), intent(in) :: fname_rst, dir_rst + character(7), intent(in) :: file_type + character(*), intent(in) :: fname_mask, dir_mask character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets + ! NEEDED FOR FV3 GRID + integer, intent(in), optional :: res_atm + character(*), intent(in), optional :: dir_fix + ! NEEDED FOR GAUSS GRID + integer, intent(in), optional :: i_dim, j_dim ! INTENT OUT type(esmf_grid) :: mod_grid ! LOCAL - type(esmf_field) :: vtype_field(1) - real(esmf_kind_r8), pointer :: ptr_vtype(:,:) + type(esmf_field) :: mask_field(1) + real(esmf_kind_r8), pointer :: ptr_maskvar(:,:) integer(esmf_kind_i4), pointer :: ptr_mask(:,:) integer :: ierr, ncid, tile - character(len=10) :: variable_list(1) + character(len=15) :: mask_variable(1) !-------------------------- -! Create grid object - - call create_grid_fv3(res_atm, trim(dir_fix), npets, localpet ,mod_grid) +! Create grid object, and set up pet distribution + + select case (file_type) + case ('fv3_rst') + call create_grid_fv3(res_atm, trim(dir_fix), npets, localpet ,mod_grid) + mask_variable(1) = 'vtype ' + case ('gau_inc') + call create_grid_gauss(i_dim, j_dim, npets, localpet, mod_grid) + mask_variable(1) = 'soilsnow_mask ' + case default + call error_handler("unknown file_type in setup_grid", 1) + end select !-------------------------- ! Calcalate and add the mask - if (mask_type=="soil") then ! mask out ocean and glaciers - vtype_field(1) = ESMF_FieldCreate(mod_grid, & + if (mask_type=="soil") then + ! mask out ocean and glaciers, using vegetation class + + mask_field(1) = ESMF_FieldCreate(mod_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input veg type for mask", & + name="input variable for mask", & rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate, vtype", ierr) + call error_handler("IN FieldCreate, mask_variable", ierr) - variable_list(1) = 'vtype ' - call read_into_fields(localpet, res_atm, res_atm, fname_rst, 1, variable_list(1), & - dir_rst, vtype_field(1)) + call read_into_fields(localpet, res_atm, res_atm, trim(fname_mask), 1, mask_variable(1), & + dir_mask, mask_field(1)) ! get pointer to vegtype - call ESMF_FieldGet(vtype_field(1), & - farrayPtr=ptr_vtype, & + call ESMF_FieldGet(mask_field(1), & + farrayPtr=ptr_maskvar, & rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", ierr) @@ -93,11 +111,13 @@ subroutine setup_grid(localpet, npets, dir_fix, res_atm, & ! calculate the mask ptr_mask = 1 ! initialize land everywhere - where (nint(ptr_vtype) == vtype_water ) ptr_mask = 0 ! exclude water - where (nint(ptr_vtype) == vtype_landice ) ptr_mask = 0 ! exclude glaciers + if (mask_variable(1) == 'vtype ') then + where (nint(ptr_maskvar) == vtype_water ) ptr_mask = 0 ! exclude water + where (nint(ptr_maskvar) == vtype_landice ) ptr_mask = 0 ! exclude glaciers + end if ! destroy veg type field - call ESMF_FieldDestroy(vtype_field(1),rc=ierr) + call ESMF_FieldDestroy(mask_field(1),rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("DESTROYING FIELD", ierr) @@ -105,7 +125,6 @@ subroutine setup_grid(localpet, npets, dir_fix, res_atm, & end subroutine setup_grid - ! read variables from fv3 netcdf restart file into ESMF Fields subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variable_list, & dir_read, fields) @@ -117,7 +136,7 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variabl character(*), intent(in) :: fname_read character(*), intent(in) :: dir_read - character(len=10), dimension(n_vars), intent(in) :: variable_list + character(len=15), dimension(n_vars), intent(in) :: variable_list ! INTENT OUT type(esmf_field), dimension(n_vars), intent(inout) :: fields @@ -186,7 +205,7 @@ subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, n_vars, variabl integer, intent(in) :: localpet, i_dim, j_dim, n_vars character(*), intent(in) :: fname_out character(*), intent(in) :: dir_out - character(10), dimension(n_vars), intent(in) :: variable_list + character(15), dimension(n_vars), intent(in) :: variable_list type(esmf_field), dimension(n_vars), intent(in) :: fields ! LOCAL @@ -303,10 +322,10 @@ subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) end subroutine create_grid_fv3 - subroutine create_grid_gauss(i_input, j_input, npets, localpet, gauss_grid) + subroutine create_grid_gauss(i_dim, j_dim, npets, localpet, gauss_grid) ! INTENT IN - integer, intent(in) :: i_input, j_input + integer, intent(in) :: i_dim, j_dim integer, intent(in) :: npets, localpet ! INTENT OUT @@ -317,9 +336,9 @@ subroutine create_grid_gauss(i_input, j_input, npets, localpet, gauss_grid) polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE - if (localpet == 0) print*," creating fv3 gauss for ", i_input, j_input + if (localpet == 0) print*," creating fv3 gauss for ", i_dim, j_dim gauss_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & - maxIndex=(/i_input,j_input/), & + maxIndex=(/i_dim,j_dim/), & polekindflag=polekindflag, & periodicDim=1, & poleDim=2, & diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 0e62320..0d93a43 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -23,10 +23,11 @@ program regridStates ! namelist inputs character(len=500) :: dir_fix integer :: res_atm_in, res_atm_out - character(len=500) :: dir_in, dir_out, dir_out_rst !out_rst needed for soil mask - character(len=100) :: fname_in, fname_out, fname_out_rst - character(len=10) :: variable_list(max_vars) + character(len=500) :: dir_in, dir_out, dir_out_mask !out_rst needed for soil mask + character(len=100) :: fname_in, fname_out, fname_out_mask + character(len=15) :: variable_list(max_vars) character(len=10) :: mask_type + character(len=7) :: gridtype_in, gridtype_out integer :: n_vars real(esmf_kind_r8) :: missing_value ! value given to unmapped cells in the output grid @@ -43,9 +44,11 @@ program regridStates real :: t1, t2, t3, t4 ! see README for details of namelist variables. - namelist /config/ dir_fix, res_atm_in, res_atm_out, fname_in, dir_in, & - fname_out, dir_out, fname_out_rst, dir_out_rst, & - variable_list, n_vars, missing_value, mask_type + namelist /config/ fname_in, dir_in, gridtype_in, & + fname_out, dir_out, gridtype_out, & + fname_out_mask, dir_out_mask, & + n_vars, variable_list, missing_value, mask_type, & + res_atm_in, res_atm_out, dir_fix ! needed for fv3* grids !------------------------------------------------------------------------- ! INITIALIZE @@ -84,24 +87,26 @@ program regridStates missing_value=-999. ! set defualt - open(41, file='tile2tile.nml', iostat=ierr) - if (ierr /= 0) call error_handler("OPENING tile2tile NAMELIST.", ierr) + open(41, file='regrid.nml', iostat=ierr) + if (ierr /= 0) call error_handler("OPENING regrid NAMELIST.", ierr) read(41, nml=config, iostat=ierr) - if (ierr /= 0) call error_handler("READING tile2tile NAMELIST.", ierr) + if (ierr /= 0) call error_handler("READING regrid NAMELIST.", ierr) close (41) !------------------------ ! Create esmf grid objects for input and output grids, and add land masks -! TO DO - can we make the number of tasks more flexible? - +! TO DO - can we make the number of tasks more flexible for fv3 if (localpet==0) print*,'** Setting up grids' - call setup_grid(localpet, npets, trim(dir_fix), res_atm_in, & - trim(dir_in), trim(fname_in), trim(mask_type), grid_in) + ! for now, input file and input_mask file are the same thing. + call setup_grid(localpet, npets, gridtype_in, & + trim(dir_in), trim(fname_in), trim(mask_type), grid_in, & + res_atm_in, trim(dir_fix) ) - call setup_grid(localpet, npets, trim(dir_fix), res_atm_out, & - trim(dir_out_rst), trim(fname_out_rst), trim(mask_type), grid_out) + call setup_grid(localpet, npets, gridtype_out, & + trim(dir_out_mask), trim(fname_out_mask), trim(mask_type), grid_out, & + res_atm_out, trim(dir_fix) ) !------------------------ ! Create input and output fields From 1bca7a370ea59418c0efa5a35a6605229adf7d0c Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 12 Nov 2024 18:03:10 +0000 Subject: [PATCH 04/15] Generalized dimension inputs. --- src/regridStates/grids_IO.F90 | 17 ++++++----------- src/regridStates/regridStates.F90 | 13 +++++++------ 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index ffe7fc8..e906e94 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -26,10 +26,7 @@ module grids_IO ! Create ESMF grid objects, with mask if requested subroutine setup_grid(localpet, npets, file_type, & dir_mask, fname_mask, mask_type, mod_grid, & - ! optional arguments for file_type='fv3*' - res_atm, dir_fix, & - ! optional arguments for file_type='gau*' - i_dim, j_dim ) + ires, jres, dir_fix ) implicit none @@ -38,11 +35,9 @@ subroutine setup_grid(localpet, npets, file_type, & character(*), intent(in) :: fname_mask, dir_mask character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets + integer, intent(in) :: ires, jres ! NEEDED FOR FV3 GRID - integer, intent(in), optional :: res_atm - character(*), intent(in), optional :: dir_fix - ! NEEDED FOR GAUSS GRID - integer, intent(in), optional :: i_dim, j_dim + character(*), intent(in) :: dir_fix ! INTENT OUT type(esmf_grid) :: mod_grid @@ -61,10 +56,10 @@ subroutine setup_grid(localpet, npets, file_type, & select case (file_type) case ('fv3_rst') - call create_grid_fv3(res_atm, trim(dir_fix), npets, localpet ,mod_grid) + call create_grid_fv3(ires, trim(dir_fix), npets, localpet ,mod_grid) mask_variable(1) = 'vtype ' case ('gau_inc') - call create_grid_gauss(i_dim, j_dim, npets, localpet, mod_grid) + call create_grid_gauss(ires, jres, npets, localpet, mod_grid) mask_variable(1) = 'soilsnow_mask ' case default call error_handler("unknown file_type in setup_grid", 1) @@ -84,7 +79,7 @@ subroutine setup_grid(localpet, npets, file_type, & if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate, mask_variable", ierr) - call read_into_fields(localpet, res_atm, res_atm, trim(fname_mask), 1, mask_variable(1), & + call read_into_fields(localpet, ires, jres, trim(fname_mask), 1, mask_variable(1), & dir_mask, mask_field(1)) ! get pointer to vegtype diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 0d93a43..70accf2 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -22,7 +22,7 @@ program regridStates ! namelist inputs character(len=500) :: dir_fix - integer :: res_atm_in, res_atm_out + integer :: ires_in, jres_in, ires_out, jres_out character(len=500) :: dir_in, dir_out, dir_out_mask !out_rst needed for soil mask character(len=100) :: fname_in, fname_out, fname_out_mask character(len=15) :: variable_list(max_vars) @@ -48,7 +48,8 @@ program regridStates fname_out, dir_out, gridtype_out, & fname_out_mask, dir_out_mask, & n_vars, variable_list, missing_value, mask_type, & - res_atm_in, res_atm_out, dir_fix ! needed for fv3* grids + ires_in, jres_in, ires_out, jres_out, & + dir_fix ! only needed for fv3* grids !------------------------------------------------------------------------- ! INITIALIZE @@ -102,11 +103,11 @@ program regridStates ! for now, input file and input_mask file are the same thing. call setup_grid(localpet, npets, gridtype_in, & trim(dir_in), trim(fname_in), trim(mask_type), grid_in, & - res_atm_in, trim(dir_fix) ) + ires_in, jres_in, trim(dir_fix) ) call setup_grid(localpet, npets, gridtype_out, & trim(dir_out_mask), trim(fname_out_mask), trim(mask_type), grid_out, & - res_atm_out, trim(dir_fix) ) + ires_out, jres_out, trim(dir_fix) ) !------------------------ ! Create input and output fields @@ -158,7 +159,7 @@ program regridStates !------------------------ ! read data into input fields - call read_into_fields(localpet, res_atm_in, res_atm_in, trim(fname_in), n_vars, & + call read_into_fields(localpet, ires_in, jres_in, trim(fname_in), n_vars, & variable_list(1:n_vars), trim(dir_in), fields_in) call cpu_time(t2) @@ -208,7 +209,7 @@ program regridStates if (localpet==0) print*,'** Writing out regridded fields' - call write_from_fields(localpet, res_atm_out, res_atm_out, trim(fname_out), & + call write_from_fields(localpet, ires_out, jres_out, trim(fname_out), & n_vars, variable_list(1:n_vars), trim(dir_out), fields_out) From af1b4fdbce6789d9bd43cc30dcd48c7055cdc9bc Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 13 Nov 2024 16:42:14 +0000 Subject: [PATCH 05/15] pre-adding derived type --- src/regridStates/grids_IO.F90 | 49 +++++++++++++++++++++---------- src/regridStates/regridStates.F90 | 6 ++-- 2 files changed, 38 insertions(+), 17 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index e906e94..718a28a 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -15,7 +15,8 @@ module grids_IO integer, public, parameter :: n_tiles=6 ! number tiles in fv3 grid integer, public, parameter :: vtype_water=0, & ! TO DO - which veg classification is this? vtype_landice=15 ! used for soil mask - + !type, public :: + public :: setup_grid, & read_into_fields, & write_from_fields @@ -24,14 +25,14 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested - subroutine setup_grid(localpet, npets, file_type, & + subroutine setup_grid(localpet, npets, grid_type, & dir_mask, fname_mask, mask_type, mod_grid, & ires, jres, dir_fix ) implicit none ! INTENT IN - character(7), intent(in) :: file_type + character(7), intent(in) :: grid_type character(*), intent(in) :: fname_mask, dir_mask character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets @@ -54,7 +55,7 @@ subroutine setup_grid(localpet, npets, file_type, & !-------------------------- ! Create grid object, and set up pet distribution - select case (file_type) + select case (grid_type) case ('fv3_rst') call create_grid_fv3(ires, trim(dir_fix), npets, localpet ,mod_grid) mask_variable(1) = 'vtype ' @@ -62,7 +63,7 @@ subroutine setup_grid(localpet, npets, file_type, & call create_grid_gauss(ires, jres, npets, localpet, mod_grid) mask_variable(1) = 'soilsnow_mask ' case default - call error_handler("unknown file_type in setup_grid", 1) + call error_handler("unknown grid_type in setup_grid", 1) end select !-------------------------- @@ -79,7 +80,7 @@ subroutine setup_grid(localpet, npets, file_type, & if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate, mask_variable", ierr) - call read_into_fields(localpet, ires, jres, trim(fname_mask), 1, mask_variable(1), & + call read_into_fields(localpet, ires, jres, trim(fname_mask), grid_type, 1, mask_variable(1), & dir_mask, mask_field(1)) ! get pointer to vegtype @@ -121,7 +122,7 @@ subroutine setup_grid(localpet, npets, file_type, & end subroutine setup_grid ! read variables from fv3 netcdf restart file into ESMF Fields - subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variable_list, & + subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, grid_type, n_vars, variable_list, & dir_read, fields) implicit none @@ -130,6 +131,7 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variabl integer, intent(in) :: localpet, i_dim, j_dim, n_vars character(*), intent(in) :: fname_read character(*), intent(in) :: dir_read + character(7), intent(in) :: grid_type character(len=15), dimension(n_vars), intent(in) :: variable_list @@ -138,6 +140,7 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variabl ! LOCAL integer :: tt, id_var, ncid, ierr, v + integer :: n_files character(len=1) :: tchar character(len=500) :: fname real(esmf_kind_r8), allocatable :: array2D(:,:) @@ -151,19 +154,34 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, n_vars, variabl allocate(array2D(0,0)) end if - do tt = 1, n_tiles + select case (grid_type) + case ('fv3_rst') + n_files=n_tiles + case ('gau_inc') + n_files=1 + case default + call error_handler("unknown grid_type in read into fields", 1) + end select + + do tt = 1, n_files ! read from restart if (localpet == 0) then - write(tchar,'(i1)') tt - fname = dir_read//"/"//fname_read//".tile"//tchar//".nc" + if ( n_files > 1) then + write(tchar,'(i1)') tt + fname = dir_read//"/"//fname_read//".tile"//tchar//".nc" + else + fname = dir_read//"/"//fname_read + endif + + print *, 'Reading ', trim(fname) ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) call netcdf_err(ierr, 'opening: '//trim(fname) ) do v =1, n_vars - print *, 'Reading ', trim(variable_list(v)), ' into field, tile', tchar + print *, 'Reading ', trim(variable_list(v)) ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var) call netcdf_err(ierr, 'reading variable id' ) @@ -298,11 +316,12 @@ subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) enddo ! mosaic file - write(rchar,'(i3.3)') res_atm + write(rchar,'(i3)') res_atm + !write(rchar,'(i)') res_atm - dir_fix_res = dir_fix//"/C"//trim(rchar)//"/" + dir_fix_res = dir_fix//"/C"//trim(adjustl(rchar))//"/" - fname = trim(dir_fix_res)//"/C"//trim(rchar)// "_mosaic.nc" + fname = trim(dir_fix_res)//"/C"//trim(adjustl(rchar))// "_mosaic.nc" ! create the grid fv3_grid = ESMF_GridCreateMosaic(filename=trim(fname), & @@ -331,7 +350,7 @@ subroutine create_grid_gauss(i_dim, j_dim, npets, localpet, gauss_grid) polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE - if (localpet == 0) print*," creating fv3 gauss for ", i_dim, j_dim + if (localpet == 0) print*," creating gauss grid for ", i_dim, j_dim gauss_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/i_dim,j_dim/), & polekindflag=polekindflag, & diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 70accf2..915a701 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -41,6 +41,8 @@ program regridStates type(esmf_routehandle) :: regrid_route real(esmf_kind_r8), pointer :: ptr_out(:,:) + + real :: t1, t2, t3, t4 ! see README for details of namelist variables. @@ -159,7 +161,7 @@ program regridStates !------------------------ ! read data into input fields - call read_into_fields(localpet, ires_in, jres_in, trim(fname_in), n_vars, & + call read_into_fields(localpet, ires_in, jres_in, trim(fname_in), gridtype_in, n_vars, & variable_list(1:n_vars), trim(dir_in), fields_in) call cpu_time(t2) @@ -175,7 +177,7 @@ program regridStates dstField=fields_out(1), dstMaskValues=(/0/), & ! allow unmapped grid cells, without returning error unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & - polemethod=ESMF_POLEMETHOD_ALLAVG, & + ! polemethod=ESMF_POLEMETHOD_ALLAVG, & ! fill un-mapped grid cells with a neighbour extrapMethod=ESMF_EXTRAPMETHOD_CREEP, & ! number of "levels" of neighbours to search for a value From 30fa44e6ee3e1585727c0718163c230fb384b91c Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 13 Nov 2024 21:26:21 +0000 Subject: [PATCH 06/15] Moved grid setup details into a derived type. --- src/regridStates/grids_IO.F90 | 52 ++++++++++++++-------------- src/regridStates/regridStates.F90 | 56 ++++++++++++++++++++++--------- 2 files changed, 68 insertions(+), 40 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 718a28a..a876b26 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -15,7 +15,15 @@ module grids_IO integer, public, parameter :: n_tiles=6 ! number tiles in fv3 grid integer, public, parameter :: vtype_water=0, & ! TO DO - which veg classification is this? vtype_landice=15 ! used for soil mask - !type, public :: + type, public :: grid_setup_type + character(7) :: descriptor + character(100) :: fname + character(100) :: dir + character(100) :: fname_mask + character(100) :: dir_mask + character(100) :: fname_coord + character(100) :: dir_coord + end type public :: setup_grid, & read_into_fields, & @@ -25,20 +33,17 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested - subroutine setup_grid(localpet, npets, grid_type, & - dir_mask, fname_mask, mask_type, mod_grid, & - ires, jres, dir_fix ) + subroutine setup_grid(localpet, npets, grid_setup, & + mask_type, mod_grid, & + ires, jres) implicit none ! INTENT IN - character(7), intent(in) :: grid_type - character(*), intent(in) :: fname_mask, dir_mask + type(grid_setup_type), intent(in) :: grid_setup character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets integer, intent(in) :: ires, jres - ! NEEDED FOR FV3 GRID - character(*), intent(in) :: dir_fix ! INTENT OUT type(esmf_grid) :: mod_grid @@ -51,19 +56,18 @@ subroutine setup_grid(localpet, npets, grid_type, & integer :: ierr, ncid, tile character(len=15) :: mask_variable(1) - !-------------------------- ! Create grid object, and set up pet distribution - select case (grid_type) + select case (grid_setup%descriptor) case ('fv3_rst') - call create_grid_fv3(ires, trim(dir_fix), npets, localpet ,mod_grid) + call create_grid_fv3(ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) mask_variable(1) = 'vtype ' case ('gau_inc') call create_grid_gauss(ires, jres, npets, localpet, mod_grid) mask_variable(1) = 'soilsnow_mask ' case default - call error_handler("unknown grid_type in setup_grid", 1) + call error_handler("unknown grid_setup%descriptor in setup_grid", 1) end select !-------------------------- @@ -80,8 +84,9 @@ subroutine setup_grid(localpet, npets, grid_type, & if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate, mask_variable", ierr) - call read_into_fields(localpet, ires, jres, trim(fname_mask), grid_type, 1, mask_variable(1), & - dir_mask, mask_field(1)) + call read_into_fields(localpet, ires, jres, trim(grid_setup%fname_mask), & + trim(grid_setup%dir_mask), grid_setup, 1, & + mask_variable(1), mask_field(1)) ! get pointer to vegtype call ESMF_FieldGet(mask_field(1), & @@ -122,8 +127,8 @@ subroutine setup_grid(localpet, npets, grid_type, & end subroutine setup_grid ! read variables from fv3 netcdf restart file into ESMF Fields - subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, grid_type, n_vars, variable_list, & - dir_read, fields) + subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & + grid_setup, n_vars, variable_list, fields) implicit none @@ -131,7 +136,7 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, grid_type, n_va integer, intent(in) :: localpet, i_dim, j_dim, n_vars character(*), intent(in) :: fname_read character(*), intent(in) :: dir_read - character(7), intent(in) :: grid_type + type(grid_setup_type), intent(in) :: grid_setup character(len=15), dimension(n_vars), intent(in) :: variable_list @@ -154,13 +159,13 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, grid_type, n_va allocate(array2D(0,0)) end if - select case (grid_type) + select case (grid_setup%descriptor) case ('fv3_rst') n_files=n_tiles case ('gau_inc') n_files=1 case default - call error_handler("unknown grid_type in read into fields", 1) + call error_handler("unknown grid_setup%descriptor in read into fields", 1) end select do tt = 1, n_files @@ -209,8 +214,8 @@ end subroutine read_into_fields ! write variables from ESMF Fields into netcdf restart-like file - subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, n_vars, variable_list, & - dir_out, fields) + subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, dir_out, & + n_vars, variable_list, fields) implicit none @@ -305,7 +310,7 @@ subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) integer :: decomptile(2,n_tiles) character(len=5) :: rchar - character(len=500) :: fname, dir_fix_res + character(len=200) :: fname, dir_fix_res if (localpet == 0) print*," creating fv3 grid for ", res_atm @@ -317,10 +322,7 @@ subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) ! mosaic file write(rchar,'(i3)') res_atm - !write(rchar,'(i)') res_atm - dir_fix_res = dir_fix//"/C"//trim(adjustl(rchar))//"/" - fname = trim(dir_fix_res)//"/C"//trim(adjustl(rchar))// "_mosaic.nc" ! create the grid diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 915a701..6460486 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -12,7 +12,8 @@ program regridStates use grids_IO, only : setup_grid, & write_from_fields, & read_into_fields, & - n_tiles + n_tiles, & + grid_setup_type use utilities, only : error_handler @@ -21,9 +22,9 @@ program regridStates integer, parameter :: max_vars = 10 ! increase if wish to specify more variables ! namelist inputs - character(len=500) :: dir_fix + character(len=100) :: dir_fix integer :: ires_in, jres_in, ires_out, jres_out - character(len=500) :: dir_in, dir_out, dir_out_mask !out_rst needed for soil mask + character(len=100) :: dir_in, dir_out, dir_out_mask character(len=100) :: fname_in, fname_out, fname_out_mask character(len=15) :: variable_list(max_vars) character(len=10) :: mask_type @@ -31,6 +32,8 @@ program regridStates integer :: n_vars real(esmf_kind_r8) :: missing_value ! value given to unmapped cells in the output grid + type(grid_setup_type) :: grid_setup_in, grid_setup_out + integer :: ierr, localpet, npets integer :: v, SRCTERM @@ -96,20 +99,41 @@ program regridStates if (ierr /= 0) call error_handler("READING regrid NAMELIST.", ierr) close (41) + + grid_setup_in%descriptor = gridtype_in + + ! if in + grid_setup_in%dir_mask = dir_in ! for fv3 and gauss use input file for mask + grid_setup_in%fname_mask = fname_in + grid_setup_in%dir = dir_in ! for fv3 and gauss use input file for mask + grid_setup_in%fname = fname_in + ! if gau_inc + !grid_setup_in%dir_coord = dir_in + grid_setup_in%dir_coord = dir_fix + !grid_setup_in%fname_coord = fname_in + + grid_setup_out%descriptor = gridtype_out + ! if out + grid_setup_out%dir_mask = dir_out_mask + grid_setup_out%fname_mask = fname_out_mask + grid_setup_out%dir = dir_out + grid_setup_out%fname = fname_out + ! if fv3 (will construct mosaic fname later) + grid_setup_out%dir_coord = dir_fix + !------------------------ ! Create esmf grid objects for input and output grids, and add land masks ! TO DO - can we make the number of tasks more flexible for fv3 if (localpet==0) print*,'** Setting up grids' - ! for now, input file and input_mask file are the same thing. - call setup_grid(localpet, npets, gridtype_in, & - trim(dir_in), trim(fname_in), trim(mask_type), grid_in, & - ires_in, jres_in, trim(dir_fix) ) + call setup_grid(localpet, npets, grid_setup_in, & + trim(mask_type), grid_in, & + ires_in, jres_in) - call setup_grid(localpet, npets, gridtype_out, & - trim(dir_out_mask), trim(fname_out_mask), trim(mask_type), grid_out, & - ires_out, jres_out, trim(dir_fix) ) + call setup_grid(localpet, npets, grid_setup_out, & + trim(mask_type), grid_out, & + ires_out, jres_out) !------------------------ ! Create input and output fields @@ -161,8 +185,9 @@ program regridStates !------------------------ ! read data into input fields - call read_into_fields(localpet, ires_in, jres_in, trim(fname_in), gridtype_in, n_vars, & - variable_list(1:n_vars), trim(dir_in), fields_in) + call read_into_fields(localpet, ires_in, jres_in, trim(grid_setup_in%fname), & + trim(grid_setup_in%dir), grid_setup_in, n_vars, & + variable_list(1:n_vars), fields_in) call cpu_time(t2) !------------------------ @@ -177,7 +202,7 @@ program regridStates dstField=fields_out(1), dstMaskValues=(/0/), & ! allow unmapped grid cells, without returning error unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & - ! polemethod=ESMF_POLEMETHOD_ALLAVG, & + polemethod=ESMF_POLEMETHOD_ALLAVG, & ! fill un-mapped grid cells with a neighbour extrapMethod=ESMF_EXTRAPMETHOD_CREEP, & ! number of "levels" of neighbours to search for a value @@ -211,8 +236,9 @@ program regridStates if (localpet==0) print*,'** Writing out regridded fields' - call write_from_fields(localpet, ires_out, jres_out, trim(fname_out), & - n_vars, variable_list(1:n_vars), trim(dir_out), fields_out) + call write_from_fields(localpet, ires_out, jres_out, trim(grid_setup_out%fname), & + trim(grid_setup_out%dir), n_vars, variable_list(1:n_vars), & + fields_out) ! clean up From 03d788e7a51fff176ebead7fd2ac07d4533695f9 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 14 Nov 2024 20:05:41 +0000 Subject: [PATCH 07/15] Fixed up new data type. --- src/regridStates/grids_IO.F90 | 12 +-- src/regridStates/regridStates.F90 | 148 ++++++++++++++++++++---------- 2 files changed, 104 insertions(+), 56 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index a876b26..05a1e57 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -23,6 +23,8 @@ module grids_IO character(100) :: dir_mask character(100) :: fname_coord character(100) :: dir_coord + integer :: ires + integer :: jres end type public :: setup_grid, & @@ -34,8 +36,7 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested subroutine setup_grid(localpet, npets, grid_setup, & - mask_type, mod_grid, & - ires, jres) + mask_type, mod_grid ) implicit none @@ -43,7 +44,6 @@ subroutine setup_grid(localpet, npets, grid_setup, & type(grid_setup_type), intent(in) :: grid_setup character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets - integer, intent(in) :: ires, jres ! INTENT OUT type(esmf_grid) :: mod_grid @@ -61,10 +61,10 @@ subroutine setup_grid(localpet, npets, grid_setup, & select case (grid_setup%descriptor) case ('fv3_rst') - call create_grid_fv3(ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) + call create_grid_fv3(grid_setup%ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) mask_variable(1) = 'vtype ' case ('gau_inc') - call create_grid_gauss(ires, jres, npets, localpet, mod_grid) + call create_grid_gauss(grid_setup%ires, grid_setup%jres, npets, localpet, mod_grid) mask_variable(1) = 'soilsnow_mask ' case default call error_handler("unknown grid_setup%descriptor in setup_grid", 1) @@ -84,7 +84,7 @@ subroutine setup_grid(localpet, npets, grid_setup, & if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate, mask_variable", ierr) - call read_into_fields(localpet, ires, jres, trim(grid_setup%fname_mask), & + call read_into_fields(localpet, grid_setup%ires, grid_setup%jres, trim(grid_setup%fname_mask), & trim(grid_setup%dir_mask), grid_setup, 1, & mask_variable(1), mask_field(1)) diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 6460486..ff3850a 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -22,15 +22,10 @@ program regridStates integer, parameter :: max_vars = 10 ! increase if wish to specify more variables ! namelist inputs - character(len=100) :: dir_fix - integer :: ires_in, jres_in, ires_out, jres_out - character(len=100) :: dir_in, dir_out, dir_out_mask - character(len=100) :: fname_in, fname_out, fname_out_mask character(len=15) :: variable_list(max_vars) - character(len=10) :: mask_type - character(len=7) :: gridtype_in, gridtype_out integer :: n_vars real(esmf_kind_r8) :: missing_value ! value given to unmapped cells in the output grid + character(len=10) :: mask_type type(grid_setup_type) :: grid_setup_in, grid_setup_out @@ -43,20 +38,14 @@ program regridStates type(esmf_field), allocatable :: fields_out(:) type(esmf_routehandle) :: regrid_route real(esmf_kind_r8), pointer :: ptr_out(:,:) - - + + integer :: ut real :: t1, t2, t3, t4 ! see README for details of namelist variables. - namelist /config/ fname_in, dir_in, gridtype_in, & - fname_out, dir_out, gridtype_out, & - fname_out_mask, dir_out_mask, & - n_vars, variable_list, missing_value, mask_type, & - ires_in, jres_in, ires_out, jres_out, & - dir_fix ! only needed for fv3* grids + namelist /config/ n_vars, variable_list, missing_value, mask_type -!------------------------------------------------------------------------- ! INITIALIZE !------------------------------------------------------------------------- @@ -93,33 +82,14 @@ program regridStates missing_value=-999. ! set defualt - open(41, file='regrid.nml', iostat=ierr) + open(newunit=ut, file='regrid.nml', iostat=ierr) if (ierr /= 0) call error_handler("OPENING regrid NAMELIST.", ierr) - read(41, nml=config, iostat=ierr) - if (ierr /= 0) call error_handler("READING regrid NAMELIST.", ierr) - close (41) - + read(ut, nml=config, iostat=ierr) + if (ierr /= 0) call error_handler("OPENING config NAMELIST.", ierr) + call readin_setup(ut,"input",grid_setup_in) + call readin_setup(ut,"output",grid_setup_out) + close (ut) - grid_setup_in%descriptor = gridtype_in - - ! if in - grid_setup_in%dir_mask = dir_in ! for fv3 and gauss use input file for mask - grid_setup_in%fname_mask = fname_in - grid_setup_in%dir = dir_in ! for fv3 and gauss use input file for mask - grid_setup_in%fname = fname_in - ! if gau_inc - !grid_setup_in%dir_coord = dir_in - grid_setup_in%dir_coord = dir_fix - !grid_setup_in%fname_coord = fname_in - - grid_setup_out%descriptor = gridtype_out - ! if out - grid_setup_out%dir_mask = dir_out_mask - grid_setup_out%fname_mask = fname_out_mask - grid_setup_out%dir = dir_out - grid_setup_out%fname = fname_out - ! if fv3 (will construct mosaic fname later) - grid_setup_out%dir_coord = dir_fix !------------------------ ! Create esmf grid objects for input and output grids, and add land masks @@ -128,12 +98,10 @@ program regridStates if (localpet==0) print*,'** Setting up grids' call setup_grid(localpet, npets, grid_setup_in, & - trim(mask_type), grid_in, & - ires_in, jres_in) + trim(mask_type), grid_in ) call setup_grid(localpet, npets, grid_setup_out, & - trim(mask_type), grid_out, & - ires_out, jres_out) + trim(mask_type), grid_out ) !------------------------ ! Create input and output fields @@ -185,9 +153,9 @@ program regridStates !------------------------ ! read data into input fields - call read_into_fields(localpet, ires_in, jres_in, trim(grid_setup_in%fname), & - trim(grid_setup_in%dir), grid_setup_in, n_vars, & - variable_list(1:n_vars), fields_in) + call read_into_fields(localpet, grid_setup_in%ires, grid_setup_in%jres, & + trim(grid_setup_in%fname), trim(grid_setup_in%dir), & + grid_setup_in, n_vars, variable_list(1:n_vars), fields_in) call cpu_time(t2) !------------------------ @@ -236,9 +204,9 @@ program regridStates if (localpet==0) print*,'** Writing out regridded fields' - call write_from_fields(localpet, ires_out, jres_out, trim(grid_setup_out%fname), & - trim(grid_setup_out%dir), n_vars, variable_list(1:n_vars), & - fields_out) + call write_from_fields(localpet, grid_setup_out%ires, grid_setup_out%jres, & + trim(grid_setup_out%fname), trim(grid_setup_out%dir), & + n_vars, variable_list(1:n_vars), fields_out) ! clean up @@ -279,3 +247,83 @@ program regridStates print*,"** DONE.", localpet end program regridStates + + subroutine readin_setup(unt,namel,grid_setup) +! subroutine to read in namelists, and convert +! values into setupgrid. +! Also fills in some values, and tests have all +! needed vals, according to the selected grid type. + + use grids_IO, only : grid_setup_type + use utilities, only : error_handler + + implicit none + + ! INPUTS + integer, intent(in) :: unt + character(*), intent(in) :: namel + ! OUTPUTS + type(grid_setup_type), intent(out) :: grid_setup + + character(len=7) :: gridtype + character(len=100) :: fname, fname_mask + character(len=100) :: dir, dir_mask, dir_fix + character(len=4) :: default_str="NULL" + integer :: ires, jres + integer :: ierr + + namelist /input/ fname, dir, gridtype, & + fname_mask, dir_mask, & + dir_fix, ires, jres + + namelist /output/ fname, dir, gridtype, & + fname_mask, dir_mask, & + dir_fix, ires, jres + + ! set defaults + fname = default_str + dir = default_str + fname_mask = default_str + dir_mask = default_str + ires = 0 + jres = 0 + + if (namel=="input") then + read(unt, nml=input, iostat=ierr) + if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) + elseif (namel=="output") then + read(unt, nml=output, iostat=ierr) + if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) + endif + + grid_setup%descriptor = gridtype + grid_setup%dir = dir ! for fv3 and gauss use input file for mask + grid_setup%fname = fname + + select case (namel) + case ("input") + grid_setup%dir_mask = dir ! for fv3 and gauss use input file for mask + grid_setup%fname_mask = fname + case ("output") + grid_setup%dir_mask = dir_mask ! for fv3 and gauss use input file for mask + grid_setup%fname_mask = fname_mask + case default + call error_handler("unknown namel in readin_setup", 1) + end select + + select case ( grid_setup%descriptor ) + case ("fv3_rst") + grid_setup%dir_coord = dir_fix + grid_setup%ires = ires + grid_setup%jres = ires ! always same for fv3 + case ("inc_gau") + grid_setup%dir_coord = dir + grid_setup%fname_coord = fname + grid_setup%ires = ires + grid_setup%jres = jres + case default + call error_handler("unknown grid_setup%descriptor in readin_setup", 1) + end select + + end subroutine +!------------------------------------------------------------------------- From b6a6819ca5bb35663aaee463e0ea8877a34997a3 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 5 Dec 2024 16:30:10 +0000 Subject: [PATCH 08/15] Sorted out Coord bugs. Started generalizing reading in lat/lon for Gauss grid. --- src/regridStates/grids_IO.F90 | 232 ++++++++++++++++++++++++++++-- src/regridStates/regridStates.F90 | 2 +- 2 files changed, 224 insertions(+), 10 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 05a1e57..4a6935a 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -35,6 +35,7 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested + subroutine setup_grid(localpet, npets, grid_setup, & mask_type, mod_grid ) @@ -64,14 +65,14 @@ subroutine setup_grid(localpet, npets, grid_setup, & call create_grid_fv3(grid_setup%ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) mask_variable(1) = 'vtype ' case ('gau_inc') - call create_grid_gauss(grid_setup%ires, grid_setup%jres, npets, localpet, mod_grid) + call create_grid_gauss(grid_setup, npets, localpet, mod_grid) mask_variable(1) = 'soilsnow_mask ' case default call error_handler("unknown grid_setup%descriptor in setup_grid", 1) end select !-------------------------- -! Calcalate and add the mask +! Calculate and add the mask if (mask_type=="soil") then ! mask out ocean and glaciers, using vegetation class @@ -180,13 +181,13 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & fname = dir_read//"/"//fname_read endif - print *, 'Reading ', trim(fname) + !print *, 'Reading ', trim(fname) ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) call netcdf_err(ierr, 'opening: '//trim(fname) ) do v =1, n_vars - print *, 'Reading ', trim(variable_list(v)) + !print *, 'Reading ', trim(variable_list(v)) ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var) call netcdf_err(ierr, 'reading variable id' ) @@ -212,6 +213,107 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & end subroutine read_into_fields + ! read lat and lon from history file (Guassian grid), and put into 2D fields + ! needed for setting up the Gauusian grid. + + subroutine latlon_read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & + grid_setup, variable_list, fields) + + implicit none + + ! INTENT IN + integer, intent(in) :: localpet, i_dim, j_dim + character(*), intent(in) :: fname_read + character(*), intent(in) :: dir_read + type(grid_setup_type), intent(in) :: grid_setup + + character(len=15), dimension(2), intent(in) :: variable_list + + ! INTENT OUT + type(esmf_field), dimension(2), intent(out) :: fields + + ! LOCAL + integer :: id_var, id_lon, ncid, ierr + integer :: i,j + character(len=500) :: fname + real(esmf_kind_r8) :: vec_lat(j_dim), vec_lon(i_dim) + real(esmf_kind_r8), allocatable :: array_lat(:,:), array_lon(:,:) + character(len=30) :: vname + + ! Create the fields + do v=1,2 + vname="gauss_grid_" // trim(var_list(v) + fields(v) = ESMF_FieldCreate(gauss_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name=vname, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", ierr) + + enddo + + if (localpet==0) then + allocate(array_lat(i_dim, j_dim)) + allocate(array_lon(i_dim, j_dim)) + else + allocate(array_lat(0,0)) + allocate(array_lon(0,0)) + end if + + ! read from restart + if (localpet == 0) then + + fname = dir_read//"/"//fname_read + + print *, 'Reading ', trim(fname) + + ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) + call netcdf_err(ierr, 'opening: '//trim(fname) ) + + ! read latitude + print *, 'Reading ', trim(variable_list(1)) + ierr=nf90_inq_varid(ncid, trim(variable_list(1)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + + ierr=nf90_get_var(ncid, id_var, vec_lat) + call netcdf_err(ierr, 'reading variable' ) + + do i = 1, i_dim + array_lat(i,:) =vec_lat + enddo + + ! read longitude + print *, 'Reading ', trim(variable_list(2)) + ierr=nf90_inq_varid(ncid, trim(variable_list(2)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + + ierr=nf90_get_var(ncid, id_var, vec_lon) + call netcdf_err(ierr, 'reading variable' ) + ierr = nf90_close(ncid) + + do i = j, j_dim + array_lon(:,j) =vec_lon + enddo + + endif + + ! scatter latitude + if (localpet==0) print *, 'scattering lat' + call ESMF_FieldScatter(fields(1), array_lat, rootpet=0, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", ierr) + + ! scatter longitude + if (localpet==0) print *, 'scattering lon' + call ESMF_FieldScatter(fields(2), array_lon, rootpet=0, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", ierr) + + ! clean up + deallocate(array_lat) + deallocate(array_lon) + + end subroutine latlon_read_into_fields ! write variables from ESMF Fields into netcdf restart-like file subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, dir_out, & @@ -338,20 +440,32 @@ subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid) end subroutine create_grid_fv3 - subroutine create_grid_gauss(i_dim, j_dim, npets, localpet, gauss_grid) + subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) ! INTENT IN - integer, intent(in) :: i_dim, j_dim + type(grid_setup_type), intent(in) :: grid_setup integer, intent(in) :: npets, localpet ! INTENT OUT - type(esmf_grid) :: gauss_grid + type(esmf_grid) :: gauss_grid - type(esmf_polekind_flag) :: polekindflag(2) - integer :: ierr + type(esmf_field) :: latlon_fields(2) + type(esmf_polekind_flag) :: polekindflag(2) + real(esmf_kind_r8), pointer :: lon_ptr_field(:,:), lon_ptr_coord(:,:) + real(esmf_kind_r8), pointer :: lat_ptr_field(:,:), lat_ptr_coord(:,:) + character(len=15), dimension(2) :: latlon_variables + + integer :: ierr, i_dim, j_dim polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + ! Gaussian increment files lat/lon name + latlon_variables(1) = 'latitude' + latlon_variables(2) = 'longitude' + + i_dim = grid_setup%ires + j_dim = grid_setup%jres + if (localpet == 0) print*," creating gauss grid for ", i_dim, j_dim gauss_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/i_dim,j_dim/), & @@ -364,5 +478,105 @@ subroutine create_grid_gauss(i_dim, j_dim, npets, localpet, gauss_grid) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN GridCreate1PeriDim", ierr) + ! read lat lon coordinates into fields + call latlon_read_into_fields(localpet, i_dim, j_dim, & + trim(grid_setup%fname_coord), trim(grid_setup%dir_coord), & + grid_setup, latlon_variables, latlon_fields) + + ! add coordinates to the grid + call ESMF_GridAddCoord(gauss_grid, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", ierr) + + nullify(lon_ptr_coord) + nullify(lon_ptr_field) + + ! get pointer to lon coord + !if (localpet == 0) print*," getting GridCoord for long" + print*," getting GridCoord for long", localpet + call ESMF_GridGetCoord(gauss_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=1, & + farrayPtr=lon_ptr_coord, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord longitude ", ierr) + + ! fetch lon from field into pointer + call ESMF_FieldGet(latlon_fields(2), & + farrayPtr=lon_ptr_field, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet longitude", ierr) + + ! set coord pointe to field pointer + lon_ptr_coord = lon_ptr_field + + nullify(lat_ptr_coord) + nullify(lat_ptr_field) + + ! get pointer to lat coord. Need bounds? + call ESMF_GridGetCoord(gauss_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=2, & + !computationalLBound=clb, & + !computationalUBound=cub, & + farrayPtr=lat_ptr_coord, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord latitude", ierr) + + ! fetch lat from field into pointer + call ESMF_FieldGet(latlon_fields(1), & + farrayPtr=lat_ptr_field, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet latitude", ierr) + + ! set lat coord to fiel values + lat_ptr_coord = lat_ptr_field + + ! assign the corners + ! call ESMF_GridAddCoord(input_grid, & + ! staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) + !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + ! call error_handler("IN GridAddCoord", rc) + + !print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + !nullify(lon_corner_src_ptr) + !call ESMF_GridGetCoord(input_grid, & + ! staggerLoc=ESMF_STAGGERLOC_CORNER, & + ! coordDim=1, & + ! farrayPtr=lon_corner_src_ptr, rc=rc) + !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + ! call error_handler("IN GridGetCoord", rc) + + !print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + !nullify(lat_corner_src_ptr) + !call ESMF_GridGetCoord(input_grid, & + ! staggerLoc=ESMF_STAGGERLOC_CORNER, & + ! coordDim=2, & + ! computationalLBound=clb, & + ! computationalUBound=cub, & + ! farrayPtr=lat_corner_src_ptr, rc=rc) + !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + ! call error_handler("IN GridGetCoord", rc) + + !do j = clb(2), cub(2) + ! do i = clb(1), cub(1) + ! lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon) + ! if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8 + ! if (j == 1) then + ! lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8 + ! cycle + ! endif + ! if (j == jp1_input) then + ! lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8 + ! cycle + ! endif + ! lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j)) + ! enddo + !enddo + end subroutine create_grid_gauss + end module grids_IO diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index ff3850a..f9d3a09 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -316,7 +316,7 @@ subroutine readin_setup(unt,namel,grid_setup) grid_setup%dir_coord = dir_fix grid_setup%ires = ires grid_setup%jres = ires ! always same for fv3 - case ("inc_gau") + case ("gau_inc") grid_setup%dir_coord = dir grid_setup%fname_coord = fname grid_setup%ires = ires From e4561ceccc0e7de570f031c196e30b9fabed2e45 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 6 Dec 2024 20:10:08 +0000 Subject: [PATCH 09/15] Finished genaalizing lat,lon input. Switched to using weightgen file for Gauss grid coords. --- src/regridStates/grids_IO.F90 | 152 ++++++++++++++++++------------ src/regridStates/regridStates.F90 | 34 +++---- 2 files changed, 107 insertions(+), 79 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 4a6935a..22ec7fd 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -216,49 +216,62 @@ end subroutine read_into_fields ! read lat and lon from history file (Guassian grid), and put into 2D fields ! needed for setting up the Gauusian grid. - subroutine latlon_read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & - grid_setup, variable_list, fields) + subroutine lonlat_read_into_fields(localpet, i_dim, j_dim, fname_read, dir_read, & + grid_setup, gauss_grid, n_vars, lonvar_list, & + latvar_list, lon_fields, lat_fields) implicit none ! INTENT IN - integer, intent(in) :: localpet, i_dim, j_dim + integer, intent(in) :: localpet, i_dim, j_dim, n_vars character(*), intent(in) :: fname_read character(*), intent(in) :: dir_read type(grid_setup_type), intent(in) :: grid_setup + type(esmf_grid), intent(in) :: gauss_grid - character(len=15), dimension(2), intent(in) :: variable_list + character(len=15), dimension(n_vars), intent(in) :: lonvar_list + character(len=15), dimension(n_vars), intent(in) :: latvar_list ! INTENT OUT - type(esmf_field), dimension(2), intent(out) :: fields + type(esmf_field), dimension(n_vars), intent(out) :: lon_fields + type(esmf_field), dimension(n_vars), intent(out) :: lat_fields ! LOCAL - integer :: id_var, id_lon, ncid, ierr - integer :: i,j + integer :: id_var, ncid, ierr + integer :: i,j,v character(len=500) :: fname - real(esmf_kind_r8) :: vec_lat(j_dim), vec_lon(i_dim) - real(esmf_kind_r8), allocatable :: array_lat(:,:), array_lon(:,:) + real(esmf_kind_r8) :: vec_lat(i_dim*j_dim), vec_lon(i_dim*j_dim) + real(esmf_kind_r8), allocatable :: array_lat(:,:,:), array_lon(:,:,:), array2d(:,:) character(len=30) :: vname ! Create the fields - do v=1,2 - vname="gauss_grid_" // trim(var_list(v) - fields(v) = ESMF_FieldCreate(gauss_grid, & + do v=1,n_vars + vname="gauss_grid_" // trim(lonvar_list(v)) + lon_fields(v) = ESMF_FieldCreate(gauss_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & name=vname, rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate", ierr) + vname="gauss_grid_" // trim(latvar_list(v)) + lat_fields(v) = ESMF_FieldCreate(gauss_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name=vname, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", ierr) enddo - + if (localpet==0) then - allocate(array_lat(i_dim, j_dim)) - allocate(array_lon(i_dim, j_dim)) + allocate(array_lat(n_vars, i_dim, j_dim)) + allocate(array_lon(n_vars, i_dim, j_dim)) + allocate(array2d( i_dim, j_dim)) else - allocate(array_lat(0,0)) - allocate(array_lon(0,0)) - end if + allocate(array_lat(n_vars,0,0)) + allocate(array_lon(n_vars,0,0)) + allocate(array2d(0,0)) + endif ! read from restart if (localpet == 0) then @@ -270,50 +283,64 @@ subroutine latlon_read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) call netcdf_err(ierr, 'opening: '//trim(fname) ) - ! read latitude - print *, 'Reading ', trim(variable_list(1)) - ierr=nf90_inq_varid(ncid, trim(variable_list(1)), id_var) - call netcdf_err(ierr, 'reading variable id' ) + do v=1, n_vars - ierr=nf90_get_var(ncid, id_var, vec_lat) - call netcdf_err(ierr, 'reading variable' ) + print *, 'Reading ', trim(lonvar_list(v)) + ierr=nf90_inq_varid(ncid, trim(lonvar_list(v)), id_var) + call netcdf_err(ierr, 'reading variable id' ) - do i = 1, i_dim - array_lat(i,:) =vec_lat - enddo + ierr=nf90_get_var(ncid, id_var, vec_lon) + call netcdf_err(ierr, 'reading variable' ) - ! read longitude - print *, 'Reading ', trim(variable_list(2)) - ierr=nf90_inq_varid(ncid, trim(variable_list(2)), id_var) - call netcdf_err(ierr, 'reading variable id' ) + array_lon(v,:,:) = reshape(vec_lon,(/i_dim,j_dim/)) - ierr=nf90_get_var(ncid, id_var, vec_lon) - call netcdf_err(ierr, 'reading variable' ) - ierr = nf90_close(ncid) + !do j = 1, j_dim + ! array_lon(v,:,j) =vec_lon + !enddo + print *, 'CSD lons', array_lon(v,1,1), array_lon(v,i_dim,1) + print *, 'CSD lons', array_lon(v,1,j_dim), array_lon(v,i_dim,j_dim) - do i = j, j_dim - array_lon(:,j) =vec_lon - enddo + print *, 'Reading ', trim(latvar_list(v)) + ierr=nf90_inq_varid(ncid, trim(latvar_list(v)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + ierr=nf90_get_var(ncid, id_var, vec_lat) + call netcdf_err(ierr, 'reading variable' ) + ierr = nf90_close(ncid) + + array_lat(v,:,:) = reshape(vec_lat,(/i_dim,j_dim/)) + !do i = 1, i_dim + ! array_lat(v,i,:) =vec_lat + !enddo + print *, 'CSD lats', array_lat(v,1,1), array_lat(v,i_dim,1) + print *, 'CSD lats', array_lat(v,1,j_dim), array_lat(v,i_dim,j_dim) + enddo + ierr = nf90_close(ncid) endif - ! scatter latitude - if (localpet==0) print *, 'scattering lat' - call ESMF_FieldScatter(fields(1), array_lat, rootpet=0, rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", ierr) - - ! scatter longitude - if (localpet==0) print *, 'scattering lon' - call ESMF_FieldScatter(fields(2), array_lon, rootpet=0, rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", ierr) - + do v =1,n_vars + + ! scatter longitudes + if (localpet==0) print *, 'scattering lon', lonvar_list(v) + array2d=array_lon(v,:,:) + call ESMF_FieldScatter(lon_fields(v), array2d, rootpet=0, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", ierr) + + ! scatter latitude + if (localpet==0) print *, 'scattering lat', latvar_list(v) + array2d=array_lat(v,:,:) + call ESMF_FieldScatter(lat_fields(v), array2d, rootpet=0, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", ierr) + enddo + ! clean up deallocate(array_lat) deallocate(array_lon) + deallocate(array2d) - end subroutine latlon_read_into_fields + end subroutine lonlat_read_into_fields ! write variables from ESMF Fields into netcdf restart-like file subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, dir_out, & @@ -449,19 +476,23 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) ! INTENT OUT type(esmf_grid) :: gauss_grid - type(esmf_field) :: latlon_fields(2) type(esmf_polekind_flag) :: polekindflag(2) + + type(esmf_field) :: lon_fields(1) + type(esmf_field) :: lat_fields(1) real(esmf_kind_r8), pointer :: lon_ptr_field(:,:), lon_ptr_coord(:,:) real(esmf_kind_r8), pointer :: lat_ptr_field(:,:), lat_ptr_coord(:,:) - character(len=15), dimension(2) :: latlon_variables + character(len=15), dimension(1) :: lon_variables, lat_variables integer :: ierr, i_dim, j_dim polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE ! Gaussian increment files lat/lon name - latlon_variables(1) = 'latitude' - latlon_variables(2) = 'longitude' + lon_variables(1) = 'grid_center_lon' + lat_variables(1) = 'grid_center_lat' + !lon_variables(1) = 'longitude' + !lat_variables(1) = 'latitude' i_dim = grid_setup%ires j_dim = grid_setup%jres @@ -479,9 +510,10 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) call error_handler("IN GridCreate1PeriDim", ierr) ! read lat lon coordinates into fields - call latlon_read_into_fields(localpet, i_dim, j_dim, & + call lonlat_read_into_fields(localpet, i_dim, j_dim, & trim(grid_setup%fname_coord), trim(grid_setup%dir_coord), & - grid_setup, latlon_variables, latlon_fields) + grid_setup, gauss_grid, 1, lon_variables, lat_variables, & + lon_fields, lat_fields) ! add coordinates to the grid call ESMF_GridAddCoord(gauss_grid, & @@ -503,7 +535,7 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) call error_handler("IN GridGetCoord longitude ", ierr) ! fetch lon from field into pointer - call ESMF_FieldGet(latlon_fields(2), & + call ESMF_FieldGet(lon_fields(1), & farrayPtr=lon_ptr_field, & rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -519,14 +551,12 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) call ESMF_GridGetCoord(gauss_grid, & staggerLoc=ESMF_STAGGERLOC_CENTER, & coordDim=2, & - !computationalLBound=clb, & - !computationalUBound=cub, & farrayPtr=lat_ptr_coord, rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN GridGetCoord latitude", ierr) ! fetch lat from field into pointer - call ESMF_FieldGet(latlon_fields(1), & + call ESMF_FieldGet(lat_fields(1), & farrayPtr=lat_ptr_field, & rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -577,6 +607,8 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) ! enddo !enddo + ! TO DO destrot the fields + end subroutine create_grid_gauss end module grids_IO diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index f9d3a09..2976a97 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -266,25 +266,29 @@ subroutine readin_setup(unt,namel,grid_setup) type(grid_setup_type), intent(out) :: grid_setup character(len=7) :: gridtype - character(len=100) :: fname, fname_mask - character(len=100) :: dir, dir_mask, dir_fix + character(len=100) :: fname, fname_mask, fname_coord + character(len=100) :: dir, dir_mask, dir_coord character(len=4) :: default_str="NULL" integer :: ires, jres integer :: ierr namelist /input/ fname, dir, gridtype, & fname_mask, dir_mask, & - dir_fix, ires, jres + fname_coord, dir_coord, & + ires, jres namelist /output/ fname, dir, gridtype, & fname_mask, dir_mask, & - dir_fix, ires, jres + fname_coord, dir_coord,& + ires, jres ! set defaults fname = default_str dir = default_str fname_mask = default_str dir_mask = default_str + fname_coord = default_str + dir_coord = default_str ires = 0 jres = 0 @@ -300,30 +304,22 @@ subroutine readin_setup(unt,namel,grid_setup) grid_setup%dir = dir ! for fv3 and gauss use input file for mask grid_setup%fname = fname + ! to-do, add routine to check if present. select case (namel) case ("input") - grid_setup%dir_mask = dir ! for fv3 and gauss use input file for mask + grid_setup%dir_mask = dir ! use input file for mask grid_setup%fname_mask = fname case ("output") - grid_setup%dir_mask = dir_mask ! for fv3 and gauss use input file for mask + grid_setup%dir_mask = dir_mask ! need a file on output grid for mask grid_setup%fname_mask = fname_mask case default call error_handler("unknown namel in readin_setup", 1) end select - select case ( grid_setup%descriptor ) - case ("fv3_rst") - grid_setup%dir_coord = dir_fix - grid_setup%ires = ires - grid_setup%jres = ires ! always same for fv3 - case ("gau_inc") - grid_setup%dir_coord = dir - grid_setup%fname_coord = fname - grid_setup%ires = ires - grid_setup%jres = jres - case default - call error_handler("unknown grid_setup%descriptor in readin_setup", 1) - end select + grid_setup%dir_coord = dir_coord + grid_setup%fname_coord = fname_coord + grid_setup%ires = ires + grid_setup%jres = jres end subroutine !------------------------------------------------------------------------- From 3578f93535d2c7866a6a2d8aac300d90ef8b1056 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 6 Dec 2024 23:10:25 +0000 Subject: [PATCH 10/15] Switched to setting up gauss grid from SCRIP file. Runs, results appear reasonable. --- src/regridStates/grids_IO.F90 | 134 +++++++++++++++++++++------------- 1 file changed, 82 insertions(+), 52 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 22ec7fd..5a9b648 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -204,6 +204,7 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & call ESMF_FieldScatter(fields(v), array2D, rootpet=0, tile=tt, rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", ierr) + enddo enddo @@ -213,39 +214,46 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & end subroutine read_into_fields - ! read lat and lon from history file (Guassian grid), and put into 2D fields - ! needed for setting up the Gauusian grid. + ! read lat and lon from SCRIP file, for use in Gaussian grid subroutine lonlat_read_into_fields(localpet, i_dim, j_dim, fname_read, dir_read, & - grid_setup, gauss_grid, n_vars, lonvar_list, & - latvar_list, lon_fields, lat_fields) + grid_setup, gauss_grid, lon_fields, lat_fields) implicit none ! INTENT IN - integer, intent(in) :: localpet, i_dim, j_dim, n_vars + integer, intent(in) :: localpet, i_dim, j_dim character(*), intent(in) :: fname_read character(*), intent(in) :: dir_read type(grid_setup_type), intent(in) :: grid_setup type(esmf_grid), intent(in) :: gauss_grid - character(len=15), dimension(n_vars), intent(in) :: lonvar_list - character(len=15), dimension(n_vars), intent(in) :: latvar_list ! INTENT OUT - type(esmf_field), dimension(n_vars), intent(out) :: lon_fields - type(esmf_field), dimension(n_vars), intent(out) :: lat_fields + type(esmf_field), dimension(2), intent(out) :: lon_fields + type(esmf_field), dimension(2), intent(out) :: lat_fields ! LOCAL integer :: id_var, ncid, ierr integer :: i,j,v character(len=500) :: fname - real(esmf_kind_r8) :: vec_lat(i_dim*j_dim), vec_lon(i_dim*j_dim) + real(esmf_kind_r8) :: vec1(i_dim*j_dim), vec4(4,i_dim*j_dim) real(esmf_kind_r8), allocatable :: array_lat(:,:,:), array_lon(:,:,:), array2d(:,:) + real(esmf_kind_r8) :: array_tmp(i_dim,j_dim) character(len=30) :: vname + character(len=15), dimension(2) :: lonvar_list + character(len=15), dimension(2) :: latvar_list + + ! center coord names + lonvar_list(1) = 'grid_center_lon' + latvar_list(1) = 'grid_center_lat' + + ! corner coord names + lonvar_list(2) = 'grid_corner_lon' + latvar_list(2) = 'grid_corner_lat' ! Create the fields - do v=1,n_vars + do v=1,2 vname="gauss_grid_" // trim(lonvar_list(v)) lon_fields(v) = ESMF_FieldCreate(gauss_grid, & typekind=ESMF_TYPEKIND_R8, & @@ -264,16 +272,16 @@ subroutine lonlat_read_into_fields(localpet, i_dim, j_dim, fname_read, dir_read, enddo if (localpet==0) then - allocate(array_lat(n_vars, i_dim, j_dim)) - allocate(array_lon(n_vars, i_dim, j_dim)) + allocate(array_lat(2, i_dim, j_dim)) + allocate(array_lon(2, i_dim, j_dim)) allocate(array2d( i_dim, j_dim)) else - allocate(array_lat(n_vars,0,0)) - allocate(array_lon(n_vars,0,0)) + allocate(array_lat(2,0,0)) + allocate(array_lon(2,0,0)) allocate(array2d(0,0)) endif - ! read from restart + ! read coordinates from restart if (localpet == 0) then fname = dir_read//"/"//fname_read @@ -283,42 +291,71 @@ subroutine lonlat_read_into_fields(localpet, i_dim, j_dim, fname_read, dir_read, ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) call netcdf_err(ierr, 'opening: '//trim(fname) ) - do v=1, n_vars + ! reading the center points - print *, 'Reading ', trim(lonvar_list(v)) - ierr=nf90_inq_varid(ncid, trim(lonvar_list(v)), id_var) - call netcdf_err(ierr, 'reading variable id' ) + print *, 'Reading ', trim(lonvar_list(1)) + ierr=nf90_inq_varid(ncid, trim(lonvar_list(1)), id_var) + call netcdf_err(ierr, 'reading variable id' ) - ierr=nf90_get_var(ncid, id_var, vec_lon) - call netcdf_err(ierr, 'reading variable' ) + ierr=nf90_get_var(ncid, id_var, vec1) + call netcdf_err(ierr, 'reading variable' ) - array_lon(v,:,:) = reshape(vec_lon,(/i_dim,j_dim/)) + array_lon(1,:,:) = reshape(vec1,(/i_dim,j_dim/)) - !do j = 1, j_dim - ! array_lon(v,:,j) =vec_lon - !enddo - print *, 'CSD lons', array_lon(v,1,1), array_lon(v,i_dim,1) - print *, 'CSD lons', array_lon(v,1,j_dim), array_lon(v,i_dim,j_dim) + print *, 'Reading ', trim(latvar_list(1)) + ierr=nf90_inq_varid(ncid, trim(latvar_list(1)), id_var) + call netcdf_err(ierr, 'reading variable id' ) - print *, 'Reading ', trim(latvar_list(v)) - ierr=nf90_inq_varid(ncid, trim(latvar_list(v)), id_var) - call netcdf_err(ierr, 'reading variable id' ) + ierr=nf90_get_var(ncid, id_var, vec1) + call netcdf_err(ierr, 'reading variable' ) - ierr=nf90_get_var(ncid, id_var, vec_lat) - call netcdf_err(ierr, 'reading variable' ) - ierr = nf90_close(ncid) - array_lat(v,:,:) = reshape(vec_lat,(/i_dim,j_dim/)) - !do i = 1, i_dim - ! array_lat(v,i,:) =vec_lat - !enddo - print *, 'CSD lats', array_lat(v,1,1), array_lat(v,i_dim,1) - print *, 'CSD lats', array_lat(v,1,j_dim), array_lat(v,i_dim,j_dim) - enddo + ! increment files are S->N, need to flip input lats from N->S + if ( grid_setup%descriptor == 'gau_inc') then + do j =1,j_dim + array_tmp = reshape(vec1,(/i_dim,j_dim/)) + array_lat(1,:,j) = array_tmp(:,j_dim+1-j) + enddo + else + array_lat(1,:,:) = reshape(vec1,(/i_dim,j_dim/)) + endif + + ! read in the corners + + print *, 'Reading ', trim(lonvar_list(2)) + ierr=nf90_inq_varid(ncid, trim(lonvar_list(2)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + + ierr=nf90_get_var(ncid, id_var, vec4) + call netcdf_err(ierr, 'reading variable' ) + + ! 2nd element of vec4 is the NW corner + array_lon(2,:,:) = reshape(vec4(2,:),(/i_dim,j_dim/)) + + print *, 'Reading ', trim(latvar_list(2)) + ierr=nf90_inq_varid(ncid, trim(latvar_list(2)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + + ierr=nf90_get_var(ncid, id_var, vec4) + call netcdf_err(ierr, 'reading variable' ) + + ! increment files are S->N, need to flip input lats from N->S + if ( grid_setup%descriptor == 'gau_inc') then + do j =1,j_dim + ! 2nd element of vec4 is the NW corner + array_tmp = reshape(vec4(2,:),(/i_dim,j_dim/)) + array_lat(2,:,j) = array_tmp(:,j_dim+1-j) + enddo + else + ! 2nd element of vec4 is the NW corner + array_lat(2,:,:) = reshape(vec4(2,:),(/i_dim,j_dim/)) + endif + ierr = nf90_close(ncid) endif - do v =1,n_vars + ! scatter the arrays into the fields + do v =1,2 ! scatter longitudes if (localpet==0) print *, 'scattering lon', lonvar_list(v) @@ -478,21 +515,15 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) type(esmf_polekind_flag) :: polekindflag(2) - type(esmf_field) :: lon_fields(1) - type(esmf_field) :: lat_fields(1) + type(esmf_field) :: lon_fields(2) + type(esmf_field) :: lat_fields(2) real(esmf_kind_r8), pointer :: lon_ptr_field(:,:), lon_ptr_coord(:,:) real(esmf_kind_r8), pointer :: lat_ptr_field(:,:), lat_ptr_coord(:,:) - character(len=15), dimension(1) :: lon_variables, lat_variables integer :: ierr, i_dim, j_dim polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE - ! Gaussian increment files lat/lon name - lon_variables(1) = 'grid_center_lon' - lat_variables(1) = 'grid_center_lat' - !lon_variables(1) = 'longitude' - !lat_variables(1) = 'latitude' i_dim = grid_setup%ires j_dim = grid_setup%jres @@ -512,8 +543,7 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) ! read lat lon coordinates into fields call lonlat_read_into_fields(localpet, i_dim, j_dim, & trim(grid_setup%fname_coord), trim(grid_setup%dir_coord), & - grid_setup, gauss_grid, 1, lon_variables, lat_variables, & - lon_fields, lat_fields) + grid_setup, gauss_grid, lon_fields, lat_fields) ! add coordinates to the grid call ESMF_GridAddCoord(gauss_grid, & From acd9943a22e9dec45f690cdbe2b4523d21ea0cf1 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Mon, 9 Dec 2024 21:02:23 +0000 Subject: [PATCH 11/15] Added soil-snow mask for gsi increment files. --- src/regridStates/grids_IO.F90 | 121 +++++++++++++++--------------- src/regridStates/regridStates.F90 | 39 ++++++---- 2 files changed, 86 insertions(+), 74 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 5a9b648..19b9b36 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -15,10 +15,14 @@ module grids_IO integer, public, parameter :: n_tiles=6 ! number tiles in fv3 grid integer, public, parameter :: vtype_water=0, & ! TO DO - which veg classification is this? vtype_landice=15 ! used for soil mask + ! mask values for soilsnow_mask calculated in the GSI EnKF + integer, public, parameter :: mtype_water=0, & + mtype_snow=2 type, public :: grid_setup_type character(7) :: descriptor character(100) :: fname character(100) :: dir + character(15) :: mask_variable(1) character(100) :: fname_mask character(100) :: dir_mask character(100) :: fname_coord @@ -36,14 +40,12 @@ module grids_IO !----------------------------------- ! Create ESMF grid objects, with mask if requested - subroutine setup_grid(localpet, npets, grid_setup, & - mask_type, mod_grid ) + subroutine setup_grid(localpet, npets, grid_setup, mod_grid ) implicit none ! INTENT IN type(grid_setup_type), intent(in) :: grid_setup - character(*), intent(in) :: mask_type integer, intent(in) :: localpet, npets ! INTENT OUT @@ -55,7 +57,6 @@ subroutine setup_grid(localpet, npets, grid_setup, & integer(esmf_kind_i4), pointer :: ptr_mask(:,:) integer :: ierr, ncid, tile - character(len=15) :: mask_variable(1) !-------------------------- ! Create grid object, and set up pet distribution @@ -63,10 +64,8 @@ subroutine setup_grid(localpet, npets, grid_setup, & select case (grid_setup%descriptor) case ('fv3_rst') call create_grid_fv3(grid_setup%ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) - mask_variable(1) = 'vtype ' case ('gau_inc') call create_grid_gauss(grid_setup, npets, localpet, mod_grid) - mask_variable(1) = 'soilsnow_mask ' case default call error_handler("unknown grid_setup%descriptor in setup_grid", 1) end select @@ -74,56 +73,59 @@ subroutine setup_grid(localpet, npets, grid_setup, & !-------------------------- ! Calculate and add the mask - if (mask_type=="soil") then ! mask out ocean and glaciers, using vegetation class - mask_field(1) = ESMF_FieldCreate(mod_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input variable for mask", & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate, mask_variable", ierr) - - call read_into_fields(localpet, grid_setup%ires, grid_setup%jres, trim(grid_setup%fname_mask), & - trim(grid_setup%dir_mask), grid_setup, 1, & - mask_variable(1), mask_field(1)) - - ! get pointer to vegtype - call ESMF_FieldGet(mask_field(1), & - farrayPtr=ptr_maskvar, & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldGet", ierr) - - ! create and get pointer to the mask - call ESMF_GridAddItem(mod_grid, & - itemflag=ESMF_GRIDITEM_MASK, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("in GridAddItem mask", ierr) - - call ESMF_GridGetItem(mod_grid, & - itemflag=ESMF_GRIDITEM_MASK, & - farrayPtr=ptr_mask, & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("in GridGetItem mask", ierr) - - ! calculate the mask - ptr_mask = 1 ! initialize land everywhere - if (mask_variable(1) == 'vtype ') then - where (nint(ptr_maskvar) == vtype_water ) ptr_mask = 0 ! exclude water - where (nint(ptr_maskvar) == vtype_landice ) ptr_mask = 0 ! exclude glaciers - end if - - ! destroy veg type field - call ESMF_FieldDestroy(mask_field(1),rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("DESTROYING FIELD", ierr) - - end if ! mask = soil + mask_field(1) = ESMF_FieldCreate(mod_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input variable for mask", & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate, mask_variable", ierr) + + call read_into_fields(localpet, grid_setup%ires, grid_setup%jres, trim(grid_setup%fname_mask), & + trim(grid_setup%dir_mask), grid_setup, 1, & + grid_setup%mask_variable(1), mask_field(1)) + +! get pointer to mask + call ESMF_FieldGet(mask_field(1), & + farrayPtr=ptr_maskvar, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", ierr) + +! create and get pointer to the mask + call ESMF_GridAddItem(mod_grid, & + itemflag=ESMF_GRIDITEM_MASK, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("in GridAddItem mask", ierr) + + call ESMF_GridGetItem(mod_grid, & + itemflag=ESMF_GRIDITEM_MASK, & + farrayPtr=ptr_mask, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("in GridGetItem mask", ierr) + +! calculate the mask + ptr_mask = 1 ! initialize land everywhere + select case (trim(grid_setup%mask_variable(1))) + case("vtype") ! removing non-land and glaciers using veg class + where (nint(ptr_maskvar) == vtype_water ) ptr_mask = 0 ! exclude water + where (nint(ptr_maskvar) == vtype_landice ) ptr_mask = 0 ! exclude glaciers + case("soilsnow_mask") ! removing snow and non-land using pre-computed mask + where (nint(ptr_maskvar) == mtype_water ) ptr_mask = 0 ! exclude non-soil + where (nint(ptr_maskvar) == mtype_snow ) ptr_mask = 0 ! exclude snow + case default + call error_handler("unknown mask_variable", 1) + end select + +! destroy veg type field + call ESMF_FieldDestroy(mask_field(1),rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("DESTROYING FIELD", ierr) end subroutine setup_grid @@ -181,13 +183,13 @@ subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, & fname = dir_read//"/"//fname_read endif - !print *, 'Reading ', trim(fname) + print *, 'Reading ', trim(fname) ierr=nf90_open(trim(fname),NF90_NOWRITE,ncid) call netcdf_err(ierr, 'opening: '//trim(fname) ) do v =1, n_vars - !print *, 'Reading ', trim(variable_list(v)) + print *, 'Reading ', trim(variable_list(v)) ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var) call netcdf_err(ierr, 'reading variable id' ) @@ -555,8 +557,7 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) nullify(lon_ptr_field) ! get pointer to lon coord - !if (localpet == 0) print*," getting GridCoord for long" - print*," getting GridCoord for long", localpet + if (localpet == 0) print*," getting GridCoord for long" call ESMF_GridGetCoord(gauss_grid, & staggerLoc=ESMF_STAGGERLOC_CENTER, & coordDim=1, & @@ -592,7 +593,7 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet latitude", ierr) - ! set lat coord to fiel values + ! set lat coord to field values lat_ptr_coord = lat_ptr_field ! assign the corners @@ -637,7 +638,7 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) ! enddo !enddo - ! TO DO destrot the fields + ! TO DO destroy the fields end subroutine create_grid_gauss diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 2976a97..6ad1916 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -25,7 +25,6 @@ program regridStates character(len=15) :: variable_list(max_vars) integer :: n_vars real(esmf_kind_r8) :: missing_value ! value given to unmapped cells in the output grid - character(len=10) :: mask_type type(grid_setup_type) :: grid_setup_in, grid_setup_out @@ -44,7 +43,7 @@ program regridStates real :: t1, t2, t3, t4 ! see README for details of namelist variables. - namelist /config/ n_vars, variable_list, missing_value, mask_type + namelist /config/ n_vars, variable_list, missing_value ! INITIALIZE !------------------------------------------------------------------------- @@ -97,11 +96,9 @@ program regridStates ! TO DO - can we make the number of tasks more flexible for fv3 if (localpet==0) print*,'** Setting up grids' - call setup_grid(localpet, npets, grid_setup_in, & - trim(mask_type), grid_in ) + call setup_grid(localpet, npets, grid_setup_in, grid_in ) - call setup_grid(localpet, npets, grid_setup_out, & - trim(mask_type), grid_out ) + call setup_grid(localpet, npets, grid_setup_out, grid_out ) !------------------------ ! Create input and output fields @@ -272,12 +269,14 @@ subroutine readin_setup(unt,namel,grid_setup) integer :: ires, jres integer :: ierr - namelist /input/ fname, dir, gridtype, & + namelist /input/ fname, dir, & + gridtype, & fname_mask, dir_mask, & fname_coord, dir_coord, & ires, jres - namelist /output/ fname, dir, gridtype, & + namelist /output/ fname, dir, & + gridtype, & fname_mask, dir_mask, & fname_coord, dir_coord,& ires, jres @@ -292,16 +291,20 @@ subroutine readin_setup(unt,namel,grid_setup) ires = 0 jres = 0 - if (namel=="input") then + select case (namel) + case ("input") read(unt, nml=input, iostat=ierr) if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) - elseif (namel=="output") then + case ("output") read(unt, nml=output, iostat=ierr) if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) - endif + case default + call error_handler("unknown namel in readin_setup", 1) + end select grid_setup%descriptor = gridtype - grid_setup%dir = dir ! for fv3 and gauss use input file for mask + + grid_setup%dir = dir grid_setup%fname = fname ! to-do, add routine to check if present. @@ -312,10 +315,18 @@ subroutine readin_setup(unt,namel,grid_setup) case ("output") grid_setup%dir_mask = dir_mask ! need a file on output grid for mask grid_setup%fname_mask = fname_mask - case default - call error_handler("unknown namel in readin_setup", 1) end select + ! set-up mask details, based on file type + select case (gridtype) + case ("fv3_rst") ! for history file and restarts, use veg type + grid_setup%mask_variable(1) = "vtype " + case ("gau_inc") ! gsi-output incr files only, use calculated mask + grid_setup%mask_variable(1) = "soilsnow_mask " + case default + call error_handler("unknown gridtype in readin_setup", 1) + end select + grid_setup%dir_coord = dir_coord grid_setup%fname_coord = fname_coord grid_setup%ires = ires From 7a01a1656255d153400430dc858202f3429e9911 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 10 Dec 2024 19:27:47 +0000 Subject: [PATCH 12/15] Destroyed all fields, added corner locations to gauss grid. --- src/regridStates/grids_IO.F90 | 154 ++++++++++++++-------------------- 1 file changed, 62 insertions(+), 92 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 19b9b36..b22b69a 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -516,13 +516,14 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) type(esmf_grid) :: gauss_grid type(esmf_polekind_flag) :: polekindflag(2) + type(esmf_staggerloc) :: staggerloc(2) type(esmf_field) :: lon_fields(2) type(esmf_field) :: lat_fields(2) real(esmf_kind_r8), pointer :: lon_ptr_field(:,:), lon_ptr_coord(:,:) real(esmf_kind_r8), pointer :: lat_ptr_field(:,:), lat_ptr_coord(:,:) - integer :: ierr, i_dim, j_dim + integer :: ierr, i_dim, j_dim, v polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE @@ -546,100 +547,69 @@ subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid) call lonlat_read_into_fields(localpet, i_dim, j_dim, & trim(grid_setup%fname_coord), trim(grid_setup%dir_coord), & grid_setup, gauss_grid, lon_fields, lat_fields) - - ! add coordinates to the grid - call ESMF_GridAddCoord(gauss_grid, & - staggerloc=ESMF_STAGGERLOC_CENTER, rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridAddCoord", ierr) - - nullify(lon_ptr_coord) - nullify(lon_ptr_field) - - ! get pointer to lon coord - if (localpet == 0) print*," getting GridCoord for long" - call ESMF_GridGetCoord(gauss_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=1, & - farrayPtr=lon_ptr_coord, rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord longitude ", ierr) - - ! fetch lon from field into pointer - call ESMF_FieldGet(lon_fields(1), & - farrayPtr=lon_ptr_field, & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldGet longitude", ierr) - - ! set coord pointe to field pointer - lon_ptr_coord = lon_ptr_field - - nullify(lat_ptr_coord) - nullify(lat_ptr_field) - ! get pointer to lat coord. Need bounds? - call ESMF_GridGetCoord(gauss_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=2, & - farrayPtr=lat_ptr_coord, rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord latitude", ierr) + staggerloc = (/ ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER /) + do v = 1,2 + ! add coordinates to the grid + call ESMF_GridAddCoord(gauss_grid, & + staggerloc=staggerloc(v), rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", ierr) + + nullify(lon_ptr_coord) + nullify(lon_ptr_field) + + ! get pointer to lon coord + if (localpet == 0) print*," getting GridCoord for long" + call ESMF_GridGetCoord(gauss_grid, & + staggerLoc=staggerloc(v), & + coordDim=1, & + farrayPtr=lon_ptr_coord, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord longitude ", ierr) + + ! fetch lon from field into pointer + call ESMF_FieldGet(lon_fields(v), & + farrayPtr=lon_ptr_field, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet longitude", ierr) + + ! set coord pointe to field pointer + lon_ptr_coord = lon_ptr_field - ! fetch lat from field into pointer - call ESMF_FieldGet(lat_fields(1), & - farrayPtr=lat_ptr_field, & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldGet latitude", ierr) - - ! set lat coord to field values - lat_ptr_coord = lat_ptr_field - - ! assign the corners - ! call ESMF_GridAddCoord(input_grid, & - ! staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) - !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - ! call error_handler("IN GridAddCoord", rc) - - !print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." - !nullify(lon_corner_src_ptr) - !call ESMF_GridGetCoord(input_grid, & - ! staggerLoc=ESMF_STAGGERLOC_CORNER, & - ! coordDim=1, & - ! farrayPtr=lon_corner_src_ptr, rc=rc) - !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - ! call error_handler("IN GridGetCoord", rc) - - !print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." - !nullify(lat_corner_src_ptr) - !call ESMF_GridGetCoord(input_grid, & - ! staggerLoc=ESMF_STAGGERLOC_CORNER, & - ! coordDim=2, & - ! computationalLBound=clb, & - ! computationalUBound=cub, & - ! farrayPtr=lat_corner_src_ptr, rc=rc) - !if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - ! call error_handler("IN GridGetCoord", rc) - - !do j = clb(2), cub(2) - ! do i = clb(1), cub(1) - ! lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon) - ! if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8 - ! if (j == 1) then - ! lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8 - ! cycle - ! endif - ! if (j == jp1_input) then - ! lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8 - ! cycle - ! endif - ! lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j)) - ! enddo - !enddo - - ! TO DO destroy the fields + nullify(lat_ptr_coord) + nullify(lat_ptr_field) + + ! get pointer to lat coord. Need bounds? + call ESMF_GridGetCoord(gauss_grid, & + staggerLoc=staggerloc(v), & + coordDim=2, & + farrayPtr=lat_ptr_coord, rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord latitude", ierr) + + ! fetch lat from field into pointer + call ESMF_FieldGet(lat_fields(v), & + farrayPtr=lat_ptr_field, & + rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet latitude", ierr) + + ! set lat coord to field values + lat_ptr_coord = lat_ptr_field + enddo + ! clean up + do v = 1,2 + call ESMF_FieldDestroy(lat_fields(v),rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("DESTROYING FIELD", ierr) + call ESMF_FieldDestroy(lon_fields(v),rc=ierr) + if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("DESTROYING FIELD", ierr) + enddo + end subroutine create_grid_gauss end module grids_IO From f81758f874c20b2121cbd7618dc58ba8707eb4f2 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Sun, 15 Dec 2024 23:04:10 +0000 Subject: [PATCH 13/15] Switched to reading in vegtype from fix files instead of restart. --- src/regridStates/grids_IO.F90 | 2 +- src/regridStates/regridStates.F90 | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index b22b69a..872872f 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -112,7 +112,7 @@ subroutine setup_grid(localpet, npets, grid_setup, mod_grid ) ! calculate the mask ptr_mask = 1 ! initialize land everywhere select case (trim(grid_setup%mask_variable(1))) - case("vtype") ! removing non-land and glaciers using veg class + case("vegetation_type") ! removing non-land and glaciers using veg class where (nint(ptr_maskvar) == vtype_water ) ptr_mask = 0 ! exclude water where (nint(ptr_maskvar) == vtype_landice ) ptr_mask = 0 ! exclude glaciers case("soilsnow_mask") ! removing snow and non-land using pre-computed mask diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 6ad1916..6f8ece1 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -320,7 +320,8 @@ subroutine readin_setup(unt,namel,grid_setup) ! set-up mask details, based on file type select case (gridtype) case ("fv3_rst") ! for history file and restarts, use veg type - grid_setup%mask_variable(1) = "vtype " + !grid_setup%mask_variable(1) = "vtype " ! if getting from a restart + grid_setup%mask_variable(1) = "vegetation_type" ! if getting from fix file case ("gau_inc") ! gsi-output incr files only, use calculated mask grid_setup%mask_variable(1) = "soilsnow_mask " case default From 4408d9ea48d17625c03de72926d4ee881549e307 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 17 Dec 2024 21:04:06 +0000 Subject: [PATCH 14/15] Tidy up README, update example namelist. --- src/regridStates/README | 32 +++++++++++------------------ src/regridStates/gauIncr2native.nml | 31 ++++++++++++++++++++++++++++ src/regridStates/grids_IO.F90 | 4 +--- src/regridStates/tile2tile.nml | 14 ------------- 4 files changed, 44 insertions(+), 37 deletions(-) create mode 100644 src/regridStates/gauIncr2native.nml delete mode 100644 src/regridStates/tile2tile.nml diff --git a/src/regridStates/README b/src/regridStates/README index 993551e..3f5ae28 100644 --- a/src/regridStates/README +++ b/src/regridStates/README @@ -1,28 +1,20 @@ This program is designed to convert individual variables between resolutions, for DA applications (specificially re-gridding for use in re-centering, and regridding increments). -Clara Draper, Aug 13, 2024. +Clara Draper, Aug 17, 2024. * Uses bi-linear interpolation, with masking of the input and output. -* Only mask option right now is "soil" (no water, no land-ice). Could add others easily. -* Does not guarantee all output grid cells with have mapped values, by design. +* Current masking options are all soil (land, no glaciers) and snow-free soil (input increment only) +* Does not guarantee all output grid cells with have mapped values, by design. CANNOT BY USED FOR RESTARTS (could be changed though). * Output grid cells that do not have a value after the interpolation step are filled with the nearest neighbouring output grid cell. If there are no mapped neighbours within two layers of surrounding neighbours, the output grid cell will remain unmapped. In this way, islands do not get filled with distant values. -* Currently, only reads in / writes out fv3 files (TO-DO, add Gaussian grid option for GSI increments). +* Reads in / writes out fv3 files and Gaussian grid (gsi output increment) files. +* Handling Gaussian files requires a scrip file with the grid details. This can be generated using the ufs_utils weigh_gen program. +* has separate namelist for input and output. Required options depend on the grid type and whether it's in/out. Look in readin_setup routine to check what is needed. +* example nml provided for input Gaussian increment -> fv3 grid. + +FEATURES TO ADD: * Can only process 2D variables (no vertical dimension; TO-DO fix this) +* Add time dimension to put all gsi increments in a single file (for reading in by IAU) +* Re-think parallelization so can do multiple ensemble members at once +* Add code to check required options are available for requested set-up grid -Options are controlled by the namelist: - dir_fix -> base directory for the orog files (one layer below C$RES dir) - res_atm_in -> atmos input resolution - res_atm_out -> atmos output resolution - fname_in -> start of the input restart file name (everything before ".tile?.nc" ) - dir_in -> directory where input restarts are - fname_out -> start of the output restart file name (everything before ".tile?.nc" ) - dir_out -> directory where output restarts will go -! the next two are needed to read in veg_type for calculating the output mask - fname_out_rst -> start of file name for example restart at output res (everything before ".tile?.nc" ) - dir_out_rst -> directory where example restarts at output res are - variable_list(1) -> list of variables names to remap. Right-pad with spaces to length of 10 - variable_list(2) -> - n_vars -> number of variables to remap (current max. set to 10). - missing_value -> value used for grid cells that remain unmapped - mask_type -> only options for now are "none", and "soil" (no water, no land-ice). diff --git a/src/regridStates/gauIncr2native.nml b/src/regridStates/gauIncr2native.nml new file mode 100644 index 0000000..b1ebeba --- /dev/null +++ b/src/regridStates/gauIncr2native.nml @@ -0,0 +1,31 @@ +&config + n_vars=4, ! number of vars in list below + variable_list(1)="soilt1_inc ", + variable_list(2)="soilt2_inc ", + variable_list(3)="slc1_inc ", + variable_list(4)="slc2_inc ", + missing_value=0., ! any un-filled grid cell takes this value. + ! use 0. for increments +/ + +! required fields depend on gridtype, and whether input or output +&input + gridtype="gau_inc", ! "gau_inc" for GSI/Gaussian increment file + ires=192, + jres=96, + fname="enkfgdas.t12z.sfci003.ensmean.nc", + dir="./input/", + fname_coord="gaussian.192.96.nc", ! this is the scrip file generated by ufs_utils weight_gen + dir_coord="./input/" +/ + +&output + gridtype="fv3_rst", ! "fv3_rst" any native model grid file (increment or restart) + ires=48, + jres=48, + fname="inc2ens.sfc_data", + dir="./output/", + fname_mask="20211221.090000.sfc_data" + dir_mask="./output_rst/", + dir_coord="./fixdir/", ! orog dir in fix +/ diff --git a/src/regridStates/grids_IO.F90 b/src/regridStates/grids_IO.F90 index 872872f..d6ccca5 100644 --- a/src/regridStates/grids_IO.F90 +++ b/src/regridStates/grids_IO.F90 @@ -73,8 +73,6 @@ subroutine setup_grid(localpet, npets, grid_setup, mod_grid ) !-------------------------- ! Calculate and add the mask - ! mask out ocean and glaciers, using vegetation class - mask_field(1) = ESMF_FieldCreate(mod_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & @@ -122,7 +120,7 @@ subroutine setup_grid(localpet, npets, grid_setup, mod_grid ) call error_handler("unknown mask_variable", 1) end select -! destroy veg type field +! destroy mask field call ESMF_FieldDestroy(mask_field(1),rc=ierr) if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("DESTROYING FIELD", ierr) diff --git a/src/regridStates/tile2tile.nml b/src/regridStates/tile2tile.nml deleted file mode 100644 index 8a8c041..0000000 --- a/src/regridStates/tile2tile.nml +++ /dev/null @@ -1,14 +0,0 @@ -&config - dir_fix="/scratch1/NCEPDEV/global/glopara/fix/orog/20231027/", - res_atm_in=768, - res_atm_out=384, - fname_in="20240120.030000.sfc_data", - dir_in="./input/", - fname_out="opsDet2Ens.sfc_data", - dir_out="./output/", - fname_out_rst="ensmean.20240120.030000.sfc_data", - dir_out_rst="./output_rst/", - variable_list(1)="snwdph ", - n_vars=1, - missing_value=-999.9 -/ From e4a4f88df43ab43d55a69355f046daab316715da Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 18 Dec 2024 17:32:57 +0000 Subject: [PATCH 15/15] Fixes for Tseganeh's review. --- README | 1 - src/regridStates/README | 2 +- src/regridStates/gauIncr2native.nml | 32 +++++++++++++++++++++-------- src/regridStates/regridStates.F90 | 2 +- 4 files changed, 25 insertions(+), 12 deletions(-) diff --git a/README b/README index 07a5b55..ec7726f 100644 --- a/README +++ b/README @@ -2,7 +2,6 @@ To compile regridding tool in stand-alone mode: load GDASApp modules - load ESMF module (module load esmf/8.6.0) > cd DA-utils > mkdir build > ecbuild ../bundle diff --git a/src/regridStates/README b/src/regridStates/README index 3f5ae28..a4f16b6 100644 --- a/src/regridStates/README +++ b/src/regridStates/README @@ -4,7 +4,7 @@ Clara Draper, Aug 17, 2024. * Uses bi-linear interpolation, with masking of the input and output. * Current masking options are all soil (land, no glaciers) and snow-free soil (input increment only) -* Does not guarantee all output grid cells with have mapped values, by design. CANNOT BY USED FOR RESTARTS (could be changed though). +* Does not guarantee all output grid cells will have mapped values, by design. CANNOT BY USED FOR RESTARTS (could be changed though). * Output grid cells that do not have a value after the interpolation step are filled with the nearest neighbouring output grid cell. If there are no mapped neighbours within two layers of surrounding neighbours, the output grid cell will remain unmapped. In this way, islands do not get filled with distant values. * Reads in / writes out fv3 files and Gaussian grid (gsi output increment) files. * Handling Gaussian files requires a scrip file with the grid details. This can be generated using the ufs_utils weigh_gen program. diff --git a/src/regridStates/gauIncr2native.nml b/src/regridStates/gauIncr2native.nml index b1ebeba..470457b 100644 --- a/src/regridStates/gauIncr2native.nml +++ b/src/regridStates/gauIncr2native.nml @@ -1,31 +1,45 @@ &config n_vars=4, ! number of vars in list below - variable_list(1)="soilt1_inc ", - variable_list(2)="soilt2_inc ", - variable_list(3)="slc1_inc ", - variable_list(4)="slc2_inc ", + variable_list="soilt1_inc", "soilt2_inc", "slc1_inc", "slc2_inc ", ! vars to be regridded missing_value=0., ! any un-filled grid cell takes this value. ! use 0. for increments / -! required fields depend on gridtype, and whether input or output +! required fields depend on gridtype, and whether input or output. See below for description. &input - gridtype="gau_inc", ! "gau_inc" for GSI/Gaussian increment file + gridtype="gau_inc", ires=192, jres=96, fname="enkfgdas.t12z.sfci003.ensmean.nc", dir="./input/", - fname_coord="gaussian.192.96.nc", ! this is the scrip file generated by ufs_utils weight_gen + fname_coord="gaussian.192.96.nc", dir_coord="./input/" / &output - gridtype="fv3_rst", ! "fv3_rst" any native model grid file (increment or restart) + gridtype="fv3_rst", ires=48, jres=48, fname="inc2ens.sfc_data", dir="./output/", - fname_mask="20211221.090000.sfc_data" + fname_mask="vegetation_type" dir_mask="./output_rst/", dir_coord="./fixdir/", ! orog dir in fix / + +!input and output namelist options: +! Note: designed for files to be linked prior to executing (arguments aren't long enough for full paths) +! +!fname ! filename to be interpolated to/from +!dir ! directory of above file +!gridtype ! "gau_inc" for GSI/Gaussian increment file , +! ! "fv3_rst" any native model grid file (increment or restart) +!fname_mask ! filename with mask information. Only needed for fv3 files. Set to the +! ! vegetation_type files in ${FIXorog}/${CRES}/sfc/ +!dir_mask ! directory of above file +!fname_coord ! File name with lat/lon info. Only needed for Gaussian files. Set to the +! ! scrip file generated by ufs_utils weight_gen +!dir_coord ! Gaussian - directory with above file +! ! FV3 - set to $FIXorog +!ires ! longitudinal dimension +!jres ! latitudinal dimension diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index 6f8ece1..c3c8589 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -297,7 +297,7 @@ subroutine readin_setup(unt,namel,grid_setup) if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) case ("output") read(unt, nml=output, iostat=ierr) - if (ierr /= 0) call error_handler("READING input NAMELIST.", ierr) + if (ierr /= 0) call error_handler("READING output NAMELIST.", ierr) case default call error_handler("unknown namel in readin_setup", 1) end select