Changeset 501 for trunk/SRC/Grid/romsdepth.pro

Timestamp:

09/20/16 20:15:00 (8 years ago)

Author:

smasson

Message:

set of bugfixes...

File:

• trunk/SRC/Grid/romsdepth.pro (modified) (3 diffs)

trunk/SRC/Grid/romsdepth.pro

@@ -25 +25 @@
 ;
 ;-
-FUNCTION romsdepth
+FUNCTION romsdepth, new_s_coord = new_s_coord
 ;
   compile_opt idl2, strictarrsubs
@@ -46 +46 @@
   grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
   hroms = hroms[firstx:lastx, firsty:lasty]
-  nt = n_elements(zeta)/nx/ny
-;
-  cff1 = 1./sinh(theta_s)
-  cff2 = 0.5/tanh(0.5*theta_s)
+  nt = n_elements(zeta[firstx:lastx, firsty:lasty, *])/nx/ny
 ;
   IF type EQ 'W' THEN BEGIN
@@ -59 +56 @@
   ENDELSE
 ;
-  cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5)
-  cff = hc*(sc-cs)
-  cff1 = cs
+  IF keyword_set(new_s_coord) THEN BEGIN
+
+    IF theta_s GT 0. THEN BEGIN
+      csrf = ( 1. - cosh(theta_s*sc) ) / ( cosh(theta_s) - 1. )
+    ENDIF ELSE BEGIN
+      csrf = - sc*sc
+    ENDELSE 
+    if theta_b GT 0. THEN BEGIN 
+      cs = ( exp(theta_b*csrf) - 1. ) / ( 1. - exp(-theta_b) )  
+    ENDIF ELSE BEGIN
+      cs = csrf
+    ENDELSE
+    
+    hh = (hroms[*]#replicate(1., jpk))
+    z0 = hh * ( replicate(1., nx*ny)#(hc*sc)[*] + hh * replicate(1., nx*ny)#cs[*] ) / ( hc + hh )
+    
+    
+  ENDIF ELSE BEGIN
+    
+    cff1 = 1./sinh(theta_s)
+    cff2 = 0.5/tanh(0.5*theta_s)
+
+    cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5)
+    cff = hc*(sc-cs)
+    cff1 = cs
+
+    z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))
+
+  ENDELSE
+
 ;
   hinv = 1./hroms
   hinv = hinv[*]#replicate(1., jpk)
 ; put a z dimension to zeta
-  zeta = transpose(temporary(zeta))
+  zeta = transpose((temporary(zeta))[firstx:lastx, firsty:lasty, *])
   zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite)
   zeta = transpose(temporary(zeta), [2, 1, 3, 0])
 
-  z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))
   z = (z0[*]#replicate(1., nt)) + temporary(zeta) * ((1.+z0*hinv)[*]#replicate(1., nt))
   z = reform(z, nx, ny, jpk, nt, /overwrite)