source: trunk/SRC/Computation/msf.pro @ 445

Last change on this file since 445 was 445, checked in by smasson, 13 years ago

improve the implementation of the partial steps + bugfix in bsf

  • Property svn:keywords set to Id
File size: 5.1 KB
Line 
1;+
2;
3; @file_comments
4; compute the meridional (along the model grid) stream function in Sv using
5; the following equation: total(total(e1v*e3v*vn),1), 2, /cumulative)/1.e6
6;
7; @categories
8; diagnostics
9;
10; @param z3d {in} {type=3D xyz array}
11; meridional velocity
12;
13; @param mask2d {in} {type=2D xy array} {optional} {default=1 everywhere}
14; land/sea (0/1) mask to be applied (in addtion to the 3D mask
15; of the model) to mask some parts of the domain (and keep for example
16; the northend atlantic).
17;
18; @keyword NOSTRUCTURE {default=0}{type=scalar: 0 or 1}
19; activate if you do not want that msf returns a structure
20; but only the array referring to the field.
21;
22; @keyword INDEXBOXZOOM
23; Set this keyword to a named variable in which msf will return the x
24; and y indexes definig the zoom-box that should be used by pltz to do
25; the plot (see example)
26;
27; @keyword MASKOUT
28; Set this keyword to a named variable in which msf will return the 2d
29; array (yz) should be used by pltz to do plot the land/sea mask(see example)
30;
31; @returns {type=2D yz array}
32; meridional stream function in Sv
33;
34; @uses
35; <pro>cm_4mesh</pro>
36; <pro>cm_4data</pro>
37; <pro>grille</pro>
38; <pro>litchamp</pro>
39; <pro>fitintobox</pro>
40;
41; @restrictions
42; - T and V boxzoom parameters must be the same
43; - to be valid, mask must be equal to 0 at the bottom and on each
44;   side of the domain along x direction
45; - define the common variables (of cm_4data)
46;       varname = 'MSF'
47;       varunit = 'Sv'
48;       vargrid = 'W' ; -> should introduce WF point...
49;
50; @examples
51;
52;   IDL> initorca05
53;   IDL> file = '/Volumes/pim/F31/oce/1y/F31_1y_0040_0040_grid_V.nc'
54;   IDL> vn = read_ncdf('vomecrty', 0, 0, /timestep, file = file)
55;   IDL> mm = msf(vn, mskatl, indexboxzoom = ind, maskout = ma)
56;   IDL> pltz, mm, -20, 20, int = 1, boxzoom = [ind, 0, 5200], /xindex  $
57;            , /yindex, zoom = 1000, maskdta = ma, /no_partial, style = 'so0so'
58;
59; @history
60; Gurvan Madec, Christian Ethe, Claude Talandier and Sebastien Masson, Dec 2008
61;
62; @version
63; $Id$
64;
65;-
66;
67FUNCTION msf, z3d, mask2d, INDEXBOXZOOM = indexboxzoom, MASKOUT = maskout, NOSTRUCTURE = nostructure
68;
69  compile_opt idl2, strictarrsubs
70;
71@cm_4mesh
72@cm_4data
73;
74  CASE 1 OF
75    firstxt NE firstxv OR lastxt NE lastxv: return, report('T and V box must be defined over the same x indexes, use /memeindices when calling domdef')
76    firstyt NE firstyv OR lastyt NE lastyv: return, report('T and V box must be defined over the same y indexes, use /memeindices when calling domdef')
77    ELSE:
78  ENDCASE
79;
80  grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, grid = 'T'
81;
82  vn = litchamp(z3d)
83  IF size(vn, /n_dimensions) NE 3 THEN return, report('input data must be a 3D array')
84  vn = fitintobox(temporary(vn))
85;
86  IF n_elements(mask2d) EQ 0 THEN mask2d = replicate(1b, nx, ny)
87  IF size(mask2d, /n_dimensions) NE 2 THEN return, report('input mask2d must be a 2D array')
88  msk = fitintobox(mask2d)
89  msk = (temporary(msk))[*]#replicate(1b, nz)
90  msk = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] * temporary(msk)
91;
92  IF total(msk[*, *, nz-1]) NE 0 THEN return, report('The bottom of the box (defined by lastzt) must be set to land value (=0)')
93  IF (total(msk[0, *, *]) NE 0 OR total(msk[nx-1, *, *]) NE 0) AND NOT (keyword_set(key_periodic) AND nx EQ jpi) THEN $
94     return, report('for non-periodic domains, eastern and western bondaries must be land point')
95  IF keyword_set(key_partialstep) AND total(msk[*, ny-1, *]) NE 0 THEN flagdata = 1
96;
97  e13 = e3v_3d(/e1)
98;
99; mask the array
100  vn = temporary(vn) * msk
101  IF arg_present(maskout) THEN BEGIN
102    msk = total(temporary(msk), 1)
103    maskout = temporary(msk) NE 0
104    maskout = shift(maskout, 0, 1) ; bottom value of W should be taken into account
105    maskout[*, 0] = maskout[*, 1]
106    maskout = reverse(temporary(maskout), 2)
107  ENDIF
108
109; defaut computation:
110;  msf = 1.e-6 * total(total(temporary(e13) * temporary(vn), 1), 2, /cumulative)
111; but we force the computation from bottom to up.
112  msf = total(temporary(vn) * temporary(e13), 1) ; -> yz array
113  msf = reverse(temporary(msf), 2)
114  msf = 1.e-6 * total(temporary(msf), 2, /cumulative)
115  msf = -1.*reverse(temporary(msf), 2)
116  IF keyword_set(flagdata) THEN msf[ny-1, *] = !values.f_nan
117
118; set msf to 0 in the largest continent... no done...
119
120; update data informations
121  varname = 'MSF'
122  varunit = 'Sv'
123  vargrid = 'W' ; -> should introduce WF point...
124;
125; we must define a domain related to the section
126;
127; are we using an ORCA grid???
128
129  IF (jpiglo EQ 182 AND jpjglo EQ 149) OR (jpiglo EQ 722 AND jpjglo EQ 511) $
130     OR (jpiglo EQ 1442 AND jpjglo EQ 1021) OR (jpiglo EQ 4322 AND jpjglo EQ 3059) THEN BEGIN
131; find the i index with the northest latitude
132    ind = (where(gphit[firstx:lastx, firsty:lasty] EQ max(gphit[firstx:lastx, firsty:lasty])))[0] mod nx
133  ENDIF ELSE ind = 0
134;
135  indexboxzoom = [ind, ind, firsty, lasty]
136;
137  IF keyword_set(nostructure) THEN return, msf
138
139  IF keyword_set(key_forgetold) THEN BEGIN
140    return, {arr:temporary(msf), grid:vargrid, unit:varunit, experiment:varexp, name:varname}
141  ENDIF ELSE BEGIN
142    return, {tab:temporary(msf), grille:vargrid, unite:varunit, experience:varexp, nom:varname}
143  ENDELSE
144
145END
Note: See TracBrowser for help on using the repository browser.