forcing interpolation not working

[d.kobashi]

Hi,

Can anyone tell me if the interpolation of forcing works? I mean that the interpolation is the interpolation done internally by ROMS.
I am trying to create forcing file (parent and child grids); however, due to the large file size, I thought it should be better for ROMS to take care of interpolation rather than creating forcing file interpolated into ROMS by me.

However, when I ran the model, all forcing data were zeros in history files and the maximum and minimum values of input data (e.g. wind velocities, heat fluxes, precipitation etc.) were 1e+35 and -1e+35 in the log, respectively (see below). These values come from Fmin and Fmax in interpolate.F.
I also checked that the values of Fout/Finp in interpolate.F and it seems that those values were NaN (did not show the values). That's why I think min and max values were 1e+35. I double-checked the values and there are no 1e+35 values in any grids of the forcing file. Also, the domain of the forcing is larger than the model domain.

Code: Select all

    GET_2DFLD   - surface u-wind component,                        2014-01-01 01:00:00.00
                   (Grid=01, Rec=0000002, Index=2, File: nwgom.frc.2014.nc)              
                   (Tmin=      16071.0000 Tmax=      16435.9583)      t =      16071.0417
                   (Min =  1.00000000E+35 Max = -1.00000000E+35)      regrid = T  

I wonder if anyone had a similar problem.
I attach the metadata of the uninterpolated forcing file below. Does the metadata look wrong?
I use COWAST ver 3.5 (ROMS 3.9?).

Code: Select all

netcdf nwgom.frc.2014 {
dimensions:
	time = UNLIMITED ; // (8760 currently)
	lon = 81 ;
	lat = 61 ;
variables:
	double time(time) ;
		time:long_name = "bulk formulation atmospheric forcing time" ;
		time:units = "days since 1970-01-01 00:00:00" ;
	double lon(lon) ;
		lon:long_name = "Longitude" ;
		lon:units = "degrees_east" ;
		lon:axis = "X" ;
	double lat(lat) ;
		lat:long_name = "Latitude" ;
		lat:units = "degrees_north" ;
		lat:axis = "Y" ;
	double Tair(time, lat, lon) ;
		Tair:long_name = "surface air temperature" ;
		Tair:units = "Celsius" ;
		Tair:coordinates = "lon lat" ;
		Tair:time = "time" ;
	double Pair(time, lat, lon) ;
		Pair:long_name = "surface pressure" ;
		Pair:units = "millibar" ;
		Pair:coordinates = "lon lat" ;
		Pair:time = "time" ;
	double Qair(time, lat, lon) ;
		Qair:long_name = "relative humidity" ;
		Qair:units = "percentage" ;
		Qair:coordinates = "lon lat" ;
		Qair:time = "time" ;
	double rain(time, lat, lon) ;
		rain:long_name = "precipitation rate" ;
		rain:units = "kg m-2 s-1" ;
		rain:coordinates = "lon lat" ;
		rain:time = "time" ;
	double swflux(time, lat, lon) ;
		swflux:long_name = "net surface freshwater flux" ;
		swflux:units = "cm day-1" ;
		swflux:coordinates = "lon lat" ;
		swflux:time = "time" ;
		swflux:positive_value = "downward flux, freshening" ;
		swflux:negative_value = "upward flux, salting" ;
	double lwrad(time, lat, lon) ;
		lwrad:long_name = "net longwave radiation flux" ;
		lwrad:units = "Watts metre-2" ;
		lwrad:coordinates = "lon lat" ;
		lwrad:time = "time" ;
		lwrad:positive_value = "downward flux, heating" ;
		lwrad:negative_value = "upward flux, cooling" ;
	double lwrad_down(time, lat, lon) ;
		lwrad_down:long_name = "downward longwave radiation" ;
		lwrad_down:units = "Watts metre-2" ;
		lwrad_down:coordinates = "lon lat" ;
		lwrad_down:time = "time" ;
	double swrad(time, lat, lon) ;
		swrad:long_name = "shortwave radiation" ;
		swrad:units = "Watts metre-2" ;
		swrad:coordinates = "lon lat" ;
		swrad:time = "time" ;
		swrad:positive_value = "downward flux, heating" ;
		swrad:negative_value = "upward flux, cooling" ;
	double shflux(time, lat, lon) ;
		shflux:long_name = "net surface heat flux" ;
		shflux:units = "Watts metre-2" ;
		shflux:coordinates = "lon lat" ;
		shflux:time = "time" ;
		shflux:positive_value = "downward flux, heating" ;
		shflux:negative_value = "upward flux, cooling" ;
	double Uwind(time, lat, lon) ;
		Uwind:long_name = "u-wind" ;
		Uwind:units = "metre second-1" ;
		Uwind:coordinates = "lon lat" ;
		Uwind:time = "time" ;
	double Vwind(time, lat, lon) ;
		Vwind:long_name = "v-wind" ;
		Vwind:units = "metre second-1" ;
		Vwind:coordinates = "lon lat" ;
		Vwind:time = "time" ;
	double cloud(time, lat, lon) ;
		cloud:long_name = "cloud fraction" ;
		cloud:units = "nondimensional" ;
		cloud:coordinates = "lon lat" ;
		cloud:time = "time" ;
	double wspd(time, lat, lon) ;
		wspd:long_name = "wind speed 10m" ;
		wspd:units = "metre second-1" ;
		wspd:coordinates = "lon lat" ;
		wspd:time = "time" ;
	double sustr(time, lat, lon) ;
		sustr:long_name = "surface u-momentum stress" ;
		sustr:units = "Newton metre-2" ;
		sustr:coordinates = "lon lat" ;
		sustr:time = "time" ;
	double svstr(time, lat, lon) ;
		svstr:long_name = "surface v-momentum stress" ;
		svstr:units = "Newton metre-2" ;
		svstr:coordinates = "lon lat" ;
		svstr:time = "time" ;

// global attributes:
		:type = "ROMS Bulk forcing file" ;
		:title = "ROMS Bulk forcing file" ;
		:author = "DJ,,,, wegener.tamu.edu" ;
		:dataset = "ECMWF ERA-5" ;
		:NCO = "\"4.6.0\"" ;
		:nco_openmp_thread_number = 1 ;
}

Thanks in advance.

-DJ

[jcwarner]

i just ran the Sandy test case with grid refinement 2 grids roms +swan and it worked fine. it reads in met data for both roms grids.

The output you show:
Min = 1.00000000E+35 Max = -1.00000000E+35)
is a pure report from roms of what it reads from your netcdf forcing file. do an ncview of that forcing file.

Yes, if the file is big you can just list it twice and let roms do the interpolation to each grid:
NFFILES == 1
FRCNAME == Projects/Sandy/romsforc_NARR_Sandy2012.nc \
Projects/Sandy/romsforc_NARR_Sandy2012.nc

so i list 1 forc file, probably should have that as
NFFILES == 1 1
then i have the same force file name repeated, once for the parent, once for the child.

-j

[d.kobashi]

Thanks, John.

It looks like that the forcing file you used for Sandy test case, i.e. romsforc_NARR_Sandy2012.nc, was interpolated into your parent grids already. Am I correct? Then you used the same file for the child grids without any modifications?

I am not sure if this is correct, but I assume that ROMS takes care of the interpolation internally for both parent and child grids. The file I created has the same resolution as the native forcing grid originated from the source (in this case ECMWF's ERA5). The resolution of ERA5 is 30 something kilometers and my grid resolution is 1.5 km and 500 m so I even don't want to interpolate the original forcing (~1.5 GB in total) into the parent grids (250+ GB in total) as well as the child grids (500+ GB in total).

So for the ROMS function for the interpolation, do we need to interpolate the external forcing into the parent grids at least?


On second thought, I've started to realize that the ROMS online nesting whether it is one way or two way, is not necessarily practical for my applications (1-year simulation) due to enormous computational cost. I start working on offline nesting to see if it is reasonable in terms of computational cost and accuracy.

Thanks,

-DJ

[jivica]

ROMS should interpolate atmo fields for you (unless you have exactly the same dimensions as your grid).
I am using the approach but have defined lon/lat as 2d fileds. This is simple, just repmat and it cost nothing as doesn't change in time.
Problem I had was if you are in the south hemisphere then location of first index does matter (1,1)

Cheers
I
p.s. for nesting we are using 1-way downscale with "stations". If interested search for the keywords.

[d.kobashi]

Thanks for your input, Ivica.

I actually did change the coordinate variables from 1D to 2D and the result was the same.
The problem about 1e+35 values in the log file comes from interpolate.F below. For Fmin and Fmax on the lines 133 and 134, Finp was empty. So I suspect ROMS did not read the data from forcing netCDF files correctly. Don't know why. I need to dig in further.

Do you have the metadata of your forcing files you used for nesting (log you get when typing ncdump)? The COAWST's Sandy test case John mentioned appears to use the forcing already interpolated into the model grids, looking at the metadata structure.

Code: Select all

 111       Fmin=1.0E+35_r8
 112       Fmax=-1.0E+35_r8
 113       write(*,*) 'linear interpolation'
 114       DO j=Jstr,Jend
 115         DO i=Istr,Iend
 116           i1=INT(Iout(i,j))
 117           i2=i1+1
 118           j1=INT(Jout(i,j))
 119           j2=j1+1
 120           IF (((LBx.le.i1).and.(i1.le.UBx)).and.                        &
 121      &        ((LBy.le.j1).and.(j1.le.UBy))) THEN
 122             p2=REAL(i2-i1,r8)*(Iout(i,j)-REAL(i1,r8))
 123             q2=REAL(j2-j1,r8)*(Jout(i,j)-REAL(j1,r8))
 124             p1=1.0_r8-p2
 125             q1=1.0_r8-q2
 126             write(*,*) p1, q1, p2, q2
 127             write(*,*) Finp(i1,j1), Finp(i2,j2)
 128             Fout(i,j)=p1*q1*Finp(i1,j1)+                                &
 129      &                p2*q1*Finp(i2,j1)+                                &
 130      &                p2*q2*Finp(i2,j2)+                                &
 131      &                p1*q2*Finp(i1,j2)
 132             write(*,*) Fout(i,j)
 133             Fmin=MIN(Fmin,Fout(i,j))
 134             Fmax=MAX(Fmax,Fout(i,j))
 135           END IF
 136         END DO
 137       END DO

Also, I am curious to know more about getting boundary conditions using station files. I found your paper using the method.

Janekovic and Powell 2012. Analysis of imposing tidal dynamics to nested numerical models. CSR.

I just skimmed the paper and if I understand correctly, this method only creates boundary file not boundary and climatology files. Am I correct? So the climatological nudging used for RadNud cannot be used. Your paper said the method uses clamped BC for all open boundaries, which does not requires clim nudging. Does the clamped BC work well for offline nesting for simulating baroclinic processes? because I have not had much success using it. I am currently trying to run non-tidal simulation.

Thanks,

-DJ

[jivica]

You have cdl example file in your ROMS/Data/ROMS/CDL/
and this is my ncdump -h forcing.nc

Code: Select all

netcdf forcing {
dimensions:
	time = UNLIMITED ; // (56 currently)
	dim1 = 221 ;
	dim2 = 341 ;
variables:
	double time(time) ;
		time:long_name = "time" ;
		time:field = "time, scalar, series" ;
		time:units = "days since 2000-01-01 00:00:00" ;
	double Uwind(time, dim1, dim2) ;
		Uwind:long_name = "U-wind componente at 10m" ;
		Uwind:time = "wind_time" ;
		Uwind:field = "Uwind, scalar, series" ;
		Uwind:units = "m s-1" ;
		Uwind:coordinates = "lon lat" ;
		Uwind:least_significant_digit = 3 ;
	double Vwind(time, dim1, dim2) ;
		Vwind:long_name = "V-wind componente at 10m" ;
		Vwind:time = "wind_time" ;
		Vwind:field = "Vwind, scalar, series" ;
		Vwind:units = "m s-1" ;
		Vwind:coordinates = "lon lat" ;
		Vwind:least_significant_digit = 3 ;
	double Qair(time, dim1, dim2) ;
		Qair:long_name = "Surface air relative humidity" ;
		Qair:time = "qair_time" ;
		Qair:field = "Qair, scalar, series" ;
		Qair:units = "[0-100]" ;
		Qair:coordinates = "lon lat" ;
		Qair:least_significant_digit = 3 ;
......

As you can see, there is nothing special (I am using nco --ppc default=.3 for additional compression as I don't need wind in mm/s etc.). Be sure that you don't have missing values, and from the output you attached I have feeling you have. Just try with ncdump and make sure there are no _missing or _filled values. Interpolating inside ROMS is super simple and there are not many places it can go wrong, so I believe your data is to blame. If you can't find the problem, ncks first record ROMS is complaining about and send it through (one time record, and use compression to make it smaller).

Regarding our paper; I would say depending what you are doing. In our case we had small nested domain which was quickly washed with boundary conditions so nudging to climatology didn't make sense and what we did seemed appropriate. If you really need to nudge it doesn't have to be at the same time stepping as boundary (i.e. you can nudge to monthly state and compute those fields from the big grid by using interpolation) and this is usually how we are using nudge time scales....

Cheers
Ivica

[jcwarner]

i just checked my forcing for the Sandy test case and the atm forcing is on a separate grid that is not identical to either of the roms grids.
I remember changing this a while ago to make the test case forcing files a lot smaller, it is just a test case.
So for the Sandy case, ROMS does read correctly the forcing on a different grid and interpolate to roms parent and child grids.

[arango]

I also have been using this capability for years and frequently and I never have any problems. Our group at Rutgers uses the same strategy on a daily basis in our operational system. Something is wrong in your files. When I do nesting, I only specify the forcing files for the first grid, and ROMS will take care of the rest. I don't need to repeat the files for every nested grid. For example, in our US East coast three nested grid application I just have:

Code: Select all

  NFFILES == 8                          ! number of unique forcing file

  FRCNAME == ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/Pair_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/Pair_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/Qair_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/Qair_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/rain_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/rain_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/Tair_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/Tair_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/Uwind_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/Uwind_nam_3hourly_MAB_and_GoM_2015.nc \

             ../om/Vwind_nam_3hourly_MAB_and_GoM_2014.nc |
             ../om/Vwind_nam_3hourly_MAB_and_GoM_2015.nc

I always suggest to the users to look at ROMS code when having problems and do not use ROMS as a black box. The routines are well documented. For example, check regrid.F and get_varcoords.F. In get_varcoords.F, you will see:

Code: Select all

#include "cppdefs.h"
      SUBROUTINE get_varcoords (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, Nx, Ny,               &
     &                          Xmin, Xmax, X, Ymin, Ymax, Y,           &
     &                          rectangular)
!
!svn $Id: get_varcoords.F 1821 2020-01-10 03:54:15Z arango $
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This routine reads the spatial positions of any NetCDF variable     !
!  having the attribute  "coordinates",  as specified by CF rules.     !
!  For example, in CDL syntax:                                         !
!                                                                      !
!       float my_var(time, lat, lon) ;                                 !
!               my_var:long_name = "my variable long name" ;           !
!               my_var:units = "my variable units" ;                   !
!               my_var:coordinates = "lon lat time" ;                  !
!               my_var:time = "time" ;                                 !
!                                                                      !
!  The following "coordinates" attribute is also allowed:              !
!                                                                      !
!               my_var:coordinates = "lon lat" ;                       !
!                                                                      !
!  That is, the time variable "time" is missing in the "coordinates"   !
!  attribute.                                                          !
!                                                                      !
!  Notice that the associated coordinate names "lon" and "lat" are     !
!  separated by a single blank space.  Both "lon" and "lat" can be     !
!  1D or 2D arrays. If 1D array, the positions are rectangular and     !
!  and full 2D arrays are filled with the same values.                 !
!                                                                      !
!  It also determines the rectangular switch  which indicates that     !
!  the spatial positions have a plaid distribution.                    !
!                                                                      !
!=======================================================================

This information is somewhere in this forum and in trac. You need to follow the metadata model, Otherwise, the interpolation WILL NOT WORK.

[mma]

Dear DJ,

When dealing with ocean modelling it is normal to use the forum before checking all the code. That is the purpose of this forum, I believe. Suggesting that a model user is lazy and uses ROMS as a black box is a very sad thing to say and unfortunately very frequent in this forum.

I found the problem and is has nothing to do with laziness, except for some parts of the code poorly documented. The problem is in subroutine hindices, inside interpolate.F. And this is the kind of elementary auxiliary tools that an user will hardly check, simply because people trust such tools, just like they trust the Matlab griddata function.

The problem is that these lines:

IF ((Ygrd(1,j ).le.Ypos(mp,np)).and. &
& (Ygrd(1,j+1).gt.Ypos(mp,np))) THEN
Jmin=j
foundj=.TRUE.
EXIT J_LOOP
END IF

assume that latitude must be increasing! If you perform a griddata in Matlab or Python you can use a decreasing rectangular latitude. But here you cannot. That is not mentioned in the code and unfortunately no error is raised when foundi or foundj are false.

If you check your ERA5 data you will see that the latitudes decrease. So, I suggest you to flip all the ERA5 variables when creating the forcing file and try again. You may also be interested in changing the interpolation method which is currently linear by default (check InterpFlag in mod_scalars.F).

best wishes
mma