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 | ; |
---|
67 | FUNCTION 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 | ; |
---|
98 | IF keyword_set(key_partialstep) THEN BEGIN |
---|
99 | |
---|
100 | ; Rebuild the T-point 3D partial steps array from 2D bottom e3t_ps |
---|
101 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
---|
102 | ; level of the bottom of the ocean |
---|
103 | bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) |
---|
104 | bottom = 0L > ( long(temporary(bottom)) - 1L ) |
---|
105 | ; the bottom of the ocean in 3D index is: |
---|
106 | bottom = lindgen(nx*ny) + temporary(bottom)*nx*ny |
---|
107 | ; 3D e3t array |
---|
108 | e3t_3D = replicate(1, nx*ny) # e3t[firstz:lastz] |
---|
109 | ; apply e3t_ps to e3t_3D at the bottom of the ocean |
---|
110 | e3t_3D[temporary(bottom)] = e3t_ps[firstx:lastx, firsty:lasty] |
---|
111 | ; 3D e3t_3D array |
---|
112 | e3t_3D = reform(e3t_3D, nx, ny, nz, /overwrite) |
---|
113 | |
---|
114 | ; Rebuild the V-point 3D partial steps array from T-point 3D e3t_3D array |
---|
115 | tmp = shift(e3t_3D, 0, -1, 0) |
---|
116 | e3v_3D = [ [ (temporary(e3t_3D))[*] ], [ (temporary(tmp))[*] ] ] |
---|
117 | e3v_3D = min(temporary(e3v_3D), dimension = 2) |
---|
118 | |
---|
119 | ; 3D e1v array |
---|
120 | e1_3D = (e1v[firstx:lastx, firsty:lasty])[*] # replicate(1., nz) |
---|
121 | ; e1v*e3v |
---|
122 | e13 = temporary(e1_3D) * temporary(e3v_3D) |
---|
123 | ENDIF ELSE BEGIN |
---|
124 | e13 = (e1v[firstx:lastx, firsty:lasty])[*] # e3t[firstz:lastz] |
---|
125 | ENDELSE |
---|
126 | ; |
---|
127 | ; mask the array |
---|
128 | vn = temporary(vn) * msk |
---|
129 | IF arg_present(maskout) THEN BEGIN |
---|
130 | msk = total(temporary(msk), 1) |
---|
131 | maskout = temporary(msk) NE 0 |
---|
132 | maskout = shift(maskout, 0, 1) ; bottom value of W should be taken into account |
---|
133 | maskout[*, 0] = maskout[*, 1] |
---|
134 | maskout = reverse(temporary(maskout), 2) |
---|
135 | ENDIF |
---|
136 | ; |
---|
137 | ; defaut computation: |
---|
138 | ; msf = 1.e-6 * total(total(temporary(e13) * temporary(vn), 1), 2, /cumulative) |
---|
139 | ; but we force the computation from bottom to up. |
---|
140 | msf = total(temporary(vn) * temporary(e13), 1) ; -> yz array |
---|
141 | msf = reverse(temporary(msf), 2) |
---|
142 | msf = 1.e-6 * total(temporary(msf), 2, /cumulative) |
---|
143 | msf = -1.*reverse(temporary(msf), 2) |
---|
144 | IF keyword_set(flagdata) THEN msf[ny-1, *] = !values.f_nan |
---|
145 | |
---|
146 | ; set msf to 0 in the largest continent... no done... |
---|
147 | |
---|
148 | ; update data informations |
---|
149 | varname = 'MSF' |
---|
150 | varunit = 'Sv' |
---|
151 | vargrid = 'W' ; -> should introduce WF point... |
---|
152 | ; |
---|
153 | ; we must define a domain related to the section |
---|
154 | ; |
---|
155 | ; are we using an ORCA grid??? |
---|
156 | |
---|
157 | IF (jpiglo EQ 182 AND jpjglo EQ 149) OR (jpiglo EQ 722 AND jpjglo EQ 511) $ |
---|
158 | OR (jpiglo EQ 1442 AND jpjglo EQ 1021) OR (jpiglo EQ 4322 AND jpjglo EQ 3059) THEN BEGIN |
---|
159 | ; find the i index with the northest latitude |
---|
160 | ind = (where(gphit[firstx:lastx, firsty:lasty] EQ max(gphit[firstx:lastx, firsty:lasty])))[0] mod nx |
---|
161 | ENDIF ELSE ind = 0 |
---|
162 | ; |
---|
163 | indexboxzoom = [ind, ind, firsty, lasty] |
---|
164 | ; |
---|
165 | IF keyword_set(nostructure) THEN return, msf |
---|
166 | |
---|
167 | IF keyword_set(key_forgetold) THEN BEGIN |
---|
168 | return, {arr:temporary(msf), grid:vargrid, unit:varunit, experiment:varexp, name:varname} |
---|
169 | ENDIF ELSE BEGIN |
---|
170 | return, {tab:temporary(msf), grille:vargrid, unite:varunit, experience:varexp, nom:varname} |
---|
171 | ENDELSE |
---|
172 | |
---|
173 | END |
---|