1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ; compute depth of ROMS outputs according to ROMS parameters |
---|
5 | ; stored in the common variable (cm_4mesh) romszinfos |
---|
6 | ; |
---|
7 | ; @categories |
---|
8 | ; Grid |
---|
9 | ; |
---|
10 | ; @returns |
---|
11 | ; the depth of the points (or -1 if error) |
---|
12 | ; |
---|
13 | ; @uses |
---|
14 | ; cm_4mesh |
---|
15 | ; cm_4data |
---|
16 | ; |
---|
17 | ; @restrictions |
---|
18 | ; common variable (cm_4mesh) romszinfos must be correctly defined |
---|
19 | ; |
---|
20 | ; @history |
---|
21 | ; Sept 2006 Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
22 | ; |
---|
23 | ; @version |
---|
24 | ; $Id$ |
---|
25 | ; |
---|
26 | ;- |
---|
27 | ; |
---|
28 | FUNCTION romsdepth |
---|
29 | ; |
---|
30 | compile_opt idl2, strictarrsubs |
---|
31 | ; |
---|
32 | @cm_4mesh |
---|
33 | @cm_4data |
---|
34 | ; |
---|
35 | theta_s = romszinfos.theta_s |
---|
36 | IF theta_s EQ -1 THEN return, -1 |
---|
37 | theta_b = romszinfos.theta_b |
---|
38 | IF theta_b EQ -1 THEN return, -1 |
---|
39 | hc = romszinfos.hc |
---|
40 | IF hc EQ -1 THEN return, -1 |
---|
41 | hroms = romszinfos.h |
---|
42 | IF hroms[0] EQ -1 THEN return, -1 |
---|
43 | zeta = romszinfos.zeta |
---|
44 | IF zeta[0] EQ -1 THEN return, -1 |
---|
45 | type = vargrid |
---|
46 | ; |
---|
47 | grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz |
---|
48 | hroms = hroms[firstx:lastx, firsty:lasty] |
---|
49 | nt = n_elements(zeta)/nx/ny |
---|
50 | ; |
---|
51 | cff1 = 1./sinh(theta_s) |
---|
52 | cff2 = 0.5/tanh(0.5*theta_s) |
---|
53 | ; |
---|
54 | IF type EQ 'W' THEN BEGIN |
---|
55 | sc = (findgen(jpk)-jpk)/jpk |
---|
56 | ; sc = (dindgen(jpk+1)-jpk)/jpk |
---|
57 | ; jpk = jpk+1 |
---|
58 | ENDIF ELSE BEGIN |
---|
59 | sc = (findgen(jpk)-jpk-0.5+1)/jpk |
---|
60 | ENDELSE |
---|
61 | ; |
---|
62 | cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5) |
---|
63 | cff = hc*(sc-cs) |
---|
64 | cff1 = cs |
---|
65 | ; |
---|
66 | hinv = 1./hroms |
---|
67 | hinv = hinv[*]#replicate(1., jpk) |
---|
68 | ; put a z dimension to zeta |
---|
69 | zeta = transpose(temporary(zeta)) |
---|
70 | zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite) |
---|
71 | zeta = transpose(temporary(zeta), [2, 1, 3, 0]) |
---|
72 | |
---|
73 | z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk)) |
---|
74 | z = (z0[*]#replicate(1., nt)) + temporary(zeta) * ((1.+z0*hinv)[*]#replicate(1., nt)) |
---|
75 | z = reform(z, nx, ny, jpk, nt, /overwrite) |
---|
76 | z = -1*reverse(temporary(z), 3) |
---|
77 | |
---|
78 | return, z |
---|
79 | end |
---|