ROMS
Loading...
Searching...
No Matches
rp_set_depth.F
Go to the documentation of this file.
1#include "cppdefs.h"
3#ifdef TL_IOMS
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes the time evolving depths of the model grid !
13! and its associated vertical transformation metric (thickness). !
14! !
15! Currently, two vertical coordinate transformations are available !
16! with various possible vertical stretching, C(s), functions, (see !
17! routine "set_scoord.F" for details). !
18! !
19! BASIC STATE variables needed: NONE !
20! Independent Variables: tl_Hz, tl_z_r, tl_z_w !
21! !
22!=======================================================================
23!
24 implicit none
25!
26 PRIVATE
27# ifdef SOLVE3D
29# ifdef ADJUST_BOUNDARY
30 PUBLIC :: rp_set_depth_bry
31# endif
32# endif
33 PUBLIC :: rp_bath, rp_bath_tile
34!
35 CONTAINS
36
37# ifdef SOLVE3D
38!
39!***********************************************************************
40 SUBROUTINE rp_set_depth (ng, tile, model)
41!***********************************************************************
42!
43 USE mod_param
44 USE mod_coupling
45 USE mod_grid
46 USE mod_ocean
47 USE mod_stepping
48!
49! Imported variable declarations.
50!
51 integer, intent(in) :: ng, tile, model
52!
53! Local variable declarations.
54!
55 character (len=*), parameter :: myfile = &
56 & __FILE__
57!
58# include "tile.h"
59!
60# ifdef PROFILE
61 CALL wclock_on (ng, model, 12, __line__, myfile)
62# endif
63 CALL rp_set_depth_tile (ng, tile, model, &
64 & lbi, ubi, lbj, ubj, &
65 & imins, imaxs, jmins, jmaxs, &
66 & nstp(ng), nnew(ng), &
67 & grid(ng) % h, &
68 & grid(ng) % tl_h, &
69# ifdef ICESHELF
70 & grid(ng) % zice, &
71# endif
72 & coupling(ng) % Zt_avg1, &
73 & coupling(ng) % tl_Zt_avg1, &
74 & grid(ng) % tl_Hz, &
75 & grid(ng) % tl_z_r, &
76 & grid(ng) % tl_z_w)
77# ifdef PROFILE
78 CALL wclock_off (ng, model, 12, __line__, myfile)
79# endif
80!
81 RETURN
82 END SUBROUTINE rp_set_depth
83!
84!***********************************************************************
85 SUBROUTINE rp_set_depth_tile (ng, tile, model, &
86 & LBi, UBi, LBj, UBj, &
87 & IminS, ImaxS, JminS, JmaxS, &
88 & nstp, nnew, &
89 & h, tl_h, &
90# ifdef ICESHELF
91 & zice, &
92# endif
93 & Zt_avg1, tl_Zt_avg1, &
94 & tl_Hz, tl_z_r, tl_z_w)
95!***********************************************************************
96!
97 USE mod_param
98 USE mod_scalars
99!
102# ifdef DISTRIBUTE
104# endif
105!
106! Imported variable declarations.
107!
108 integer, intent(in) :: ng, tile, model
109 integer, intent(in) :: lbi, ubi, lbj, ubj
110 integer, intent(in) :: imins, imaxs, jmins, jmaxs
111 integer, intent(in) :: nstp, nnew
112!
113# ifdef ASSUMED_SHAPE
114 real(r8), intent(in) :: h(lbi:,lbj:)
115# ifdef ICESHELF
116 real(r8), intent(in) :: zice(lbi:,lbj:)
117# endif
118 real(r8), intent(in) :: zt_avg1(lbi:,lbj:)
119 real(r8), intent(in) :: tl_zt_avg1(lbi:,lbj:)
120 real(r8), intent(inout) :: tl_h(lbi:,lbj:)
121 real(r8), intent(out) :: tl_hz(lbi:,lbj:,:)
122 real(r8), intent(out) :: tl_z_r(lbi:,lbj:,:)
123 real(r8), intent(out) :: tl_z_w(lbi:,lbj:,0:)
124# else
125 real(r8), intent(in) :: h(lbi:ubi,lbj:ubj)
126# ifdef ICESHELF
127 real(r8), intent(in) :: zice(lbi:ubi,lbj:ubj)
128# endif
129 real(r8), intent(in) :: zt_avg1(lbi:ubi,lbj:ubj)
130 real(r8), intent(in) :: tl_zt_avg1(lbi:ubi,lbj:ubj)
131 real(r8), intent(inout) :: tl_h(lbi:ubi,lbj:ubj)
132 real(r8), intent(out) :: tl_hz(lbi:ubi,lbj:ubj,n(ng))
133 real(r8), intent(out) :: tl_z_r(lbi:ubi,lbj:ubj,n(ng))
134 real(r8), intent(out) :: tl_z_w(lbi:ubi,lbj:ubj,0:n(ng))
135# endif
136!
137! Local variable declarations.
138!
139 integer :: i, j, k
140
141 real(r8) :: cff, cff_r, cff1_r, cff2_r, cff_w, cff1_w, cff2_w
142 real(r8) :: hinv, hwater, z_r0, z_w0
143 real(r8) :: tl_cff2_r, tl_cff2_w
144 real(r8) :: tl_hinv, tl_hwater, tl_z_r0, tl_z_w0
145
146# ifdef WET_DRY
147 real(r8), parameter :: eps = 1.0e-14_r8
148# endif
149
150# include "set_bounds.h"
151!
152!-----------------------------------------------------------------------
153! Original formulation: Compute vertical depths (meters, negative) at
154! RHO- and W-points, and vertical grid
155! thicknesses. Various stretching functions are possible.
156!
157! z_w(x,y,s,t) = Zo_w + zeta(x,y,t) * [1.0 + Zo_w / h(x,y)]
158!
159! Zo_w = hc * [s(k) - C(k)] + C(k) * h(x,y)
160!
161!-----------------------------------------------------------------------
162!
163 IF (vtransform(ng).eq.1) THEN
164 DO j=jstrt,jendt
165 DO i=istrt,iendt
166!^ z_w(i,j,0)=-h(i,j)
167!^
168 tl_z_w(i,j,0)=-tl_h(i,j)
169 END DO
170 DO k=1,n(ng)
171 cff_r=hc(ng)*(scalars(ng)%sc_r(k)-scalars(ng)%Cs_r(k))
172 cff_w=hc(ng)*(scalars(ng)%sc_w(k)-scalars(ng)%Cs_w(k))
173 cff1_r=scalars(ng)%Cs_r(k)
174 cff1_w=scalars(ng)%Cs_w(k)
175 DO i=istrt,iendt
176 hwater=h(i,j)
177# ifdef ICESHELF
178 hwater=hwater-abs(zice(i,j))
179# endif
180 tl_hwater=tl_h(i,j)
181 hinv=1.0_r8/hwater
182# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
183 tl_hinv=-hinv*hinv*tl_hwater+ &
184# ifdef TL_IOMS
185 & 2.0_r8*hinv
186# endif
187# endif
188 z_w0=cff_w+cff1_w*hwater
189# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
190 tl_z_w0=cff1_w*tl_hwater+ &
191# ifdef TL_IOMS
192 & cff_w
193# endif
194# endif
195!^ z_w(i,j,k)=z_w0+Zt_avg1(i,j)*(1.0_r8+z_w0*hinv)
196!^
197# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
198 tl_z_w(i,j,k)=tl_z_w0+ &
199 & tl_zt_avg1(i,j)*(1.0_r8+z_w0*hinv)+ &
200 & zt_avg1(i,j)*z_w0*tl_hinv+ &
201 & zt_avg1(i,j)*tl_z_w0*hinv- &
202# ifdef TL_IOMS
203 & 2.0_r8*zt_avg1(i,j)*z_w0*hinv
204# endif
205# else
206 tl_z_w(i,j,k)=tl_zt_avg1(i,j)*(1.0_r8+z_w0*hinv)+ &
207# ifdef TL_IOMS
208 & z_w0
209# endif
210# endif
211 z_r0=cff_r+cff1_r*hwater
212# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
213 tl_z_r0=cff1_r*tl_hwater+ &
214# ifdef TL_IOMS
215 & cff_r
216# endif
217# endif
218!^ z_r(i,j,k)=z_r0+Zt_avg1(i,j)*(1.0_r8+z_r0*hinv)
219!^
220# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
221 tl_z_r(i,j,k)=tl_z_r0+ &
222 & tl_zt_avg1(i,j)*(1.0_r8+z_r0*hinv)+ &
223 & zt_avg1(i,j)*z_r0*tl_hinv+ &
224 & zt_avg1(i,j)*tl_z_r0*hinv- &
225# ifdef TL_IOMS
226 & 2.0_r8*zt_avg1(i,j)*z_r0*hinv
227# endif
228# else
229 tl_z_r(i,j,k)=tl_zt_avg1(i,j)*(1.0_r8+z_r0*hinv)+ &
230# ifdef TL_IOMS
231 & z_r0
232# endif
233# endif
234# ifdef ICESHELF
235!^ z_w(i,j,k)=z_w(i,j,k)-ABS(zice(i,j))
236!^ z_r(i,j,k)=z_r(i,j,k)-ABS(zice(i,j))
237# endif
238!^ Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1)
239!^
240 tl_hz(i,j,k)=tl_z_w(i,j,k)-tl_z_w(i,j,k-1)
241 END DO
242 END DO
243 END DO
244!
245!-----------------------------------------------------------------------
246! New formulation: Compute vertical depths (meters, negative) at
247! RHO- and W-points, and vertical grid thicknesses.
248! Various stretching functions are possible.
249!
250! z_w(x,y,s,t) = zeta(x,y,t) + [zeta(x,y,t)+ h(x,y)] * Zo_w
251!
252! Zo_w = [hc * s(k) + C(k) * h(x,y)] / [hc + h(x,y)]
253!
254!-----------------------------------------------------------------------
255!
256 ELSE IF (vtransform(ng).eq.2) THEN
257 DO j=jstrt,jendt
258 DO i=istrt,iendt
259!^ z_w(i,j,0)=-h(i,j)
260!^
261 tl_z_w(i,j,0)=-tl_h(i,j)
262 END DO
263 DO k=1,n(ng)
264 cff_r=hc(ng)*scalars(ng)%sc_r(k)
265 cff_w=hc(ng)*scalars(ng)%sc_w(k)
266 cff1_r=scalars(ng)%Cs_r(k)
267 cff1_w=scalars(ng)%Cs_w(k)
268 DO i=istrt,iendt
269 hwater=h(i,j)
270# ifdef ICESHELF
271 hwater=hwater-abs(zice(i,j))
272# endif
273 tl_hwater=tl_h(i,j)
274 hinv=1.0_r8/(hc(ng)+hwater)
275# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
276 tl_hinv=-hinv*hinv*tl_hwater+ &
277# ifdef TL_IOMS
278 & 2.0_r8*hinv
279# endif
280# endif
281
282 cff2_w=(cff_w+cff1_w*hwater)*hinv
283# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
284 tl_cff2_w=cff1_w*tl_hwater*hinv+ &
285 & (cff_w+cff1_w*hwater)*tl_hinv- &
286# ifdef TL_IOMS
287 & cff1_w*hwater*hinv
288# endif
289# endif
290
291!^ z_w(i,j,k)=Zt_avg1(i,j)+ &
292!^ & (Zt_avg1(i,j)+hwater)*cff2_w
293!^
294# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
295 tl_z_w(i,j,k)=tl_zt_avg1(i,j)+ &
296 & (tl_zt_avg1(i,j)+tl_hwater)*cff2_w+ &
297 & (zt_avg1(i,j)+hwater)*tl_cff2_w- &
298# ifdef TL_IOMS
299 & (zt_avg1(i,j)+hwater)*cff2_w
300# endif
301# else
302 tl_z_w(i,j,k)=tl_zt_avg1(i,j)+ &
303 & tl_zt_avg1(i,j)*cff2_w+ &
304# ifdef TL_IOMS
305 & hwater*cff2_w
306# endif
307# endif
308
309 cff2_r=(cff_r+cff1_r*hwater)*hinv
310# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
311 tl_cff2_r=cff1_r*tl_hwater*hinv+ &
312 & (cff_r+cff1_r*hwater)*tl_hinv- &
313# ifdef TL_IOMS
314 & hwater*hinv
315# endif
316# endif
317
318!^ z_r(i,j,k)=Zt_avg1(i,j)+ &
319!^ & (Zt_avg1(i,j)+hwater)*cff2_r
320!^
321# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
322 tl_z_r(i,j,k)=tl_zt_avg1(i,j)+ &
323 & (tl_zt_avg1(i,j)+tl_hwater)*cff2_r+ &
324 & (zt_avg1(i,j)+hwater)*tl_cff2_r- &
325# ifdef TL_IOMS
326 & (zt_avg1(i,j)+hwater)*cff2_r
327# endif
328# else
329 tl_z_r(i,j,k)=tl_zt_avg1(i,j)+ &
330 & tl_zt_avg1(i,j)*cff2_r+ &
331# ifdef TL_IOMS
332 & hwater*cff2_r
333# endif
334# endif
335# ifdef ICESHELF
336!^ z_w(i,j,k)=z_w(i,j,k)-ABS(zice(i,j))
337!^ z_r(i,j,k)=z_r(i,j,k)-ABS(zice(i,j))
338# endif
339!^ Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1)
340!^
341 tl_hz(i,j,k)=tl_z_w(i,j,k)-tl_z_w(i,j,k-1)
342 END DO
343 END DO
344 END DO
345 END IF
346!
347!-----------------------------------------------------------------------
348! Exchange boundary information.
349!-----------------------------------------------------------------------
350!
351 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
352!^ CALL exchange_r2d_tile (ng, tile, &
353!^ & LBi, UBi, LBj, UBj, &
354!^ & h)
355!^
356 CALL exchange_r2d_tile (ng, tile, &
357 & lbi, ubi, lbj, ubj, &
358 & tl_h)
359!^ CALL exchange_w3d_tile (ng, tile, &
360!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
361!^ & z_w)
362!^
363 CALL exchange_w3d_tile (ng, tile, &
364 & lbi, ubi, lbj, ubj, 0, n(ng), &
365 & tl_z_w)
366!^ CALL exchange_r3d_tile (ng, tile, &
367!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
368!^ & z_r)
369!^
370 CALL exchange_r3d_tile (ng, tile, &
371 & lbi, ubi, lbj, ubj, 1, n(ng), &
372 & tl_z_r)
373!^ CALL exchange_r3d_tile (ng, tile, &
374!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
375!^ & Hz)
376!^
377 CALL exchange_r3d_tile (ng, tile, &
378 & lbi, ubi, lbj, ubj, 1, n(ng), &
379 & tl_hz)
380 END IF
381
382# ifdef DISTRIBUTE
383!^ CALL mp_exchange2d (ng, tile, model, 1, &
384!^ & LBi, UBi, LBj, UBj, &
385!^ & NghostPoints, &
386!^ & EWperiodic(ng), NSperiodic(ng), &
387!^ & h)
388!^
389 CALL mp_exchange2d (ng, tile, model, 1, &
390 & lbi, ubi, lbj, ubj, &
391 & nghostpoints, &
392 & ewperiodic(ng), nsperiodic(ng), &
393 & tl_h)
394!^ CALL mp_exchange3d (ng, tile, model, 1, &
395!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
396!^ & NghostPoints, &
397!^ & EWperiodic(ng), NSperiodic(ng), &
398!^ & z_w)
399!^
400 CALL mp_exchange3d (ng, tile, model, 1, &
401 & lbi, ubi, lbj, ubj, 0, n(ng), &
402 & nghostpoints, &
403 & ewperiodic(ng), nsperiodic(ng), &
404 & tl_z_w)
405!^ CALL mp_exchange3d (ng, tile, model, 2, &
406!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
407!^ & NghostPoints, &
408!^ & EWperiodic(ng), NSperiodic(ng), &
409!^ & z_r, Hz)
410!^
411 CALL mp_exchange3d (ng, tile, model, 2, &
412 & lbi, ubi, lbj, ubj, 1, n(ng), &
413 & nghostpoints, &
414 & ewperiodic(ng), nsperiodic(ng), &
415 & tl_z_r, tl_hz)
416# endif
417!
418 RETURN
419 END SUBROUTINE rp_set_depth_tile
420
421# ifdef ADJUST_BOUNDARY
422!
423!***********************************************************************
424 SUBROUTINE rp_set_depth_bry (ng, tile, model)
425!***********************************************************************
426!
427 USE mod_param
428 USE mod_grid
429!
430! Imported variable declarations.
431!
432 integer, intent(in) :: ng, tile, model
433!
434! Local variable declarations.
435!
436 character (len=*), parameter :: myfile = &
437 & __FILE__//", rp_set_depth_bry"
438!
439# include "tile.h"
440!
441# ifdef PROFILE
442 CALL wclock_on (ng, model, 12, __line__, myfile)
443# endif
444 CALL rp_set_depth_bry_tile (ng, tile, model, &
445 & lbi, ubi, lbj, ubj, lbij, ubij, &
446 & imins, imaxs, jmins, jmaxs, &
447 & grid(ng) % h, &
448 & grid(ng) % tl_h, &
449# ifdef ICESHELF
450 & grid(ng) % zice, &
451# endif
452 & grid(ng) % tl_Hz_bry)
453# ifdef PROFILE
454 CALL wclock_off (ng, model, 12, __line__, myfile)
455# endif
456!
457 RETURN
458 END SUBROUTINE rp_set_depth_bry
459!
460!***********************************************************************
461 SUBROUTINE rp_set_depth_bry_tile (ng, tile, model, &
462 & LBi, UBi, LBj, UBj, LBij, UBij, &
463 & IminS, ImaxS, JminS, JmaxS, &
464 & h, tl_h, &
465# ifdef ICESHELF
466 & zice, &
467# endif
468 & tl_Hz_bry)
469!***********************************************************************
470!
471 USE mod_param
472 USE mod_boundary
473 USE mod_ncparam
474 USE mod_scalars
475!
476# ifdef DISTRIBUTE
478# endif
479!
480! Imported variable declarations.
481!
482 integer, intent(in) :: ng, tile, model
483 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
484 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
485!
486# ifdef ASSUMED_SHAPE
487 real(r8), intent(in) :: h(LBi:,LBj:)
488 real(r8), intent(in) :: tl_h(LBi:,LBj:)
489# ifdef ICESHELF
490 real(r8), intent(in) :: zice(LBi:,LBj:)
491# endif
492 real(r8), intent(out) :: tl_Hz_bry(LBij:,:,:)
493# else
494 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
495 real(r8), intent(in) :: tl_h(LBi:UBi,LBj:UBj)
496# ifdef ICESHELF
497 real(r8), intent(in) :: zice(LBi:UBi,LBj:UBj)
498# endif
499 real(r8), intent(out) :: tl_Hz_bry(LBij:UBij,N(ng),4)
500# endif
501!
502! Local variable declarations.
503!
504 integer :: i, ibry, j, k
505
506 real(r8) :: cff_w, cff1_w, cff2_w
507 real(r8) :: hinv, hwater, z_w0
508 real(r8) :: tl_cff2_w, tl_hinv, tl_hwater, tl_z_w0
509
510 real(r8), dimension(0:N(ng)) :: tl_Zw
511
512# include "set_bounds.h"
513!
514!-----------------------------------------------------------------------
515! This routine is special and equivalent to "tl_set_bry_tile". The
516! additional terms that are added to RPM transformation are not needed
517! here because of the way that the adjustments to the boundary are
518! computed.
519!-----------------------------------------------------------------------
520!
521!-----------------------------------------------------------------------
522! Original formulation: Compute vertical depths (meters, negative) at
523! RHO- and W-points, and vertical grid
524! thicknesses. Various stretching functions are possible.
525!
526! z_w(x,y,s,t) = Zo_w + zeta(x,y,t) * [1.0 + Zo_w / h(x,y)]
527!
528! Zo_w = hc * [s(k) - C(k)] + C(k) * h(x,y)
529!
530!-----------------------------------------------------------------------
531!
532 IF (vtransform(ng).eq.1) THEN
533
534 IF (tl_lbc(iwest,isfsur,ng)%acquire.and. &
535 & domain(ng)%Western_Edge(tile)) THEN
536 i=bounds(ng)%edge(iwest,r2dvar)
537 DO j=jstrt,jendt
538 hwater=h(i,j)
539# ifdef ICESHELF
540 hwater=hwater-abs(zice(i,j))
541# endif
542 tl_hwater=tl_h(i,j)
543 hinv=1.0_r8/hwater
544# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
545 tl_hinv=-hinv*hinv*tl_hwater
546# endif
547!^ Zw(0)=-h(i,j)
548!^
549 tl_zw(0)=-tl_h(i,j)
550 DO k=1,n(ng)
551 cff_w=hc(ng)*(scalars(ng)%sc_w(k)-scalars(ng)%Cs_w(k))
552 cff1_w=scalars(ng)%Cs_w(k)
553 z_w0=cff_w+cff1_w*hwater
554# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
555 tl_z_w0=cff1_w*tl_hwate
556# endif
557!^ Zw(k)=z_w0+BOUNDARY(ng)%zeta_west(j)*(1.0_r8+z_w0*hinv)
558!^
559# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
560 tl_zw(k)=tl_z_w0+ &
561 & boundary(ng)%tl_zeta_west(j)* &
562 & (1.0_r8+z_w0*hinv)+ &
563 & boundary(ng)%zeta_west(j)*z_w0*tl_hinv+ &
564 & boundary(ng)%zeta_west(j)*tl_z_w0*hinv
565# else
566 tl_zw(k)=boundary(ng)%tl_zeta_west(j)*(1.0_r8+z_w0*hinv)
567# endif
568# ifdef ICESHELF
569!^ Zw(k)=Zw(k)-ABS(zice(i,j))
570# endif
571!^ Hz_bry(j,k,iwest)=Zw(k)-Zw(k-1)
572!^
573 tl_hz_bry(j,k,iwest)=tl_zw(k)-tl_zw(k-1)
574 END DO
575 END DO
576 END IF
577
578 IF (tl_lbc(ieast,isfsur,ng)%acquire.and. &
579 & domain(ng)%Eastern_Edge(tile)) THEN
580 i=bounds(ng)%edge(ieast,r2dvar)
581 DO j=jstrt,jendt
582 hwater=h(i,j)
583# ifdef ICESHELF
584 hwater=hwater-abs(zice(i,j))
585# endif
586 tl_hwater=tl_h(i,j)
587 hinv=1.0_r8/hwater
588# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
589 tl_hinv=-hinv*hinv*tl_hwater
590# endif
591!^ Zw(0)=-h(i,j)
592!^
593 tl_zw(0)=-tl_h(i,j)
594 DO k=1,n(ng)
595 cff_w=hc(ng)*(scalars(ng)%sc_w(k)-scalars(ng)%Cs_w(k))
596 cff1_w=scalars(ng)%Cs_w(k)
597 z_w0=cff_w+cff1_w*hwater
598# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
599 tl_z_w0=cff1_w*tl_hwater
600# endif
601!^ Zw(k)=z_w0+BOUNDARY(ng)%zeta_east(j)*(1.0_r8+z_w0*hinv)
602!^
603# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
604 tl_zw(k)=tl_z_w0+ &
605 & boundary(ng)%tl_zeta_east(j)* &
606 & (1.0_r8+z_w0*hinv)+ &
607 & boundary(ng)%zeta_east(j)*z_w0*tl_hinv+ &
608 & boundary(ng)%zeta_east(j)*tl_z_w0*hinv
609# else
610 tl_zw(k)=boundary(ng)%tl_zeta_east(j)*(1.0_r8+z_w0*hinv)
611# endif
612# ifdef ICESHELF
613!^ Zw(k)=Zw(k)-ABS(zice(i,j))
614# endif
615!^ Hz_bry(j,k,ieast)=Zw(k)-Zw(k-1)
616!^
617 tl_hz_bry(j,k,ieast)=tl_zw(k)-tl_zw(k-1)
618 END DO
619 END DO
620 END IF
621
622 IF (tl_lbc(isouth,isfsur,ng)%acquire.and. &
623 & domain(ng)%Southern_Edge(tile)) THEN
624 j=bounds(ng)%edge(isouth,r2dvar)
625 DO i=istrt,iendt
626 hwater=h(i,j)
627# ifdef ICESHELF
628 hwater=hwater-abs(zice(i,j))
629# endif
630 tl_hwater=tl_h(i,j)
631 hinv=1.0_r8/hwater
632# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
633 tl_hinv=-hinv*hinv*tl_hwater
634# endif
635!^ Zw(0)=-h(i,j)
636!^
637 tl_zw(0)=-tl_h(i,j)
638 DO k=1,n(ng)
639 cff_w=hc(ng)*(scalars(ng)%sc_w(k)-scalars(ng)%Cs_w(k))
640 cff1_w=scalars(ng)%Cs_w(k)
641 z_w0=cff_w+cff1_w*hwater
642# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
643 tl_z_w0=cff1_w*tl_hwater
644# endif
645!^ Zw(k)=z_w0+BOUNDARY(ng)%zeta_south(i)*(1.0_r8+z_w0*hinv)
646!^
647# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
648 tl_zw(k)=tl_z_w0+ &
649 & boundary(ng)%tl_zeta_south(i)* &
650 & (1.0_r8+z_w0*hinv)+ &
651 & boundary(ng)%zeta_south(i)*z_w0*tl_hinv+ &
652 & boundary(ng)%zeta_south(i)*tl_z_w0*hinv
653# else
654 tl_zw(k)=boundary(ng)%tl_zeta_south(i)*(1.0_r8+z_w0*hinv)
655# endif
656# ifdef ICESHELF
657!^ Zw(k)=Zw(k)-ABS(zice(i,j))
658# endif
659!^ Hz_bry(i,k,isouth)=Zw(k)-Zw(k-1)
660!^
661 tl_hz_bry(i,k,isouth)=tl_zw(k)-tl_zw(k-1)
662 END DO
663 END DO
664 END IF
665
666 IF (tl_lbc(inorth,isfsur,ng)%acquire.and. &
667 & domain(ng)%Northern_Edge(tile)) THEN
668 j=bounds(ng)%edge(inorth,r2dvar)
669 DO i=istrt,iendt
670 hwater=h(i,j)
671# ifdef ICESHELF
672 hwater=hwater-abs(zice(i,j))
673# endif
674 tl_hwater=tl_h(i,j)
675 hinv=1.0_r8/hwater
676# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
677 tl_hinv=-hinv*hinv*tl_hwater
678# endif
679!^ Zw(0)=-h(i,j)
680!^
681 tl_zw(0)=-tl_h(i,j)
682 DO k=1,n(ng)
683 cff_w=hc(ng)*(scalars(ng)%sc_w(k)-scalars(ng)%Cs_w(k))
684 cff1_w=scalars(ng)%Cs_w(k)
685 z_w0=cff_w+cff1_w*hwater
686# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
687 tl_z_w0=cff1_w*tl_hwater
688# endif
689!^ Zw(k)=z_w0+BOUNDARY(ng)%zeta_north(i)*(1.0_r8+z_w0*hinv)
690!^
691# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
692 tl_zw(k)=tl_z_w0+ &
693 & boundary(ng)%tl_zeta_north(i)* &
694 & (1.0_r8+z_w0*hinv)+ &
695 & boundary(ng)%zeta_north(i)*z_w0*tl_hinv+ &
696 & boundary(ng)%zeta_north(i)*tl_z_w0*hinv
697# else
698 tl_zw(k)=boundary(ng)%tl_zeta_north(i)*(1.0_r8+z_w0*hinv)
699# endif
700# ifdef ICESHELF
701!^ Zw(k)=Zw(k)-ABS(zice(i,j))
702# endif
703!^ Hz_bry(i,k,inorth)=Zw(k)-Zw(k-1)
704!^
705 tl_hz_bry(i,k,inorth)=tl_zw(k)-tl_zw(k-1)
706 END DO
707 END DO
708 END IF
709!
710!-----------------------------------------------------------------------
711! New formulation: Compute vertical depths (meters, negative) at
712! RHO- and W-points, and vertical grid thicknesses.
713! Various stretching functions are possible.
714!
715! z_w(x,y,s,t) = zeta(x,y,t) + [zeta(x,y,t)+ h(x,y)] * Zo_w
716!
717! Zo_w = [hc * s(k) + C(k) * h(x,y)] / [hc + h(x,y)]
718!
719!-----------------------------------------------------------------------
720!
721 ELSE IF (vtransform(ng).eq.2) THEN
722
723 IF (tl_lbc(iwest,isfsur,ng)%acquire.and. &
724 & domain(ng)%Western_Edge(tile)) THEN
725 i=bounds(ng)%edge(iwest,r2dvar)
726 DO j=jstrt,jendt
727 hwater=h(i,j)
728# ifdef ICESHELF
729 hwater=hwater-abs(zice(i,j))
730# endif
731 tl_hwater=tl_h(i,j)
732 hinv=1.0_r8/(hc(ng)+hwater)
733# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
734 tl_hinv=-hinv*hinv*tl_hwater
735# endif
736!^ Zw(0)=-h(i,j)
737!^
738 tl_zw(0)=-tl_h(i,j)
739 DO k=1,n(ng)
740 cff_w=hc(ng)*scalars(ng)%sc_w(k)
741 cff1_w=scalars(ng)%Cs_w(k)
742 cff2_w=(cff_w+cff1_w*hwater)*hinv
743# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
744 tl_cff2_w=cff1_w*tl_hwater*hinv+ &
745 & (cff_w+cff1_w*hwater)*tl_hinv
746# endif
747!^ Zw(k)=BOUNDARY(ng)%zeta_west(j)+
748!^ & (BOUNDARY(ng)%zeta_west(j)+hwater)*cff2_w
749!^
750# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
751 tl_zw(k)=boundary(ng)%tl_zeta_west(j)+ &
752 & (boundary(ng)%tl_zeta_west(j)+ &
753 & tl_hwater)*cff2_w+ &
754 & (boundary(ng)%zeta_west(j)+ &
755 & hwater)*tl_cff2_w
756# else
757 tl_zw(k)=boundary(ng)%tl_zeta_west(j)+ &
758 & boundary(ng)%tl_zeta_west(j)*cff2_w
759# endif
760# ifdef ICESHELF
761!^ Zw(k)=Zw(k)-ABS(zice(i,j))
762# endif
763!^ Hz_bry(j,k,iwest)=Zw(k)-Zw(k-1)
764!^
765 tl_hz_bry(j,k,iwest)=tl_zw(k)-tl_zw(k-1)
766 END DO
767 END DO
768 END IF
769
770 IF (tl_lbc(ieast,isfsur,ng)%acquire.and. &
771 & domain(ng)%Eastern_Edge(tile)) THEN
772 i=bounds(ng)%edge(ieast,r2dvar)
773 DO j=jstrt,jendt
774 hwater=h(i,j)
775# ifdef ICESHELF
776 hwater=hwater-abs(zice(i,j))
777# endif
778 tl_hwater=tl_h(i,j)
779 hinv=1.0_r8/(hc(ng)+hwater)
780# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
781 tl_hinv=-hinv*hinv*tl_hwater
782# endif
783!^ Zw(0)=-h(i,j)
784!^
785 tl_zw(0)=-tl_h(i,j)
786 DO k=1,n(ng)
787 cff_w=hc(ng)*scalars(ng)%sc_w(k)
788 cff1_w=scalars(ng)%Cs_w(k)
789 cff2_w=(cff_w+cff1_w*hwater)*hinv
790# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
791 tl_cff2_w=cff1_w*tl_hwater*hinv+ &
792 & (cff_w+cff1_w*hwater)*tl_hinv
793# endif
794!^ Zw(k)=BOUNDARY(ng)%zeta_east(j)+ &
795!^ & (BOUNDARY(ng)%zeta_east(j)+hwater)*cff2_w
796!^
797# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
798 tl_zw(k)=boundary(ng)%tl_zeta_east(j)+ &
799 & (boundary(ng)%tl_zeta_east(j)+ &
800 & tl_hwater)*cff2_w+ &
801 & (boundary(ng)%zeta_east(j)+ &
802 & hwater)*tl_cff2_w
803# else
804 tl_zw(k)=boundary(ng)%tl_zeta_east(j)+ &
805 & boundary(ng)%tl_zeta_east(j)*cff2_w
806# endif
807# ifdef ICESHELF
808!^ Zw(k)=Zw(k)-ABS(zice(i,j))
809# endif
810!^ Hz_bry(j,k,ieast)=Zw(k)-Zw(k-1)
811!^
812 tl_hz_bry(j,k,ieast)=tl_zw(k)-tl_zw(k-1)
813 END DO
814 END DO
815 END IF
816
817 IF (tl_lbc(isouth,isfsur,ng)%acquire.and. &
818 & domain(ng)%Southern_Edge(tile)) THEN
819 j=bounds(ng)%edge(isouth,r2dvar)
820 DO i=istrt,iendt
821 hwater=h(i,j)
822# ifdef ICESHELF
823 hwater=hwater-abs(zice(i,j))
824# endif
825 tl_hwater=tl_h(i,j)
826 hinv=1.0_r8/(hc(ng)+hwater)
827# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
828 tl_hinv=-hinv*hinv*tl_hwater
829# endif
830!^ Zw(0)=-h(i,j)
831!^
832 tl_zw(0)=-tl_h(i,j)
833 DO k=1,n(ng)
834 cff_w=hc(ng)*scalars(ng)%sc_w(k)
835 cff1_w=scalars(ng)%Cs_w(k)
836 cff2_w=(cff_w+cff1_w*hwater)*hinv
837# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
838 tl_cff2_w=cff1_w*tl_hwater*hinv+ &
839 & (cff_w+cff1_w*hwater)*tl_hinv
840# endif
841!^ Zw(k)=BOUNDARY(ng)%zeta_south(i)+ &
842!^ & (BOUNDARY(ng)%zeta_south(i)+hwater)*cff2_w
843!^
844# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
845 tl_zw(k)=boundary(ng)%tl_zeta_south(i)+ &
846 & (boundary(ng)%tl_zeta_south(i)+ &
847 & tl_hwater)*cff2_w+ &
848 & (boundary(ng)%zeta_south(i)+ &
849 & hwater)*tl_cff2_w
850# else
851 tl_zw(k)=boundary(ng)%tl_zeta_south(i)+ &
852 & boundary(ng)%tl_zeta_south(i)*cff2_w
853# endif
854# ifdef ICESHELF
855!^ Zw(k)=Zw(k)-ABS(zice(i,j))
856# endif
857!^ Hz_bry(i,k,isouth)=Zw(k)-Zw(k-1)
858!^
859 tl_hz_bry(i,k,isouth)=tl_zw(k)-tl_zw(k-1)
860 END DO
861 END DO
862 END IF
863
864 IF (tl_lbc(inorth,isfsur,ng)%acquire.and. &
865 & domain(ng)%Northern_Edge(tile)) THEN
866 j=bounds(ng)%edge(inorth,r2dvar)
867 DO i=istrt,iendt
868 hwater=h(i,j)
869# ifdef ICESHELF
870 hwater=hwater-abs(zice(i,j))
871# endif
872 tl_hwater=tl_h(i,j)
873 hinv=1.0_r8/(hc(ng)+hwater)
874# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
875 tl_hinv=-hinv*hinv*tl_hwater
876# endif
877!^ Zw(0)=-h(i,j)
878!^
879 tl_zw(0)=-tl_h(i,j)
880 DO k=1,n(ng)
881 cff_w=hc(ng)*scalars(ng)%sc_w(k)
882 cff1_w=scalars(ng)%Cs_w(k)
883 cff2_w=(cff_w+cff1_w*hwater)*hinv
884# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
885 tl_cff2_w=cff1_w*tl_hwater*hinv+ &
886 & (cff_w+cff1_w*hwater)*tl_hinv
887# endif
888!^ Zw(k)=BOUNDARY(ng)%zeta_north(i)+ &
889!^ & (BOUNDARY(ng)%zeta_north(i)+hwater)*cff2_w
890!^
891# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
892 tl_zw(k)=boundary(ng)%tl_zeta_north(i)+ &
893 & (boundary(ng)%tl_zeta_north(i)+ &
894 & tl_hwater)*cff2_w+ &
895 & (boundary(ng)%zeta_north(i)+ &
896 & hwater)*tl_cff2_w
897# else
898 tl_zw(k)=boundary(ng)%tl_zeta_north(i)+ &
899 & boundary(ng)%tl_zeta_north(i)*cff2_w
900# endif
901# ifdef ICESHELF
902!^ Zw(k)=Zw(k)-ABS(zice(i,j))
903# endif
904!^ Hz_bry(i,k,inorth)=Zw(k)-Zw(k-1)
905!^
906 tl_hz_bry(i,k,inorth)=tl_zw(k)-tl_zw(k-1)
907 END DO
908 END DO
909 END IF
910 END IF
911
912# ifdef DISTRIBUTE
913!
914!-----------------------------------------------------------------------
915! Exchange boundary information.
916!-----------------------------------------------------------------------
917!
918 DO ibry=1,4
919!^ CALL mp_exchange3d_bry (ng, tile, model, 1, ibry, &
920!^ & LBij, UBij, 1, N(ng), &
921!^ & NghostPoints, &
922!^ & EWperiodic(ng), NSperiodic(ng), &
923!^ & Hz_bry(:,:,ibry))
924!^
925 CALL mp_exchange3d_bry (ng, tile, model, 1, ibry, &
926 & lbij, ubij, 1, n(ng), &
927 & nghostpoints, &
928 & ewperiodic(ng), nsperiodic(ng), &
929 & tl_hz_bry(:,:,ibry))
930 END DO
931# endif
932!
933 RETURN
934 END SUBROUTINE rp_set_depth_bry_tile
935# endif
936# endif
937!
938!***********************************************************************
939 SUBROUTINE rp_bath (ng, tile)
940!***********************************************************************
941!
942 USE mod_param
943 USE mod_grid
944!
945! Imported variable declarations.
946!
947 integer, intent(in) :: ng, tile
948!
949! Local variable declarations.
950!
951 character (len=*), parameter :: myfile = &
952 & __FILE__//", rp_bath"
953!
954# include "tile.h"
955!
956# ifdef PROFILE
957 CALL wclock_on (ng, irpm, 12, __line__, myfile)
958# endif
959 CALL rp_bath_tile (ng, tile, &
960 & lbi, ubi, lbj, ubj, &
961 & imins, imaxs, jmins, jmaxs, &
962 & grid(ng) % h, &
963 & grid(ng) % tl_h)
964# ifdef PROFILE
965 CALL wclock_off (ng, irpm, 12, __line__, myfile)
966# endif
967!
968 RETURN
969 END SUBROUTINE rp_bath
970!
971!***********************************************************************
972 SUBROUTINE rp_bath_tile (ng, tile, &
973 & LBi, UBi, LBj, UBj, &
974 & IminS, ImaxS, JminS, JmaxS, &
975 & h, tl_h)
976!***********************************************************************
977!
978 USE mod_param
979 USE mod_scalars
980!
982# ifdef DISTRIBUTE
984# endif
985!
986! Imported variable declarations.
987!
988 integer, intent(in) :: ng, tile
989 integer, intent(in) :: lbi, ubi, lbj, ubj
990 integer, intent(in) :: imins, imaxs, jmins, jmaxs
991!
992# ifdef ASSUMED_SHAPE
993 real(r8), intent(in) :: h(lbi:,lbj:)
994 real(r8), intent(out) :: tl_h(lbi:,lbj:)
995# else
996 real(r8), intent(in) :: h(lbi:ubi,lbj:ubj)
997 real(r8), intent(out) :: tl_h(lbi:ubi,lbj:ubj)
998# endif
999!
1000! Local variable declarations.
1001!
1002 integer :: i, j
1003
1004# include "set_bounds.h"
1005!
1006!-----------------------------------------------------------------------
1007! Initialize tangent linear bathymetry tl_h(i,j) to h(i,j) so some of
1008! the terms are cancelled in the barotropic pressure gradient.
1009!-----------------------------------------------------------------------
1010!
1011 DO j=jstrt,jendt
1012 DO i=istrt,iendt
1013 tl_h(i,j)=h(i,j)
1014 END DO
1015 END DO
1016
1017 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1018 CALL exchange_r2d_tile (ng, tile, &
1019 & lbi, ubi, lbj, ubj, &
1020 & tl_h)
1021 END IF
1022
1023# ifdef DISTRIBUTE
1024 CALL mp_exchange2d (ng, tile, irpm, 1, &
1025 & lbi, ubi, lbj, ubj, &
1026 & nghostpoints, &
1027 & ewperiodic(ng), nsperiodic(ng), &
1028 & tl_h)
1029# endif
1030!
1031 RETURN
1032 END SUBROUTINE rp_bath_tile
1033#endif
1034 END MODULE rp_set_depth_mod
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
subroutine exchange_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
type(t_boundary), dimension(:), allocatable boundary
type(t_coupling), dimension(:), allocatable coupling
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
integer isfsur
integer, parameter irpm
Definition mod_param.F:664
integer, dimension(:), allocatable n
Definition mod_param.F:479
type(t_bounds), dimension(:), allocatable bounds
Definition mod_param.F:232
type(t_lbc), dimension(:,:,:), allocatable tl_lbc
Definition mod_param.F:379
integer nghostpoints
Definition mod_param.F:710
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer, parameter r2dvar
Definition mod_param.F:717
logical, dimension(:), allocatable ewperiodic
integer, parameter iwest
logical, dimension(:), allocatable nsperiodic
real(dp), dimension(:), allocatable hc
integer, parameter isouth
type(t_scalars), dimension(:), allocatable scalars
Definition mod_scalars.F:65
integer, parameter ieast
integer, parameter inorth
integer, dimension(:), allocatable vtransform
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable nstp
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine mp_exchange3d_bry(ng, tile, model, nvar, boundary, lbij, ubij, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public rp_set_depth_bry(ng, tile, model)
subroutine, public rp_set_depth(ng, tile, model)
subroutine rp_set_depth_bry_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, h, tl_h, zice, tl_hz_bry)
subroutine, public rp_bath_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, h, tl_h)
subroutine, public rp_set_depth_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, h, tl_h, zice, zt_avg1, tl_zt_avg1, tl_hz, tl_z_r, tl_z_w)
subroutine, public rp_bath(ng, tile)
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3