4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
144
145 implicit none
146
147
148
149 integer, intent(in) :: ng
150
151
152
153 integer :: k
154
155 real(dp) :: Aweight, Bweight, Cweight, Cbot, Csur, Hscale
156 real(dp) :: ds, exp_bot, exp_sur, rk, rN, sc_r, sc_w
157 real(dp) :: cff, cff1, cff2
158 real(dp) :: zhc, z1, z2, z3
159
160
161
162
163
164
165
166
167
168
169
171# if defined WET_DRY
173# else
175# endif
178 END IF
179
180
181
182
183
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202 IF (
theta_s(ng).ne.0.0_dp)
THEN
204 cff2=0.5_dp/tanh(0.5_dp*
theta_s(ng))
205 END IF
208 ds=1.0_dp/real(
n(ng),dp)
210 scalars(ng)%sc_w(k)=ds*real(k-
n(ng),dp)
211 scalars(ng)%sc_r(k)=ds*(real(k-
n(ng),dp)-0.5_dp)
212 IF (
theta_s(ng).ne.0.0_dp)
THEN
219 & 0.5_dp))- &
220 & 0.5_dp)
227 & 0.5_dp))- &
228 & 0.5_dp)
229 ELSE
232 END IF
233 END DO
234
235
236
237
238
239
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267 aweight=1.0_dp
268 bweight=1.0_dp
269 ds=1.0_dp/real(
n(ng),dp)
270
274 sc_w=ds*real(k-
n(ng),dp)
276 IF (
theta_s(ng).gt.0.0_dp)
THEN
277 csur=(1.0_dp-cosh(
theta_s(ng)*sc_w))/ &
279 IF (
theta_b(ng).gt.0.0_dp)
THEN
280 cbot=sinh(
theta_b(ng)*(sc_w+1.0_dp))/ &
282 cweight=(sc_w+1.0_dp)**aweight* &
283 & (1.0_dp+(aweight/bweight)* &
284 & (1.0_dp-(sc_w+1.0_dp)**bweight))
285 scalars(ng)%Cs_w(k)=cweight*csur+(1.0_dp-cweight)*cbot
286 ELSE
288 END IF
289 ELSE
291 END IF
292 END DO
295
297 sc_r=ds*(real(k-
n(ng),dp)-0.5_dp)
299 IF (
theta_s(ng).gt.0.0_dp)
THEN
300 csur=(1.0_dp-cosh(
theta_s(ng)*sc_r))/ &
302 IF (
theta_b(ng).gt.0.0_dp)
THEN
303 cbot=sinh(
theta_b(ng)*(sc_r+1.0_dp))/ &
305 cweight=(sc_r+1.0_dp)**aweight* &
306 & (1.0_dp+(aweight/bweight)* &
307 & (1.0_dp-(sc_r+1.0_dp)**bweight))
308 scalars(ng)%Cs_r(k)=cweight*csur+(1.0_dp-cweight)*cbot
309 ELSE
311 END IF
312 ELSE
314 END IF
315 END DO
316
317
318
319
320
321
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
359 hscale=3.0_dp
360 ds=1.0_dp/real(
n(ng),dp)
361
365 sc_w=ds*real(k-
n(ng),dp)
367 cbot= log(cosh(hscale*(sc_w+1.0_dp)**exp_bot))/ &
368 & log(cosh(hscale))-1.0_dp
369 csur=-log(cosh(hscale*abs(sc_w)**exp_sur))/ &
370 & log(cosh(hscale))
371 cweight=0.5_dp*(1.0_dp-tanh(hscale*(sc_w+0.5_dp)))
372 scalars(ng)%Cs_w(k)=cweight*cbot+(1.0_dp-cweight)*csur
373 END DO
376
378 sc_r=ds*(real(k-
n(ng),dp)-0.5_dp)
380 cbot= log(cosh(hscale*(sc_r+1.0_dp)**exp_bot))/ &
381 & log(cosh(hscale))-1.0_dp
382 csur=-log(cosh(hscale*abs(sc_r)**exp_sur))/ &
383 & log(cosh(hscale))
384 cweight=0.5_dp*(1.0_dp-tanh(hscale*(sc_r+0.5_dp)))
385 scalars(ng)%Cs_r(k)=cweight*cbot+(1.0_dp-cweight)*csur
386 END DO
387
388
389
390
391
392
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433 ds=1.0_dp/real(
n(ng),dp)
434
438 sc_w=ds*real(k-
n(ng),dp)
440 IF (
theta_s(ng).gt.0.0_dp)
THEN
441 csur=(1.0_dp-cosh(
theta_s(ng)*sc_w))/ &
443 ELSE
444 csur=-sc_w**2
445 END IF
446 IF (
theta_b(ng).gt.0.0_dp)
THEN
447 cbot=(exp(
theta_b(ng)*csur)-1.0_dp)/ &
450 ELSE
452 END IF
453 END DO
456
458 sc_r=ds*(real(k-
n(ng),dp)-0.5_dp)
460 IF (
theta_s(ng).gt.0.0_dp)
THEN
461 csur=(1.0_dp-cosh(
theta_s(ng)*sc_r))/ &
463 ELSE
464 csur=-sc_r**2
465 END IF
466 IF (
theta_b(ng).gt.0.0_dp)
THEN
467 cbot=(exp(
theta_b(ng)*csur)-1.0_dp)/ &
470 ELSE
472 END IF
473 END DO
474
475
476
477
478
479
480
481
482
483
484
485
490 rk=real(k,dp)
492 sc_w=-(rk*rk - 2.0_dp*rk*rn + rk + rn*rn - rn)/(rn*rn - rn)- &
493 0.01_dp*(rk*rk - rk*rn)/(1.0_dp - rn)
495 IF (
theta_s(ng).gt.0.0_dp)
THEN
496 csur=(1.0_dp-cosh(
theta_s(ng)*sc_w))/ &
498 ELSE
499 csur=-sc_w**2
500 END IF
501 IF (
theta_b(ng).gt.0.0_dp)
THEN
502 cbot=(exp(
theta_b(ng)*csur)-1.0_dp)/ &
505 ELSE
507 END IF
508 END DO
511
513 rk=real(k,dp)-0.5_dp
515 sc_r=-(rk*rk - 2.0_dp*rk*rn + rk + rn*rn - rn)/(rn*rn - rn)- &
516 0.01_dp*(rk*rk - rk*rn)/(1.0_dp - rn)
518 IF (
theta_s(ng).gt.0.0_dp)
THEN
519 csur=(1.0_dp-cosh(
theta_s(ng)*sc_r))/ &
521 ELSE
522 csur=-sc_r**2
523 END IF
524 IF (
theta_b(ng).gt.0.0_dp)
THEN
525 cbot=(exp(
theta_b(ng)*csur)-1.0_dp)/ &
528 ELSE
530 END IF
531 END DO
532 END IF
533
534
535
536
537
554 IF (
hc(ng).gt.
hmax(ng))
THEN
555 zhc=z3
556 ELSE
557 zhc=0.5_dp*
hc(ng)*(
scalars(ng)%sc_w(k)+ &
559 END IF
560 END IF
562 &
scalars(ng)%Cs_w(k), z1, zhc, z2, z3
563 END DO
564 END IF
565
566 10 FORMAT (/,' Vertical S-coordinate System, Grid ',i2.2,':'/,/, &
567 & ' level S-coord Cs-curve Z',3x, &
568 & 'at hmin at hc half way at hmax',/)
569 20 FORMAT (i6,2f12.7,1x,4f12.3)
570
571 RETURN
integer, dimension(:), allocatable n
real(dp), dimension(:), allocatable hmin
real(dp), dimension(:), allocatable theta_s
real(dp), dimension(:), allocatable tcline
real(dp), dimension(:), allocatable theta_b
logical, dimension(:), allocatable lwrtinfo
real(dp), dimension(:), allocatable hc
type(t_scalars), dimension(:), allocatable scalars
real(dp), dimension(:), allocatable hmax
integer, dimension(:), allocatable vstretching
integer, dimension(:), allocatable vtransform