[172] | 1 | ;+ |
---|
[231] | 2 | ; |
---|
[172] | 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 |
---|
[238] | 14 | ; cm_4mesh |
---|
| 15 | ; cm_4data |
---|
[172] | 16 | ; |
---|
| 17 | ; @restrictions |
---|
| 18 | ; common variable (cm_4mesh) romszinfos must be correctly defined |
---|
| 19 | ; |
---|
[226] | 20 | ; @history |
---|
[172] | 21 | ; Sept 2006 Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
| 22 | ; |
---|
| 23 | ; @version |
---|
| 24 | ; $Id$ |
---|
[238] | 25 | ; |
---|
[172] | 26 | ;- |
---|
[231] | 27 | ; |
---|
[172] | 28 | FUNCTION romsdepth |
---|
| 29 | ; |
---|
[238] | 30 | compile_opt idl2, strictarrsubs |
---|
| 31 | ; |
---|
[172] | 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 | ; |
---|
[192] | 47 | grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz |
---|
| 48 | hroms = hroms[firstx:lastx, firsty:lasty] |
---|
[226] | 49 | nt = n_elements(zeta)/nx/ny |
---|
[172] | 50 | ; |
---|
| 51 | cff1 = 1./sinh(theta_s) |
---|
| 52 | cff2 = 0.5/tanh(0.5*theta_s) |
---|
| 53 | ; |
---|
| 54 | IF type EQ 'W' THEN BEGIN |
---|
[226] | 55 | sc = (findgen(jpk)-jpk)/jpk |
---|
| 56 | ; sc = (dindgen(jpk+1)-jpk)/jpk |
---|
[172] | 57 | ; jpk = jpk+1 |
---|
| 58 | ENDIF ELSE BEGIN |
---|
[226] | 59 | sc = (findgen(jpk)-jpk-0.5+1)/jpk |
---|
[172] | 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) |
---|
[226] | 68 | ; put a z dimension to zeta |
---|
[172] | 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 |
---|