Changeset 247 for trunk/SRC/Interpolation

Timestamp:

05/29/07 10:58:52 (17 years ago)

Author:

smasson

Message:

add functionalities

Location:

trunk/SRC/Interpolation

Files:

• angle.pro (modified) (12 diffs)
• lbcorca.pro (added)

trunk/SRC/Interpolation/angle.pro

@@ +1 @@
+
 ;+
 ;
@@ -50 +51 @@
 ; @param fileocemesh {in}{required}{type=scalar string}
 ; a netcdf file that contains (at least) the following variables:
+;        glamt, gphit: longitudes and latitudes at T-points
 ;        glamu, gphiu: longitudes and latitudes at U-points
 ;        glamv, gphiv: longitudes and latitudes at V-points
@@ -72 +74 @@
 ; sinus of the angle between grid lines at T points and direction of the North pole
 ;
+; @param gcosf {out}{type=2d array}
+; cosinus of the angle between grid lines at F points and direction of the North pole
+;
+; @param gsinf {out}{type=2d array}
+; sinus of the angle between grid lines at F points and direction of the North pole
+;
 ; @keyword IODIRECTORY {type=scalar string}{default=''}
 ; the directory path where is located fileocemesh
 ;
-; @keyword DOUBLE {type=1 ou 2}{default=0}
+; @keyword DOUBLE {type=1 or 0}{default=0}
 ; put 1 to use double precision (default is float)
 ;
@@ -83 +91 @@
 ;       =>    gcosu[*, 0] = gcosu[*, 1]
 ;       =>    gsinu[*, 0] = gsinu[*, 1]
-;     (2) the grid follows OPA x periodicity rule, first column is
-;     equal to the next to last column
-;       =>    gcosv[0, *] = gcosv[jpj-2, *]
-;       =>    gsinv[0, *] = gsinv[jpj-2, *]
+;     (2) the grid follows OPA x and north pole periodicity rule (see <pro>lbcorca.pro</pro>)
 ;
 ; @history
@@ -92 +97 @@
 ;                   98-06 (G. Madec)
 ;       Feb 2005: IDL adaptation S. Masson
+;       May 2007: F points + call to <pro>lbcorca.pro</pro>
 ;
 ; @version
@@ -98 +104 @@
 ;-
 ;
-PRO angle, fileocemesh, gcosu, gsinu, gcosv, gsinv, gcost, gsint $
+PRO angle, fileocemesh, gcosu, gsinu, gcosv, gsinv, gcost, gsint, gcosf, gsinf $
            , IODIRECTORY = iodirectory, DOUBLE = double
 ;
@@ -137 +143 @@
   gphiv = reform(gphiv, /over)
   gphif = reform(gphif, /over)
-  jpj = (size(glamf, /dimension))[1]
+;
+  sz = size(glamf, /dimension)
+  jpi = sz[0]
+  jpj = sz[1]
 ;
 ; I. Compute the cosinus and sinus
@@ -151 +160 @@
 ;   ... north pole direction & modulous (at u-point)
   znpu = fsnspp( glamu, gphiu, DOUBLE = double )
+  znpu00 = znpu
+  znpu0i = fsnspp( shift(glamu, 0, -1), shift(gphiu, 0, -1), DOUBLE = double )
   glamu = -1 & gphiu = -1; free memory
   znpu.x = - znpu.x
@@ -163 +174 @@
   znpv.y = - znpv.y
   znnpv = znpv.x*znpv.x + znpv.y*znpv.y
-;   ... f-point
-  znpf00 = fsnspp( glamf, gphif, DOUBLE = double )
+;   ... north pole direction & modulous (at f-point)
+  znpf = fsnspp( glamf, gphif, DOUBLE = double )
+  znpf00 = znpf
   znpf01 = fsnspp( shift(glamf, 0, 1), shift(gphif, 0, 1), DOUBLE = double )
   znpf10 = fsnspp( shift(glamf, 1, 0), shift(gphif, 1, 0), DOUBLE = double )
   glamf = -1 & gphif = -1; free memory
+  znpf.x = - znpf.x
+  znpf.y = - znpf.y
+  znnpf = znpf.x*znpf.x + znpf.y*znpf.y
+;
+
 ;   ... j-direction: v-point segment direction (t-point)
   zxvvt = znpv00.x - znpv01.x
   zyvvt = znpv00.y - znpv01.y
   zmnpvt = sqrt ( temporary(znnpt) * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
-  znpv00 = -1; free memory
-  znpv01 = -1; free memory
+  znpv00 = -1 & znpv01 = -1; free memory
   IF keyword_set(double) THEN zmnpvt = 1.e-14 > zmnpvt $
   ELSE zmnpvt = 1.e-6 > zmnpvt
@@ -186 +202 @@
   zxffv = znpf00.x - znpf10.x
   zyffv = znpf00.y - znpf10.y
-  znpf00 = -1 &  znpf10 = -1; free memory
+  znpf00 = -1 & znpf10 = -1 ; free memory
   zmnpfv = sqrt ( temporary(znnpv) * ( zxffv*zxffv + zyffv*zyffv )  )
   IF keyword_set(double) THEN zmnpfv = 1.e-14 > zmnpfv $
   ELSE zmnpfv = 1.e-6 > zmnpfv
+;   ... j-direction: u-point segment direction (f-point)
+  zxuuf = znpu0i.x - znpu00.x
+  zyuuf = znpu0i.y - znpu00.y
+  zmnpuf = sqrt ( temporary(znnpf) * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
+  znpu00 = -1 & znpu0i = -1; free memory
+  IF keyword_set(double) THEN zmnpuf = 1.e-14 > zmnpuf $
+  ELSE zmnpuf = 1.e-6 > zmnpuf
+
 ;   ... cosinus and sinus using scalar and vectorial products
   gsint = ( znpt.x*zyvvt - znpt.y*zxvvt ) / zmnpvt
@@ -200 +224 @@
   gsinv =  ( znpv.x*zxffv + znpv.y*zyffv ) / zmnpfv
   gcosv = -( znpv.x*zyffv - znpv.y*zxffv ) / zmnpfv
+;   ... cosinus and sinus using scalar and vectorial products
+  gsinf = ( znpf.x*zyuuf - znpf.y*zxuuf ) / zmnpuf
+  gcosf = ( znpf.x*zxuuf + znpf.y*zyuuf ) / zmnpuf
 ;
 ; II. Geographic mesh
@@ -222 +249 @@
   gcosu[*, 0] = gcosu[*, 1]
   gsinu[*, 0] = gsinu[*, 1]
-  gcosv[0, *] = gcosv[jpj-2, *]
-  gsinv[0, *] = gsinv[jpj-2, *]
-;
+;
+  IF keyword_set(double) THEN sgn = 1.d ELSE sgn = 1.
+  dummy = lbcorca(gcost, 'T',  sgn, /correction)
+  dummy = lbcorca(gsint, 'T', -sgn, /correction)
+  dummy = lbcorca(gcosu, 'U',  sgn, /correction)
+  dummy = lbcorca(gsinu, 'U', -sgn, /correction)
+  dummy = lbcorca(gcosv, 'V',  sgn, /correction)
+  dummy = lbcorca(gsinv, 'V', -sgn, /correction)
+  dummy = lbcorca(gcosf, 'F',  sgn, /correction)
+  dummy = lbcorca(gsinf, 'F', -sgn, /correction)
+;
+
   RETURN
 END