Commit 08faa0f4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

new operator mrotuv

parent 6c6025fd
......@@ -184,6 +184,7 @@ src/Merge.c -text
src/Mergegrid.c -text
src/Mergetime.c -text
src/Merstat.c -text
src/Mrotuv.c -text
src/Ninfo.c -text
src/Nmltest.c -text
src/Output.c -text
......
/*
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2006 Uwe Schulzweida, schulzweida@dkrz.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
*/
/*
This module contains the following operators:
Mrotuv mrotuvb Backward rotation for MPIOM data
*/
#include <string.h>
#include <ctype.h>
#include "cdi.h"
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
/*
!----------------------------------------------------------------------
!
! rotation of vectors: in ocean models with rotated grids velocity
! vectors are given in the direction of grid lines and rows. they
! have to be rotated in latitudinal and longitudinal direction.
!
! note: this routine assumes positive meridional flow for a flow
! from grid point(i,j) to grid point(i,j+1) and positive
! zonal flow for a flow from grid point(i,j) to point(i+1,j).
! this is not the case for mpi-om!
!
! if this routine is used to rotate data of mpi-om, the
! logical change_sign_v needs to be true.
!j. jungclaus: 22.01.04:
!note here for the coupling fields u-i,v_j are on the non-verlapping
! (ie-2=ix) grid, furthermore, the velocity fields were previously
! interpolated onto the scalar points !
!
!h.haak: 07.10.2005 vectorisation and omp directives
!malte: use outside mpiom 02.06.2006
!----------------------------------------------------------------------
*/
void rotate_uv2(double *u_i, double *v_j, int ix, int iy,
double *lon, double *lat, double *u_lon, double *v_lat)
{
/*
real,intent(in) :: u_i(ix,iy,iz),v_j(ix,iy,iz) ! vector components in i-j-direction
real,intent(out) :: u_lon(ix,iy,iz),v_lat(ix,iy,iz) ! vector components in lon-lat direction
real,intent(in) :: lat(ix,iy),lon(ix,iy) ! latitudes and longitudes
*/
double dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j;
double lat_factor;
double absold, absnew; /* velocity vector lengths */
int i, j, ip1, im1, jp1, jm1;
int change_sign_u, change_sign_v;
double pi = 3.14159265359;
/* specification whether change in sign is needed for the input arrays */
change_sign_u = FALSE;
change_sign_v = TRUE;
/* initialization */
for ( i = 0; i < ix*iy; i++ )
{
v_lat[i] = 0;
u_lon[i] = 0;
}
/* change sign */
if ( change_sign_u )
for ( i = 0; i < ix*iy; i++ ) u_i[i] *= -1;
if ( change_sign_v )
for ( i = 0; i < ix*iy; i++ ) v_j[i] *= -1;
/* rotation */
for ( j = 0; j < iy; j++ )
for ( i = 0; i < ix; i++ )
{
ip1 = i + 1;
im1 = i - 1;
jp1 = j + 1;
jm1 = j - 1;
if ( ip1 >= ix ) ip1 = 0; /* the 0-meridian */
if ( im1 < 0 ) im1 = ix-1;
if ( jp1 >= iy ) jp1 = j; /* treatment of the last.. */
if ( jm1 < 0 ) jm1 = j; /* .. and the fist grid-row */
/* difference in latitudes */
dlat_i = lat[IX2D(j,ip1,ix)] - lat[IX2D(j,im1,ix)];
dlat_j = lat[IX2D(jp1,i,ix)] - lat[IX2D(jm1,i,ix)];
/* difference in longitudes */
dlon_i = lon[IX2D(j,ip1,ix)] - lon[IX2D(j,im1,ix)];
if ( dlon_i > pi ) dlon_i -= 2*pi;
if ( dlon_i < (-pi) ) dlon_i += 2*pi;
dlon_j = lon[IX2D(jp1,i,ix)] - lon[IX2D(jm1,i,ix)];
if ( dlon_j > pi ) dlon_j -= 2*pi;
if ( dlon_j < (-pi) ) dlon_j += 2*pi;
lat_factor = cos(lat[IX2D(j,i,ix)]);
dlon_i = dlon_i * lat_factor;
dlon_j = dlon_j * lat_factor;
/* projection by scalar product */
u_lon[IX2D(j,i,ix)] = u_i[IX2D(j,i,ix)]*dlon_i + v_j[IX2D(j,i,ix)]*dlat_i;
v_lat[IX2D(j,i,ix)] = u_i[IX2D(j,i,ix)]*dlon_j + v_j[IX2D(j,i,ix)]*dlat_j;
dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i);
dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j);
if ( fabs(dist_i) > 0 && fabs(dist_j) > 0 )
{
u_lon[IX2D(j,i,ix)] /= dist_i;
v_lat[IX2D(j,i,ix)] /= dist_j;
}
else
{
u_lon[IX2D(j,i,ix)] = 0.0;
v_lat[IX2D(j,i,ix)] = 0.0;
}
absold = sqrt(u_i[IX2D(j,i,ix)]*u_i[IX2D(j,i,ix)] + v_j[IX2D(j,i,ix)]*v_j[IX2D(j,i,ix)]);
absnew = sqrt(u_lon[IX2D(j,i,ix)]*u_lon[IX2D(j,i,ix)] + v_lat[IX2D(j,i,ix)]*v_lat[IX2D(j,i,ix)]);
/*
printf("(absold,absnew) %d %d %g %g %g %g %g %g\n", j+1, i+1, absold, absnew, u_i[IX2D(j,i,ix)], v_j[IX2D(j,i,ix)], u_lon[IX2D(j,i,ix)], v_lat[IX2D(j,i,ix)]);
*/
/* test orthogonality */
/*
if ( (dlon_i*dlon_j + dlat_j*dlat_i) > 0.1 )
fprintf(stderr, "orthogonal? %d %d %g\n", j, i, (dlon_i*dlon_j + dlat_j*dlat_i));
*/
}
}
void *Mrotuv(void *argument)
{
static char func[] = "Mrotuv";
int streamID1, streamID2, streamID3;
int nrecs, nrecs2;
int tsID, recID, levelID;
int varID1, varID2;
int nvars;
int gridID1, gridID2, gridID3;
int gridsize, gridsizex;
int nlon, nlat;
int vlistID1, vlistID2, vlistID3;
int i, j;
int taxisID1, taxisID3;
int nmiss1, nmiss2;
int code1, code2;
double missval1, missval2;
double *ufield = NULL, *vfield = NULL;
double *urfield = NULL, *vrfield = NULL;
double *uhelp = NULL, *vhelp = NULL;
double *grid1x = NULL, *grid3x = NULL, *gxhelp = NULL;
double *grid1y = NULL, *grid3y = NULL, *gyhelp = NULL;
cdoInitialize(argument);
streamID1 = streamOpenRead(cdoStreamName(0));
if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
streamID2 = streamOpenRead(cdoStreamName(1));
if ( streamID2 < 0 ) cdiError(streamID2, "Open failed on %s", cdoStreamName(1));
vlistID1 = streamInqVlist(streamID1);
vlistID2 = streamInqVlist(streamID2);
nvars = vlistNvars(vlistID1);
if ( nvars > 1 ) cdoAbort("More than one variable found in %s", cdoStreamName(0));
nvars = vlistNvars(vlistID2);
if ( nvars > 1 ) cdoAbort("More than one variable found in %s", cdoStreamName(1));
gridID1 = vlistGrid(vlistID1, 0);
gridID2 = vlistGrid(vlistID2, 0);
gridsize = gridInqSize(gridID1);
if ( gridID1 == gridID2 ) cdoAbort("Input grids are the same!");
if ( gridsize != gridInqSize(gridID2) ) cdoAbort("Grids have different size!");
if ( gridInqType(gridID1) != GRID_LONLAT &&
gridInqType(gridID1) != GRID_GAUSSIAN &&
gridInqType(gridID1) != GRID_CURVILINEAR )
cdoAbort("Grid %s unsupported!", gridNamePtr(gridInqType(gridID1)));
if ( gridInqType(gridID1) != GRID_CURVILINEAR )
gridID1 = gridToCurvilinear(gridID1);
if ( gridsize != gridInqSize(gridID1) ) cdoAbort("Internal problem: gridsize changed!");
if ( gridInqType(gridID2) != GRID_CURVILINEAR )
gridID2 = gridToCurvilinear(gridID2);
if ( gridsize != gridInqSize(gridID2) ) cdoAbort("Internal problem: gridsize changed!");
nlon = gridInqXsize(gridID1);
nlat = gridInqYsize(gridID1);
grid1x = (double *) malloc(gridsize*sizeof(double));
grid1y = (double *) malloc(gridsize*sizeof(double));
grid3x = (double *) malloc(gridsize*sizeof(double));
grid3y = (double *) malloc(gridsize*sizeof(double));
gridsizex = (nlon+2)*nlat;
gxhelp = (double *) malloc(gridsizex*sizeof(double));
gyhelp = (double *) malloc(gridsizex*sizeof(double));
gridInqXvals(gridID1, grid1x);
gridInqYvals(gridID1, grid1y);
/* load to a help field */
for ( j = 0; j < nlat; j++ )
for ( i = 0; i < nlon; i++ )
{
gxhelp[IX2D(j,i+1,nlon+2)] = grid1x[IX2D(j,i,nlon)];
gyhelp[IX2D(j,i+1,nlon+2)] = grid1y[IX2D(j,i,nlon)];
}
/* make help field cyclic */
for ( j = 0; j < nlat; j++ )
{
gxhelp[IX2D(j,0,nlon+2)] = gxhelp[IX2D(j,nlon,nlon+2)];
gxhelp[IX2D(j,nlon+1,nlon+2)] = gxhelp[IX2D(j,1,nlon+2)];
gyhelp[IX2D(j,0,nlon+2)] = gyhelp[IX2D(j,nlon,nlon+2)];
gyhelp[IX2D(j,nlon+1,nlon+2)] = gyhelp[IX2D(j,1,nlon+2)];
}
/* interpolate on pressure points */
for ( j = 0; j < nlat; j++ )
for ( i = 0; i < nlon; i++ )
{
grid3x[IX2D(j,i,nlon)] = (gxhelp[IX2D(j,i,nlon+2)]+gxhelp[IX2D(j,i+1,nlon+2)])*0.5;
if ( (gxhelp[IX2D(j,i,nlon+2)] > 340 && gxhelp[IX2D(j,i+1,nlon+2)] < 20) ||
(gxhelp[IX2D(j,i,nlon+2)] < 20 && gxhelp[IX2D(j,i+1,nlon+2)] > 340) )
{
if ( grid3x[IX2D(j,i,nlon)] < 180 )
grid3x[IX2D(j,i,nlon)] += 180;
else
grid3x[IX2D(j,i,nlon)] -= 180;
}
grid3y[IX2D(j,i,nlon)] = (gyhelp[IX2D(j,i,nlon+2)]+gyhelp[IX2D(j,i+1,nlon+2)])*0.5;
}
if ( grid1x ) free(grid1x);
if ( grid1y ) free(grid1y);
if ( gxhelp ) free(gxhelp);
if ( gyhelp ) free(gyhelp);
/*
for ( j = 0; j < nlat; j++ )
for ( i = 0; i < nlon; i++ )
printf(" %3d %3d %7.2f %7.2f\n", j+1,i+1,grid3x[IX2D(j,i,nlon)], grid3y[IX2D(j,i,nlon)]);
*/
gridID3 = gridCreate(GRID_CURVILINEAR, gridsize);
gridDefXsize(gridID3, nlon);
gridDefYsize(gridID3, nlat);
gridDefXvals(gridID3, grid3x);
gridDefYvals(gridID3, grid3y);
for ( i = 0; i < gridsize; i++ )
{
grid3x[i] *= DEG2RAD;
grid3y[i] *= DEG2RAD;
}
vlistID3 = vlistCreate();
vlistCopy(vlistID3, vlistID1);
vlistCat(vlistID3, vlistID2);
code1 = vlistInqVarCode(vlistID1, 0);
code2 = vlistInqVarCode(vlistID2, 0);
if ( code1 == code2 ) vlistDefVarCode(vlistID3, 1, code1+1);
vlistChangeGrid(vlistID3, gridID1, gridID3);
vlistChangeGrid(vlistID3, gridID2, gridID3);
taxisID1 = vlistInqTaxis(vlistID1);
taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID3);
if ( cdoVerbose ) vlistPrint(vlistID3);
streamID3 = streamOpenWrite(cdoStreamName(2), cdoFiletype());
if ( streamID3 < 0 ) cdiError(streamID3, "Open failed on %s", cdoStreamName(2));
streamDefVlist(streamID3, vlistID3);
missval1 = vlistInqVarMissval(vlistID1, 0);
missval2 = vlistInqVarMissval(vlistID2, 0);
ufield = (double *) malloc(gridsize*sizeof(double));
vfield = (double *) malloc(gridsize*sizeof(double));
urfield = (double *) malloc(gridsize*sizeof(double));
vrfield = (double *) malloc(gridsize*sizeof(double));
uhelp = (double *) malloc(gridsizex*sizeof(double));
vhelp = (double *) malloc(gridsizex*sizeof(double));
tsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
taxisCopyTimestep(taxisID3, taxisID1);
streamDefTimestep(streamID3, tsID);
nrecs2 = streamInqTimestep(streamID2, tsID);
if ( nrecs != nrecs2 ) cdoAbort("Input streams have different number of levels!");
for ( recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID1, &varID1, &levelID);
streamInqRecord(streamID2, &varID2, &levelID);
streamReadRecord(streamID1, ufield, &nmiss1);
streamReadRecord(streamID2, vfield, &nmiss2);
/* remove missing values */
if ( nmiss1 || nmiss2 )
{
for ( i = 0; i < gridsize; i++ )
{
if ( DBL_IS_EQUAL(ufield[i], missval1) ) ufield[i] = 0;
if ( DBL_IS_EQUAL(vfield[i], missval2) ) vfield[i] = 0;
}
}
/* load to a help field */
for ( j = 0; j < nlat; j++ )
for ( i = 0; i < nlon; i++ )
{
uhelp[IX2D(j,i+1,nlon+2)] = ufield[IX2D(j,i,nlon)];
vhelp[IX2D(j,i+1,nlon+2)] = vfield[IX2D(j,i,nlon)];
}
/* make help field cyclic */
for ( j = 0; j < nlat; j++ )
{
uhelp[IX2D(j,0,nlon+2)] = uhelp[IX2D(j,nlon,nlon+2)];
uhelp[IX2D(j,nlon+1,nlon+2)] = uhelp[IX2D(j,1,nlon+2)];
vhelp[IX2D(j,0,nlon+2)] = vhelp[IX2D(j,nlon,nlon+2)];
vhelp[IX2D(j,nlon+1,nlon+2)] = vhelp[IX2D(j,1,nlon+2)];
}
/* interpolate on pressure points */
for ( j = 1; j < nlat; j++ )
for ( i = 0; i < nlon; i++ )
{
ufield[IX2D(j,i,nlon)] = (uhelp[IX2D(j,i,nlon+2)]+uhelp[IX2D(j,i+1,nlon+2)])*0.5;
vfield[IX2D(j,i,nlon)] = (vhelp[IX2D(j-1,i+1,nlon+2)]+vhelp[IX2D(j,i+1,nlon+2)])*0.5;
}
for ( i = 0; i < nlon; i++ )
{
ufield[IX2D(0,i,nlon)] = 0;
vfield[IX2D(0,i,nlon)] = 0;
}
/* rotate*/
rotate_uv2(ufield, vfield, nlon, nlat, grid3x, grid3y, urfield, vrfield);
/* calc lat, lon, Auv and alpha */
/*
{
double lat, lon, auv, alpha;
for ( j = 1; j < nlat-1; j += 3 )
for ( i = 0; i < nlon; i += 3 )
{
lat = grid3y[IX2D(j,i,nlon)]*RAD2DEG;
lon = grid3x[IX2D(j,i,nlon)]*RAD2DEG;
auv = sqrt(urfield[IX2D(j,i,nlon)]*urfield[IX2D(j,i,nlon)] +
vrfield[IX2D(j,i,nlon)]*vrfield[IX2D(j,i,nlon)]);
alpha = atan2(vrfield[IX2D(j,i,nlon)], urfield[IX2D(j,i,nlon)]);
alpha = 90. - alpha*RAD2DEG;
if ( alpha < 0 ) alpha += 360.;
if ( alpha > 360 ) alpha -= 360.;
printf("%g %g %g %g\n", lon, lat, alpha, auv);
}
}
*/
streamDefRecord(streamID3, 0, levelID);
streamWriteRecord(streamID3, urfield, nmiss1);
streamDefRecord(streamID3, 1, levelID);
streamWriteRecord(streamID3, vrfield, nmiss2);
}
tsID++;
}
streamClose(streamID3);
streamClose(streamID2);
streamClose(streamID1);
if ( ufield ) free(ufield);
if ( vfield ) free(vfield);
if ( urfield ) free(urfield);
if ( vrfield ) free(vrfield);
if ( uhelp ) free(uhelp);
if ( vhelp ) free(vhelp);
if ( grid3x ) free(grid3x);
if ( grid3y ) free(grid3y);
cdoFinish();
return (0);
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment