ROMS
Loading...
Searching...
No Matches
prsgrd32.h
Go to the documentation of this file.
1 MODULE prsgrd_mod
2!
3!git $Id$
4!=======================================================================
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md Hernan G. Arango !
8!========================================== Alexander F. Shchepetkin ===
9! !
10! This subroutine evaluates the nonlinear baroclinic, hydrostatic !
11! pressure gradient term using a nonconservative Density-Jacobian !
12! scheme, based on cubic polynomial fits for "rho" and "z_r" as !
13! functions of nondimensional coordinates (XI,ETA,s), that is, its !
14! respective array indices. The cubic polynomials are monotonized !
15! by using harmonic mean instead of linear averages to interpolate !
16! slopes. This scheme retains exact anti-symmetry: !
17! !
18! J(rho,z_r)=-J(z_r,rho). !
19! !
20! If parameter OneFifth (below) is set to zero, the scheme becomes !
21! identical to standard Jacobian. !
22! !
23! Reference: !
24! !
25! Shchepetkin A.F and J.C. McWilliams, 2003: A method for !
26! computing horizontal pressure gradient force in an ocean !
27! model with non-aligned vertical coordinate, JGR, 108, !
28! 1-34. !
29! !
30!=======================================================================
31!
32 implicit none
33!
34 PRIVATE
35 PUBLIC :: prsgrd
36!
37 CONTAINS
38!
39!***********************************************************************
40 SUBROUTINE prsgrd (ng, tile)
41!***********************************************************************
42!
43 USE mod_param
44#ifdef DIAGNOSTICS
45 USE mod_diags
46#endif
47#ifdef ATM_PRESS
48 USE mod_forces
49#endif
50 USE mod_grid
51 USE mod_ocean
52 USE mod_stepping
53!
54! Imported variable declarations.
55!
56 integer, intent(in) :: ng, tile
57!
58! Local variable declarations.
59!
60 character (len=*), parameter :: MyFile = &
61 & __FILE__
62!
63#include "tile.h"
64!
65#ifdef PROFILE
66 CALL wclock_on (ng, inlm, 23, __line__, myfile)
67#endif
68 CALL prsgrd32_tile (ng, tile, &
69 & lbi, ubi, lbj, ubj, &
70 & imins, imaxs, jmins, jmaxs, &
71 & nrhs(ng), &
72#ifdef MASKING
73 & grid(ng) % umask, &
74 & grid(ng) % vmask, &
75#endif
76#ifdef WET_DRY
77 & grid(ng)%umask_wet, &
78 & grid(ng)%vmask_wet, &
79#endif
80 & grid(ng) % om_v, &
81 & grid(ng) % on_u, &
82 & grid(ng) % Hz, &
83 & grid(ng) % z_r, &
84 & grid(ng) % z_w, &
85 & ocean(ng) % rho, &
86#ifdef TIDE_GENERATING_FORCES
87 & ocean(ng) % eq_tide, &
88#endif
89#ifdef WEC_VF
90 & ocean(ng) % zetat, &
91#endif
92#ifdef ATM_PRESS
93 & forces(ng) % Pair, &
94#endif
95#ifdef DIAGNOSTICS_UV
96 & diags(ng) % DiaRU, &
97 & diags(ng) % DiaRV, &
98#endif
99 & ocean(ng) % ru, &
100 & ocean(ng) % rv)
101#ifdef PROFILE
102 CALL wclock_off (ng, inlm, 23, __line__, myfile)
103#endif
104!
105 RETURN
106 END SUBROUTINE prsgrd
107!
108!***********************************************************************
109 SUBROUTINE prsgrd32_tile (ng, tile, &
110 & LBi, UBi, LBj, UBj, &
111 & IminS, ImaxS, JminS, JmaxS, &
112 & nrhs, &
113#ifdef MASKING
114 & umask, vmask, &
115#endif
116#ifdef WET_DRY
117 & umask_wet, vmask_wet, &
118#endif
119 & om_v, on_u, &
120 & Hz, z_r, z_w, &
121 & rho, &
122#ifdef TIDE_GENERATING_FORCES
123 & eq_tide, &
124#endif
125#ifdef WEC_VF
126 & zetat, &
127#endif
128#ifdef ATM_PRESS
129 & Pair, &
130#endif
131#ifdef DIAGNOSTICS_UV
132 & DiaRU, DiaRV, &
133#endif
134 & ru, rv)
135!***********************************************************************
136!
137 USE mod_param
138 USE mod_scalars
139!
140! Imported variable declarations.
141!
142 integer, intent(in) :: ng, tile
143 integer, intent(in) :: LBi, UBi, LBj, UBj
144 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
145 integer, intent(in) :: nrhs
146
147#ifdef ASSUMED_SHAPE
148# ifdef MASKING
149 real(r8), intent(in) :: umask(LBi:,LBj:)
150 real(r8), intent(in) :: vmask(LBi:,LBj:)
151# endif
152# ifdef WET_DRY
153 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
154 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
155# endif
156 real(r8), intent(in) :: om_v(LBi:,LBj:)
157 real(r8), intent(in) :: on_u(LBi:,LBj:)
158 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
159 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
160 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
161 real(r8), intent(in) :: rho(LBi:,LBj:,:)
162# ifdef TIDE_GENERATING_FORCES
163 real(r8), intent(in) :: eq_tide(LBi:,LBj:)
164# endif
165# ifdef WEC_VF
166 real(r8), intent(in) :: zetat(LBi:,LBj:)
167# endif
168# ifdef ATM_PRESS
169 real(r8), intent(in) :: Pair(LBi:,LBj:)
170# endif
171# ifdef DIAGNOSTICS_UV
172 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
173 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
174# endif
175 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
176 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
177#else
178# ifdef MASKING
179 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
180 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
181# endif
182# ifdef WET_DRY
183 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
184 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
185# endif
186 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
187 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
188 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
189 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
190 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
191 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
192# ifdef TIDE_GENERATING_FORCES
193 real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
194# endif
195# ifdef WEC_VF
196 real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
197# endif
198# ifdef ATM_PRESS
199 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
200# endif
201# ifdef DIAGNOSTICS_UV
202 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
203 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
204# endif
205 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
206 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
207#endif
208!
209! Local variable declarations.
210!
211 integer :: i, j, k
212
213 real(r8), parameter :: OneFifth = 0.2_r8
214 real(r8), parameter :: OneTwelfth = 1.0_r8/12.0_r8
215 real(r8), parameter :: eps = 1.0e-10_r8
216
217 real(r8) :: GRho, GRho0, HalfGRho
218 real(r8) :: cff, cff1, cff2
219#ifdef ATM_PRESS
220 real(r8) :: OneAtm, fac
221#endif
222 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: P
223
224 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dR
225 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dZ
226
227 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FC
228 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: aux
229 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dRx
230 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZx
231!
232#include "set_bounds.h"
233!
234!-----------------------------------------------------------------------
235! Preliminary step (same for XI- and ETA-components):
236!-----------------------------------------------------------------------
237!
238 grho=g/rho0
239 grho0=1000.0_r8*grho
240 halfgrho=0.5_r8*grho
241#ifdef ATM_PRESS
242 oneatm=1013.25_r8 ! 1 atm = 1013.25 mb
243 fac=100.0_r8/rho0
244#endif
245!
246! Compute kinematic pressure: P/rho0 (m2/s2).
247!
248 DO j=jstrv-1,jend
249 DO k=1,n(ng)-1
250 DO i=istru-1,iend
251 dr(i,k)=rho(i,j,k+1)-rho(i,j,k)
252 dz(i,k)=z_r(i,j,k+1)-z_r(i,j,k)
253 END DO
254 END DO
255 DO i=istru-1,iend
256 dr(i,n(ng))=dr(i,n(ng)-1)
257 dz(i,n(ng))=dz(i,n(ng)-1)
258 dr(i,0)=dr(i,1)
259 dz(i,0)=dz(i,1)
260 END DO
261 DO k=n(ng),1,-1
262 DO i=istru-1,iend
263 cff=2.0_r8*dr(i,k)*dr(i,k-1)
264 IF (cff.gt.eps) THEN
265 dr(i,k)=cff/(dr(i,k)+dr(i,k-1))
266 ELSE
267 dr(i,k)=0.0_r8
268 END IF
269 dz(i,k)=2.0_r8*dz(i,k)*dz(i,k-1)/(dz(i,k)+dz(i,k-1))
270 END DO
271 END DO
272 DO i=istru-1,iend
273 cff1=1.0_r8/(z_r(i,j,n(ng))-z_r(i,j,n(ng)-1))
274 cff2=0.5_r8*(rho(i,j,n(ng))-rho(i,j,n(ng)-1))* &
275 & (z_w(i,j,n(ng))-z_r(i,j,n(ng)))*cff1
276 p(i,j,n(ng))=g*z_w(i,j,n(ng))+ &
277#ifdef WEC_VF
278 & zetat(i,j)+ &
279#endif
280#ifdef ATM_PRESS
281 & fac*(pair(i,j)-oneatm)+ &
282#endif
283 & grho*(rho(i,j,n(ng))+cff2)* &
284 & (z_w(i,j,n(ng))-z_r(i,j,n(ng)))
285#ifdef TIDE_GENERATING_FORCES
286 p(i,j,n(ng))=p(i,j,n(ng))-g*eq_tide(i,j)
287#endif
288 END DO
289 DO k=n(ng)-1,1,-1
290 DO i=istru-1,iend
291 p(i,j,k)=p(i,j,k+1)+ &
292 & halfgrho*((rho(i,j,k+1)+rho(i,j,k))* &
293 & (z_r(i,j,k+1)-z_r(i,j,k))- &
294 & onefifth* &
295 & ((dr(i,k+1)-dr(i,k))* &
296 & (z_r(i,j,k+1)-z_r(i,j,k)- &
297 & onetwelfth* &
298 & (dz(i,k+1)+dz(i,k)))- &
299 & (dz(i,k+1)-dz(i,k))* &
300 & (rho(i,j,k+1)-rho(i,j,k)- &
301 & onetwelfth* &
302 & (dr(i,k+1)+dr(i,k)))))
303 END DO
304 END DO
305 END DO
306!
307!-----------------------------------------------------------------------
308! Compute XI-component pressure gradient term.
309!-----------------------------------------------------------------------
310!
311 DO k=n(ng),1,-1
312 DO j=jstr,jend
313 DO i=istru-1,iend+1
314 aux(i,j)=z_r(i,j,k)-z_r(i-1,j,k)
315#ifdef MASKING
316 aux(i,j)=aux(i,j)*umask(i,j)
317#endif
318 fc(i,j)=rho(i,j,k)-rho(i-1,j,k)
319#ifdef MASKING
320 fc(i,j)=fc(i,j)*umask(i,j)
321#endif
322 END DO
323 END DO
324!
325 DO j=jstr,jend
326 DO i=istru-1,iend
327 cff=2.0_r8*aux(i,j)*aux(i+1,j)
328 IF (cff.gt.eps) THEN
329 cff1=1.0_r8/(aux(i,j)+aux(i+1,j))
330 dzx(i,j)=cff*cff1
331 ELSE
332 dzx(i,j)=0.0_r8
333 END IF
334 cff1=2.0_r8*fc(i,j)*fc(i+1,j)
335 IF (cff1.gt.eps) THEN
336 cff2=1.0_r8/(fc(i,j)+fc(i+1,j))
337 drx(i,j)=cff1*cff2
338 ELSE
339 drx(i,j)=0.0_r8
340 END IF
341 END DO
342 END DO
343!
344 DO j=jstr,jend
345 DO i=istru,iend
346 ru(i,j,k,nrhs)=on_u(i,j)*0.5_r8* &
347 & (hz(i,j,k)+hz(i-1,j,k))* &
348 & (p(i-1,j,k)-p(i,j,k)- &
349 & halfgrho* &
350 & ((rho(i,j,k)+rho(i-1,j,k))* &
351 & (z_r(i,j,k)-z_r(i-1,j,k))- &
352 & onefifth* &
353 & ((drx(i,j)-drx(i-1,j))* &
354 & (z_r(i,j,k)-z_r(i-1,j,k)- &
355 & onetwelfth* &
356 & (dzx(i,j)+dzx(i-1,j)))- &
357 & (dzx(i,j)-dzx(i-1,j))* &
358 & (rho(i,j,k)-rho(i-1,j,k)- &
359 & onetwelfth* &
360 & (drx(i,j)+drx(i-1,j))))))
361#ifdef WET_DRY
362 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
363#endif
364#ifdef DIAGNOSTICS_UV
365 diaru(i,j,k,nrhs,m3pgrd)=ru(i,j,k,nrhs)
366#endif
367 END DO
368 END DO
369 END DO
370!
371!-----------------------------------------------------------------------
372! ETA-component pressure gradient term.
373!-----------------------------------------------------------------------
374!
375 DO k=n(ng),1,-1
376 DO j=jstrv-1,jend+1
377 DO i=istr,iend
378 aux(i,j)=z_r(i,j,k)-z_r(i,j-1,k)
379#ifdef MASKING
380 aux(i,j)=aux(i,j)*vmask(i,j)
381#endif
382 fc(i,j)=rho(i,j,k)-rho(i,j-1,k)
383#ifdef MASKING
384 fc(i,j)=fc(i,j)*vmask(i,j)
385#endif
386 END DO
387 END DO
388!
389 DO j=jstrv-1,jend
390 DO i=istr,iend
391 cff=2.0_r8*aux(i,j)*aux(i,j+1)
392 IF (cff.gt.eps) THEN
393 cff1=1.0_r8/(aux(i,j)+aux(i,j+1))
394 dzx(i,j)=cff*cff1
395 ELSE
396 dzx(i,j)=0.0_r8
397 END IF
398 cff1=2.0_r8*fc(i,j)*fc(i,j+1)
399 IF (cff1.gt.eps) THEN
400 cff2=1.0_r8/(fc(i,j)+fc(i,j+1))
401 drx(i,j)=cff1*cff2
402 ELSE
403 drx(i,j)=0.0_r8
404 END IF
405 END DO
406 END DO
407!
408 DO j=jstrv,jend
409 DO i=istr,iend
410 rv(i,j,k,nrhs)=om_v(i,j)*0.5_r8* &
411 & (hz(i,j,k)+hz(i,j-1,k))* &
412 & (p(i,j-1,k)-p(i,j,k)- &
413 & halfgrho* &
414 & ((rho(i,j,k)+rho(i,j-1,k))* &
415 & (z_r(i,j,k)-z_r(i,j-1,k))- &
416 & onefifth* &
417 & ((drx(i,j)-drx(i,j-1))* &
418 & (z_r(i,j,k)-z_r(i,j-1,k)- &
419 & onetwelfth* &
420 & (dzx(i,j)+dzx(i,j-1)))- &
421 & (dzx(i,j)-dzx(i,j-1))* &
422 & (rho(i,j,k)-rho(i,j-1,k)- &
423 & onetwelfth* &
424 & (drx(i,j)+drx(i,j-1))))))
425#ifdef WET_DRY
426 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
427#endif
428#ifdef DIAGNOSTICS_UV
429 diarv(i,j,k,nrhs,m3pgrd)=rv(i,j,k,nrhs)
430#endif
431 END DO
432 END DO
433 END DO
434!
435 RETURN
436 END SUBROUTINE prsgrd32_tile
437
438 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::prsgrd32_tile
subroutine prsgrd32_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, umask, vmask, umask_wet, vmask_wet, om_v, on_u, hz, z_r, z_w, rho, eq_tide, zetat, pair, diaru, diarv, ru, rv)
Definition prsgrd32.h:135
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