Changeset 25

Timestamp:

11/28/07 17:04:24 (16 years ago)

Author:

kolasinski

Message:

Simplification of mesh_gaussian routine by use of computegrid

File:

• trunk/procs/meshes/mesh_gaussian.pro (modified) (4 diffs)

trunk/procs/meshes/mesh_gaussian.pro

@@ -5 +5 @@
 ; init grid, sf, masks for gaussian grid for SINTEX
 ;
-   print,'T',trunc, ' gaussian grid inits.',format='(A5,I3,A20)'
+   print,' T',trunc, ' gaussian grid inits.',format='(A5,I3,A20)'
  
-; use key_shiftg if grid and data file are not organized the same
- ; key_shiftg(0) = longitude shift
- ; t106 for SQs is jpi/2
- ; key_shiftg(1) = -1 if latitude reversing needed
+   key_shift = 0
 
-   key_shift = 0
-   key_shiftg = [0, 0]
+; initialisation of character variables used in the execution of initncdf
+ ; (and computegrid as a consequence)
+   shift_txt =  ''
+   plain_txt =  ''
+
+   IF keyword_set(NO_SHIFT) THEN BEGIN
+      shift_txt = ', SHIFT = 0'
+   ENDIF
+   IF keyword_set(WHOLE_ARRAYS) THEN BEGIN
+      plain_txt = ', PLAIN = 1'
+      ;; Force YREVERSE = 0, ZREVERSE = 0, 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.
+   ENDIF
+
+   IF debug_w THEN print, ' key_yreverse = ', key_yreverse      
 
    CASE trunc OF
@@ -43 +54 @@
 
    glamt = 360.0*findgen(jpi)/float(jpi)
-   glamt = glamt#replicate(1, jpj)
 
 ; 2 Define latitudes
@@ -135 +145 @@
    ENDCASE 
 
-   gphit = replicate(1, jpi)#gphit
-
 ; 4. Define mask
 
@@ -168 +176 @@
       END 
    ENDCASE
-
-   IF key_shiftg(0) NE 0 THEN BEGIN 
-      glamt = shift(glamt,  key_shiftg(0), 0)
-      glamt(0:key_shiftg(0)-1, *) = glamt(0:key_shiftg(0)-1, *) - 360.0
-      gphit = shift(gphit,  key_shiftg(0), 0)
-      tmask = shift(tmask,  key_shiftg(0), 0)
-   ENDIF
-
-   IF key_shiftg(1) EQ -1 THEN BEGIN
-      glamt = reverse(glamt, 2)
-      gphit = reverse(gphit, 2)
-      tmask = reverse(tmask, 2)
-   ENDIF 
 ;
-; definition of key_shift
+; compute grid
 ;
-   if keyword_set(glamboundary) AND jpi GT 1 then BEGIN
-     xaxis = glamt[*, 0]
-     case 1 OF 
-       glamboundary[0] LT xaxis[0]:BEGIN
-         toobig = where(xaxis GT glamboundary[1])
-         if toobig[0] NE -1 then begin
-           key_shift = n_elements(toobig)
-           toobig = where(glamt GT glamboundary[1])
-           glamt[toobig] = glamt[toobig]-360
-         endif ELSE key_shift = 0
-       END
-       glamboundary[1] GT xaxis[jpi-1]:BEGIN
-         toosmall = where(xaxis LT glamboundary[0])
-         if toosmall[0] NE -1 then begin
-           key_shift = -n_elements(toosmall)
-           toosmall = where(glamt LT glamboundary[0])
-           glamt[toosmall] = glamt[toosmall]+360
-         ENDIF ELSE key_shift = 0
-       END
-       ELSE:key_shift = 0
-     endcase
-   ENDIF ELSE key_shift = 0
-
-   IF keyword_set(NO_SHIFT) THEN key_shift = 0
-   IF keyword_set(WHOLE_ARRAY) THEN key_shift = 0
-
-   print, '      key_shift =', key_shift
-
-   glamt = shift(glamt,key_shift, 0)
-   gphit = shift(gphit,key_shift, 0)
-   tmask = shift(tmask,  key_shift, 0)
-
-; 3 Define scale factors
-
-   zrad=6371229.0
-
-   ; Exact method: integration on a sphere
-
-   e1t = abs(2*!pi*zrad*cos(gphit*!pi/180.0)/jpi)
-
-   IF gphit(1, 0) GT gphit(0, 0) THEN BEGIN
-      ytop = (shift(gphit, 0, -1)+gphit)*0.5*!pi/180.0
-      ybot = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0
-   
-      ytop(*, jpj-1) = !pi/2
-      ybot(*,     0) = -!pi/2
-
-      e2t = abs(zrad*(ytop-ybot))
-   ENDIF ELSE BEGIN 
-      ybot = (shift(gphit, 0, -1)+gphit)*0.5*!pi/180.0
-      ytop = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0
-   
-      ybot(*, jpj-1) = !pi/2
-      ytop(*,     0) = -!pi/2
-
-      e2t = abs(zrad*(ybot-ytop))
-   ENDELSE 
-
-; 5. vertical grid (hPa)
-
-;   gdept = reverse(gdept)
-   gdepw = gdept
-
-   e3t = shift(gdept, 1)-gdept
-   e3t[0] = e3t[1]
-
-
-   key_periodique=1
-
-   glamu = glamt
-   gphiu = gphit
-   e1u = e1t
-   e2u = e2t
-
-   glamv = glamt
-   gphiv = gphit
-   e1v = e1t
-   e2v = e2t
-
-; used for coast plot
-
-   glamf = 0.5*(shift(glamt, -1, 0)+glamt)
-   glamf(jpi-1, *) =  glamf(jpi-2, *) + (glamf(jpi-2, *)-glamf(jpi-3, *))
-
-   IF gphit(1, 0) GT gphit(0, 0) THEN BEGIN
-      gphif = 0.5*(shift(gphit, 0, -1)+gphit)
-      gphif(*, jpj-1) =  gphif(*, jpj-2) + (gphif(*, jpj-2)-gphif(*, jpj-3))
-   ENDIF ELSE BEGIN
-      gphif = 0.5*(shift(gphit, 0, 1)+gphit)
-      gphif(*, 0) =  gphif(*, 1) + (gphif(*, 2)-gphif(*, 1))
-   ENDELSE 
-
-   e1f = e1t
-   e2f = e2t
-   
-   e3w = e3t
-
+   computegrid, xaxis = glamt, yaxis = gphit, mask = tmask, GLAMBOUNDARY = glamboundary,  /periodic
+ 
    key_offset = [0, 0, 0]
-
-   masked_data = 0
-   mesh_type =  'atm'
-;
-; indice i pour grille j moyenne zonale
-;
-   diaznl_idx = 1
-
-   tmask = reform(tmask, jpi*jpj)
-   tmask = reform(tmask#replicate(1, jpk), jpi, jpj, jpk)
+                
+   print,'    End of initialisation for gaussian grid T', trunc
 
   return