;+
;
; @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
;
; 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}{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 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.
;
; @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}
; a string containing the calling command used to
; call computegrid (this is used by xxx)
;
; @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 decreasing order (along each column if 2d array).
;
; @keyword XYINDEX activate to specify that the horizontal grid should
; be simply defined by using the index of the points
; (xaxis = findgen(nx) and yaxis = findgen(ny))
; using this keyword forces key_onearth=0
;
; @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
; Note that if XMAXMESH < 0 then ixmaxmesh is defined as ixmaxmesh = jpiglo -1 + xmaxmesh
;
; @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
; Note that if YMAXMESH < 0 then iymaxmesh is defined as iymaxmesh = jpjglo -1 + ymaxmesh
;
; @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
; Note that if ZMAXMESH < 0 then izmaxmesh is defined as izmaxmesh = jpkglo -1 + maxmesh
;
; @keyword FBASE2TBASE
; 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 V points
; - the last line of T and U points
; - if the grid is not x-periodic, the first column of F and U points
; - if the grid is not x-periodic, 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 in that case, key_gridtype = 'c_f' and not 'c' (-> used in read_ncdf)
; Note that activate FBASE2TBASE forces FULLCGRID=1
;
; @keyword UBASE2TBASE
; Activate when the model is a C grid based on a U point
; (with a U point at the bottom-left corner and a T point at the
; upper-right corner). In this case, we ignore
; - if the grid is not x-periodic, the first column of F and U points
; - if the grid is not x-periodic, 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 in that case, key_gridtype = 'c_u' and not 'c' (-> used in read_ncdf)
; Note that activate UBASE2TBASE forces FULLCGRID=1
;
; @keyword VBASE2TBASE
; Activate when the model is a C grid based on a V point
; (with a V 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 V points
; - the last line of T and U points
; => we are back to a C grid based on T point as for OPA model.
; Note that in that case, key_gridtype = 'c_v' and not 'c' (-> used in read_ncdf)
; Note that activate VBASE2TBASE forces FULLCGRID=1
;
; @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 YREVERSE {default=computed according to gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]}{type=scalar}
; Force the manual definition of the y reverse that must be apply to the data.
; The resulting value will be stored in the common (cm_4mesh) variable key_yreverse
;
; @keyword ZAXIS {type=1D}
; Specify the vertical axis. Must be sorted in the increasing or decreasing order
;
; @keyword ZREVERSE {default=computed according to gdept[0] GT gdept[1]}{type=scalar}
; Force the manual definition of the z reverse that must be apply to the data.
; The resulting value will be stored in the common (cm_4mesh) variable key_zreverse
;
; @keyword _EXTRA
; not used in the present case ...
;
; @uses
; cm_4mesh
; cm_4data
; cm_4cal
;
; @restrictions
; if the grid has x/y periodicity overlap and/or if
; the mask has 0 everywhere at the border (like a closed 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 it be forced
; to use cell_fill = 1).
;
; FUV points definition is not exact if the grid is irregular
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 2000-04-20
; Sept 2004, several bug fix to suit C grid type...
; Aug 2005, rewritte almost everything...
;
; @version
; $Id$
;
;-
PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $
, XAXIS=xaxis, YAXIS=yaxis, 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 $
, YREVERSE=yreverse, ZREVERSE=zreverse $
, FULLCGRID=fullcgrid, XYINDEX=xyindex $
, UBASE2TBASE=ubase2tbase, VBASE2TBASE=vbase2tbase $
, FBASE2TBASE=fbase2tbase $
, STRCALLING=strcalling, ROMSH=romsh, _EXTRA=ex
;
compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
ENDIF
;---------------------------------------------------------
;------------------------------------------------------------
time1 = systime(1) ; for key_performance
;------------------------------------------------------------
;
;====================================================
; Check input parameters
;====================================================
;
; xaxis related parameters
;
if n_elements(xaxis) NE 0 then BEGIN
CASE (size(xaxis))[0] OF
0:nx = 1L
1:nx = (size(xaxis))[1]
2:nx = (size(xaxis))[1]
ENDCASE
ENDIF ELSE BEGIN
IF n_elements(startx) EQ 0 THEN BEGIN
dummy = report('If xaxis is not given, startx must be defined')
return
ENDIF
CASE n_elements(stepxin) OF
0:BEGIN
dummy = report('If xaxis is not given, stepxin must be defined')
return
END
1:BEGIN
IF n_elements(nxin) EQ 0 THEN BEGIN
dummy = report('If xaxis is not given and stepxin has only one element, nx must be defined')
return
ENDIF ELSE nx = nxin
END
ELSE:nx = n_elements(stepxin)
ENDCASE
ENDELSE
;
; yaxis related parameters
;
if n_elements(yaxis) NE 0 then BEGIN
CASE (size(yaxis))[0] OF
0:ny = 1L
1:ny = (size(yaxis))[1]
2:ny = (size(yaxis))[2]
ENDCASE
ENDIF ELSE BEGIN
IF n_elements(starty) EQ 0 THEN BEGIN
dummy = report('If yaxis is not given, starty must be defined')
return
ENDIF
CASE n_elements(stepyin) OF
0:BEGIN
dummy = report('If yaxis is not given, stepyin must be defined')
return
END
1:BEGIN
IF n_elements(nyin) EQ 0 THEN BEGIN
dummy = report('If yaxis is not given and stepyin has only one element, ny must be defined')
return
ENDIF ELSE ny = nyin
END
ELSE:ny = n_elements(stepyin)
ENDCASE
ENDELSE
;
; zaxis related parameters
;
if n_elements(zaxis) NE 0 then BEGIN
CASE (size(zaxis))[0] OF
0:nz = 1L
1:nz = (size(zaxis))[1]
ELSE:BEGIN
ras = report( 'not coded')
stop
END
ENDCASE
ENDIF ELSE nz = 1L
;
;====================================================
; Others automatic definitions...
;====================================================
;
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:
;
IF keyword_set(plain) THEN BEGIN
yreverse = 0
zreverse = 0
periodic = 0
shift = 0
stride = [1, 1, 1]
ENDIF
;
; avoid basics errors...
;
jpidta = jpiglo
jpjdta = jpjglo
jpkdta = jpkglo
ixmindta = 0
ixmaxdta = jpidta-1
iymindta = 0
iymaxdta = jpjdta-1
izmindta = 0
izmaxdta = jpkdta-1
;
key_partialstep = 0
if n_elements(stride) eq 3 then key_stride = stride $
ELSE key_stride = [1, 1, 1]
;
; check xyindex and its consequences
;
if keyword_set(xyindex) then onearth = 0
;
; check onearth and its consequences
;
IF n_elements(onearth) EQ 0 THEN key_onearth = 1b $
ELSE key_onearth = keyword_set(onearth)
IF NOT key_onearth THEN BEGIN
periodic = 0
shift = 0
ENDIF
r = 6371000.
;
;====================================================
; X direction : glamt
;====================================================
;
; def of glamt
;
if n_elements(xaxis) NE 0 then BEGIN
if keyword_set(xyindex) THEN glamt = findgen(jpiglo) ELSE glamt = xaxis
ENDIF ELSE BEGIN
if keyword_set(xyindex) THEN stepx = 1. ELSE stepx = stepxin
CASE 1 OF
n_elements(stepx):glamt = startx + findgen(jpiglo)*stepx
size(stepx, /n_dimensions):glamt = startx + total(stepx, /cumulative)
ELSE:BEGIN
dummy = report('Wrong definition of stepx...')
return
END
ENDCASE
ENDELSE
;
; apply glamboundary
;
IF keyword_set(glamboundary) AND key_onearth THEN BEGIN
IF glamboundary[0] GE glamboundary[1] THEN stop
IF glamboundary[1]-glamboundary[0] GT 360 THEN stop
glamt = 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.
ENDIF
;
; force glamt to have 2 dimensions
;
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
IF jpj EQ 1 THEN glamt = reform(glamt, jpi, jpj, /over)
;
;====================================================
; Y direction : gphit
;====================================================
;
; def of gphit
;
if n_elements(yaxis) NE 0 THEN BEGIN
if keyword_set(xyindex) THEN gphit = findgen(jpjglo) ELSE gphit = yaxis
ENDIF ELSE BEGIN
if keyword_set(xyindex) THEN stepy = 1. ELSE stepy = stepyin
CASE 1 OF
n_elements(stepy):gphit = starty + findgen(jpjglo)*stepy
size(stepy, /n_dimensions):gphit = starty + total(stepy, /cumulative)
ELSE:BEGIN
dummy = report('Wrong definition of stepy...')
return
END
ENDCASE
ENDELSE
;
; force gphit to have 2 dimensions
;
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
IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over)
;
;====================================================
; check y periodicity... Only according to ORCA grid
;====================================================
; check the periodicity if iyminmesh and iymaxmesh have the default definitions...
IF NOT keyword_set(plain) AND key_onearth EQ 1 $
AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 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 periodicity if ixminmesh and ixmaxmesh have the default definitions...
IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth 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 $
, STRCALLING = strcalling $
, ROMSH = romsh, _extra = ex
return
ENDIF
;
;====================================================
; def key_yreverse
;====================================================
;
IF n_elements(yreverse) EQ 0 THEN BEGIN
IF jpj GT 1 THEN BEGIN
IF total(gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]) GT jpj/2 THEN key_yreverse = 1 ELSE key_yreverse = 0
ENDIF ELSE key_yreverse = 0
ENDIF ELSE key_yreverse = yreverse
IF keyword_set(key_yreverse) THEN BEGIN
gphit = reverse(gphit, 2)
glamt = reverse(glamt, 2)
ENDIF
;
;====================================================
; def of key_shift
;====================================================
;
; definition of key_shift by shifting the array to have the min
; values of glamt[*, 0] in glamt[0, 0]
;
IF n_elements(shift) EQ 0 THEN BEGIN
IF jpi GT 1 then BEGIN
xtest = glamt[*, 0]
key_shift = (where(xtest EQ min(xtest)))[0]
IF key_shift NE 0 THEN key_shift = jpi - key_shift
ENDIF ELSE key_shift = 0
ENDIF ELSE key_shift = shift
;
;====================================================
; def of key_periodic
;====================================================
;
IF n_elements(periodic) EQ 0 THEN BEGIN
IF jpi GT 1 THEN BEGIN
xtest = shift(glamt[*, 0], key_shift)
; check that xtest is now sorted in the increasing order
IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN
ras = report(['WARNING: we cannot sort the xaxis with a simple shift...', $
'we force key_periodic = 0 and key_shift = 0', $
'only horizontal plot may be ok...'])
key_periodic = 0
xnotsorted = 1
ENDIF ELSE BEGIN
key_periodic = (xtest[jpi-1]+2*(xtest[jpi-1]-xtest[jpi-2])) $
GE (xtest[0]+360)
ENDELSE
ENDIF ELSE key_periodic = 0
ENDIF ELSE key_periodic = keyword_set(periodic)
;
; update key_shift
;
key_shift = key_shift * (key_periodic EQ 1)
;
IF (key_gridtype EQ 'c_u' OR key_gridtype EQ 'c_f') AND NOT keyword_set(key_periodic) THEN BEGIN
ixmaxmesh = ixmaxmesh-1
jpi = jpi-1
ENDIF
;
;====================================================
; apply key_shift
;====================================================
;
if keyword_set(key_shift) then BEGIN
glamt = shift(glamt, key_shift, 0)
gphit = shift(gphit, key_shift, 0)
IF jpj EQ 1 THEN BEGIN
glamt = reform(glamt, jpi, jpj, /over)
gphit = reform(gphit, jpi, jpj, /over)
ENDIF
ENDIF
;
;====================================================
; Are we using a "regular" grid (that can be described
; with x vector and y vector)?
;====================================================
;
; to get faster, we first test the most basic cases before
; testing the full array.
;
CASE 1 OF
keyword_set(xyindex):key_irregular = 0b
jpi EQ 1 OR jpj EQ 1:key_irregular = 0b
n_elements(xaxis) EQ 0 AND n_elements(yaxis) EQ 0:key_irregular = 0b
size(reform(xaxis), /n_dimensions) EQ 1 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
n_elements(xaxis) EQ 0 AND size(reform(yaxis), /n_dimensions) EQ 1:key_irregular = 0b
n_elements(yaxis) EQ 0 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
array_equal(glamt[*, 0], glamt[*, jpj-1]) EQ 0:key_irregular = 1b
array_equal(gphit[0, *], gphit[jpi-1, *]) EQ 0:key_irregular = 1b
array_equal(glamt, glamt[*, 0]#replicate(1, jpj)) EQ 0:key_irregular = 1b
array_equal(gphit, replicate(1, jpi)#(gphit[0, *])[*]) EQ 0:key_irregular = 1b
ELSE:key_irregular = 0b
ENDCASE
;
;====================================================
; def of glamf: defined as the middle of T(i,j) T(i+1,j+1)
;====================================================
;
IF jpi GT 1 THEN BEGIN
; we must compute stepxf: x distance between T(i,j) T(i+1,j+1)
CASE 1 OF
n_elements(stepx):stepxf = stepx
size(stepx, /n_dimensions):stepxf = stepx#replicate(1, jpj)
ELSE:BEGIN
if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
stepxf = (glamt + 720) MOD 360
IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over)
stepxf = shift(stepxf, -1, -1) - stepxf
stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
stepxf = min(abs(stepxf), dimension = 3)
IF NOT keyword_set(key_periodic) THEN $
stepxf[jpi-1, *] = stepxf[jpi-2, *]
ENDIF ELSE BEGIN
stepxf = shift(glamt, -1, -1) - glamt
IF keyword_set(key_periodic) THEN $
stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
ENDELSE
IF jpj GT 1 THEN BEGIN
stepxf[*, jpj-1] = stepxf[*, jpj-2]
stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
ENDIF
END
ENDCASE
glamf = glamt + 0.5 * stepxf
ENDIF ELSE glamf = glamt + 0.5
;
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
;
IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
;
;====================================================
; def of gphif: defined as the middle of T(i,j) T(i+1,j+1)
;====================================================
;
IF jpj GT 1 THEN BEGIN
; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
CASE 1 OF
n_elements(stepy):stepyf = stepy
size(stepy, /n_dimensions):stepyf = replicate(1, jpi)#stepy
ELSE:BEGIN
stepyf = shift(gphit, -1, -1) - gphit
stepyf[*, jpj-1] = stepyf[*, jpj-2]
IF jpi GT 1 THEN BEGIN
if NOT keyword_set(key_periodic) THEN $
stepyf[jpi-1, *] = stepyf[jpi-2, *]
stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
ENDIF
END
ENDCASE
gphif = gphit + 0.5 * stepyf
ENDIF ELSE gphif = gphit + 0.5
IF key_onearth THEN gphif = -90. > gphif < 90.
;
IF jpj EQ 1 THEN gphif = reform(gphif, jpi, jpj, /over)
;
;====================================================
; e1t: x distance between U(i-1,j) and U(i,j)
;====================================================
;
; *-|-*---|---*---|
;
IF jpi GT 1 THEN BEGIN
IF n_elements(stepx) NE 1 THEN BEGIN
IF keyword_set(irregular) THEN BEGIN
; we must compute stepxu: x distance between T(i,j) T(i+1,j)
IF keyword_set(key_periodic) THEN BEGIN
stepxu = (glamt + 720) MOD 360
stepxu = shift(stepxu, -1, 0) - stepxu
stepxu = [ [[stepxu]], [[stepxu + 360]], [[stepxu - 360]] ]
stepxu = min(abs(stepxu), dimension = 3)
ENDIF ELSE BEGIN
stepxu = shift(glamt, -1, 0) - glamt
stepxu[jpi-1, *] = stepxf[jpi-2, *]
ENDELSE
ENDIF ELSE stepxu = stepxf
IF jpj EQ 1 THEN stepxu = reform(stepxu, jpi, jpj, /over)
e1t = 0.5*(stepxu+shift(stepxu, 1, 0))
IF NOT keyword_set(key_periodic) THEN $
e1t[0, *] = e1t[1, *]
ENDIF ELSE e1t = replicate(stepx, jpi, jpj)
ENDIF ELSE e1t = replicate(1b, jpi, jpj)
;
IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over)
;
;====================================================
; e2t: y distance between V(i,j-1) and V(i,j)
;====================================================
;
IF jpj GT 1 THEN BEGIN
; we must compute stepyv: y distance between T(i,j) T(i,j+1)
IF n_elements(stepy) NE 1 THEN BEGIN
IF keyword_set(key_irregular) THEN BEGIN
stepyv = shift(gphit, 0, -1) - gphit
stepyv[*, jpj-1] = stepyv[*, jpj-2]
ENDIF ELSE stepyv = stepyf
e2t = 0.5*(stepyv+shift(stepyv, 0, 1))
e2t[*, 0] = e2t[*, 1]
ENDIF ELSE e2t = replicate(stepy, jpi, jpj)
ENDIF ELSE e2t = replicate(1b, jpi, jpj)
;
IF key_onearth THEN e2t = r * !pi/180. * temporary(e2t)
;
IF jpj EQ 1 THEN e2t = reform(e2t, jpi, jpj, /over)
;
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
IF keyword_set(fullcgrid) THEN BEGIN
;
;====================================================
; def of glamu: defined as the middle of T(i,j) T(i+1,j)
;====================================================
;
IF keyword_set(irregular) THEN BEGIN
glamu = glamt + 0.5 * stepxu
IF keyword_set(glamboundary) AND key_onearth THEN $
glamu = glamboundary[0] > temporary(glamu) < glamboundary[1]
ENDIF ELSE glamu = glamf
;
IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over)
;
;====================================================
; def of gphiu: defined as the middle of T(i,j) T(i+1,j)
;====================================================
;
IF jpi GT 1 THEN BEGIN
; we must compute stepyu: y distance between T(i+1,j) T(i,j)
IF keyword_set(key_irregular) THEN BEGIN
stepyu = shift(gphit, -1, 0) - gphit
IF NOT keyword_set(key_periodic) THEN $
stepyu[jpi-1, *] = stepyu[jpi-2, *]
gphiu = gphit + 0.5 * stepyu
ENDIF ELSE gphiu = gphit
ENDIF ELSE gphiu = gphit
IF key_onearth THEN gphiu = -90. > gphiu < 90.
;
IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over)
;
;====================================================
; def of glamv: defined as the middle of T(i,j) T(i,j+1)
;====================================================
;
IF jpj GT 1 THEN BEGIN
; we must compute stepxv: x distance between T(i,j) T(i,j+1)
IF keyword_set(irregular) THEN BEGIN
IF keyword_set(key_periodic) THEN BEGIN
stepxv = (glamt + 720) MOD 360
stepxv = shift(stepxv, 0, -1) - stepxv
stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ]
stepxv = min(abs(stepxv), dimension = 3)
ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt
stepxv[*, jpj-1] = stepxv[*, jpj-2]
glamv = glamt + 0.5 * stepxv
IF keyword_set(glamboundary) AND key_onearth THEN $
glamv = glamboundary[0] > temporary(glamv) < glamboundary[1]
ENDIF ELSE glamv = glamt
ENDIF ELSE glamv = glamt
;
;====================================================
; def of gphiv: defined as the middle of T(i,j) T(i,j+1)
;====================================================
;
IF keyword_set(key_irregular) THEN $
gphiv = gphit + 0.5 * stepyv $
ELSE gphiv = gphif
IF key_onearth THEN gphiv = -90. > gphiv < 90.
;
IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over)
;
;====================================================
; e1u: x distance between T(i,j) and T(i+1,j)
;====================================================
;
IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $
e1u = stepxu ELSE e1u = e1t
;
;====================================================
; e2u: y distance between F(i,j-1) and F(i,j)
;====================================================
;
IF keyword_set(key_irregular) THEN BEGIN
e2u = gphif - shift(gphif, 0, 1)
e2u[*, 0] = e2u[*, 1]
IF key_onearth THEN e2u = r * !pi/180. * temporary(e2u)
ENDIF ELSE e2u = e2t
;
IF jpj EQ 1 THEN e2u = reform(e2u, jpi, jpj, /over)
;
;====================================================
; e1v: x distance between F(i-1,j) and F(i,j)
;====================================================
;
IF keyword_set(irregular) THEN BEGIN
IF keyword_set(key_periodic) THEN BEGIN
e1v = (glamf + 720) MOD 360
e1v = e1v - shift(e1v, 1, 0)
e1v = [ [[e1v]], [[e1v + 360]], [[e1v - 360]] ]
e1v = min(abs(e1v), dimension = 3)
ENDIF ELSE BEGIN
e1v = glamf - shift(glamf, 1, 0)
e1v[0, *] = stepxf[1, *]
ENDELSE
ENDIF ELSE e1v = e1t
;
IF jpj EQ 1 THEN e1v = reform(e1v, jpi, jpj, /over)
;
;====================================================
; e2v: y distance between T(i,j) and T(i+1,j)
;====================================================
;
IF jpj GT 1 and n_elements(stepy) NE 1 THEN BEGIN
e2v = stepyv
IF key_onearth THEN e2v = r * !pi/180. * temporary(e2v)
ENDIF ELSE e2v = e2t
;
;====================================================
; e1f: x distance between V(i,j) and V(i+1,j)
;====================================================
;
IF keyword_set(irregular) THEN BEGIN
IF keyword_set(key_periodic) THEN BEGIN
e1f = (glamv + 720) MOD 360
e1f = shift(e1f, -1, 0) - e1f
e1f = [ [[e1f]], [[e1f + 360]], [[e1f - 360]] ]
e1f = min(abs(e1f), dimension = 3)
ENDIF ELSE BEGIN
e1f = shift(glamv, -1, 0) - glamt
e1f[jpi-1, *] = stepxf[jpi-2, *]
ENDELSE
ENDIF ELSE e1f = e1u
;
IF jpj EQ 1 THEN e1f = reform(e1f, jpi, jpj, /over)
;
;====================================================
; e2f: y distance between U(i,j) and U(i,j+1)
;====================================================
;
IF keyword_set(key_irregular) THEN BEGIN
e2f = shift(gphiu, 0, -1) - gphiu
e2f[*, jpj-1] = e2f[*, jpj-2]
IF key_onearth THEN e2f = r * !pi/180. * temporary(e2f)
ENDIF ELSE e2f = e2v
;
IF jpj EQ 1 THEN e2f = reform(e2f, jpi, jpj, /over)
;
ENDIF
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;
;
;====================================================
; e1[tuvf] from degree to meters
;====================================================
;
IF keyword_set(key_onearth) THEN BEGIN
e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit)
IF keyword_set(fullcgrid) THEN BEGIN
e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu)
e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv)
e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif)
ENDIF
ENDIF
;
IF jpj EQ 1 THEN BEGIN
e1t = reform(e1t, jpi, jpj, /over)
IF keyword_set(fullcgrid) THEN BEGIN
e1u = reform(e1u, jpi, jpj, /over)
e1v = reform(e1v, jpi, jpj, /over)
e1f = reform(e1f, jpi, jpj, /over)
ENDIF
ENDIF
;
;====================================================
; if not fullcgrid: make sure we don't use glam[uv], gphi[uv], e[12][uvf]
;====================================================
;
IF NOT keyword_set(fullcgrid) THEN BEGIN
glamu = !values.f_nan & glamv = !values.f_nan
gphiu = !values.f_nan & gphiv = !values.f_nan
e1u = !values.f_nan & e1v = !values.f_nan & e1f = !values.f_nan
e2u = !values.f_nan & e2v = !values.f_nan & e2f = !values.f_nan
firstxu = !values.f_nan & lastxu = !values.f_nan & nxu = !values.f_nan
firstyu = !values.f_nan & lastyu = !values.f_nan & nyu = !values.f_nan
firstxv = !values.f_nan & lastxv = !values.f_nan & nxv = !values.f_nan
firstyv = !values.f_nan & lastyv = !values.f_nan & nyv = !values.f_nan
ENDIF
;
;====================================================
; Z direction
;====================================================
;
; z axis
;
CASE n_elements(zaxis) OF
0:BEGIN
gdept = 0.
key_zreverse = 0
END
1:BEGIN
gdept = zaxis
key_zreverse = 0
END
ELSE:BEGIN
gdept = zaxis[zmin:zmax]
IF n_elements(zreverse) EQ 0 THEN BEGIN
IF jpk GT 1 THEN BEGIN
if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0
ENDIF ELSE key_zreverse = 0
ENDIF ELSE key_zreverse = zreverse
IF keyword_set(key_zreverse) THEN gdept = reverse(gdept)
END
ENDCASE
;
if n_elements(gdept) GT 1 then BEGIN
stepz = shift(gdept, -1)-gdept
stepz[jpk-1] = stepz[jpk-2]
gdepw = 0. > (gdept-stepz/2.)
ENDIF ELSE BEGIN
stepz = 1.
gdepw = gdept
ENDELSE
IF keyword_set(romsh) THEN gdepw = gdept
;
;====================================================
; e3[tw]:
;====================================================
;
e3t = stepz
IF n_elements(stepz) GT 1 THEN BEGIN
e3w = 0.5*(stepz+shift(stepz, 1))
e3w[0] = 0.5*e3t[0]
ENDIF ELSE e3w = e3t
;
;====================================================
; Mask
;====================================================
;
; default mask eq 1
if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask
;
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 nxx*nyy AND nzz GT 1 THEN BEGIN
tmask = tmask[*]#replicate(1b, nzz)
tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
ENDIF
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
; dimension of tmask = 1, we call reform.
IF jpk EQ 1 THEN tmask = reform(tmask, /over)
IF key_yreverse EQ 1 THEN tmask = reverse(tmask, 2)
IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
IF key_zreverse EQ 1 THEN tmask = reverse(tmask, 3)
IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
IF keyword_set(fullcgrid) THEN BEGIN
IF keyword_set(key_periodic) THEN BEGIN
msk = tmask*shift(tmask, -1, 0, 0)
umaskred = msk[jpi-1, *, *]
ENDIF ELSE umaskred = tmask[jpi-1, *, *]
vmaskred = tmask[*, jpj-1, *]
ENDIF
fmaskredy = tmask[jpi-1, *, *]
fmaskredx = tmask[*, jpj-1, *]
ENDIF ELSE BEGIN
tmask = replicate(1b, jpi, jpj, jpk)
IF keyword_set(fullcgrid) THEN BEGIN
umaskred = replicate(1b, jpj, jpk)
vmaskred = replicate(1b, jpi, jpk)
ENDIF
fmaskredy = replicate(1b, jpj, jpk)
fmaskredx = replicate(1b, jpi, jpk)
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 $
, UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $
, FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
, ROMSH = romsh, _extra = ex
return
ENDIF
;
IF NOT keyword_set(fullcgrid) THEN BEGIN
umaskred = !values.f_nan
vmaskred = !values.f_nan
ENDIF
;
;====================================================
; stride...
;====================================================
;
IF total(key_stride) GT 3 THEN BEGIN
IF key_shift NE 0 THEN BEGIN
; for explanation, see header of read_ncdf_varget.pro
jpiright = key_shift
jpileft = jpi - key_shift - ( (key_stride[0]-1)-((key_shift-1) MOD key_stride[0]) )
jpi = ((jpiright-1)/key_stride[0]+1) + ((jpileft-1)/key_stride[0]+1)
ENDIF ELSE jpi = (jpi-1)/key_stride[0]+1
jpj = (jpj-1)/key_stride[1]+1
jpk = (jpk-1)/key_stride[2]+1
;
glamt = (temporary(glamt))[0:*:stride[0], 0:*:stride[1]]
gphit = (temporary(gphit))[0:*:stride[0], 0:*:stride[1]]
e1t = (temporary(e1t))[0:*:stride[0], 0:*:stride[1]]
e2t = (temporary(e2t))[0:*:stride[0], 0:*:stride[1]]
tmask = (temporary(tmask))[0:*:stride[0], 0:*:stride[1], 0:*:stride[2]]
gdept = gdept[0:*:stride[2]]
gdepw = gdepw[0:*:stride[2]]
e3t = e3t[0:*:stride[2]]
e3w = e3w[0:*:stride[2]]
; we must recompute glamf and gphif...
IF jpi GT 1 THEN BEGIN
if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
stepxf = (glamt + 720) MOD 360
stepxf = shift(stepxf, -1, -1) - stepxf
stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
stepxf = min(abs(stepxf), dimension = 3)
IF NOT keyword_set(key_periodic) THEN $
stepxf[jpi-1, *] = stepxf[jpi-2, *]
ENDIF ELSE BEGIN
stepxf = shift(glamt, -1, -1) - glamt
IF keyword_set(key_periodic) THEN $
stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
ENDELSE
IF jpj GT 1 THEN BEGIN
stepxf[*, jpj-1] = stepxf[*, jpj-2]
stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
ENDIF
glamf = glamt + 0.5 * stepxf
ENDIF ELSE glamf = glamt + 0.5
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
IF jpj GT 1 THEN BEGIN
; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
stepyf = shift(gphit, -1, -1) - gphit
stepyf[*, jpj-1] = stepyf[*, jpj-2]
IF jpi GT 1 THEN BEGIN
if NOT keyword_set(key_periodic) THEN $
stepyf[jpi-1, *] = stepyf[jpi-2, *]
stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
ENDIF
gphif = gphit + 0.5 * stepyf
ENDIF ELSE gphif = gphit + 0.5
;
IF jpj EQ 1 THEN BEGIN
glamt = reform(glamt, jpi, jpj, /over)
gphit = reform(gphit, jpi, jpj, /over)
glamf = reform(glamf, jpi, jpj, /over)
gphif = reform(gphif, jpi, jpj, /over)
e1t = reform(e1t, jpi, jpj, /over)
e2t = reform(e2t, jpi, jpj, /over)
ENDIF
;
IF keyword_set(fullcgrid) THEN BEGIN
glamu = (temporary(glamu))[0:*:stride[0], 0:*:stride[1]]
gphiu = (temporary(gphiu))[0:*:stride[0], 0:*:stride[1]]
e1u = (temporary(e1u))[0:*:stride[0], 0:*:stride[1]]
e2u = (temporary(e2u))[0:*:stride[0], 0:*:stride[1]]
glamv = (temporary(glamv))[0:*:stride[0], 0:*:stride[1]]
gphiv = (temporary(gphiv))[0:*:stride[0], 0:*:stride[1]]
e1v = (temporary(e1v))[0:*:stride[0], 0:*:stride[1]]
e2v = (temporary(e2v))[0:*:stride[0], 0:*:stride[1]]
e1f = (temporary(e1f))[0:*:stride[0], 0:*:stride[1]]
e2f = (temporary(e2f))[0:*:stride[0], 0:*:stride[1]]
umaskred = (temporary(umaskred))[0, 0:*:stride[1], 0:*:stride[2]]
vmaskred = (temporary(vmaskred))[0:*:stride[0], 0, 0:*:stride[2]]
fmaskredy = (temporary(fmaskredy))[0, 0:*:stride[1], 0:*:stride[2]]
fmaskredx = (temporary(fmaskredx))[0:*:stride[0], 0, 0:*:stride[2]]
IF jpj EQ 1 THEN BEGIN
glamu = reform(glamu, jpi, jpj, /over)
gphiu = reform(gphiu, jpi, jpj, /over)
e1u = reform(e1u, jpi, jpj, /over)
e2u = reform(e2u, jpi, jpj, /over)
glamv = reform(glamv, jpi, jpj, /over)
gphiv = reform(gphiv, jpi, jpj, /over)
e1v = reform(e1v, jpi, jpj, /over)
e2v = reform(e2v, jpi, jpj, /over)
e1f = reform(e1f, jpi, jpj, /over)
e2f = reform(e2f, jpi, jpj, /over)
ENDIF
ENDIF
ENDIF
;
;====================================================
; apply all the grid parameters
;====================================================
;
@updateold
domdef
;
;====================================================
; Triangulation
;====================================================
;
IF total(tmask) EQ jpi*jpj*jpk $
AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
ELSE BEGIN
; are we using ORCA2 ?
IF jpiglo EQ 182 AND jpi EQ 180 AND jpjglo EQ 149 AND jpj EQ 148 THEN $
triangles_list = triangule() ELSE triangles_list = triangule(/keep_cont)
ENDELSE
;
;====================================================
; time axis (default definition)
;====================================================
;
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
IF n_elements(glamt) GE 2 THEN BEGIN
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]
ENDIF ELSE BEGIN
glaminfo = glamt
gphiinfo = gphit
ENDELSE
IF keyword_set(romsh) THEN $
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}
ccmeshparameters = {filename:strcalling $
, glaminfo:float(string(glaminfo, format = '(E11.4)')) $
, gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $
, 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 $
print, 'time computegrid', systime(1)-time1
;------------------------------------------------------------
return
end