Changeset 69 for trunk/Grid/computegrid.pro

Timestamp:

05/11/06 12:35:53 (18 years ago)

Author:

smasson

Message:

debug + new xxx

File:

• trunk/Grid/computegrid.pro (modified) (10 diffs)

trunk/Grid/computegrid.pro

@@ -70 +70 @@
 ;       computed by using the first line of glamt. 
 ;
+;       /PLAIN: force PERIODIC = 0, SHIFT = 0, STRIDE = [1, 1, 1] and
+;       suppress the automatic redefinition of the domain in case of 
+;       x periodicity overlap, y periodicity overlap (ORCA type only)
+;       and mask border to 0.
+;
 ;       SHIFT = scalar to force the manual definition of key_shift. By
 ;       debault, key_shift is automaticaly computed according to
 ;       glamboundary (when defined) by using the FIRST LINE of glamt. if
 ;       key_periodic=0 then in any case key_shift = 0. 
+;
+;       STRCALLING: a string containing the calling command used to
+;       call computegrid (this is used by xxx.pro)
 ;
 ;       STRIDE = : a 3 elements vector to specify the stride in x, y, z
@@ -114 +122 @@
 ; COMMON BLOCKS: cm_4mesh cm_4data cm_4cal
 ;
-; SIDE EFFECTS:
+; SIDE EFFECTS: if the grid has x/y periodicity orverlap and/or if
+;    the mask has 0 everywhere at the border (like a close sea) and
+;    if (we did not activate /plain and xminmesh, xmaxmesh, yminmesh,
+;    ymaxmesh keywords are defined to their default values), we redefine
+;    xminmesh, xmaxmesh, yminmesh, ymaxmesh in order to reove the
+;    overlapping part and/or to open the domain (avoid ti be forced
+;    to use cell_fill = 1).
 ;
 ; RESTRICTIONS:FUV points definition...
@@ -135 +149 @@
                  , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
                  , ONEARTH = onearth, PERIODIC = periodic $
-                 , SHIFT = shift, STRIDE = stride $
+                 , PLAIN = plain, SHIFT = shift, STRIDE = stride $
                  , FULLCGRID = fullcgrid, XYINDEX = xyindex $
-                 , FBASE2TBASE = fbase2tbase, _extra = ex 
+                 , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
+                 , _extra = ex 
 ;---------------------------------------------------------
 @cm_4mesh
@@ -232 +247 @@
   jpkglo = long(nz)
 ;
+; impact of plain keyword:
+;
+  IF keyword_set(plain) THEN BEGIN
+    periodic = 0
+    shift = 0
+    stride = [1, 1, 1]
+  ENDIF
+;
   IF n_elements(xminmesh) NE 0 THEN ixminmesh = long(xminmesh[0]) ELSE ixminmesh  = 0l
   IF n_elements(xmaxmesh) NE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh  = jpiglo-1
@@ -282 +305 @@
     periodic = 0
     shift = 0
-  ENDIF;
+  ENDIF
 
   r = 6371000.
@@ -358 +381 @@
   IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over) 
 ;
+;====================================================
+; check y periodicity... Only according to ORCA grid
+;==================================================== 
+; check the peridicity if iyminmesh and iymaxmesh have the default definitions...
+  IF NOT keyword_set(plain) AND key_onearth EQ 1 AND key_stride[1] EQ 1 $
+    AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 THEN BEGIN
+
+    CASE 1 OF
+      ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
+        AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN 
+; T pivot
+        ymaxmesh = -1
+        recall = 1
+      END 
+      ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
+         AND array_equal(gphit[*, jpj-1], reverse(shift(gphit[*, jpj-3], -1))) EQ 1:BEGIN 
+; T pivot
+        ymaxmesh = -1
+        recall = 1
+      END 
+      ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
+       AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN 
+; F pivot
+        ymaxmesh = -1
+        recall = 1
+      END 
+      ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
+         AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN 
+; F pivot
+        ymaxmesh = -1
+        recall = 1
+      END 
+      ELSE:
+    ENDCASE
+  ENDIF 
+;
+;====================================================
+; check x periodicity...
+;==================================================== 
+IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic) 
+;                                                    check the peridicity if ixminmesh and ixmaxmesh have the default definitions...
+  IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth EQ 1 $
+     AND key_stride[0] EQ 1 AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 AND jpi GE 3 THEN BEGIN
+    CASE 0 OF
+      total((glamt[0, *] - glamt[jpi-2, *]) MOD 360) $
+      + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN 
+        xminmesh = 1
+        xmaxmesh = -1
+        recall = 1
+      END 
+      total((glamt[0, *] - glamt[jpi-2, *]) MOD 360):BEGIN 
+        xminmesh = 1
+        recall = 1
+      END 
+      total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN 
+        xmaxmesh = -1
+        recall = 1
+      END
+      ELSE:
+    ENDCASE
+  ENDIF 
+;====================================================
+; recall computegrid if needed...
+;==================================================== 
+  IF keyword_set(recall) THEN BEGIN
+    computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
+                 , MASK = mask, GLAMBOUNDARY = glamboundary $
+                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
+                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
+                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
+                 , PERIODIC = periodic, SHIFT = shift, STRIDE = stride $
+                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $
+                 , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
+                 , _extra = ex 
+    return
+  ENDIF
 ;====================================================
 ; def of key_shift
@@ -484 +583 @@
     IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
   ENDIF ELSE glamf = glamt + 0.5
-  IF keyword_set(glamboundary) AND key_onearth THEN $
-    glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
+;
+  IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN
+    IF NOT keyword_set(glamboundary) THEN BEGIN 
+      bigger = where(glamf GE min(glamt)+360)
+      glamf[bigger] = glamf[bigger]-360.
+    ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
+  ENDIF
 ;
 ;====================================================
@@ -795 +899 @@
     ENDIF 
   ENDELSE
+;
+  IF jpi GT 2 AND jpj GT 2 AND NOT keyword_set(plain) $
+     AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
+     AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 $
+     AND total(tmask[*, 0, *]) EQ 0 AND total(tmask[*, jpj-1, *]) EQ 0 $
+     AND total(tmask[0, *, *]) EQ 0 AND total(tmask[jpi-1, *, *]) EQ 0 THEN BEGIN
+        xminmesh = 1
+        xmaxmesh = -1
+        yminmesh = 1
+        ymaxmesh = -1
+        computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
+                     , MASK = mask, GLAMBOUNDARY = glamboundary $
+                     , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
+                     , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
+                     , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
+                     , ONEARTH = onearth, PERIODIC = periodic $
+                     , PLAIN = plain, SHIFT = shift, STRIDE = stride $
+                     , FULLCGRID = fullcgrid, XYINDEX = xyindex $
+                     , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
+                     , _extra = ex 
+        return
+  ENDIF 
 ;
   IF NOT keyword_set(fullcgrid) THEN BEGIN
@@ -890 +1016 @@
 ;====================================================
 ;
-  IF total(mask) EQ jpi*jpj*jpk $
+  IF total(tmask) EQ jpi*jpj*jpk $
     AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
-  ELSE triangles_list = triangule(/keep_cont)
+  ELSE BEGIN
+; are we using ORCA2 ?
+    IF jpiglo EQ 182 AND jpi EQ 181 AND jpjglo EQ 149 AND jpj EQ 148 THEN $
+       triangles_list = triangule() ELSE triangles_list = triangule(/keep_cont)
+  ENDELSE 
 ;
 ;====================================================
@@ -898 +1028 @@
 ;====================================================
 ;
-  jpt = 1
-  time = 0
+  IF n_elements(time) EQ 0 OR n_elements(jpt) EQ 0 THEN BEGIN 
+    jpt = 1
+    time = 0
+  ENDIF
 ;
   IF NOT keyword_set(key_forgetold) THEN BEGIN
 @updateold
   ENDIF
+;====================================================
+; grid parameters used by xxx
+;====================================================
+;
+  IF NOT keyword_set(strcalling) THEN BEGIN 
+    IF n_elements(ccmeshparameters) EQ 0 THEN strcalling = 'computegrid' $
+    ELSE strcalling = ccmeshparameters.filename
+  ENDIF
+  glaminfo = moment(glamt)
+  IF finite(glaminfo[2]) EQ 0 THEN glaminfo = glaminfo[0:1]
+  gphiinfo = moment(gphit)
+  IF finite(gphiinfo[2]) EQ 0 THEN gphiinfo = gphiinfo[0:1]
+  ccmeshparameters = {filename:strcalling  $
+          , glaminfo:glaminfo $
+          , gphiinfo:gphiinfo $
+          , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
+          , jpi:jpi, jpj:jpj, jpk:jpk $
+          , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
+          , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
+          , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
+          , key_shift:key_shift, key_periodic:key_periodic $
+          , key_stride:key_stride, key_gridtype:key_gridtype $
+          , key_yreverse:key_yreverse, key_zreverse:key_zreverse $
+          , key_partialstep:key_partialstep, key_onearth:key_onearth}
+
+  ccreadparameters = {funclec_name:'read_ncdf' $
+          , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $
+          , ixmindta:ixmindta, ixmaxdta:ixmaxdta $
+          , iymindta:iymindta, iymaxdta:iymaxdta $
+          , izmindta:izmindta, izmaxdta:izmaxdta}
 ;------------------------------------------------------------
   IF keyword_set(key_performance) EQ 1 THEN $