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