Regridding of winds and their rotation within ROMS

[lanerolle]

It is my belief that the surface forcing (meteorological/atmospheric forcing) in ROMS can be done in two ways:

(1) prescribe the two wind stress components and the net heat flux directly via values read from a NetCDF file or

(2) use the Bulk Flux formulation in ROMS which needs air temperature, pressure, relative humidity, two wind speed components, net shortwave radiation and if available, downward longwave radiation (otherwise use Berliand formulation).

In both the above scenarios, we need to provide all of the quantities/variables such that:

(a) they are all at rho-points in the horizontal staggered grid,
(b) they have been interpolated on to the ROMS grid and,
(c) the wind speeds/stresses have been rotated to be locally orthogonal within the ROMS grid.

For long simulations, even if the different surface forcing variables/quantities are stored in separate NetCDF files, these files can get very large. Hence, a good alternative is to store the variables/quantities in NetCDF files not on the ROMS grid (after interpolation) but on the native grid associated with these variables which may or may not be gridded. If the horizontal coverage of these points exceeds the limits of the ROMS grid then, ROMS will automatically interpolate these variables on to the grid internally (and bi-linear, bi-cubic choices can be selected via the "linear"/"cubic" parameters in the mod_scalars.F file).

I have a question (sorry if it is silly and if I have overlooked something) : if we do the spatial interpolation internally within ROMS, do the wind speeds/stresses get appropriately rotated before/during/after the interpolation procedure? I could not find any code segments to do this in the routines contained in the ROMS/Utility directory.

[kate]

Yes, the code is there. Look for CosAngler in set_data.F.

[arango]

I have a question (sorry if it is silly and if I have overlooked something) : if we do the spatial interpolation internally within ROMS, do the wind speeds/stresses get appropriately rotated before/during/after the interpolation procedure? I could not find any code segments to do this in the routines contained in the ROMS/Utility directory.

Yes, ROMS will rotate the vector components to (xi,eta) coordinates. Check set_data.F:

Code: Select all

!
!  If input point surface winds or interpolated from coarse data, rotate
!  to curvilinear grid.
!
      IF (.not.Linfo(1,idUair,ng).or.                                   &
     &    (Iinfo(5,idUair,ng).ne.Lm(ng)+2).or.                          &
     &    (Iinfo(6,idUair,ng).ne.Mm(ng)+2)) THEN
        DO j=JstrR,JendR
          DO i=IstrR,IendR
            cff1=FORCES(ng)%Uwind(i,j)*GRID(ng)%CosAngler(i,j)+         &
     &           FORCES(ng)%Vwind(i,j)*GRID(ng)%SinAngler(i,j)
            cff2=FORCES(ng)%Vwind(i,j)*GRID(ng)%CosAngler(i,j)-         &
     &           FORCES(ng)%Uwind(i,j)*GRID(ng)%SinAngler(i,j)
            FORCES(ng)%Uwind(i,j)=cff1
            FORCES(ng)%Vwind(i,j)=cff2
          END DO
        END DO
        ...

The interpolation assumes that the 2D input data in the forcing NetCDF file is grided and vectors components are true North and East. When I coded regrid.F I did thought about the vector fields having a different rotation but it was kind of dangeorus because angles are not difined in the same way. We needed a standard NetCDF attribute for this. So I set the input rotation angle to zero:

Code: Select all

!
!  Set input gridded data rotation angle.
!
      DO i=1,Nx
        DO j=1,Ny
          angle(i,j)=0.0_r8
        END DO
      END DO

until I figure out a way to do this rotation from original atmospheric grid to true geographical North and East. Anyway grided vectors are usually like this. The rotation can be accounted in model coupling with ESMF or MCT but this is a different story.

Therefore, if you are using a coarser grided forcing you need to make sure that vector fields are oriented with respect to true geographical North and East.

[lanerolle]

Thanks very much for your informative reply.

So, just to confirm, are you saying that in addition to the winds being in true Northings and Eastings (which seems very reasonable to me) they also need to be provided in a NetCDF file in a 2D array format as opposed a 1D array/vector format?

What we normally do is to get say gridded NARR (North American Regional Reanalysis) met products and then blend the winds and the atmospheric pressure together with observed values (with say 5-10 observation stations scattered within the computational domain) to increase spatio-temporal resolution. This can be easily done with OA or in Matlab but I guess it can not be done in a straightforward manner via interpolation within ROMS because by including the discrete observations, I am destroying the 2D gridded array structure of the met forcings and the only way to include say NARR + obs is to collapse them both in to a single 1D vector which does not work with the ROMS interpolation. I guess I will have to find a workaround to achieve this!

[kate]

So, just to confirm, are you saying that in addition to the winds being in true Northings and Eastings (which seems very reasonable to me) they also need to be provided in a NetCDF file in a 2D array format as opposed a 1D array/vector format?

Yes, that's true. ROMS is currently only handling the simplest of cases. The intent here was to save disk space since some gridded fields are much coarser than our current ROMS grid. If you've got denser information, then that savings doesn't hold. You might as well do the blending onto the ROMS grid as best you can outside of ROMS or else you will have to add that functionality to ROMS itself.

[wilkin]

If the wind forcing is provided as a 1-D vector (in time) with no spatial coordinates the rotation to ROMS coordinates is ALWAYS done according to the angle in the grid file, i.e. the 1-D winds (u10 or stress) are assumed in east-north coordinates. You can find the code in set_data.F. This has to be the case because the ROMS grid is assumed curvilinear and therefore the wind must be rotated differently in different regions of the grid.

If you are blending 2-D spatially resolved NARR wind with 1-D point time series observations (though be aware that NCEP may have already assimilated those obs for you!) then you'll be producing a modified 2-D spatially resolved product. So just keep the same format as the NARR data. All it requires is that wind components have the necessary 'coordinates' attributes to allow the ROMS regrid function to work. Keep the output of your OA in east-north. If you're combining NARR and point obs to create a point time series that is assumed spatially uniform (not sure why you would do this if you have NARR) then the same applies - the 1-D time series HAS to be in north-east if your grid is curvilinear.

[azhang]

I read all messages about this topic, but I am lost at the end. For my understanding, there are two kinds of wind forcing formats for ROMS: one is single point wind time series (1-D) and the other is spatial variation wind forcing (2-D). For 1-D wind, the u and v are true east and north winds, ROMS will rotate u and v to ROMS curvilinear local coordinates. For the 2-D wind, the u and v in ROMS surface wind forcing file have to be rotated to ROMS curvilinear local coordinates outside ROMS by modeler before running ROMS. From the John's reply, it seems that ROMS can take coarser grid wind forcings from an atmospheric model (e.g. ETA, NAM, NARR), if this is true, are u and v from the coarser grid model are true east and true north wind, and then ROMS rotates u and v into ROMS curvilinear local coordinates inside ROMS? How does ROMS know the u and v are on ROMS local curvilunear coordinates or on earth coordinates (true east and north)?
Thank you very much for your clarifying.

AJ

[jcwarner]

Wouldn't it be nice to put this information on the wiki and then the users could just read that. Then it could be updated with clarification if needed.

Isn't that the purpose of the wiki???

[kate]

Wouldn't it be nice to put this information on the wiki and then the users could just read that. Then it could be updated with clarification if needed.

Isn't that the purpose of the wiki???

Yes, that would be the purpose - or you could go add it and then post the linky here (3.1). A wiki is not built by one, but by many.

[lanerolle]

I have done some experiments with ROMS where:

(1) I created the fields necessary for the bulk fluxes via Matlab using bicubic interpolation on the ROMS grid and rotated the winds to be locally orthogonal on the grid and also

(2) I created the fields necessary for bulk fluxes on a coarser, bigger grid (5km x 5km) with winds in the true Northings and Eastings

The scenario (2) causes ROMS to carry out the spatial interpolation of the surface met forcings within itself. In both cases, the time interpolation of the data (linear) is carried out within ROMS.

Then, I ran (1) and (2) with ROMS (with cubic=1 in mod_scalars.F so that it will use its bicubic interpolation scheme) and I examined the ensuing wind stresses (the two components) and the surface net heat flux in a ROMS history file.

I am happy to announce that the two sets of winds stresses and net heat fluxes resemble each other (respectively) closely!!! So it appears that the ROMS internal bicubic interpolation is working as expected and the wind rotations are also done correctly.

[lolhsson]

So, on the topic of regridding winds from NARR data... here's a question that has more to do with basic netCDF processing than ROMS-specifics. (Sorry for that.)

I've noticed that all of my data files from NCEP's NARR have a time dimension, and then x and y dimensions -- longitude and latitude are stored as variables lon(x,y) and lat(x,y) instead of as dimensions.

Like so:

ncdump -h unpackedcutair.sfc.2005.nc
netcdf unpackedcutair.sfc.2005 {
dimensions:
tair_time = UNLIMITED ; // (2920 currently)
y = 36 ;
x = 32 ;
variables:
int Lambert_Conformal ;
Lambert_Conformal:grid_mapping_name = "lambert_conformal_conic" ;
Lambert_Conformal:standard_parallel = 50., 50. ;
Lambert_Conformal:longitude_of_central_meridian = -107. ;
Lambert_Conformal:latitude_of_projection_origin = 50. ;
Lambert_Conformal:false_easting = 5632642.22547 ;
Lambert_Conformal:false_northing = 4612545.65137 ;
float Tair(tair_time, y, x) ;
Tair:units = "degK" ;
Tair:long_name = "3-hourly air temperature at Surface" ;
Tair:unpacked_valid_range = 180.f, 365.f ;
Tair:precision = 2s ;
Tair:actual_range = 208.25f, 338.02f ;
Tair:missing_value = 32766s ;
Tair:valid_range = -32766s, -14266s ;
Tair:_FillValue = -32767.f ;
Tair:GRIB_name = "TMP" ;
Tair:GRIB_id = 11 ;
Tair:var_desc = "air temperature" ;
Tair:standard_name = "air_temperature" ;
Tair:level_desc = "Surface" ;
Tair:dataset = "NARR 3-hourly" ;
Tair:statistic = "Individual Obs" ;
Tair:parent_stat = "Other" ;
Tair:grid_mapping = "Lambert_Conformal" ;
Tair:coordinates = "lat lon" ;
float lat(y, x) ;
lat:long_name = "latitude coordinate" ;
lat:units = "degrees_north" ;
lat:axis = "Y" ;
lat:coordinate_defines = "point" ;
lat:standard_name = "latitude" ;
float lon(y, x) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude coordinate" ;
lon:axis = "X" ;
lon:coordinate_defines = "point" ;
lon:standard_name = "longitude" ;
double tair_time(tair_time) ;
tair_time:units = "hours since 1800-1-1 00:00:0.0" ;
tair_time:long_name = "analysis time" ;
tair_time:axis = "T" ;
tair_time:standard_name = "time" ;
tair_time:coordinate_defines = "start" ;
tair_time:delta_t = "0000-00-00 03:00:00" ;
tair_time:actual_range = 1797000., 1805757. ;
float x(x) ;
x:long_name = "eastward distance from southwest corner of domain in projection coordinates" ;
x:units = "m" ;
x:standard_name = "projection_x_coordinate" ;
float y(y) ;
y:long_name = "northward distance from southwest corner of domain in projection coordinates" ;
y:units = "m" ;
y:standard_name = "projection_y_coordinate" ;

// global attributes:
:standardpar1 = 50. ;
:standardpar2 = 50.000001 ;
:centerlon = -107. ;
:centerlat = 50. ;
:latcorners = 1.000001f, 0.897945f, 46.3544f, 46.63433f ;
:loncorners = -145.5f, -68.32005f, -2.569891f, 148.6418f ;
:stream = "s1" ;
:title = "8x Daily NARR" ;
:Conventions = "CF-1.0" ;
:history = "Mon May 2 13:00:43 2011: ncrename -v air,Tair unpackedcutair.sfc.2005.nc\n",
"Mon May 2 13:00:02 2011: ncrename -d time,tair_time unpackedcutair.sfc.2005.nc\n",
"Mon May 2 12:59:28 2011: ncrename -v time,tair_time unpackedcutair.sfc.2005.nc\n",
"Mon May 2 12:49:34 2011: ncpdq -U cutair.sfc.2005.nc unpackedcutair.sfc.2005.nc\n",
"Mon May 2 12:46:23 2011: ncks -d x,109,140 -d y,106,141 air.sfc.2005.nc cutair.sfc.2005.nc\n",
"created 2006/07 by data management at NOAA/ESRL/PSD" ;
:institution = "National Centers for Environmental Prediction" ;
:platform = "Model" ;
:references = "http://wwwt.emc.ncep.noaa.gov/mmb/rreanl/index.html\n";,
"http://www.esrl.noaa.gov/psd/data/gridd ... .narr.html"; ;
:nco_openmp_thread_number = 1 ;

I would really like these files to have dimensions of lon and lat, cutting out the middlemen of x and y, over which my variable Tair could be found -- I assume this is necessary for ROMS to use it as a forcing file that will interpolate onto my grid. I dislike MATLAB with the fiery passion of a thousand suns, and decided it was long past time for me to learn NCO. I thought this would be easy to do (and it probably is, for the experts) but after floundering a bit in the ncpdq and ncap2 sections of the NCO manual I have come to acknowledge that I have no idea how to even start.

Do we have any NCO fans in here in the ROMS forums who've had to do this sort of thing before and could offer some guidance? Or should I give up on this tactic and figure out which of the ROMS MATLAB tools might be able to help me out?

Thanks for your time and patience!

-Liz

[kate]

I assume the NARR output is in x,y like that because the grid is not on a regular lat/lon grid. ROMS is not set up to regrid from a generic WRF type grid (or other random square plopped onto the sphere). You need to regrid to either a lat/lon grid or to your ROMS grid.

I have no idea how to do this in NCO. I could probably figure out how to do it in NCL, but we have instead gone down the road of doing these operations in Python, having it call scrip. Let me know if that sounds like more fun than Matlab to you.

[lolhsson]

Thanks for the prompt reply!

Before I go further, I guess I need to make sure -- does the lat/lon grid that ROMS expects have to be regular (in the 'evenly spaced lat/lon points' sense), thus requiring interpolation from the very irregular lat/lon grid in my dataset before I show it to ROMS at all, or can I 'simply' invert from x,y to those irregular lat,lon points, and then the ROMS core will take care of the interpolation later on?

[kate]

ROMS needs the input to be either (a) on the ROMS grid or (b) on a grid that's regular in lat and in lon, such as a global 2 degree grid. For (b), it will expect dimension variables lat(lat) and lon(lon), none of this lat(x,y) business.

[lolhsson]

Interpolation it is, then. I think I will try to use the myriad MATLAB tools available through myroms.org - there are so many, and they're so well-documented, that I might as well take advantage. (To those of you who've developed them, thanks!) In the long term I'm planning on three one-way nested grids to improve spatial resolution around a region of interest (ala the Coastal Gulf of Alaska team's work) -- is it worth interpolating in pre-processing forcing files for all three ROMS grids, or just to the coarsest of the three, which is then also used to force the other two (essentially, option a) for the first one, then b) for the finer two)? The major reason I might guess it is worth gridding all three in advance is to save on computation time later on, but on the other hand that's a bunch of extra space taken up by two extra sets of forcing files, and my scratch directory is far from infinite, so maybe not. Is the computation savings significant for a) over b), have people found?

Thanks again!

[ilgar]

I run ROMS with wind stress based on NARR wind fields. Although I can run with wind stress manually interpolated and rotated to my nested grids, in order to avoid diskspace issues and rotation several times, I would like to be able to run directly with the wind stresses I interpolate from NARR data to a large native grid (with lat&lon) covering all my nested grids. Wiki, and several postings here and within the forum reveal that it is possible. I generated a windstress file, oriented wrt true North and East on the large native grid, with variables and dimensions ROMS would have expected if this large grid and ROMS grid had been overlapping (with lat_u, lon_u, lat_v, lon_v ; and sustr and svstr on these lat_u&lon_u and lat_v&lon_v pairs, respectively) which might be the issue if the code expects data on different variables (just “lat” and “lon” maybe, as Kate said above):

Code: Select all

narr_windstress {
dimensions:
	sms_time = 113 ;
	xi_psi = 980 ;
	xi_rho = 981 ;
	xi_u = 980 ;
	xi_v = 981 ;
	eta_psi = 760 ;
	eta_rho = 761 ;
	eta_u = 761 ;
	eta_v = 760 ;
variables:
	double sms_time(sms_time) ;
		sms_time:long_name = "surface stress time" ;
		sms_time:units = "days since 1858-11-17 00:00:00 UTC" ;
		sms_time:field = "scalar" ;
	double lat_u(eta_u, xi_u) ;
		lat_u:long_name = "latitude of U-points" ;
		lat_u:units = "degree_north" ;
	double lon_u(eta_u, xi_u) ;
		lon_u:long_name = "longitude of U-points" ;
		lon_u:units = "degree_east" ;
	double lat_v(eta_v, xi_v) ;
		lat_v:long_name = "latitude of V-points" ;
		lat_v:units = "degree_north" ;
	double lon_v(eta_v, xi_v) ;
		lon_v:long_name = "longitude of V-points" ;
		lon_v:units = "degree_east" ;
	double sustr(sms_time, eta_u, xi_u) ;
		sustr:long_name = "surface stress in xi dir" ;
		sustr:units = "N/m^2" ;
		sustr:coordinates = "lon_u lat_u" ;
	double svstr(sms_time, eta_v, xi_v) ;
		svstr:long_name = "surface stress in eta dir" ;
		svstr:units = "N/m^2" ;
		svstr:coordinates = "lon_v lat_v" ;

// global attributes:
		:type = "Gridpak file" ;
		:gridid = "combined grid" ;
		:history = "Created by \"\" on 07-May-2012 20:24:19" ;
		:CPP_options = "DCOMPLEX, DBLEPREC, NCARG_32, PLOTS," ;

I am a bit confused with the 1-D, 2-D discussion above, so could someone please clarify which variables should be in the windstress file if it is defined on a large native grid and the user wants ROMS to interpolate and rotate the windstress to ROMS grids? What exactly is meant by the 'regular lat/lon grid' for the large native grid ROMS will interpolate from? Because when I try to run with the wind stresses given at the large native grid (the file shown above), I am getting the following error --cropped for simplicity-- indicating a dimension inconsistency between my ROMS grid and windstress grid:

Code: Select all

NETCDF_CHECK_DIM - inconsistent size of dimension: xi_psi    980  696
..
..
..
ROMS/TOMS - Input error ............. exit_flag:   2


ERROR: Abnormal termination: NetCDF INPUT.

where 980 and 696 indicate the large native grid and ROMS grid dimensions, respectively. I might be missing something basic, as I am seeing that several users run their cases the way I want. If it is not the content and format of my windstress file that is causing this issue, could you please let me know whether there are certain CPP options that need to be defined or undefined in order to make ROMS interpolate (and rotate) winds from a large native wind field?

[ilgar]

Given the following netcdf setup of windstresses (on a large native grid covering the ROMS grid(s)) and possibly with several similar setups, ROMS does interpolate and rotate winds to its own grid(s). The main thing seems to be to introduce variable dimensions different than those in regular ROMS files.

Code: Select all

netcdf narr_windstress {
dimensions:
	xpsi = 980 ;
	xrho = 981 ;
	xu = 980 ;
	xv = 981 ;
	epsi = 760 ;
	erho = 761 ;
	eu = 761 ;
	ev = 760 ;
	time = 113 ;
variables:
	double time(time) ;
		time:long_name = "field time" ;
		time:units = "days since 1858-11-17 00:00:00 UTC" ;
		time:field = "time,scalar,series" ;
	float lon(erho, xrho) ;
		lon:long_name = "lon" ;
		lon:units = "degrees_east" ;
		lon:_FillValue = 9.999e+20 ;
		lon:missing_value = 9.999e+20 ;
		lon:field = "xp,scalar,series" ;
	float lat(erho, xrho) ;
		lat:long_name = "lat" ;
		lat:units = "degrees_north" ;
		lat:missing_value = 9.999e+20 ;
		lat:_FillValue = 9.999e+20 ;
		lat:field = "yp,scalar,series" ;
	double sustr(time, erho, xrho) ;
		sustr:long_name = "sustr" ;
		sustr:units = "N/m^2" ;
		sustr:field = "sustr,scalar,series" ;
		sustr:time = "time" ;
		sustr:coordinates = "lon lat" ;
		sustr:missing_value = 9.999e+20 ;
		sustr:_FillValue = 9.999e+20 ;
	double svstr(time, erho, xrho) ;
		svstr:long_name = "svstr" ;
		svstr:units = "N/m^2" ;
		svstr:field = "svstr,scalar,series" ;
		svstr:time = "time" ;
		svstr:coordinates = "lon lat" ;
		svstr:missing_value = 9.999e+20 ;
		svstr:_FillValue = 9.999e+20 ;

// global attributes:
		:history = "Created by \"narr_windstress_rotation\" on 15-May-2012 19:14:10" ;
		:type = "surface forcing file from COAWST_2" ;
}