Changeset 172 for trunk/SRC/Grid

Timestamp:

09/11/06 09:11:26 (18 years ago)

Author:

smasson

Message:

bugfix + manage roms outputs

Location:

trunk/SRC/Grid

Files:

• computegrid.pro (modified) (11 diffs)
• ncdf_meshread.pro (modified) (5 diffs)
• ncdf_meshroms.pro (added)
• romsdepth.pro (added)

trunk/SRC/Grid/computegrid.pro

@@ -4 +4 @@
 ;+
 ;
-; @file_comments compute the grid parameters from cm_4mesh common:
-;
-;       computegrid, startx, starty, stepx, stepy, nx, ny
-;       computegrid, startx, starty, stepx, stepy
-;       computegrid, xaxis = xaxis, yaxis = yaxis
-;       or a suitable mix...
-;
-; glamt
-; glamf
-; gphit
-; gphit
-; e1t
-; e2t
-; horizontal parameters
-;
-; glamu {in}
-; glamv {in}
-; gphiu {in}
-; gphiv {in}
-; e1u {in}
-; e1v {in}
-; e1f {in}
-; e2u {in}
-; e2v {in}
-; e2f {in}
-; horizontal parameters if FULLCGRID keyword is defined
-;
-; gdept
-; gdepw
-; e3t
-; e3w
-; verticals parameters
-;
+; @file_comments 
+;
+; compute the grid parameters (cm_4mesh) common
+;
+; domains sizes:
+;---------------
+; jpi, jpj, jpk, jpiglo, jpjglo, jpkglo, jpidta, jpjdta, jpkdta
+;
+; domains positions regarding to the original grid:
+;--------------------------------------------------
+; ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh
+; ixmindta, ixmaxdta, iymindta, iymaxdta, izmindta, izmaxdta
+;
+; horizontal parameters:
+; ----------------------
+; glamt, glamf, gphit, gphit, e1t, e2t
+;
+; additional horizontal parameters if FULLCGRID keyword is defined:
+; -----------------------------------------------------------------
+; glamu, glamv, gphiu, gphiv, e1u, e1v, e1f, e2u, e2v, e2f
+;
+; verticals parameters:
+; ---------------------
+; gdept, gdepw, e3t, e3w
+;
+; masks:
+; ------
 ; tmask
-; masks
-;
-; umaskred {in}
-; vmaskred {in}
-; fmaskredx {in}
-; fmaskredy {in}
-; masks if FULLCGRID keyword is defined
-;
-; triangles_list
+;
+; additional masks if FULLCGRID keyword is defined:
+; -------------------------------------------------
+; umaskred, vmaskred, fmaskredx, fmaskredy
+;
+; triangles_list:
+; ---------------
 ; triangulation
+;
+; key_* variables:
+; ----------------
+; key_onearth, key_periodic, key_shift, key_stride, key_partialstep,
+; key_yreverse, key_zreverse, key_gridtype
+;
+; xxx related variables:
+; ----------------------
+; ccmeshparameters, ccreadparameters
 ;
 ; @categories 
 ; Grid
 ;
-; @param startx {in}{required} scalar, x starting point
-; @param starty {in}{required} scalar, y starting point
-; @param stepxin {in}{required} scalar or vector: x direction step, must be > 0
-;             if vector nx is not used
-; @param stepyin {in}{required} scalar or vector: y direction step,
-;             could be > 0 (south to north) or < 0 (north to south)
-;             if vector ny is not used
-; @param nxin {in}{required} scalar, number of points in x direction
-; @param nyin {in}{required} scalar, number of points in y direction
-;
-; @keyword FULLCGRID activate to specify that you want to compute
-;       all the parameters of a C grid. Computation of glam[uv],
-;       gphi[uv], e1[uvf], e2[uvf], [uv]maskred and fmaskred[xy]
-;       will be add to the default computations
-;
-; @keyword GLAMBOUNDARY a 2 elements vector, [lon1,lon2], the longitude
-;       boundaries that should be used to visualize the data.
-;       we must have lon2 > lon1 and lon2 - lon1 le 360
-;       key_shift will be defined automaticaly computed according to
-;       glamboundary by using the FIRST LINE of glamt but
-;       key_shift will /= 0 only if key_periodic = 1
-;
-; @keyword MASK to specify the mask with a 2 or 3 dimension array
-;
-; @keyword ONEARTH {default= key_onearth=1}
-;       = 0 or 1 to force the manual definition of
-;       key_onearth (to specify if the data are on earth -> use longitude
-;       /latitude etc...).
-;       note that ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0,
-;       and is cancelling GLAMBOUNDARY
-;
-; @keyword PERIODIC {default=key_periodic is automatically computed by using the first line of glamt}
-;       = 0 or 1 to force the manual definition of key_periodic.
-;
-; @keyword PLAIN force PERIODIC = 0, SHIFT = 0, STRIDE = [1, 1, 1] and
+; @param startx {in}{optional}{type=scalar}
+;       x starting point, optional if [XY]AXIS keyword is used
+;
+; @param starty {in}{optional}{type=scalar}
+;       y starting point, optional if [XY]AXIS keyword is used
+;
+; @param stepxin {in}{optional}{type=scalar or vector}
+;       x direction step, optional if [XY]AXIS keyword is used, must be > 0
+;       if stepxin is a vector nx is not used
+;
+; @param stepyin {in}{optional}{type=scalar or vector}
+;       y direction step, optional if [XY]AXIS keyword is used,
+;       could be > 0 (south to north) or < 0 (north to south)
+;       if stepyin is a vector ny is not used
+;
+; @param nxin {in}{optional}{type=scalar}
+;       number of points in x direction,
+;       optional if [XY]AXIS keyword is used or stepxin is a vector
+;
+; @param nyin {in}{optional}{type=scalar} 
+;       number of points in y direction,
+;       optional if [XY]AXIS keyword is used or stepyin is a vector
+;
+; @keyword FULLCGRID {default=0}{type=scalar: 0 or 1}
+;       Activate to specify that you want to compute all the C grid parameters:
+;       definition of glam[uv], gphi[uv], e1[uvf], e2[uvf], [uv]maskred and
+;       fmaskred[xy] will be add to the default computations
+;
+; @keyword GLAMBOUNDARY {default=those defined in the file}{type=2 elements vector}
+;       Longitude boundaries that should be used to visualize the data.
+;         lon2 > lon1
+;         lon2 - lon1 le 360
+;       By default, the common (cm_4mesh) variable key_shift will be automatically
+;       defined according to GLAMBOUNDARY.
+;
+; @keyword MASK {default=array of 1}{type=2D or 3D array}
+;       Specify the land(0)/sea(1) mask
+;
+; @keyword ONEARTH {default=1}{type=scalar: 0 or 1}
+;       Force the manual definition of data localization on the earth or not
+;          0) if the data are not on the earth
+;          1) if the data are on earth (in that case we can for example use
+;             the labels 'longitude', 'latitude' in plots).
+;       The resulting value will be stored in the common (cm_4mesh) variable key_onearth
+;       ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0 and is cancelling GLAMBOUNDARY
+;
+; @keyword PERIODIC {default=computed by using the first line of glamt}{type=scalar: 0 or 1}
+;       Force the manual definition of the grid zonal periodicity.
+;       The resulting value will be stored in the common (cm_4mesh) variable key_periodic
+;       PERIODIC = 0 forces SHIFT = 0
+;
+; @keyword PLAIN {default=0}{type=scalar: 0 or 1}
+;       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.
 ;
-; @keyword SHIFT = scalar to force the manual definition of key_shift. By
-;       default, key_shift is automatically computed according to
-;       glamboundary (when defined) by using the FIRST LINE of glamt. if
-;       key_periodic=0 then in any case key_shift = 0.
+; @keyword SHIFT {default=computed according to glamboundary}{type=scalar}
+;       Force the manual definition of the zonal shift that must be apply to the data. 
+;       The resulting value will be stored in the common (cm_4mesh) variable key_shift
+;       Note that if key_periodic=0 then in any case key_shift = 0.
 ;
 ; @keyword STRCALLING {type=string}
@@ -100 +119 @@
 ;       call computegrid (this is used by xxx.pro)
 ;
-; @keyword STRIDE {default=[1, 1, 1]} a 3 elements vector to specify the stride in x, y, z
-;       direction. The resulting value
-;       will be stored in the common (cm_4mesh) variable key_stride
-;
-; @keyword XAXIS to specify longitude1 with a 1 or 2 dimension array, in
-;       this case startx, stepx and nx are not used but could be
-;       necessary if the y axis is not defined with yaxis. It must be
-;       possible to sort the first line of xaxis in the increasing
-;       order by shifting its elements.
-;
-; @keyword YAXIS to specify latitudes with a 1 or 2 dimension array, in
-;       this case starty, stepy and ny are not used but starty and
-;       stepy could be necessary if the x axis is not defined with xaxis.
-;       It must be sorted in the increasing or deceasing order
-;       (along each column if 2d array).
+; @keyword STRIDE {default=[1, 1, 1]}{type=3 elements vector}
+;       Specify the stride in x, y and z direction. The resulting
+;       value will be stored in the common (cm_4mesh) variable key_stride
+;
+; @keyword XAXIS {type=1D or 2D array}
+;       Specify longitudes in this case startx, stepx and nx are not used but
+;       could be necessary if the y axis is not defined with yaxis. It must be
+;       possible to sort the first line of xaxis in the increasing order by
+;       shifting its elements.
+;
+; @keyword YAXIS {type=1D or 2D array}
+;       Specify latitudes in this case starty, stepy and ny are not used but
+;       starty and stepy could be necessary if the x axis is not defined with
+;       xaxis. It must be sorted in the increasing or deceasing order (along each column if 2d array).
 ;
 ; @keyword XYINDEX activate to specify that the horizontal grid should
@@ -121 +139 @@
 ;       using this keyword forces key_onearth=0
 ;
-; @keyword XMINMESH {default=0L}
-; @keyword YMINMESH {default=0L}
-; @keyword ZMINMESH {default=0L}
-;       to define the common variables i[xyz]minmesh
-;       used to define the grid only in a zoomed part of the original
-;       grid. max value is [XYZ]MAXMESH
-;
-; @keyword XMAXMESH {default=jpiglo-1}
-; @keyword YMAXMESH {default=jpjglo-1}
-; @keyword ZMAXMESH {default=jpkglo-1}
-;       to define the common variables i[xyz]maxmesh
-;       used to define the grid only in a zoomed part of the original
-;       grid. max value is jp[ijk]glo-1.
-;       if [XYZ]MAXMESH is negative, then we define i[xyz]maxmesh as
-;       jp[ijk]glo - 1 + [XYZ]MAXMESH instead of [XYZ]MAXMESH   
+; @keyword XMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables ixminmesh used to define the localization
+;       of the first point of the grid along the x direction in a zoom of the original grid
+;
+; @keyword YMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables iyminmesh used to define the localization
+;       of the first point of the grid along the y direction in a zoom of the original grid
+;
+; @keyword ZMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables izminmesh used to define the localization
+;       of the first point of the grid along the z direction in a zoom of the original grid
+;
+; @keyword XMAXMESH {default=jpiglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables ixmaxmesh used to define the localization
+;       of the last point of the grid along the x direction in a zoom of the original grid
+;
+; @keyword YMAXMESH {default=jpjglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables iymaxmesh used to define the localization
+;       of the last point of the grid along the y direction in a zoom of the original grid
+;
+; @keyword ZMAXMESH {default=jpkglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables izmaxmesh used to define the localization
+;       of the last point of the grid along the z direction in a zoom of the original grid
 ;
 ; @keyword FBASE2TBASE
-;
-; @keyword STRCALLING
-;
-; @keyword ZAXIS to specify the vertical axis with a 1 dimension
-;       array. Must be sorted in the increasing or deceasing order
-;
-; @keyword _EXTRA used to pass your keywords to the created function.
+;       Activate when the model is a C grid based on a F point
+;       (with a F point at the bottom-left corner and a T point at the
+;       upper-right corner). In this case, we ignore 
+;         - the first line of F and U points
+;         - the first column of F and V points
+;         - the last line of T and U points
+;         - the last column of T and V points.
+;      => we are back to a C grid based on T point as for OPA model.
+;
+; @keyword ROMSH {type=2D array}
+;       This array is the final bathymetry at RHO-points. It is stored in the common
+;       variable (cm_4mesh) romszinfos.h
+;       Used when the model is a ROMS C-grid with one more point
+;       in longitude for T and V grid and one more point in latitude
+;       for T and U grid. In this case, we ignore 
+;         - the last line of T and U points
+;         - the last column of T and V points.
+;      => we are back to a C grid based on T point as for OPA model.
+;       Note that activate ROMSH forces FULLCGRID=1
+;
+; @keyword STRCALLING {type=scalar string}
+;       Used by xxx...
+;
+; @keyword ZAXIS {type=1D}
+;       Specify the vertical axis. Must be sorted in the increasing or deceasing order
+;
+; @keyword _EXTRA 
+;    (not used in the present case...)
 ;
 ; @uses cm_4mesh cm_4data cm_4cal
@@ -156 +203 @@
 ;    to use cell_fill = 1).
 ;
-; @restrictions FUV points definition...
+; @restrictions FUV points definition is not exact if the grid is irregular
 ;
 ; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
@@ -179 +226 @@
                  , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                  , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                 , _extra = ex
+                 , ROMSH = romsh, _extra = ex
 ;---------------------------------------------------------
 ;
@@ -276 +323 @@
   jpjglo = long(ny)
   jpkglo = long(nz)
+  IF keyword_set(romsh) THEN BEGIN
+    jpiglo = jpiglo - 1
+    jpjglo = jpjglo - 1
+    fullcgrid = 1
+  ENDIF
 ;
 ; impact of plain keyword:
@@ -484 +536 @@
                  , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                  , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                 , _extra = ex
+                 , ROMSH = romsh, _extra = ex
     return
   ENDIF
@@ -910 +962 @@
       gdepw = gdept
     ENDELSE
+    IF keyword_set(romsh) THEN gdepw = gdept 
 ;
 ;====================================================
@@ -926 +979 @@
 ;
 ; default mask eq 1
-  if NOT keyword_set(mask) then mask = -1
-;
-  if mask[0] NE -1 then BEGIN
-    tmask = byte(mask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh])
+  if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask
+  IF keyword_set(romsh) AND tmask[0] NE -1 THEN tmask = tmask[0:jpiglo-1, 0:jpjglo-1]
+  IF array_equal(size(tmask, /dimensions), [jpiglo, jpjglo]) EQ 1 AND jpkglo GT 1 THEN BEGIN
+    tmask = tmask[*]#replicate(1, jpkglo)
+    tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite)
+  ENDIF
+;
+  if tmask[0] NE -1 then BEGIN
+    tmask = byte(tmask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh])
     tmask = reform(tmask, jpi, jpj, jpk, /over)
     if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0)
@@ -976 +1034 @@
                      , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                      , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                     , _extra = ex
+                     , ROMSH = romsh, _extra = ex
         return
   ENDIF
@@ -1131 +1189 @@
     gphiinfo = gphit
   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} $
+  ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
+
   ccmeshparameters = {filename:strcalling  $
           , glaminfo:float(string(glaminfo, format = '(E11.4)')) $

trunk/SRC/Grid/ncdf_meshread.pro

@@ -12 +12 @@
 ; IDL> ncdf_meshread [,' filename']
 ;
-; @param filename {in}{optional}{default=meshmask.nc} 
-;    the name of the meshmask file to read. 
-;    if this name does not contain any "/" and if
-;    iodirectory keyword is not specify, then the common variable
+; @param filename {in}{optional}{default='meshmask.nc'}{type=scalar string}
+;    Name of the meshmask file to read. If this name does not contain any "/"
+;    and if iodirectory keyword is not specify, then the common variable
 ;    iodir will be use to define the mesh file path.
 ;
-; @keyword GLAMBOUNDARY a 2 elements vector, [lon1,lon2], the longitude
-;    boundaries that should be used to visualize the data.
+; @keyword GLAMBOUNDARY {default=those defined in the file}{type=2 elements vector}
+;    Longitude boundaries that should be used to visualize the data.
 ;      lon2 > lon1
 ;      lon2 - lon1 le 360
-;    key_shift will be automatically defined according to GLAMBOUNDARY.
+;    By default, the common (cm_4mesh) variable key_shift will be automatically
+;    defined according to GLAMBOUNDARY.
 ;
 ; @keyword CHECKDAT Suppressed. Use micromeshmask.pro to create an
 ;    appropriate meshmask.
 ;
-; @keyword ONEARTH {default=key_onearth=1}
-;    = 0 or 1 to force the manual definition of
-;    key_onearth (to specify if the data are on earth -> use longitude
-;    /latitude etc...).
-;    note that ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0 and is
-;    cancelling GLAMBOUNDARY
-;
-; @keyword GETDIMENSIONS
+; @keyword ONEARTH {default=1}{type=scalar: 0 or 1}
+;    Force the manual definition of data localization on the earth or not
+;       0) if the data are not on the earth
+;       1) if the data are on earth (in that case we can for example use
+;          the labels 'longitude', 'latitude' in plots).
+;    The resulting value will be stored in the common (cm_4mesh) variable key_onearth
+;    ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0 and is cancelling GLAMBOUNDARY
+;
+; @keyword GETDIMENSIONS {default=0}{type=scalar: 0 or 1}
 ;    Activate this keywords if you only want to know the dimension
 ;    of the domain stored in the mesh file. This dimensions will be
 ;    defined in jpiglo, jpjglo, jpkglo (cm_4mesh common variables)
 ;
-; @keyword PERIODIC {default=key_periodic is automatically
-;    computed by using the first line of glamt}
-;    = 0 or 1 to force the manual definition of
-;    key_periodic.
-;
-; @keyword SHIFT {default= key_shift is automatically computed according to the glamboundary}
-;    to force the manual definition of key_shift. By
-;    default, key_shift is automatically computed according to the
-;    glamboundary (when defined) by using the first line of glamt. if
-;    key_periodic=0 then in any case key_shift = 0.
-;
-; @keyword STRCALLING {type=string}
-;    a string containing the calling command used to call
-;    computegrid (this is used by xxx.pro)
-;
-; @keyword STRIDE {default=key_stride} 
-;    a 3 elements vector to specify the stride in x, y, z direction. 
-;    The resulting value will be stored in the common (cm_4mesh) variable 
-;    key_stride
-;
-; @keyword _EXTRA used to pass your keywords to the created function.
+; @keyword PERIODIC {default=computed by using the first line of glamt}{type=scalar: 0 or 1}
+;    Force the manual definition of the grid zonal periodicity.
+;    The resulting value will be stored in the common (cm_4mesh) variable key_periodic
+;    PERIODIC = 0 forces SHIFT = 0
+;
+; @keyword SHIFT {default=computed according to glamboundary}{type=scalar}
+;    Force the manual definition of the zonal shift that must be apply to the data. 
+;    The resulting value will be stored in the common (cm_4mesh) variable key_shift
+;    Note that if key_periodic=0 then in any case key_shift = 0.
+;
+; @keyword STRCALLING {type=scalar string}
+;    the calling command used to call computegrid (this is used by xxx.pro)
+;
+; @keyword STRIDE {default=[1, 1, 1]}{type=3 elements vector}
+;    Specify the stride in x, y and z direction. The resulting
+;    value will be stored in the common (cm_4mesh) variable key_stride
+;
+; @keyword _EXTRA used to pass isafile keywords
 ;
 ; @uses cm_4mesh cm_4data cm_4cal
-;
-; @restrictions define and/or use common variables from
-;               cm_4mesh, cm_4data, cm_4cal
 ;
 ; @restrictions
@@ -160 +155 @@
   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
 ; definition of jpi,jpj,jpj
-  jpi    = long(ixmaxmesh-ixminmesh+1)
-  jpj    = long(iymaxmesh-iyminmesh+1)
-  jpk    = long(izmaxmesh-izminmesh+1)
+  jpi = long(ixmaxmesh-ixminmesh+1)
+  jpj = long(iymaxmesh-iyminmesh+1)
+  jpk = long(izmaxmesh-izminmesh+1)
 ;-------------------------------------------------------
 ; check onearth and its consequences
@@ -315 +310 @@
     name = varcontient.name
 @read_ncdf_varget
-    command = namevar[i]+'=float(res)'
+    command = namevar[i]+'=float(temporary(res))'
     nothing = execute(command)
   ENDFOR
@@ -514 +509 @@
   if keyword_set(glamboundary) AND key_onearth then BEGIN
     if glamboundary[0] NE glamboundary[1] then BEGIN
-      glamt = glamt MOD 360
+      glamt = temporary(glamt) MOD 360
       smaller = where(glamt LT glamboundary[0])
       if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
       bigger = where(glamt GE glamboundary[1])
       if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
-      glamu = glamu MOD 360
+      glamu = temporary(glamu) MOD 360
       smaller = where(glamu LT glamboundary[0])
       if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
       bigger = where(glamu GE glamboundary[1])
       if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
-      glamv = glamv MOD 360
+      glamv = temporary(glamv) MOD 360
       smaller = where(glamv LT glamboundary[0])
       if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
       bigger = where(glamv GE glamboundary[1])
       if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
-      glamf = glamf MOD 360
+      glamf = temporary(glamf) MOD 360
       smaller = where(glamf LT glamboundary[0])
       if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
@@ -595 +590 @@
     gphiinfo = gphit
   ENDELSE
+  romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
   ccmeshparameters = {filename:strcalling  $
           , glaminfo:float(string(glaminfo, format = '(E11.4)')) $