Curvilinear grid, (Uwind,Vwind) and (sustr,svstr) angles

[lfparejar]

I'm trying to force a shallow bay with uniform, northward wind.

When I use ana_smflux and prescribe svstr=0.1/rho0 and no eastward wind, I get this (ok) pattern:




But when I use ana_winds and prescribe Vwind=10.0_r8 and no eastward wind, I get vectors that seem to get curved with the grid:



I wonder what is happening here an how can I fix Vwind to follow true north?

[jcwarner]

do you have curvgrid defined?
i think the ana_wind is doing the right thing because it looks like your grid is curved. So a 'north' wind on a curved grid is in the eta direction. The set_data routine does not seem to perform any rotation for calls to either ana_wind or ana_smflux. So i think the issue is in your ana_smflux. Do you have a section of code that does something like:
Ewind=......
Nwind=..
DO j=JstrT,JendT
DO i=IstrP,IendT
val1=0.5_r8*(angler(i-1,j)+angler(i,j))
sustr(i,j)=Ewind*COS(val1)+Nwind*SIN(val1)

..... Is there a rotation going on in the section of ana_smflux that you are using????
-j

[lfparejar]

I have no curvgrid on the project header file.

In ana_winds.h I have this (no rotation):

Code: Select all

   Wmag=13.0_r8
      DO j=JstrT,JendT
        DO i=IstrT,IendT
          Uwind(i,j)=0.0_r8
          Vwind(i,j)=10.0_r8
        END DO
      END DO

In ana_smflux.h I have:

Code: Select all

DO j=JstrT,JendT
        DO i=IstrP,IendT
          sustr(i,j)=0.0_r8
          svstr(i,j)=0.1_r8/rho0
        END DO
      END DO

I did a grep search in /Build to see if there was any interaction between winds/stresses and 'angler' but didn't find anything...

[wilkin]

You must #define CURVGRID or you won't have the correct nonlinear momentum terms that arise from non-Cartesian coordinates.

Without knowing how you made those Matlab vector plots and whether you have showing output sustr,svstr vectors or output Uwind,Vwind, it's hard to be sure where the issue lies, but to see what ROMS does with Uwind a handy trick (posted by Kate on the forum a while ago) is, at the unix $prompt ...

$ grep -r Uwind .

which turns up many things - but ignore Tangent, Adjoint etc. and focus on Utility and Nonlinear and you quickly see an entry for

Code: Select all

./Nonlinear/set_data.F:     &                     FORCES(ng)%UwindG,                           &
./Nonlinear/set_data.F:     &                     FORCES(ng)%Uwind,                            &
./Nonlinear/set_data.F:            cff1=FORCES(ng)%Uwind(i,j)*GRID(ng)%CosAngler(i,j)+         &
./Nonlinear/set_data.F:     &           FORCES(ng)%Uwind(i,j)*GRID(ng)%SinAngler(i,j)
./Nonlinear/set_data.F:            FORCES(ng)%Uwind(i,j)=cff1

If you view that code you'll see that when winds are read in from netcdf they are always rotated to ROMS coordinates, whereas if winds are set in ana_winds.h then it's on the user to do the rotation.

[lfparejar]

thanks

I defined CURVGRID and tried rotating the winds in ana_winds:

I took a look at set_data.f and found how the clockwise rotation matrix is applied. Looks something like this:

Code: Select all

cff1=Uwind*CosAngler + Vwind*SinAngler
cff2=Vwind*CosAngler - Uwind*SinAngler
            
Uwind(i,j)=cff1
Vwind(i,j)=cff2

----
So I modified ana_winds.h to call for 'angler' and wrote:

Code: Select all

u_wind=0.0_r8
      v_wind=10.0_r8
          DO j=JstrT,JendT
             DO i=IstrP,IendT
                   Uwind(i,j)=u_wind*COS(angler(i,j))+v_wind*SIN(angler(i,j))
                   Vwind(i,j)=v_wind*COS(angler(i,j))-v_wind*SIN(angler(i,j))
             END DO
          END DO

Figure from history file drawn as roms_quivergrd((Uwind(:,:)),(Vwind(:,:)),grd,20,0.01);





Magnitude increasing landward?

----

For reference, this is a figure of the grid angle from the history output (grid 'grd' obtained with roms_get_grid)
Colorbar shows angles in degrees.

[han,data]=roms_quivergrd(cos(grd.angle),sin(grd.angle),grd,20,0.1)

[wilkin]

roms_quivergrd is not the function for what you are trying to do.

roms_quivergrd plots vectors that are already in ROMS xi,eta coordinates and separated out onto their respective u,v points of the staggered Arakawa-C grid.

You are asking it to plot vectors that are centered on the common rho-points of the grid. I would have expected it to give an error in this usage, but I see in your call to roms_quivergrd you opted to decimate all vectors by a factor of 20, which made the component arrays the same dimension.

If you plot all vectors (dec = 1) the function throws the error:

Matrix dimensions must agree.
Error in roms_quivergrd (line 96)
uveitheta = (u+sqrt(-1)*v).*exp(sqrt(-1)*angle);

I guess I need to add a check of the input u,v that tests if they match the ROMS staggered dimensions for velocity.

I see what you're trying to do in your test plots, but this is not the function for it.

Did you make this grid and can say how the 'angle' was computed?

[lfparejar]

Ok, I see. The pcolor is of Umag wind speed magnitude, though.

I didn't make this grid, don't know how the angle was computed. I got the grid .nc file and configured everything else from scratch.

[wilkin]

If you have saved the output sustr,svstr you can plot those with roms_quivergrd - they are on the staggered grid and are in ROMS coordinates.

What not just give it a simple netCDF file with uniform winds - it can be a 1-D (function of time only) file and the regrid step will automatically propagate that to all the horizontal grid.

But even so, you might want to check whether the grid file angle makes sense.

[lfparejar]

I made a netdcf file from a .cdl template like this below. They're on rho points.

Not sure how to treat 'time' here. What should I do so that I can use the forcing file akin to the analytical option (without having to change the forcing file if I change NTIMES in ocean.in)
---

Code: Select all

NetCDF-3 Classic atmos_forc.nc {

dimensions:
	xrho = 1276 ;
	yrho = 184 ;
	time = 337 ;
	wind_time = 337 ;

variables:
	// Preference 'PRESERVE_FVD':  false,
	// dimensions consistent with ncBrowse, not with native MATLAB netcdf package.
	double time(time), shape = [337]
		time:long_name = "atmospheric forcing time" 
		time:units = "days" 
		time:field = "time, scalar, series" 
	double lon(yrho,xrho), shape = [184 1276]
		lon:long_name = "longitude" 
		lon:units = "degrees_east" 
		lon:field = "xp, scalar, series" 
	double lat(yrho,xrho), shape = [184 1276]
		lat:long_name = "latitude" 
		lat:units = "degrees_north" 
		lat:field = "yp, scalar, series" 
	double wind_time(wind_time), shape = [337]
		wind_time:long_name = "wind_time" 
		wind_time:units = "days" 
		wind_time:field = "Uwind_time, scalar, series" 
	double Uwind(wind_time,yrho,xrho), shape = [337 184 1276]
		Uwind:long_name = "surface u-wind component" 
		Uwind:units = "meter second-1" 
		Uwind:field = "Uwind, scalar, series" 
		Uwind:coordinates = "lon lat" 
		Uwind:time = "wind_time" 
	double Vwind(wind_time,yrho,xrho), shape = [337 184 1276]
		Vwind:long_name = "surface v-wind component" 
		Vwind:units = "meter second-1" 
		Vwind:field = "Vwind, scalar, series" 
		Vwind:coordinates = "lon lat" 
		Vwind:time = "wind_time" 

//global Attributes:
		:history = "Created by create_roms_forcings on 26-Jan-2018 10:45:20" 
		:type = "bulk fluxes forcing file"