ROMS
Loading...
Searching...
No Matches
prsgrd42.h
Go to the documentation of this file.
1#undef NEUMANN
2 MODULE prsgrd_mod
3!
4!git $Id$
5!=======================================================================
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md Hernan G. Arango !
9!========================================== Alexander F. Shchepetkin ===
10! !
11! This subroutine evaluates the baroclinic, hydrostatic pressure !
12! gradient term using a finite-volume pressure Jacobian scheme. !
13! The scheme is based on local, conservative, limited-oscillation !
14! vertical quartic polynomial reconstruction of density field and !
15! subsequent projection fits of the derivatives of density into !
16! isosurface of vertical coordinate. The monotonicity constraint !
17! uses a PPM-style limitting algorithm. !
18! !
19! The pressure gradient terms (m4/s2) are loaded into right-hand- !
20! side arrays "ru" and "rv". !
21! !
22! Reference: !
23! !
24! Shchepetkin A.F and J.C. McWilliams, 2003: A method for !
25! computing horizontal pressure gradient force in an ocean !
26! model with non-aligned vertical coordinate, JGR, 108, !
27! 1-34. !
28! !
29!=======================================================================
30!
31 implicit none
32!
33 PRIVATE
34 PUBLIC :: prsgrd
35!
36 CONTAINS
37!
38!***********************************************************************
39 SUBROUTINE prsgrd (ng, tile)
40!***********************************************************************
41!
42 USE mod_param
43#ifdef DIAGNOSTICS
44 USE mod_diags
45#endif
46#ifdef ATM_PRESS
47 USE mod_forces
48#endif
49 USE mod_grid
50 USE mod_ocean
51 USE mod_stepping
52!
53! Imported variable declarations.
54!
55 integer, intent(in) :: ng, tile
56!
57! Local variable declarations.
58!
59 character (len=*), parameter :: MyFile = &
60 & __FILE__
61!
62#include "tile.h"
63!
64#ifdef PROFILE
65 CALL wclock_on (ng, inlm, 23, __line__, myfile)
66#endif
67 CALL prsgrd42_tile (ng, tile, &
68 & lbi, ubi, lbj, ubj, &
69 & imins, imaxs, jmins, jmaxs, &
70 & nrhs(ng), &
71#ifdef MASKING
72 & grid(ng) % umask, &
73 & grid(ng) % vmask, &
74#endif
75#ifdef WET_DRY
76 & grid(ng)%umask_wet, &
77 & grid(ng)%vmask_wet, &
78#endif
79 & grid(ng) % Hz, &
80 & grid(ng) % om_v, &
81 & grid(ng) % on_u, &
82 & grid(ng) % z_w, &
83 & ocean(ng) % rho, &
84#ifdef TIDE_GENERATING_FORCES
85 & ocean(ng) % eq_tide, &
86#endif
87#ifdef WEC_VF
88 & ocean(ng) % zetat, &
89#endif
90#ifdef ATM_PRESS
91 & forces(ng) % Pair, &
92#endif
93#ifdef DIAGNOSTICS_UV
94 & diags(ng) % DiaRU, &
95 & diags(ng) % DiaRV, &
96#endif
97 & ocean(ng) % ru, &
98 & ocean(ng) % rv)
99#ifdef PROFILE
100 CALL wclock_off (ng, inlm, 23, __line__, myfile)
101#endif
102!
103 RETURN
104 END SUBROUTINE prsgrd
105!
106!***********************************************************************
107 SUBROUTINE prsgrd42_tile (ng, tile, &
108 & LBi, UBi, LBj, UBj, &
109 & IminS, ImaxS, JminS, JmaxS, &
110 & nrhs, &
111#ifdef MASKING
112 & umask, vmask, &
113#endif
114#ifdef WET_DRY
115 & umask_wet, vmask_wet, &
116#endif
117 & Hz, om_v, on_u, z_w, &
118 & rho, &
119#ifdef TIDE_GENERATING_FORCES
120 & eq_tide, &
121#endif
122#ifdef WEC_VF
123 & zetat, &
124#endif
125#ifdef ATM_PRESS
126 & Pair, &
127#endif
128#ifdef DIAGNOSTICS_UV
129 & DiaRU, DiaRV, &
130#endif
131 & ru, rv)
132!***********************************************************************
133!
134 USE mod_param
135 USE mod_scalars
136!
137! Imported variable declarations.
138!
139 integer, intent(in) :: ng, tile
140 integer, intent(in) :: LBi, UBi, LBj, UBj
141 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
142 integer, intent(in) :: nrhs
143!
144#ifdef ASSUMED_SHAPE
145# ifdef MASKING
146 real(r8), intent(in) :: umask(LBi:,LBj:)
147 real(r8), intent(in) :: vmask(LBi:,LBj:)
148# endif
149# ifdef WET_DRY
150 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
151 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
152# endif
153 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
154 real(r8), intent(in) :: om_v(LBi:,LBj:)
155 real(r8), intent(in) :: on_u(LBi:,LBj:)
156 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
157 real(r8), intent(in) :: rho(LBi:,LBj:,:)
158# ifdef TIDE_GENERATING_FORCES
159 real(r8), intent(in) :: eq_tide(LBi:,LBj:)
160# endif
161# ifdef WEC_VF
162 real(r8), intent(in) :: zetat(LBi:,LBj:)
163# endif
164# ifdef ATM_PRESS
165 real(r8), intent(in) :: Pair(LBi:,LBj:)
166# endif
167# ifdef DIAGNOSTICS_UV
168 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
169 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
170# endif
171 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
172 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
173#else
174# ifdef MASKING
175 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
176 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
177# endif
178# ifdef WET_DRY
179 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
180 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
181# endif
182 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
183 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
184 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
185 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
186 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
187# ifdef TIDE_GENERATING_FORCES
188 real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
189# endif
190# ifdef WEC_VF
191 real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
192# endif
193# ifdef ATM_PRESS
194 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
195# endif
196# ifdef DIAGNOSTICS_UV
197 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
198 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
199# endif
200 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
201 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
202#endif
203!
204! Local variable declarations.
205!
206 integer :: i, j, k
207
208 real(r8), parameter :: eps = 1.0e-8_r8
209
210 real(r8) :: cff, cff1, cff2, cffL, cffR
211 real(r8) :: deltaL, deltaR, dh, delP, rr
212#ifdef ATM_PRESS
213 real(r8) :: OneAtm, fac
214#endif
215 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FX
216 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: P
217 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: r
218
219 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
220 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aL
221 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aR
222 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dL
223 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dR
224
225#include "set_bounds.h"
226!
227!---------------------------------------------------------------------
228! Finite-volume pressure gradient force algorithm.
229!---------------------------------------------------------------------
230!
231#ifdef ATM_PRESS
232 oneatm=1013.25_r8 ! 1 atm = 1013.25 mb
233 fac=100.0_r8/g
234#endif
235 cff2=1.0_r8/6.0_r8
236 DO j=jstrv-2,jend+1
237 DO k=n(ng)-1,1,-1
238 DO i=istru-2,iend+1
239 fc(i,k)=(rho(i,j,k+1)-rho(i,j,k))/(hz(i,j,k+1)+hz(i,j,k))
240 END DO
241 END DO
242!
243! Parabolic WENO reconstruction of density field. Compute left and
244! right side limits aL and aR for the density assuming monotonized
245! parabolic distributions within each grid box. Also compute dL and
246! dR which are used as a measure of quadratic variation during
247! subsquent WENO reconciliation of side limits.
248!
249 DO k=2,n(ng)-1
250 DO i=istru-2,iend+1
251 deltar=hz(i,j,k)*fc(i,k)
252 deltal=hz(i,j,k)*fc(i,k-1)
253 IF ((deltar*deltal).lt.0.0_r8) THEN
254 deltar=0.0_r8
255 deltal=0.0_r8
256 END IF
257 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
258 cffr=cff*fc(i,k)
259 cffl=cff*fc(i,k-1)
260 IF (abs(deltar).gt.abs(cffl)) deltar=cffl
261 IF (abs(deltal).gt.abs(cffr)) deltal=cffr
262 cff=(deltar-deltal)/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
263 deltar=deltar-cff*hz(i,j,k+1)
264 deltal=deltal+cff*hz(i,j,k-1)
265 ar(i,k)=rho(i,j,k)+deltar
266 al(i,k)=rho(i,j,k)-deltal
267 dr(i,k)=(2.0_r8*deltar-deltal)**2
268 dl(i,k)=(2.0_r8*deltal-deltar)**2
269 END DO
270 END DO
271!
272 DO i=istru-2,iend+1
273 al(i,n(ng))=ar(i,n(ng)-1)
274 ar(i,n(ng))=2.0_r8*rho(i,j,n(ng))-al(i,n(ng))
275 dr(i,n(ng))=(2.0_r8*ar(i,n(ng))+al(i,n(ng))- &
276 & 3.0_r8*rho(i,j,n(ng)))**2
277 dl(i,n(ng))=(3.0_r8*rho(i,j,n(ng))- &
278 & 2.0_r8*al(i,n(ng))-ar(i,n(ng)))**2
279 ar(i,1)=al(i,2)
280 al(i,1)=2.0_r8*rho(i,j,1)-ar(i,1)
281 dr(i,1)=(2.0_r8*ar(i,1)+al(i,1)-3.0_r8*rho(i,j,1))**2
282 dl(i,1)=(3.0_r8*rho(i,j,1)-2.0_r8*al(i,1)-ar(i,1))**2
283 END DO
284!
285 DO k=1,n(ng)-1
286 DO i=istru-2,iend+1
287 deltal=max(dl(i,k ),eps)
288 deltar=max(dr(i,k+1),eps)
289 r(i,j,k)=(deltar*ar(i,k)+deltal*al(i,k+1))/ &
290 & (deltar+deltal)
291 END DO
292 END DO
293!
294 DO i=istru-2,iend+1
295#ifdef NEUMANN
296 r(i,j,n(ng))=1.5_r8*rho(i,j,n(ng))-0.5_r8*r(i,j,n(ng)-1)
297 r(i,j,0)=1.5_r8*rho(i,j,1)-0.5_r8*r(i,j,1 )
298#else
299 r(i,j,n(ng))=2.0_r8*rho(i,j,n(ng))-r(i,j,n(ng)-1)
300 r(i,j,0)=2.0_r8*rho(i,j,1)-r(i,j,1 )
301#endif
302 END DO
303!
304! Compute pressure (P) and lateral pressure force (FX). Initialize
305! pressure at the free-surface as zero
306!
307 DO i=istru-2,iend+1
308 p(i,j,n(ng))=0.0_r8
309#ifdef WEC_VF
310 p(i,j,n(ng))=p(i,j,n(ng))+zetat(i,j)
311#endif
312#ifdef ATM_PRESS
313 p(i,j,n(ng))=p(i,j,n(ng))+fac*(pair(i,j)-oneatm)
314#endif
315#ifdef TIDE_GENERATING_FORCES
316 p(i,j,n(ng))=p(i,j,n(ng))-g*eq_tide(i,j)
317#endif
318 END DO
319 DO k=n(ng),1,-1
320 DO i=istru-2,iend+1
321 p(i,j,k-1)=p(i,j,k)+hz(i,j,k)*rho(i,j,k)
322 deltar=r(i,j,k)-rho(i,j,k)
323 deltal=rho(i,j,k)-r(i,j,k-1)
324 IF ((deltar*deltal).lt.0.0_r8) THEN
325 rr=0.0_r8
326 ELSE IF (abs(deltar).gt.(2.0_r8*abs(deltal))) THEN
327 rr=3.0_r8*deltal
328 ELSE IF (abs(deltal).gt.(2.0_r8*abs(deltar))) THEN
329 rr=3.0_r8*deltar
330 ELSE
331 rr=deltar+deltal
332 END IF
333 fx(i,j,k)=0.5_r8*hz(i,j,k)* &
334 & (p(i,j,k)+p(i,j,k-1)+cff2*rr*hz(i,j,k))
335 END DO
336 END DO
337!
338! Compute net pressure gradient forces in the XI-directions.
339! Set pressure at free-surface as zero.
340!
341 IF ((j.ge.jstr).and.(j.le.jend)) THEN
342 DO i=istru-1,iend+1
343 fc(i,n(ng))=0.0_r8
344 END DO
345 DO k=n(ng),1,-1
346 DO i=istru-1,iend+1
347 delp=p(i-1,j,k-1)-p(i,j,k-1)
348 dh=z_w(i,j,k-1)-z_w(i-1,j,k-1)
349 deltar=dh*r(i,j,k-1)-delp
350 deltal=delp-dh*r(i-1,j,k-1)
351 IF ((deltar*deltal).lt.0.0_r8) THEN
352 rr=0.0_r8
353 ELSE IF (abs(deltar).gt.(2.0_r8*abs(deltal))) THEN
354 rr=3.0_r8*deltal
355 ELSE IF (abs(deltal).gt.(2.0_r8*abs(deltar))) THEN
356 rr=3.0_r8*deltar
357 ELSE
358 rr=deltar+deltal
359 END IF
360 fc(i,k-1)=0.5_r8*dh*(p(i,j,k-1)+p(i-1,j,k-1)+cff2*rr)
361 ru(i,j,k,nrhs)=2.0_r8*(fx(i-1,j,k)-fx(i,j,k)+ &
362 & fc(i,k)-fc(i,k-1))/ &
363 & (hz(i-1,j,k)+hz(i,j,k))
364#ifdef MASKING
365 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask(i,j)
366#endif
367#ifdef WET_DRY
368 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
369#endif
370 END DO
371 END DO
372 END IF
373!
374! Compute net pressure gradient forces in the ETA-directions.
375! Set pressure at free-surface as zero.
376!
377 IF (j.ge.jstrv-1) THEN
378 DO i=istr,iend
379 fc(i,n(ng))=0.0_r8
380 END DO
381 DO k=n(ng),1,-1
382 DO i=istr,iend
383 delp=p(i,j-1,k-1)-p(i,j,k-1)
384 dh=z_w(i,j,k-1)-z_w(i,j-1,k-1)
385 deltar=dh*r(i,j,k-1)-delp
386 deltal=delp-dh*r(i,j-1,k-1)
387 IF ((deltar*deltal).lt.0.0_r8) THEN
388 rr=0.0_r8
389 ELSE IF (abs(deltar).gt.(2.0_r8*abs(deltal))) THEN
390 rr=3.0_r8*deltal
391 ELSE IF (abs(deltal).gt.(2.0_r8*abs(deltar))) THEN
392 rr=3.0_r8*deltar
393 ELSE
394 rr=deltar+deltal
395 END IF
396 fc(i,k-1)=0.5_r8*dh*(p(i,j,k-1)+p(i,j-1,k-1)+cff2*rr)
397 rv(i,j,k,nrhs)=2.0_r8*(fx(i,j-1,k)-fx(i,j,k)+ &
398 & fc(i,k)-fc(i,k-1))/ &
399 & (hz(i,j-1,k)+hz(i,j,k))
400#ifdef MASKING
401 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask(i,j)
402#endif
403#ifdef WET_DRY
404 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
405#endif
406 END DO
407 END DO
408 END IF
409 END DO
410!
411 rr=g/(24.0_r8*rho0)
412 cff=0.5_r8*g
413 cff1=0.5_r8*g/rho0
414 DO j=jstr,jend
415 DO k=n(ng)-1,1,-1
416 DO i=istru,iend
417 dh=rr*(z_w(i,j,k)-z_w(i-1,j,k))
418 fc(i,k)=max(dh,0.0_r8)* &
419 & (ru(i,j,k+1,nrhs)+ru(i+1,j,k ,nrhs)- &
420 & ru(i,j,k ,nrhs)-ru(i-1,j,k+1,nrhs))+ &
421 & min(dh,0.0_r8)* &
422 & (ru(i,j,k ,nrhs)+ru(i+1,j,k+1,nrhs)- &
423 & ru(i,j,k+1,nrhs)-ru(i-1,j,k ,nrhs))
424 END DO
425 END DO
426 DO i=istru,iend
427 fc(i,n(ng))=0.0_r8
428 dh=rr*(z_w(i,j,0)-z_w(i-1,j,0))
429 fc(i,0)=max(dh,0.0_r8)* &
430 & (ru(i ,j,1,nrhs)-ru(i-1,j,1,nrhs))+ &
431 & min(dh,0.0_r8)* &
432 & (ru(i+1,j,1,nrhs)-ru(i ,j,1,nrhs))
433 END DO
434 DO k=1,n(ng)
435 DO i=istru,iend
436 ru(i,j,k,nrhs)=(cff*(z_w(i-1,j,n(ng))-z_w(i,j,n(ng)))+ &
437 & cff1*ru(i,j,k,nrhs))* &
438 & (hz(i-1,j,k)+hz(i,j,k))*on_u(i,j)+ &
439 & (fc(i,k)-fc(i,k-1))*on_u(i,j)
440#ifdef DIAGNOSTICS_UV
441 diaru(i,j,k,nrhs,m3pgrd)=ru(i,j,k,nrhs)
442#endif
443 END DO
444 END DO
445 END DO
446!
447 DO j=jstrv,jend
448 DO k=n(ng)-1,1,-1
449 DO i=istr,iend
450 dh=rr*(z_w(i,j,k)-z_w(i,j-1,k))
451 fx(i,j,k)=max(dh,0.0_r8)* &
452 & (rv(i,j,k+1,nrhs)+rv(i+1,j ,k ,nrhs)- &
453 & rv(i,j,k ,nrhs)-rv(i ,j-1,k+1,nrhs))+ &
454 & min(dh,0.0_r8)* &
455 & (rv(i,j,k ,nrhs)+rv(i+1,j ,k+1,nrhs)- &
456 & rv(i,j,k+1,nrhs)-rv(i ,j-1,k ,nrhs))
457 END DO
458 END DO
459 DO i=istr,iend
460 fx(i,j,n(ng))=0.0_r8
461 dh=rr*(z_w(i,j,0)-z_w(i,j-1,0))
462 fx(i,j,0)=max(dh,0.0_r8)* &
463 & (rv(i ,j,1,nrhs)-rv(i,j-1,1,nrhs))+ &
464 & min(dh,0.0_r8)* &
465 & (rv(i+1,j,1,nrhs)-rv(i,j ,1,nrhs))
466 END DO
467 END DO
468 DO j=jstrv,jend
469 DO k=1,n(ng)
470 DO i=istr,iend
471 rv(i,j,k,nrhs)=(cff*(z_w(i,j-1,n(ng))-z_w(i,j,n(ng)))+ &
472 & cff1*rv(i,j,k,nrhs))* &
473 & (hz(i,j-1,k)+hz(i,j,k))*om_v(i,j)+ &
474 & (fx(i,j,k)-fx(i,j,k-1))*om_v(i,j)
475#ifdef DIAGNOSTICS_UV
476 diarv(i,j,k,nrhs,m3pgrd)=rv(i,j,k,nrhs)
477#endif
478 END DO
479 END DO
480 END DO
481!
482 RETURN
483 END SUBROUTINE prsgrd42_tile
484
485 END MODULE prsgrd_mod
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_scalars
Definition mod_scalars.F:2
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::m3pgrd
integer m3pgrd
Definition mod_scalars.F:190
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
prsgrd_mod
Definition prsgrd31.h:1
prsgrd_mod::prsgrd42_tile
subroutine prsgrd42_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, umask, vmask, umask_wet, vmask_wet, hz, om_v, on_u, z_w, rho, eq_tide, zetat, pair, diaru, diarv, ru, rv)
Definition prsgrd42.h:132
prsgrd_mod::prsgrd
subroutine, public prsgrd(ng, tile)
Definition prsgrd31.h:36
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3