Changeset 440 for trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro

Timestamp:

12/21/10 14:26:05 (13 years ago)

Author:

smasson

Message:

take into account the ssh in moyenne and grossemoyenne

File:

• trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro (modified) (8 diffs)

trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro

@@ -61 +61 @@
 ; mean even if NAN is activated (see SPATIALFIRST explanations...)
 ;
+; @keyword SSH
+; array giving the sea level elenation (in meter) that must be added
+; to the thickness of the first level when computing vertical mean
+;
 ; @keyword WDEPTH
 ; to specify that the field is at W depth instead of T
@@ -89 +93 @@
 ;-
 FUNCTION grossemoyenne, tab, direc, BOXZOOM=boxzoom, INTEGRATION=integration, KEEPBOTTOM = keepbottom $
-                      , NAN=nan, NODOMDEF=nodomdef, WDEPTH=wdepth $
+                      , NAN=nan, NODOMDEF=nodomdef, WDEPTH=wdepth, ssh=ssh $
                       , SPATIALFIRST=spatialfirst, TEMPORALFIRST=temporalfirst $
                       , _EXTRA=ex
@@ -164 +168 @@
 ;---------------------------------------------------------------
   grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth
+;---------------------------------------------------------------
+; check the use of ssh keyword...
+;---------------------------------------------------------------
+  IF keyword_set(ssh) AND dirz NE 0 AND firstz EQ 0 AND taille[0] GT 2 THEN BEGIN
+    szssh = size(ssh)
+    IF szssh[1] NE taille[1] OR szssh[2] NE taille[2] THEN $
+       return, report('ssh and input arrays must have the same size')
+    sshh = ssh
+  ENDIF
 ;------------------------------------------------------------
 ; I.4) if dirt equal 1, we make the temporal average and we send it in moyenne
 ;------------------------------------------------------------
-  if dirt EQ 1 AND NOT keyword_set(spatialfirst) then begin
+  if dirt EQ 1 AND NOT keyword_set(spatialfirst) then BEGIN
+    IF keyword_set(sshh) THEN return, report('ssh and temporal mean: case not coded...')
     if dim EQ 3d then BEGIN
       case 1 of
@@ -338 +352 @@
 ;---------------------------------------------------------------
     case 1 of
-      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $
+      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt:BEGIN 
         res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *]
-      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
+        IF keyword_set(sshh) THEN sshh = sshh[firstx:lastx, firsty:lasty, *]
+      END
+      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt:BEGIN 
         res = tab[firstx:lastx, firsty:lasty, *, *]
+        IF keyword_set(sshh) THEN sshh = sshh[firstx:lastx, firsty:lasty, *]
+      END
       taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res = tab
       taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $
@@ -455 +473 @@
         echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt)
         echelle = reform(echelle, nx, ny, nz, jpt, /over)
+        IF keyword_set(sshh) THEN echelle[*, *, 0, *] = echelle[*, *, 0, *] + sshh
         if keyword_set(integration) then divi = 1 ELSE begin
           IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $
@@ -506 +525 @@
         echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt)
         echelle = reform(echelle, nx, ny, nz, jpt, /over)
+        IF keyword_set(sshh) THEN echelle[*, *, 0, *] = echelle[*, *, 0, *] + sshh*( e1[*]#replicate(1., jpt) )
         if keyword_set(integration) then divi = 1 ELSE begin
           IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $
@@ -527 +547 @@
         echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt)
         echelle = reform(echelle, nx, ny, nz, jpt, /over)
+        IF keyword_set(sshh) THEN echelle[*, *, 0, *] = echelle[*, *, 0, *] + sshh*( e2[*]#replicate(1., jpt) )
         if keyword_set(integration) then divi = 1 ELSE begin
           IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $
@@ -548 +569 @@
         echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt)
         echelle = reform(echelle, nx, ny, nz, jpt, /over)
+        IF keyword_set(sshh) THEN echelle[*, *, 0, *] = echelle[*, *, 0, *] + sshh*( (e1[*]*e2[*])#replicate(1., jpt) )
         if keyword_set(integration) then divi = 1 ELSE begin
           IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $