diff --git a/README b/README index ac0089d..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 essmf/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/README b/src/regridStates/README index 993551e..a4f16b6 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 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. -* 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/fv3_grid.F90 b/src/regridStates/fv3_grid.F90 deleted file mode 100644 index 15699f6..0000000 --- a/src/regridStates/fv3_grid.F90 +++ /dev/null @@ -1,284 +0,0 @@ - module fv3_grid - ! ESMF grid-specific routines for GFS regridding program, including IO. - ! - ! Clara Draper, Aug 2024. - - use esmf - use netcdf - use ESMF_LogPublicMod - use utilities, only : error_handler, netcdf_err - - implicit none - - private - - 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 - - public :: setup_grid, & - read_into_fields, & - write_from_fields - - contains - -!----------------------------------- -! 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) - - implicit none - - ! INTENT IN - character(*), intent(in) :: dir_fix - integer, intent(in) :: res_atm - character(*), intent(in) :: fname_rst, dir_rst - character(*), intent(in) :: mask_type - integer, intent(in) :: localpet, npets - - ! INTENT OUT - type(esmf_grid) :: fv3_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) - 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)//"/" - - 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) - -!-------------------------- -! Calcalate and add the mask - - if (mask_type=="soil") then ! mask out ocean and glaciers - vtype_field(1) = ESMF_FieldCreate(fv3_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input veg type for mask", & - rc=ierr) - if(ESMF_logFoundError(rcToCheck=ierr,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate, vtype", 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)) - - ! get pointer to vegtype - call ESMF_FieldGet(vtype_field(1), & - farrayPtr=ptr_vtype, & - 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(fv3_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, & - 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 - where (nint(ptr_vtype) == vtype_water ) ptr_mask = 0 ! exclude water - where (nint(ptr_vtype) == vtype_landice ) ptr_mask = 0 ! exclude glaciers - - ! destroy veg type field - call ESMF_FieldDestroy(vtype_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 - - 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) - - implicit none - - ! INTENT IN - integer, intent(in) :: localpet, i_dim, j_dim, n_vars - character(*), intent(in) :: fname_read - character(*), intent(in) :: dir_read - - character(len=10), dimension(n_vars), intent(in) :: variable_list - - ! INTENT OUT - type(esmf_field), dimension(n_vars), intent(inout) :: fields - - ! LOCAL - integer :: tt, id_var, ncid, ierr, v - character(len=1) :: tchar - character(len=500) :: fname - real(esmf_kind_r8), allocatable :: array2D(:,:) - real(esmf_kind_r8), allocatable :: array_in(:,:,:) - - if (localpet==0) then - allocate(array_in(n_vars,i_dim, j_dim)) - allocate(array2D(i_dim, j_dim)) - else - allocate(array_in(0,0,0)) - allocate(array2D(0,0)) - end if - - do tt = 1, n_tiles - - ! read from restart - if (localpet == 0) then - - write(tchar,'(i1)') tt - fname = dir_read//"/"//fname_read//".tile"//tchar//".nc" - - 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 - ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var) - call netcdf_err(ierr, 'reading variable id' ) - - ierr=nf90_get_var(ncid, id_var, array_in(v,:,:)) - call netcdf_err(ierr, 'reading variable' ) - enddo - ierr = nf90_close(ncid) - - endif - - ! scatter - do v =1, n_vars - array2D=array_in(v,:,:) ! scatter misbehaves if given indexed 3D array. - 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 - - ! clean up - deallocate(array_in) - - 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) - - implicit none - - ! INTENT IN - 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 - type(esmf_field), dimension(n_vars), intent(in) :: fields - - ! LOCAL - integer :: tt, id_var, ncid, ierr, & - id_x, id_y, v - character(len=1) :: tchar - character(len=500) :: fname - real(esmf_kind_r8), allocatable :: array2D(:,:) - real(esmf_kind_r8), allocatable :: array_out(:,:,:) - - do v = 1, n_vars - if (localpet == 0) print *, 'Writing ', trim(variable_list(v)), ' into field' - enddo - - if (localpet==0) then - allocate(array_out(n_vars, i_dim, j_dim)) - allocate(array2D(i_dim, j_dim)) - else - allocate(array_out(0,0,0)) - allocate(array2D(0,0)) - end if - - do tt = 1, n_tiles - - ! fetch data (all PETs) - do v=1, n_vars - call ESMF_FieldGather(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 FieldGather", ierr) - array_out(v,:,:) = array2D - enddo - - ! write to netcdf - if (localpet == 0) then - - ! open file, set dimensions - write(tchar,'(i1)') tt - fname = dir_out//"/"//fname_out//".tile"//tchar//".nc" - - ierr = nf90_create(trim(fname), NF90_NETCDF4, ncid) - call netcdf_err(ierr, 'creating file='//trim(fname) ) - - ierr = nf90_def_dim(ncid, 'xaxis_1', i_dim, id_x) - call netcdf_err(ierr, 'defining xaxis dimension' ) - - ierr = nf90_def_dim(ncid, 'yaxis_1', j_dim, id_y) - call netcdf_err(ierr, 'defining yaxis dimension' ) - - do v=1, n_vars - - ierr = nf90_def_var(ncid, trim(variable_list(v)), NF90_DOUBLE, (/id_x, id_y/), id_var) - call netcdf_err(ierr, 'defining '//variable_list(v) ) - - ierr = nf90_put_var( ncid, id_var, array_out(v,:,:) ) - call netcdf_err(ierr, 'writing '//variable_list(v) ) - - enddo - - ierr = nf90_close(ncid) - - endif - - enddo - - ! clean up - deallocate(array_out) - - end subroutine write_from_fields - -end module fv3_grid diff --git a/src/regridStates/gauIncr2native.nml b/src/regridStates/gauIncr2native.nml new file mode 100644 index 0000000..470457b --- /dev/null +++ b/src/regridStates/gauIncr2native.nml @@ -0,0 +1,45 @@ +&config + n_vars=4, ! number of vars in list below + 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. See below for description. +&input + gridtype="gau_inc", + ires=192, + jres=96, + fname="enkfgdas.t12z.sfci003.ensmean.nc", + dir="./input/", + fname_coord="gaussian.192.96.nc", + dir_coord="./input/" +/ + +&output + gridtype="fv3_rst", + ires=48, + jres=48, + fname="inc2ens.sfc_data", + dir="./output/", + 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/grids_IO.F90 b/src/regridStates/grids_IO.F90 new file mode 100644 index 0000000..d6ccca5 --- /dev/null +++ b/src/regridStates/grids_IO.F90 @@ -0,0 +1,613 @@ + module grids_IO + ! ESMF grid-specific routines for GFS regridding program, including IO. + ! + ! Clara Draper, Aug 2024. + + use esmf + use netcdf + use ESMF_LogPublicMod + use utilities, only : error_handler, netcdf_err + + implicit none + + private + + 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 + character(100) :: dir_coord + integer :: ires + integer :: jres + end type + + public :: setup_grid, & + read_into_fields, & + write_from_fields + + contains + +!----------------------------------- +! Create ESMF grid objects, with mask if requested + + subroutine setup_grid(localpet, npets, grid_setup, mod_grid ) + + implicit none + + ! INTENT IN + type(grid_setup_type), intent(in) :: grid_setup + integer, intent(in) :: localpet, npets + + ! INTENT OUT + type(esmf_grid) :: mod_grid + + ! LOCAL + type(esmf_field) :: mask_field(1) + real(esmf_kind_r8), pointer :: ptr_maskvar(:,:) + integer(esmf_kind_i4), pointer :: ptr_mask(:,:) + + integer :: ierr, ncid, tile + +!-------------------------- +! Create grid object, and set up pet distribution + + select case (grid_setup%descriptor) + case ('fv3_rst') + call create_grid_fv3(grid_setup%ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid) + case ('gau_inc') + call create_grid_gauss(grid_setup, npets, localpet, mod_grid) + case default + call error_handler("unknown grid_setup%descriptor in setup_grid", 1) + end select + +!-------------------------- +! Calculate and add the mask + + 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("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 + 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 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) + + 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, dir_read, & + grid_setup, n_vars, variable_list, fields) + + implicit none + + ! INTENT IN + 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 + + character(len=15), dimension(n_vars), intent(in) :: variable_list + + ! INTENT OUT + type(esmf_field), dimension(n_vars), intent(inout) :: fields + + ! 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(:,:) + real(esmf_kind_r8), allocatable :: array_in(:,:,:) + + if (localpet==0) then + allocate(array_in(n_vars,i_dim, j_dim)) + allocate(array2D(i_dim, j_dim)) + else + allocate(array_in(0,0,0)) + allocate(array2D(0,0)) + end if + + 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_setup%descriptor in read into fields", 1) + end select + + do tt = 1, n_files + + ! read from restart + if (localpet == 0) then + + 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)) + ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var) + call netcdf_err(ierr, 'reading variable id' ) + + ierr=nf90_get_var(ncid, id_var, array_in(v,:,:)) + call netcdf_err(ierr, 'reading variable' ) + enddo + ierr = nf90_close(ncid) + + endif + + ! scatter + do v =1, n_vars + array2D=array_in(v,:,:) ! scatter misbehaves if given indexed 3D array. + 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 + + ! clean up + deallocate(array_in) + + end subroutine read_into_fields + + ! 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, lon_fields, lat_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 + type(esmf_grid), intent(in) :: gauss_grid + + + ! INTENT OUT + 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) :: 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,2 + 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(2, i_dim, j_dim)) + allocate(array_lon(2, i_dim, j_dim)) + allocate(array2d( i_dim, j_dim)) + else + allocate(array_lat(2,0,0)) + allocate(array_lon(2,0,0)) + allocate(array2d(0,0)) + endif + + ! read coordinates 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) ) + + ! reading the center points + + 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, vec1) + call netcdf_err(ierr, 'reading variable' ) + + array_lon(1,:,:) = reshape(vec1,(/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' ) + + ierr=nf90_get_var(ncid, id_var, vec1) + 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 + 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 + + ! scatter the arrays into the fields + do v =1,2 + + ! 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 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, & + n_vars, variable_list, fields) + + implicit none + + ! INTENT IN + integer, intent(in) :: localpet, i_dim, j_dim, n_vars + character(*), intent(in) :: fname_out + character(*), intent(in) :: dir_out + character(15), dimension(n_vars), intent(in) :: variable_list + type(esmf_field), dimension(n_vars), intent(in) :: fields + + ! LOCAL + integer :: tt, id_var, ncid, ierr, & + id_x, id_y, v + character(len=1) :: tchar + character(len=500) :: fname + real(esmf_kind_r8), allocatable :: array2D(:,:) + real(esmf_kind_r8), allocatable :: array_out(:,:,:) + + do v = 1, n_vars + if (localpet == 0) print *, 'Writing ', trim(variable_list(v)), ' into field' + enddo + + if (localpet==0) then + allocate(array_out(n_vars, i_dim, j_dim)) + allocate(array2D(i_dim, j_dim)) + else + allocate(array_out(0,0,0)) + allocate(array2D(0,0)) + end if + + do tt = 1, n_tiles + + ! fetch data (all PETs) + do v=1, n_vars + call ESMF_FieldGather(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 FieldGather", ierr) + array_out(v,:,:) = array2D + enddo + + ! write to netcdf + if (localpet == 0) then + + ! open file, set dimensions + write(tchar,'(i1)') tt + fname = dir_out//"/"//fname_out//".tile"//tchar//".nc" + + ierr = nf90_create(trim(fname), NF90_NETCDF4, ncid) + call netcdf_err(ierr, 'creating file='//trim(fname) ) + + ierr = nf90_def_dim(ncid, 'xaxis_1', i_dim, id_x) + call netcdf_err(ierr, 'defining xaxis dimension' ) + + ierr = nf90_def_dim(ncid, 'yaxis_1', j_dim, id_y) + call netcdf_err(ierr, 'defining yaxis dimension' ) + + do v=1, n_vars + + ierr = nf90_def_var(ncid, trim(variable_list(v)), NF90_DOUBLE, (/id_x, id_y/), id_var) + call netcdf_err(ierr, 'defining '//variable_list(v) ) + + ierr = nf90_put_var( ncid, id_var, array_out(v,:,:) ) + call netcdf_err(ierr, 'writing '//variable_list(v) ) + + enddo + + ierr = nf90_close(ncid) + + endif + + enddo + + ! clean up + deallocate(array_out) + + 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=200) :: 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)') 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 + 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(grid_setup, npets, localpet, gauss_grid) + + ! INTENT IN + type(grid_setup_type), intent(in) :: grid_setup + integer, intent(in) :: npets, localpet + + ! INTENT OUT + 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, v + + polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + + + 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/), & + 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) + + ! 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, lon_fields, lat_fields) + + 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 + + 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 diff --git a/src/regridStates/regridStates.F90 b/src/regridStates/regridStates.F90 index f1d4973..c3c8589 100644 --- a/src/regridStates/regridStates.F90 +++ b/src/regridStates/regridStates.F90 @@ -9,10 +9,11 @@ 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 + n_tiles, & + grid_setup_type use utilities, only : error_handler @@ -21,15 +22,12 @@ program regridStates integer, parameter :: max_vars = 10 ! increase if wish to specify more variables ! 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=10) :: mask_type + 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 + type(grid_setup_type) :: grid_setup_in, grid_setup_out + integer :: ierr, localpet, npets integer :: v, SRCTERM @@ -39,15 +37,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/ 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/ n_vars, variable_list, missing_value -!------------------------------------------------------------------------- ! INITIALIZE !------------------------------------------------------------------------- @@ -84,24 +81,24 @@ program regridStates missing_value=-999. ! set defualt - open(41, file='tile2tile.nml', iostat=ierr) - if (ierr /= 0) call error_handler("OPENING tile2tile NAMELIST.", ierr) - read(41, nml=config, iostat=ierr) - if (ierr /= 0) call error_handler("READING tile2tile NAMELIST.", ierr) - close (41) + open(newunit=ut, file='regrid.nml', iostat=ierr) + if (ierr /= 0) call error_handler("OPENING regrid NAMELIST.", ierr) + 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) + !------------------------ ! 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) + call setup_grid(localpet, npets, grid_setup_in, grid_in ) - 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, grid_setup_out, grid_out ) !------------------------ ! Create input and output fields @@ -153,8 +150,9 @@ program regridStates !------------------------ ! read data into input fields - call read_into_fields(localpet, res_atm_in, res_atm_in, trim(fname_in), n_vars, & - variable_list(1:n_vars), trim(dir_in), 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) !------------------------ @@ -203,8 +201,9 @@ program regridStates if (localpet==0) print*,'** Writing out regridded fields' - call write_from_fields(localpet, res_atm_out, res_atm_out, trim(fname_out), & - n_vars, variable_list(1:n_vars), trim(dir_out), 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 @@ -245,3 +244,94 @@ 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, 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, & + fname_coord, dir_coord, & + ires, jres + + namelist /output/ fname, dir, & + gridtype, & + fname_mask, dir_mask, & + 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 + + select case (namel) + case ("input") + read(unt, nml=input, iostat=ierr) + 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 output NAMELIST.", ierr) + case default + call error_handler("unknown namel in readin_setup", 1) + end select + + grid_setup%descriptor = gridtype + + grid_setup%dir = dir + grid_setup%fname = fname + + ! to-do, add routine to check if present. + select case (namel) + case ("input") + grid_setup%dir_mask = dir ! use input file for mask + grid_setup%fname_mask = fname + case ("output") + grid_setup%dir_mask = dir_mask ! need a file on output grid for mask + grid_setup%fname_mask = fname_mask + 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 " ! 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 + 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 + grid_setup%jres = jres + + end subroutine +!------------------------------------------------------------------------- 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 -/