add fortran code for ice or ocn post #25

DeniseWorthen · 2023-10-14T17:28:51Z

fixes add fortran code to remap tripole ocean and ice output to rectilinear #24
depends on Fix weights generation for tripole->rectilinear output of ocean and ice ufs-community/UFS_UTILS#858

This code is verified to reproduce the original NCL scripts w/in ~O10-8, with the exception of the ocean field speed which is now re-calculated from the remapped SSU and SSV fields.

The code will remap either ocean or ice model history files via a run-time flag do_ocnpost set according to the ftype value in the nml file. It is assumed that the model history file is present in the run directory with the name ftype.nc, ie either ocean.nc or ice.nc.

The model history file will be remapped to three rectilinear resolutions, 1/4 deg, 1/2 deg and 1/deg. The output files will be named, for example, ftype.0p25.nc for the 1/4 deg case.

A debug flag is available in the nml file. If set, information will be sent to ftype.post.log to assist debugging. Intermediate netCDF files will also be created during the remapping process.

The required ESMF weights are obtained from the location specified in wgtsdir. The same weights are used for either the ocean or ice case. The current weights in /scratch1/NCEPDEV/climate/climpara/S2S/FIX/fix_UFSp6/fix_reg2grb2/ need to be updated with the weights generated by the UFS_UTILS cpld_gridgen utility after the PR noted above.

Two run directories are available showing how the code works at /scratch2/NCEPDEV/stmp1/Denise.Worthen/ForRahul

HenryRWinterbottom · 2023-10-14T17:33:11Z

src/ocnicepost.fd/outputvars.F90

+
+  end subroutine ocnvars_typedefine
+
+  subroutine icevars_typedefine


Nothing should be hard coded. Any options, including variable names to be used, should be read in either via a FORTRAN namelist or a CSV file. Hard-coding features limits the capabilities of the application.

HenryRWinterbottom · 2023-10-14T17:33:13Z

src/ocnicepost.fd/outputvars.F90

+    !set defaults
+    outvars(:)%input_var_name=''
+    outvars(:)%var_grid  = 'Ct'
+    outvars(:)%var_remapmethod  = 'bilinear'


What happens if a variable requires a different type of interpolation type (e.g., nearest-neighbor)? This should not be hard coded.

State variables should be remapped with linear interpolation. Fluxes could be remapped with nearest neighbor. Is this something desired for MOM6 or CICE6 remapped fields?

@DeniseWorthen What happens if an additional variable is needed in the future and requires nearest-neighbor interpolation? How do you do this without needing to rebuild the entire application?

@HenryWinterbottom-NOAA We can address the customization as an improvement activity.
For now, let us focus on this improved functionality that can be leveraged

HenryRWinterbottom · 2023-10-14T17:33:54Z

src/ocnicepost.fd/outputvars.F90

+    integer :: ii = 0
+
+    !set defaults
+    outvars(:)%input_var_name=''


This entire subroutine needs to be defined via a namelist or CSV file. Nothing should be hard coded as it limits the capabilities of the software.

I'm not sure I understand the comments regarding generality. Variables in MOM6/CICE6 output have specified names, on specified grid locations. MOM6 temperature is never written in the output file on anything other than the center grid locations.

I think @HenryWinterbottom-NOAA's point is to make this configurable so that one can add/remove variables without having to rebuild executables. It separates action (code) from options (configuration).

This can be addressed later as an improvement activity.

HenryRWinterbottom · 2023-10-14T17:36:08Z

src/ocnicepost.fd/ocnicepost.F90

@@ -0,0 +1,577 @@
+program ocnicepost


The netCDF interface should be modified to use the https://github.com/NOAA-EMC/NCEPLIBS-ncio. It is operationally supported code and the netCDF interface is redundant.

The file ocnicepost.F90 is too long and thus difficult to follow. The applications within need to be modified such that they are within standalone modules and subroutines to reduce the bloat of this code.

Re. ncio, I think we can accept code without nceplibs-ncio. It may make things easier to read, but it requires a template to write. So it doesn't fully satisfy the needs of this program.

HenryRWinterbottom · 2023-10-14T17:36:57Z

src/ocnicepost.fd/ocnicepost.F90

+  integer, parameter :: ndest = 3
+  integer, parameter, dimension(ndest) :: nxrs = (/1440, 720, 360/)
+  integer, parameter, dimension(ndest) :: nyrs = (/ 721, 361, 181/)
+  character(len=4), dimension(ndest)   :: dstgrds = (/'0p25', '0p5 ', '1p0 '/)


Hard coding of variables is not supported. The hard coded attributes in this instance should be read into the program via either a FORTRAN namelist or a CSV file.

HenryRWinterbottom · 2023-10-14T17:42:51Z

src/ocnicepost.fd/utils_mod.F90

+  !----------------------------------------------------------
+  ! write a bare netcdf file of a 2D packed field
+  !----------------------------------------------------------
+  subroutine dumpnc2d(fname, vname, dims, nflds, field)


This can be accomplished via ncio.

HenryRWinterbottom · 2023-10-14T17:43:15Z

src/ocnicepost.fd/utils_mod.F90

+  !----------------------------------------------------------
+  ! handle netcdf errors
+  !----------------------------------------------------------
+  subroutine handle_err(ierr,string)


See comments regarding ncio.

HenryRWinterbottom · 2023-10-14T17:43:20Z

src/ocnicepost.fd/utils_mod.F90

+  !----------------------------------------------------------
+  ! write a bare netcdf file of an unpacked 2D field
+  !----------------------------------------------------------
+  subroutine dumpnc1d(fname, vname, dims, field)


See comments regarding ncio.

HenryRWinterbottom · 2023-10-14T17:43:30Z

src/ocnicepost.fd/utils_mod.F90

+  !----------------------------------------------------------
+  ! write a bare netcdf file of an unpacked 3D field
+  !----------------------------------------------------------
+  subroutine dumpnc3dk(fname, vname, dims, field)


See comments regarding ncio.

HenryRWinterbottom · 2023-10-14T17:43:35Z

src/ocnicepost.fd/utils_mod.F90

+  !----------------------------------------------------------
+  ! write a bare netcdf file of a packed 3D field
+  !----------------------------------------------------------
+  subroutine dumpnc3d(fname, vname, dims, nk, nflds, field)


See comments regarding ncio.

HenryRWinterbottom · 2023-10-14T18:18:42Z

src/ocnicepost.fd/outputvars.F90

+    ii = ii + 1; outvars(ii)%input_var_name  = 'SSH'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'SST'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'SSS'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'speed'
+    !ii = ii + 1; outvars(ii)%input_var_name  = 'mld'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'ePBL'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'MLD_003'
+    ii = ii + 1; outvars(ii)%input_var_name  = 'MLD_0125'


@DeniseWorthen You are hard coding variables below (see above). Hard coding paths, variables, etc., that can be read via an external file is not permitted.

aerorahul

Thank you for the PR.
A few comments that would help with the code.

I agree w/ @HenryWinterbottom-NOAA, the main program could be better organized, but I understand and appreciate the contribution when it was not your responsibility.

aerorahul · 2023-10-17T17:16:12Z

src/ocnicepost.fd/ocnicepost.F90

+  integer, parameter :: ndest = 3
+  integer, parameter, dimension(ndest) :: nxrs = (/1440, 720, 360/)
+  integer, parameter, dimension(ndest) :: nyrs = (/ 721, 361, 181/)
+  character(len=4), dimension(ndest)   :: dstgrds = (/'0p25', '0p5 ', '1p0 '/)


Can we make this namelist configurable so that the number of and the destination grid is not hard-wired.

The number of destination grids is dependent on the source grid in the sense that you don't remap the 1/2 deg tripole to the 1/4deg rectilinear. So for 1/4 deg tripole, you have 3 destination grids, for 1/2 deg tripole there are 2 destination grids and for 1deg tripole there is 1 destination grid. We've never set up a remapping of the 5deg tripole, because there was never any intention of anyone looking at that output.

The dimensions of the destination grids can be set based on whether it is a 1/4 deg, 1/2deg or 1 deg. The grid dimensions of those destination grids are fixed though. They are no more variable than the dimensions of the 1/4deg tripole grid.

aerorahul · 2023-10-17T17:17:08Z

src/ocnicepost.fd/ocnicepost.F90

+  ! --------------------------------------------------------
+
+  nvalid = 0
+  rc = nf90_open(trim(input_file), nf90_nowrite, ncid)


Can we please wrap error checking so that the netCDF calls here and below terminate the program in the event of an error?

An example is:

subroutine nc_check(status) integer, intent ( in) :: status if(status /= nf90_noerr) then write(6, '(a)') 'FATAL ERROR: ' // trim(nf90_strerror(status)) stop status end if end subroutine nc_check

and then

call nc_check(nf90_open(trim(input_file), nf90_nowrite, ncid))

I personally do not prefer this style of checking, but I can use it if preferred here.

Understood, but we need to be able to check the return code rc to be non-zero and exit if it is.
I was merely providing a typical example used in a lot of our codes so it feels familiar.
However, feel free to use whatever you prefer.

aerorahul · 2023-10-17T17:18:35Z

src/ocnicepost.fd/ocnicepost.F90

+     rc = nf90_inq_varid(ncid, 'temp', varid)
+     rc = nf90_get_var(ncid, varid, tmp3d)


why not use getfield?

The getfield is not used so that I can retain the fill value, which I use to define the 3-D interpolation mask.

aerorahul · 2023-10-17T17:21:24Z

src/ocnicepost.fd/ocnicepost.F90

+  character(len=120) :: input_file, wgtsfile, output_file
+
+  ! source grid, tripole 1/4 deg, 40 vertical levels
+  integer, parameter :: nxt = 1440, nyt = 1080, nlevs = 40


Is the input file required to be 1/4 degree? Will this program work for MOM6/CICE6 output at any other resolution e.g. 5 degree?

The required weights are generated for 1/4 tripole->(rect 1/4deg, 1/2deg, 1deg), the 1/2 tripole -> (rect 1/2deg, 1deg) and 1deg tripole>rect 1deg. This initial code was intended as a demonstration of how the original NCL scripts could be rewritten in fortran. The NCL scripts could do any of the three tripole resolutions.

The other resolutions can all be added as variations to the fortran version. As noted above however, the 5deg ocean output was never intended as being used for anything other than plumbing. Are you suggesting that the 5deg tripole output will be used for some purpose? Or do you want it for plumbing in this case also?

The destination grid resolutions themselves are fixed. That is, the weights are pre-generated and work only for a rectilinear 1/4deg of 1440x721, for example The code is not expected to operate between completely arbitrary destination and sources.

We discussed this offline.
We will put in place checks to ensure that the target grid resolution does not exceed the input grid resolution; i.e. a 5-deg tripolar input will only produce 5-deg rectilinear output. 1/4-deg tripolar input can produce 1/4, 1/2, 1, 5-deg rectillinear outputs. This will be namelist configurable rather than looping over in the code.

aerorahul · 2023-10-17T18:06:56Z

src/ocnicepost.fd/ocnicepost.F90

@@ -0,0 +1,577 @@
+program ocnicepost


@DeniseWorthen
At a minimum, please add a few lines of what this code does, what the inputs are, etc.
Note that the weight files have hard-wired names down in this code. Please be sure to document that.
Also include the authorship, so users can know who to ask questions in the case an error is encountered.

* add csv files to set source field info * restrict destination grid to a single resolution

DeniseWorthen · 2023-10-20T16:54:36Z

The code has been restructured. Per conversation with @aerorahul, looping over destination grids has been removed. Samples of the input nml and csv files are found in staged directories /scratch2/NCEPDEV/stmp1/Denise.Worthen/ForRahul.

The code does still hard-wire the name of the time and vertical axis in the source file. The implementation of netcdf-error checking and making the time axis and vertical grid names configurable are relatively minor changes and have not been implemented at this time.

Weights are not yet available for mapping the 5-deg tripole or for remapping to a 5-deg rectilinear but there should be no restriction on the code being able to do that once weights are available.

aerorahul

Thank you.

add fortran code for ice or ocn post

a9170c2

HenryRWinterbottom self-assigned this Oct 14, 2023

HenryRWinterbottom suggested changes Oct 14, 2023

View reviewed changes

HenryRWinterbottom reviewed Oct 14, 2023

View reviewed changes

aerorahul requested changes Oct 17, 2023

View reviewed changes

aerorahul reviewed Oct 17, 2023

View reviewed changes

DeniseWorthen added 3 commits October 20, 2023 10:29

restructure code, use configurable files

6518d72

* add csv files to set source field info * restrict destination grid to a single resolution

tweaks, remove unused files

7709ef8

add program description

0a404fe

aerorahul self-requested a review October 20, 2023 17:30

aerorahul approved these changes Oct 20, 2023

View reviewed changes

aerorahul merged commit 50a06ee into NOAA-EMC:develop Oct 20, 2023

DeniseWorthen mentioned this pull request Oct 27, 2023

add ESMF weights for tripole->rectilinear post NOAA-EMC/global-workflow#1992

Closed

4 tasks

aerorahul mentioned this pull request Nov 21, 2023

Trap netcdf errors when they happen in ocnicepost.fd. #35

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add fortran code for ice or ocn post #25

add fortran code for ice or ocn post #25

DeniseWorthen commented Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

DeniseWorthen Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

aerorahul Oct 17, 2023

HenryRWinterbottom Oct 14, 2023

DeniseWorthen Oct 14, 2023

aerorahul Oct 17, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

aerorahul Oct 17, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

HenryRWinterbottom Oct 14, 2023

aerorahul left a comment

aerorahul Oct 17, 2023

DeniseWorthen Oct 17, 2023

aerorahul Oct 17, 2023

DeniseWorthen Oct 17, 2023

aerorahul Oct 17, 2023

aerorahul Oct 17, 2023

DeniseWorthen Oct 17, 2023

aerorahul Oct 17, 2023

DeniseWorthen Oct 17, 2023

aerorahul Oct 17, 2023

aerorahul Oct 17, 2023

DeniseWorthen Oct 17, 2023

DeniseWorthen commented Oct 20, 2023

aerorahul left a comment


		end subroutine ocnvars_typedefine

		subroutine icevars_typedefine

		rc = nf90_inq_varid(ncid, 'temp', varid)
		rc = nf90_get_var(ncid, varid, tmp3d)

add fortran code for ice or ocn post #25

add fortran code for ice or ocn post #25

Conversation

DeniseWorthen commented Oct 14, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aerorahul left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DeniseWorthen commented Oct 20, 2023

aerorahul left a comment

Choose a reason for hiding this comment