Changeset 399 for trunk/SRC/Grid/computegrid.pro

Timestamp:

07/02/09 10:38:54 (15 years ago)

Author:

smasson

Message:

wrong bugfix in previous changeset:398

File:

• trunk/SRC/Grid/computegrid.pro (modified) (7 diffs)

trunk/SRC/Grid/computegrid.pro

@@ -353 +353 @@
 ;====================================================
 ;
-  jpiglo = long(nx)
-  jpjglo = long(ny)
-  jpkglo = long(nz)
-  IF keyword_set(romsh) THEN BEGIN
-    jpiglo = jpiglo - 1
-    jpjglo = jpjglo - 1
-    fullcgrid = 1
-  ENDIF
+  IF keyword_set(romsh) THEN fullcgrid = 1
+  CASE 1 OF
+    keyword_set(fbase2tbase):key_gridtype = 'c_f'
+    keyword_set(ubase2tbase):key_gridtype = 'c_u'
+    keyword_set(vbase2tbase):key_gridtype = 'c_v'
+    else:key_gridtype = 'c'
+  ENDCASE
+  IF strlen(key_gridtype) EQ 3 THEN fullcgrid = 1
+;
+  IF n_elements(xminmesh) NE 0 AND n_elements(xmaxmesh) NE 0 THEN BEGIN
+    IF nx EQ jpi AND xminmesh EQ ixminmesh AND xmaxmesh EQ ixmaxmesh THEN xalreadycut = 1
+  ENDIF
+  IF keyword_set(xalreadycut) THEN BEGIN
+    xmin = 0 
+    xmax = jpi - 1 
+    nxx = jpi
+  ENDIF ELSE BEGIN
+    jpiglo = long(nx)
+    IF keyword_set(romsh) THEN jpiglo = jpiglo - 1    
+    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
+    IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh
+    ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1)
+    ixminmesh = 0 > ixminmesh < ixmaxmesh
+    jpi = ixmaxmesh-ixminmesh+1
+    xmin = ixminmesh 
+    xmax = ixmaxmesh
+    nxx = jpiglo
+  ENDELSE
+
+  IF n_elements(yminmesh) NE 0 AND n_elements(ymaxmesh) NE 0 THEN BEGIN
+    IF ny EQ jpj AND yminmesh EQ iyminmesh AND ymaxmesh EQ iymaxmesh THEN yalreadycut = 1
+  ENDIF
+  IF keyword_set(yalreadycut) THEN BEGIN
+    ymin = 0 
+    ymax = jpj - 1 
+    nyy = jpj
+  ENDIF ELSE BEGIN 
+    jpjglo = long(ny)
+    IF keyword_set(romsh) THEN jpjglo = jpjglo - 1
+    IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh  = 0l
+    IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh  = jpjglo-1
+    IF key_gridtype EQ 'c_v' OR key_gridtype EQ 'c_f' THEN iymaxmesh = iymaxmesh-1
+    IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh
+    iymaxmesh = 0 > iymaxmesh < (jpjglo-1)
+    iyminmesh = 0 > iyminmesh < iymaxmesh
+    jpj = iymaxmesh-iyminmesh+1
+    ymin = iyminmesh 
+    ymax = iymaxmesh
+    nyy = jpjglo
+  ENDELSE
+
+  IF n_elements(zminmesh) NE 0 AND n_elements(zmaxmesh) NE 0 THEN BEGIN
+    IF nz EQ jpk AND zminmesh EQ izminmesh AND zmaxmesh EQ izmaxmesh THEN zalreadycut = 1
+  ENDIF
+  IF keyword_set(zalreadycut) THEN BEGIN
+    zmin = 0 
+    zmax = jpk - 1 
+    nzz = jpk
+  ENDIF ELSE BEGIN 
+   jpkglo = long(nz)
+    IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh  = 0l
+    IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh  = jpkglo-1
+    IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh
+    izmaxmesh = 0 > izmaxmesh < (jpkglo-1)
+    izminmesh = 0 > izminmesh < izmaxmesh
+    jpk = izmaxmesh-izminmesh+1
+    zmin = izminmesh 
+    zmax = izmaxmesh
+    nzz = jpkglo
+  ENDELSE
 ;
 ; impact of plain keyword:
@@ -372 +435 @@
   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
-  IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh  = 0l
-  IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh  = jpjglo-1
-  IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh  = 0l
-  IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh  = jpkglo-1
-;
-  CASE 1 OF
-    keyword_set(fbase2tbase):key_gridtype = 'c_f'
-    keyword_set(ubase2tbase):key_gridtype = 'c_u'
-    keyword_set(vbase2tbase):key_gridtype = 'c_v'
-    else:key_gridtype = 'c'
-  ENDCASE
-  IF key_gridtype EQ 'c_v' OR key_gridtype EQ 'c_f' THEN BEGIN
-    iymaxmesh = iymaxmesh-1
-  ENDIF
-  IF strlen(key_gridtype) EQ 3 THEN fullcgrid = 1
-;
-  IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh
-  IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh
-  IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh
 ; avoid basics errors...
-  ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1)
-  ixminmesh = 0 > ixminmesh < ixmaxmesh
-  iymaxmesh = 0 > iymaxmesh < (jpjglo-1)
-  iyminmesh = 0 > iyminmesh < iymaxmesh
-  izmaxmesh = 0 > izmaxmesh < (jpkglo-1)
-  izminmesh = 0 > izminmesh < izmaxmesh
-;
-  jpi = ixmaxmesh-ixminmesh+1
-  jpj = iymaxmesh-iyminmesh+1
-  jpk = izmaxmesh-izminmesh+1
 ;
   jpidta = jpiglo
@@ -468 +500 @@
 ; force glamt to have 2 dimensions
 ;
-  CASE size(reform(glamt), /n_dimensions) OF
-    0:glamt = replicate(glamt, jpi, jpj)
-    1:glamt = glamt[ixminmesh:ixmaxmesh]#replicate(1, jpj)
-    2:glamt = glamt[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh]
+  IF n_elements(glamt) EQ nxx*nyy THEN glamt = reform(glamt, nxx, nyy, /over) $
+  ELSE glamt = reform(glamt, nxx, /over)
+  CASE size(glamt, /n_dimensions) OF
+    1:BEGIN 
+      IF n_elements(glamt) EQ 1 THEN glamt = replicate(glamt[0], jpi, jpj) $
+      ELSE glamt = glamt[xmin:xmax]#replicate(1, jpj)
+    END
+    2:glamt = glamt[xmin:xmax, ymin:ymax]
   ENDCASE
 ; keep 2d array even with degenerated dimension
@@ -498 +534 @@
 ; force gphit to have 2 dimensions
 ;
-  CASE size(reform(gphit), /n_dimensions) OF
-    0:gphit = replicate(gphit, jpi, jpj)
-    1:gphit = replicate(1, jpi)#gphit[iyminmesh:iymaxmesh]
-    2:gphit = gphit[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh]
+  IF n_elements(gphit) EQ nxx*nyy THEN gphit = reform(gphit, nxx, nyy, /over) $
+  ELSE gphit = reform(gphit, nyy, /over)
+  CASE size(gphit, /n_dimensions) OF
+    1:BEGIN 
+      IF n_elements(gphit) EQ 1 THEN gphit = replicate(gphit[0], jpi, jpj) $
+      ELSE gphit = replicate(1, jpi)#gphit[ymin:ymax]
+    END
+    2:gphit = gphit[xmin:xmax, ymin:ymax]
   ENDCASE
 ; keep 2d array even with degenerated dimension
@@ -990 +1030 @@
       END
       ELSE:BEGIN
-        gdept = zaxis[izminmesh:izmaxmesh]
+        gdept = zaxis[zmin:zmax]
         IF n_elements(zreverse) EQ 0 THEN BEGIN
           IF jpk GT 1 THEN BEGIN
@@ -1028 +1068 @@
 ;
   if tmask[0] NE -1 then BEGIN
+    tmask = byte(temporary(tmask))
     IF keyword_set(romsh) THEN tmask = tmask[0:jpiglo-1, 0:jpjglo-1]
-    IF n_elements(mask) EQ jpiglo*jpjglo AND jpkglo GT 1 THEN BEGIN
-      tmask = tmask[*]#replicate(1, jpkglo)
-      tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite)
+    IF n_elements(mask) EQ nxx*nyy AND nzz GT 1 THEN BEGIN
+      tmask = tmask[*]#replicate(1b, nzz)
+      tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
     ENDIF
-    IF jpiglo EQ 1 OR jpjglo EQ 1 THEN tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite)
-    tmask = byte(tmask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh])
-    tmask = reform(tmask, jpi, jpj, jpk, /over)
+    IF nxx EQ 1 OR nyy EQ 1 OR nzz EQ 1 THEN tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
+    tmask = tmask[xmin:xmax, ymin:ymax, zmin:zmax]
+    IF jpi EQ 1 OR jpj EQ 1 OR jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
     if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0)
 ; because tmask = reverse(tmask, 2) is not working if the 3rd
@@ -1242 +1283 @@
   ENDELSE
   IF keyword_set(romsh) THEN $
-     romszinfos = {h:romsh[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh], zeta:replicate(0., jpi, jpj), theta_s:-1, theta_b:-1, hc:-1} $
+     romszinfos = {h:romsh[xmin:xmax, ymin:ymax], zeta:replicate(0., jpi, jpj), theta_s:-1, theta_b:-1, hc:-1} $
   ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}