Ticket #313: set_scoord.F

File set_scoord.F, 20.0 KB (added by rsignell, 16 years ago)

set_scoord.F that works with WET_DRY and Vstretching==2 or 3

Line 
1#include "cppdefs.h"
2#ifdef SOLVE3D
3 SUBROUTINE set_scoord (ng)
4!
5!svn $Id: set_scoord.F 332 2009-03-13 03:38:51Z arango $
6!=======================================================================
7! Copyright (c) 2002-2009 The ROMS/TOMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.txt Hernan G. Arango !
10!========================================== Alexander F. Shchepetkin ===
11! !
12! This routine sets and initializes relevant variables associated !
13! with the vertical terrain-following coordinates transformation. !
14! !
15! Definitions: !
16! !
17! N(ng) : Number of vertical levels for each nested grid. !
18! !
19! zeta : time-varying free-surface, zeta(x,y,t), (m) !
20! !
21! h : bathymetry, h(x,y), (m, positive, maybe time-varying) !
22! !
23! hc : critical (thermocline, pycnocline) depth (m, positive) !
24! !
25! z : vertical depths, z(x,y,s,t), meters, negative !
26! z_w(x,y,0:N(ng)) at W-points (top/bottom cell) !
27! z_r(z,y,1:N(ng)) at RHO-points (cell center) !
28! !
29! z_w(x,y,0 ) = -h(x,y) !
30! z_w(x,y,N(ng)) = zeta(x,y,t) !
31! !
32! s : nondimensional stretched vertical coordinate, !
33! -1 <= s <= 0 !
34! !
35! s = 0 at the free-surface, z(x,y, 0,t) = zeta(x,y,t) !
36! s = -1 at the bottom, z(x,y,-1,t) = - h(x,y,t) !
37! !
38! sc_w(k) = (k-N(ng))/N(ng) k=0:N, W-points !
39! sc_r(k) = (k-N(ng)-0.5)/N(ng) k=1:N, RHO-points !
40! !
41! C : nondimensional vertical stretching function, C(s), !
42! -1 <= C(s) <= 0 !
43! !
44! C(s) = 0 for s = 0, at the free-surface !
45! C(s) = -1 for s = -1, at the bottom !
46! !
47! Cs_w(k) = F(s,theta_s,theta_b) k=0:N, W-points !
48! Cs_r(k) = C(s,theta_s,theta_b) k=1:N, RHO-points !
49! !
50! Zo : vertical transformation functional, Zo(x,y,s): !
51! !
52! Zo(x,y,s) = H(x,y)C(s) separable functions !
53! !
54! !
55! Two vertical transformations are supported, z => z(x,y,s,t): !
56! !
57! (1) Original transformation (Shchepetkin and McWilliams, 2005): In !
58! ROMS since 1999 (version 1.8): !
59! !
60! z(x,y,s,t) = Zo(x,y,s) + zeta(x,y,t) * [1 + Zo(x,y,s)/h(x,y)] !
61! !
62! where !
63! !
64! Zo(x,y,s) = hc * s + [h(x,y) - hc] * C(s) !
65! !
66! Zo(x,y,s) = 0 for s = 0, C(s) = 0, at the surface !
67! Zo(x,y,s) = -h(x,y) for s = -1, C(s) = -1, at the bottom !
68! !
69! (2) New transformation: In UCLA-ROMS since 2005: !
70! !
71! z(x,y,s,t) = zeta(x,y,t) + [zeta(x,y,t) + h(x,y)] * Zo(x,y,s) !
72! !
73! where !
74! !
75! Zo(x,y,s) = [hc * s(k) + h(x,y) * C(k)] / [hc + h(x,y)] !
76! !
77! Zo(x,y,s) = 0 for s = 0, C(s) = 0, at the surface !
78! Zo(x,y,s) = -1 for s = -1, C(s) = -1, at the bottom !
79! !
80! At the rest state, corresponding to zero free-surface, this !
81! transformation yields the following unperturbed depths, zhat: !
82! !
83! zhat = z(x,y,s,0) = h(x,y) * Zo(x,y,s) !
84! !
85! = h(x,y) * [hc * s(k) + h(x,y) * C(k)] / [hc + h(x,y)] !
86! !
87! and !
88! !
89! d(zhat) = ds * h(x,y) * hc / [hc + h(x,y)] !
90! !
91! As a consequence, the uppermost grid box retains very little !
92! dependency from bathymetry in the areas where hc << h(x,y), !
93! that is deep areas. For example, if hc=250 m, and h(x,y) !
94! changes from 2000 to 6000 meters, the uppermost grid box !
95! changes only by a factor of 1.08 (less than 10%). !
96! !
97! Notice that: !
98! !
99! * Regardless of the design of C(s), transformation (2) behaves !
100! like equally-spaced sigma-coordinates in shallow areas, where !
101! h(x,y) << hc. This is advantageous because high vertical !
102! resolution and associated CFL limitation is avoided in these !
103! areas. !
104! !
105! * Near-surface refinement is close to geopotential coordinates !
106! in deep areas (level thickness do not depend or weakly-depend !
107! on the bathymetry). Contrarily, near-bottom refinement is !
108! like sigma-coordinates with thicknesses roughly proportional !
109! to depth reducing high r-factors in these areas. !
110! !
111! !
112! This generic transformation design facilitates numerous vertical !
113! stretching functions, C(s). These functions are set-up in this !
114! routine in terms of several stretching parameters specified in !
115! the standard input file. !
116! !
117! C(s) vertical stretching function properties: !
118! !
119! * a nonlinear, monotonic function !
120! * a continuous differentiable function, or !
121! * a piecewise function with smooth transition and differentiable !
122! * must be constrained by -1 <= C(s) <= 0, with C(0)=0 at the !
123! free-surface and C(-1)=-1 at the bottom (bathymetry). !
124! !
125! References: !
126! !
127! Shchepetkin, A.F. and J.C. McWilliams, 2005: The regional oceanic !
128! modeling system (ROMS): a split-explicit, free-surface, !
129! topography-following-coordinate oceanic model, Ocean !
130! Modelling, 9, 347-404. !
131! !
132! Song, Y. and D. Haidvogel, 1994: A semi-implicit ocean !
133! circulation model using a generalized topography- !
134! following coordinate system, J. Comp. Physics, !
135! 115, 228-244. !
136! !
137!=======================================================================
138!
139 USE mod_param
140 USE mod_parallel
141 USE mod_grid
142 USE mod_iounits
143 USE mod_scalars
144!
145 implicit none
146!
147! Imported variable declarations.
148!
149 integer, intent(in) :: ng
150!
151! Local variable declarations.
152!
153 integer :: k
154
155 real(r8) :: Aweight, Bweight, Cweight, Cbot, Csur, Hscale
156 real(r8) :: ds, exp_bot, exp_sur, sc_r, sc_w
157 real(r8) :: cff, cff1, cff2, cff3
158!
159!-----------------------------------------------------------------------
160! Original vertical strectching function, Song and Haidvogel (1994).
161!-----------------------------------------------------------------------
162!
163# if defined WET_DRY
164 hc(ng)=MIN(MAX(hmin(ng),0.0_r8),Tcline(ng))
165# else
166 hc(ng)=MIN(hmin(ng),Tcline(ng))
167# endif
168 IF (Vstretching(ng).eq.1) THEN
169!
170! This vertical stretching function is defined as:
171!
172! C(s) = (1 - b) * [SINH(s * a) / SINH(a)] +
173!
174! b * [-0.5 + 0.5 * TANH(a * (s + 0.5)) / TANH(0.5 * a)]
175!
176! where the stretching parameters (a, b) are specify at input:
177!
178! a = theta_s 0 < theta_s <= 8
179! b = theta_b 0 <= theta_b <= 1
180!
181! If theta_b=0, the refinement is surface intensified as theta_s is
182! increased.
183! If theta_b=1, the refinement is both bottom ans surface intensified
184! as theta_s is increased.
185!
186
187 IF (theta_s(ng).ne.0.0_r8) THEN
188 cff1=1.0_r8/SINH(theta_s(ng))
189 cff2=0.5_r8/TANH(0.5_r8*theta_s(ng))
190 END IF
191 SCALARS(ng)%sc_w(0)=-1.0_r8
192 SCALARS(ng)%Cs_w(0)=-1.0_r8
193 ds=1.0_r8/REAL(N(ng),r8)
194 DO k=1,N(ng)
195 SCALARS(ng)%sc_w(k)=ds*REAL(k-N(ng),r8)
196 SCALARS(ng)%sc_r(k)=ds*(REAL(k-N(ng),r8)-0.5_r8)
197 IF (theta_s(ng).ne.0.0_r8) THEN
198 SCALARS(ng)%Cs_w(k)=(1.0_r8-theta_b(ng))* &
199 & cff1*SINH(theta_s(ng)* &
200 & SCALARS(ng)%sc_w(k))+ &
201 & theta_b(ng)* &
202 & (cff2*TANH(theta_s(ng)* &
203 & (SCALARS(ng)%sc_w(k)+ &
204 & 0.5_r8))- &
205 & 0.5_r8)
206 SCALARS(ng)%Cs_r(k)=(1.0_r8-theta_b(ng))* &
207 & cff1*SINH(theta_s(ng)* &
208 & SCALARS(ng)%sc_r(k))+ &
209 & theta_b(ng)* &
210 & (cff2*TANH(theta_s(ng)* &
211 & (SCALARS(ng)%sc_r(k)+ &
212 & 0.5_r8))- &
213 & 0.5_r8)
214 ELSE
215 SCALARS(ng)%Cs_w(k)=SCALARS(ng)%sc_w(k)
216 SCALARS(ng)%Cs_r(k)=SCALARS(ng)%sc_r(k)
217 END IF
218 END DO
219!
220!-----------------------------------------------------------------------
221! A. Shchepetkin new vertical stretching function.
222!-----------------------------------------------------------------------
223!
224 ELSE IF (Vstretching(ng).eq.2) THEN
225!
226! This vertical stretching function is defined, in the simplest form,
227! as:
228!
229! C(s) = [1.0 - COSH(theta_s * s)] / [COSH(theta_s) - 1.0]
230!
231! it is similar in meaning to the original vertical stretcing function
232! (Song and Haidvogel, 1994), but note that hyperbolic functions are
233! COSH, and not SINH.
234!
235! Note that the above definition results in
236!
237! -1 <= C(s) <= 0
238!
239! as long as
240!
241! -1 <= s <= 0
242!
243! and, unlike in any previous definition
244!
245! d[C(s)]/ds --> 0 if s --> 0
246!
247! For the purpose of bottom boundary layer C(s) is further modified
248! to allow near-bottom refinement. This is done by blending it with
249! another function.
250!
251
252 Aweight=1.0_r8
253 Bweight=1.0_r8
254 ds=1.0_r8/REAL(N(ng),r8)
255!
256 SCALARS(ng)%sc_w(N(ng))=0.0_r8
257 SCALARS(ng)%Cs_w(N(ng))=0.0_r8
258 DO k=N(ng)-1,1,-1
259 sc_w=ds*REAL(k-N(ng),r8)
260 SCALARS(ng)%sc_w(k)=sc_w
261 IF (theta_s(ng).gt.0.0_r8) THEN
262 Csur=(1.0_r8-COSH(theta_s(ng)*sc_w))/ &
263 & (COSH(theta_s(ng))-1.0_r8)
264 IF (theta_b(ng).gt.0.0_r8) THEN
265 Cbot=SINH(theta_b(ng)*(sc_w+1.0_r8))/ &
266 & SINH(theta_b(ng))-1.0_r8
267 Cweight=(sc_w+1.0_r8)**Aweight* &
268 & (1.0_r8+(Aweight/Bweight)* &
269 & (1.0_r8-(sc_w+1.0_r8)**Bweight))
270 SCALARS(ng)%Cs_w(k)=Cweight*Csur+(1.0_r8-Cweight)*Cbot
271 ELSE
272 SCALARS(ng)%Cs_w(k)=Csur
273 END IF
274 ELSE
275 SCALARS(ng)%Cs_w(k)=sc_w
276 END IF
277 END DO
278 SCALARS(ng)%sc_w(0)=-1.0_r8
279 SCALARS(ng)%Cs_w(0)=-1.0_r8
280!
281 DO k=1,N(ng)
282 sc_r=ds*(REAL(k-N(ng),r8)-0.5_r8)
283 SCALARS(ng)%sc_r(k)=sc_r
284 IF (theta_s(ng).gt.0.0_r8) THEN
285 Csur=(1.0_r8-COSH(theta_s(ng)*sc_r))/ &
286 & (COSH(theta_s(ng))-1.0_r8)
287 IF (theta_b(ng).gt.0.0_r8) THEN
288 Cbot=SINH(theta_b(ng)*(sc_r+1.0_r8))/ &
289 & SINH(theta_b(ng))-1.0_r8
290 Cweight=(sc_r+1.0_r8)**Aweight* &
291 & (1.0_r8+(Aweight/Bweight)* &
292 & (1.0_r8-(sc_r+1.0_r8)**Bweight))
293 SCALARS(ng)%Cs_r(k)=Cweight*Csur+(1.0_r8-Cweight)*Cbot
294 ELSE
295 SCALARS(ng)%Cs_r(k)=Csur
296 END IF
297 ELSE
298 SCALARS(ng)%Cs_r(k)=sc_r
299 END IF
300 END DO
301
302!
303!-----------------------------------------------------------------------
304! R. Geyer stretching function for high bottom boundary layer
305! resolution.
306!-----------------------------------------------------------------------
307!
308 ELSE IF (Vstretching(ng).eq.3) THEN
309!
310! This stretching function is intended for very shallow coastal
311! applications, like gravity sediment flows.
312!
313! At the surface, C(s=0)=0
314!
315! C(s) = - LOG(COSH(Hscale * ABS(s) ** alpha)) /
316! LOG(COSH(Hscale))
317!
318! At the bottom, C(s=-1)=-1
319!
320! C(s) = LOG(COSH(Hscale * (s + 1) ** beta)) /
321! LOG(COSH(Hscale)) - 1
322!
323! where
324!
325! Hscale : scale value for all hyperbolic functions
326! Hscale = 3.0 set internally here
327! alpha : surface stretching exponent
328! alpha = 0.65 minimal increase of surface resolution
329! 1.0 significant amplification
330! beta : bottoom stretching exponent
331! beta = 0.58 no amplification
332! 1.0 significant amplification
333! 3.0 super-high bottom resolution
334! s : stretched vertical coordinate, -1 <= s <= 0
335! s(k) = (k-N)/N k=0:N, W-points (s_w)
336! s(k) = (k-N-0.5)/N k=1:N, RHO-points (s_rho)
337!
338! The stretching exponents (alpha, beta) are specify at input:
339!
340! alpha = theta_s
341! beta = theta_b
342!
343
344 exp_sur=theta_s(ng)
345 exp_bot=theta_b(ng)
346 Hscale=3.0_r8
347 ds=1.0_r8/REAL(N(ng),r8)
348!
349 SCALARS(ng)%sc_w(N(ng))=0.0_r8
350 SCALARS(ng)%Cs_w(N(ng))=0.0_r8
351 DO k=N(ng)-1,1,-1
352 sc_w=ds*REAL(k-N(ng),r8)
353 SCALARS(ng)%sc_w(k)=sc_w
354 Cbot= LOG(COSH(Hscale*(sc_w+1.0_r8)**exp_bot))/ &
355 & LOG(COSH(Hscale))-1.0_r8
356 Csur=-LOG(COSH(Hscale*ABS(sc_w)**exp_sur))/ &
357 & LOG(COSH(Hscale))
358 Cweight=0.5_r8*(1.0_r8-TANH(Hscale*(sc_w+0.5_r8)))
359 SCALARS(ng)%Cs_w(k)=Cweight*Cbot+(1.0_r8-Cweight)*Csur
360 END DO
361 SCALARS(ng)%sc_w(0)=-1.0_r8
362 SCALARS(ng)%Cs_w(0)=-1.0_r8
363!
364 DO k=1,N(ng)
365 sc_r=ds*(REAL(k-N(ng),r8)-0.5_r8)
366 SCALARS(ng)%sc_r(k)=sc_r
367 Cbot= LOG(COSH(Hscale*(sc_r+1.0_r8)**exp_bot))/ &
368 & LOG(COSH(Hscale))-1.0_r8
369 Csur=-LOG(COSH(Hscale*ABS(sc_r)**exp_sur))/ &
370 & LOG(COSH(Hscale))
371 Cweight=0.5_r8*(1.0_r8-TANH(Hscale*(sc_r+0.5_r8)))
372 SCALARS(ng)%Cs_r(k)=Cweight*Cbot+(1.0_r8-Cweight)*Csur
373 END DO
374
375 END IF
376!
377!-----------------------------------------------------------------------
378! Report information about vertical transformation.
379!-----------------------------------------------------------------------
380!
381 IF (Master.and.LwrtInfo(ng)) THEN
382 WRITE (stdout,20)
383 DO k=N(ng),0,-1
384 IF (Vstretching(ng).eq.2) THEN
385 cff=0.5_r8*(hmax(ng)+hmin(ng))
386 cff1=hc(ng) *(SCALARS(ng)%sc_w(k)*hc(ng)+ &
387 & SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+hc(ng))
388 cff2=cff *(SCALARS(ng)%sc_w(k)*cff+ &
389 & SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+cff)
390 cff3=hmax(ng)*(SCALARS(ng)%sc_w(k)*hmax(ng)+ &
391 & SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+hmax(ng))
392 ELSE
393 cff1=SCALARS(ng)%sc_w(k)*hc(ng)+ &
394 & (hmin(ng)-hc(ng))*SCALARS(ng)%Cs_w(k)
395 cff2=SCALARS(ng)%sc_w(k)*hc(ng)+ &
396 & (0.5*(hmin(ng)+hmax(ng))-hc(ng))*SCALARS(ng)%Cs_w(k)
397 cff3=SCALARS(ng)%sc_w(k)*hc(ng)+ &
398 & (hmax(ng)-hc(ng))*SCALARS(ng)%Cs_w(k)
399 END IF
400 WRITE (stdout,30) k, SCALARS(ng)%sc_w(k), &
401 & SCALARS(ng)%Cs_w(k), cff1, cff2, cff3
402 END DO
403 END IF
404
405 10 FORMAT (/,' SET_SCOORD - Specified critical depth parameter, ', &
406 & ' hc =',f7.3,/,14x, &
407 & 'exceeds minimum unmasked bathymetry, hmin = ',f7.3)
408 20 FORMAT (/,' Vertical S-coordinate System: ',/,/, &
409 & ' level S-coord Cs-curve',10x, &
410 & 'at_hmin over_slope at_hmax',/)
411 30 FORMAT (i6,2f12.7,4x,3f12.3)
412
413 RETURN
414 END SUBROUTINE set_scoord
415#else
416 SUBROUTINE set_scoord
417 RETURN
418 END SUBROUTINE set_scoord
419#endif
420