ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
prsgrd_mod Module Reference

Functions/Subroutines

subroutine, public prsgrd (ng, tile)
 
subroutine prsgrd31_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, umask_wet, vmask_wet, hz, om_v, on_u, z_r, z_w, rho, eq_tide, zetat, pair, diaru, diarv, ru, rv)
 
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)
 
subroutine prsgrd40_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, umask_wet, vmask_wet, om_v, on_u, hz, z_w, rho, eq_tide, zetat, pair, diaru, diarv, ru, rv)
 
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)
 
subroutine prsgrd44_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, umask_wet, vmask_wet, hz, om_v, on_u, z_w, rho, eq_tide, zetat, pair, diaru, diarv, ru, rv)
 

Function/Subroutine Documentation

◆ prsgrd()

subroutine public prsgrd_mod::prsgrd ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 35 of file prsgrd31.h.

36!***********************************************************************
37!
38 USE mod_param
39#ifdef DIAGNOSTICS
40 USE mod_diags
41#endif
42#ifdef ATM_PRESS
43 USE mod_forces
44#endif
45 USE mod_grid
46 USE mod_ocean
47 USE mod_stepping
48!
49! Imported variable declarations.
50!
51 integer, intent(in) :: ng, tile
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, inlm, 23, __line__, myfile)
62#endif
63 CALL prsgrd31_tile (ng, tile, &
64 & lbi, ubi, lbj, ubj, &
65 & imins, imaxs, jmins, jmaxs, &
66 & nrhs(ng), &
67#ifdef WET_DRY
68 & grid(ng)%umask_wet, &
69 & grid(ng)%vmask_wet, &
70#endif
71 & grid(ng) % Hz, &
72 & grid(ng) % om_v, &
73 & grid(ng) % on_u, &
74 & grid(ng) % z_r, &
75 & grid(ng) % z_w, &
76 & ocean(ng) % rho, &
77#ifdef TIDE_GENERATING_FORCES
78 & ocean(ng) % eq_tide, &
79#endif
80#ifdef WEC_VF
81 & ocean(ng) % zetat, &
82#endif
83#ifdef ATM_PRESS
84 & forces(ng) % Pair, &
85#endif
86#ifdef DIAGNOSTICS_UV
87 & diags(ng) % DiaRU, &
88 & diags(ng) % DiaRV, &
89#endif
90 & ocean(ng) % ru, &
91 & ocean(ng) % rv)
92#ifdef PROFILE
93 CALL wclock_off (ng, inlm, 23, __line__, myfile)
94#endif
95!
96 RETURN
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_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
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

References mod_diags::diags, mod_forces::forces, mod_grid::grid, mod_param::inlm, mod_stepping::nrhs, mod_ocean::ocean, prsgrd31_tile(), wclock_off(), and wclock_on().

Referenced by rhs3d_mod::rhs3d().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ prsgrd31_tile()

subroutine prsgrd_mod::prsgrd31_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) eq_tide,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zetat,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pair,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diaru,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diarv,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) ru,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) rv )
private

Definition at line 100 of file prsgrd31.h.

122!***********************************************************************
123!
124 USE mod_param
125 USE mod_scalars
126!
127! Imported variable declarations.
128!
129 integer, intent(in) :: ng, tile
130 integer, intent(in) :: LBi, UBi, LBj, UBj
131 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
132 integer, intent(in) :: nrhs
133
134#ifdef ASSUMED_SHAPE
135# ifdef WET_DRY
136 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
137 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
138# endif
139 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
140 real(r8), intent(in) :: om_v(LBi:,LBj:)
141 real(r8), intent(in) :: on_u(LBi:,LBj:)
142 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
143 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
144 real(r8), intent(in) :: rho(LBi:,LBj:,:)
145# ifdef TIDE_GENERATING_FORCES
146 real(r8), intent(in) :: eq_tide(LBi:,LBj:)
147# endif
148# ifdef WEC_VF
149 real(r8), intent(in) :: zetat(LBi:,LBj:)
150# endif
151# ifdef ATM_PRESS
152 real(r8), intent(in) :: Pair(LBi:,LBj:)
153# endif
154# ifdef DIAGNOSTICS_UV
155 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
156 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
157# endif
158 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
159 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
160#else
161# ifdef WET_DRY
162 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
163 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
164# endif
165 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
166 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
167 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
168 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
169 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
170 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
171# ifdef TIDE_GENERATING_FORCES
172 real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
173# endif
174# ifdef WEC_VF
175 real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
176# endif
177# ifdef ATM_PRESS
178 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
179# endif
180# ifdef DIAGNOSTICS_UV
181 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
182 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
183# endif
184 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
185 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
186#endif
187!
188! Local variable declarations.
189!
190 integer :: i, j, k
191 real(r8) :: fac, fac1, fac2, fac3
192 real(r8) :: cff1, cff2, cff3, cff4
193#ifdef WJ_GRADP
194 real(r8) :: gamma
195#endif
196
197 real(r8), dimension(IminS:ImaxS) :: phie
198 real(r8), dimension(IminS:ImaxS) :: phix
199
200#include "set_bounds.h"
201!
202!-----------------------------------------------------------------------
203! Calculate pressure gradient in the XI-direction (m4/s2).
204!-----------------------------------------------------------------------
205!
206! Compute surface baroclinic pressure gradient.
207!
208#ifdef ATM_PRESS
209 fac=100.0_r8/rho0
210#endif
211 fac1=0.5_r8*g/rho0
212 fac2=1000.0_r8*g/rho0
213 fac3=0.25_r8*g/rho0
214
215 DO j=jstr,jend
216 DO i=istru,iend
217#ifdef TIDE_GENERATING_FORCES
218 cff1=z_w(i ,j,n(ng))-eq_tide(i ,j)-z_r(i ,j,n(ng))+ &
219 & z_w(i-1,j,n(ng))-eq_tide(i-1,j)-z_r(i-1,j,n(ng))
220#else
221 cff1=z_w(i ,j,n(ng))-z_r(i ,j,n(ng))+ &
222 & z_w(i-1,j,n(ng))-z_r(i-1,j,n(ng))
223#endif
224 phix(i)=fac1*(rho(i,j,n(ng))-rho(i-1,j,n(ng)))*cff1
225#ifdef WEC_VF
226 phix(i)=phix(i)+zetat(i,j)-zetat(i-1,j)
227#endif
228#ifdef ATM_PRESS
229 phix(i)=phix(i)+fac*(pair(i,j)-pair(i-1,j))
230#endif
231#ifdef RHO_SURF
232 phix(i)=phix(i)+ &
233 & (fac2+fac1*(rho(i,j,n(ng))+rho(i-1,j,n(ng))))* &
234 & (z_w(i,j,n(ng))-z_w(i-1,j,n(ng)))
235#endif
236 ru(i,j,n(ng),nrhs)=-0.5_r8*(hz(i,j,n(ng))+hz(i-1,j,n(ng)))* &
237 & phix(i)*on_u(i,j)
238#ifdef WET_DRY
239 ru(i,j,n(ng),nrhs)=ru(i,j,n(ng),nrhs)*umask_wet(i,j)
240#endif
241#ifdef DIAGNOSTICS_UV
242 diaru(i,j,n(ng),nrhs,m3pgrd)=ru(i,j,n(ng),nrhs)
243#endif
244 END DO
245!
246! Compute interior baroclinic pressure gradient. Differentiate and
247! then vertically integrate.
248!
249 DO k=n(ng)-1,1,-1
250 DO i=istru,iend
251#ifdef WJ_GRADP
252 cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* &
253 & (z_r(i-1,j,k+1)-z_r(i-1,j,k)))
254 cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ &
255 & z_r(i ,j,k+1)-z_r(i-1,j,k+1)
256 cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- &
257 & z_r(i-1,j,k+1)+z_r(i-1,j,k )
258 gamma=0.125_r8*cff1*cff2*cff3
259
260 cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i-1,j,k+1))+ &
261 & (1.0_r8-gamma)*(rho(i,j,k )-rho(i-1,j,k ))
262 cff2=rho(i,j,k+1)+rho(i-1,j,k+1)- &
263 & rho(i,j,k )-rho(i-1,j,k )
264 cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
265 & z_r(i,j,k )-z_r(i-1,j,k )
266 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ &
267 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k ))
268 phix(i)=phix(i)+ &
269 & fac3*(cff1*cff3-cff2*cff4)
270#else
271 cff1=rho(i,j,k+1)-rho(i-1,j,k+1)+ &
272 & rho(i,j,k )-rho(i-1,j,k )
273 cff2=rho(i,j,k+1)+rho(i-1,j,k+1)- &
274 & rho(i,j,k )-rho(i-1,j,k )
275 cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
276 & z_r(i,j,k )-z_r(i-1,j,k )
277 cff4=z_r(i,j,k+1)-z_r(i-1,j,k+1)+ &
278 & z_r(i,j,k )-z_r(i-1,j,k )
279 phix(i)=phix(i)+ &
280 & fac3*(cff1*cff3-cff2*cff4)
281#endif
282 ru(i,j,k,nrhs)=-0.5_r8*(hz(i,j,k)+hz(i-1,j,k))* &
283 & phix(i)*on_u(i,j)
284#ifdef WET_DRY
285 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
286#endif
287#ifdef DIAGNOSTICS_UV
288 diaru(i,j,k,nrhs,m3pgrd)=ru(i,j,k,nrhs)
289#endif
290 END DO
291 END DO
292!
293!-----------------------------------------------------------------------
294! Calculate pressure gradient in the ETA-direction (m4/s2).
295!-----------------------------------------------------------------------
296!
297! Compute surface baroclinic pressure gradient.
298!
299 IF (j.ge.jstrv) THEN
300 DO i=istr,iend
301#ifdef TIDE_GENERATING_FORCES
302 cff1=z_w(i,j ,n(ng))-eq_tide(i,j )-z_r(i,j ,n(ng))+ &
303 & z_w(i,j-1,n(ng))-eq_tide(i,j-1)-z_r(i,j-1,n(ng))
304#else
305 cff1=z_w(i,j ,n(ng))-z_r(i,j ,n(ng))+ &
306 & z_w(i,j-1,n(ng))-z_r(i,j-1,n(ng))
307#endif
308 phie(i)=fac1*(rho(i,j,n(ng))-rho(i,j-1,n(ng)))*cff1
309#ifdef WEC_VF
310 phie(i)=phie(i)+zetat(i,j)-zetat(i,j-1)
311#endif
312#ifdef ATM_PRESS
313 phie(i)=phie(i)+fac*(pair(i,j)-pair(i,j-1))
314#endif
315#ifdef RHO_SURF
316 phie(i)=phie(i)+ &
317 & (fac2+fac1*(rho(i,j,n(ng))+rho(i,j-1,n(ng))))* &
318 & (z_w(i,j,n(ng))-z_w(i,j-1,n(ng)))
319#endif
320 rv(i,j,n(ng),nrhs)=-0.5_r8*(hz(i,j,n(ng))+hz(i,j-1,n(ng)))* &
321 & phie(i)*om_v(i,j)
322#ifdef WET_DRY
323 rv(i,j,n(ng),nrhs)=rv(i,j,n(ng),nrhs)*vmask_wet(i,j)
324#endif
325#ifdef DIAGNOSTICS_UV
326 diarv(i,j,n(ng),nrhs,m3pgrd)=rv(i,j,n(ng),nrhs)
327#endif
328 END DO
329!
330! Compute interior baroclinic pressure gradient. Differentiate and
331! then vertically integrate.
332!
333 DO k=n(ng)-1,1,-1
334 DO i=istr,iend
335#ifdef WJ_GRADP
336 cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* &
337 & (z_r(i,j-1,k+1)-z_r(i,j-1,k)))
338 cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ &
339 & z_r(i,j ,k+1)-z_r(i,j-1,k+1)
340 cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- &
341 & z_r(i,j-1,k+1)+z_r(i,j-1,k )
342 gamma=0.125_r8*cff1*cff2*cff3
343
344 cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i,j-1,k+1))+ &
345 & (1.0_r8-gamma)*(rho(i,j,k )-rho(i,j-1,k ))
346 cff2=rho(i,j,k+1)+rho(i,j-1,k+1)- &
347 & rho(i,j,k )-rho(i,j-1,k )
348 cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
349 & z_r(i,j,k )-z_r(i,j-1,k )
350 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ &
351 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k ))
352 phie(i)=phie(i)+ &
353 & fac3*(cff1*cff3-cff2*cff4)
354#else
355 cff1=rho(i,j,k+1)-rho(i,j-1,k+1)+ &
356 & rho(i,j,k )-rho(i,j-1,k )
357 cff2=rho(i,j,k+1)+rho(i,j-1,k+1)- &
358 & rho(i,j,k )-rho(i,j-1,k )
359 cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
360 & z_r(i,j,k )-z_r(i,j-1,k )
361 cff4=z_r(i,j,k+1)-z_r(i,j-1,k+1)+ &
362 & z_r(i,j,k )-z_r(i,j-1,k )
363 phie(i)=phie(i)+ &
364 & fac3*(cff1*cff3-cff2*cff4)
365#endif
366 rv(i,j,k,nrhs)=-0.5_r8*(hz(i,j,k)+hz(i,j-1,k))* &
367 & phie(i)*om_v(i,j)
368#ifdef WET_DRY
369 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
370#endif
371#ifdef DIAGNOSTICS_UV
372 diarv(i,j,k,nrhs,m3pgrd)=rv(i,j,k,nrhs)
373#endif
374 END DO
375 END DO
376 END IF
377 END DO
378!
379 RETURN
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
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

References mod_scalars::g, mod_scalars::m3pgrd, and mod_scalars::rho0.

Referenced by prsgrd().

Here is the caller graph for this function:

◆ prsgrd32_tile()

subroutine prsgrd_mod::prsgrd32_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) eq_tide,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zetat,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pair,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diaru,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diarv,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) ru,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) rv )
private

Definition at line 109 of file prsgrd32.h.

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

References mod_scalars::g, mod_scalars::m3pgrd, and mod_scalars::rho0.

◆ prsgrd40_tile()

subroutine prsgrd_mod::prsgrd40_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
umask_wet,
vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) eq_tide,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zetat,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pair,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diaru,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diarv,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) ru,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) rv )
private

Definition at line 97 of file prsgrd40.h.

120!***********************************************************************
121!
122 USE mod_param
123 USE mod_scalars
124!
125! Imported variable declarations.
126!
127 integer, intent(in) :: ng, tile
128 integer, intent(in) :: LBi, UBi, LBj, UBj
129 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
130 integer, intent(in) :: nrhs
131
132#ifdef ASSUMED_SHAPE
133 real(r8), intent(in) :: om_v(LBi:,LBj:)
134 real(r8), intent(in) :: on_u(LBi:,LBj:)
135 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
136 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
137 real(r8), intent(in) :: rho(LBi:,LBj:,:)
138# ifdef TIDE_GENERATING_FORCES
139 real(r8), intent(in) :: eq_tide(LBi:,LBj:)
140# endif
141# ifdef WEC_VF
142 real(r8), intent(in) :: zetat(LBi:,LBj:)
143# endif
144# ifdef ATM_PRESS
145 real(r8), intent(in) :: Pair(LBi:,LBj:)
146# endif
147# ifdef DIAGNOSTICS_UV
148 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
149 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
150# endif
151 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
152 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
153#else
154 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
155 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
156 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
157 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
158 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
159# ifdef TIDE_GENERATING_FORCES
160 real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
161# endif
162# ifdef WEC_VF
163 real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
164# endif
165# ifdef ATM_PRESS
166 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
167# endif
168# ifdef DIAGNOSTICS_UV
169 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
170 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
171# endif
172 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
173 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
174#endif
175!
176! Local variable declarations.
177!
178 integer :: i, j, k
179
180 real(r8) :: cff, cff1, dh
181#ifdef ATM_PRESS
182 real(r8) :: OneAtm, fac
183#endif
184#ifdef WEC_VF
185 real(r8) :: fac1
186#endif
187 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
188
189 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: FX
190
191 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: P
192
193#include "set_bounds.h"
194!
195!-----------------------------------------------------------------------
196! Finite Volume pressure gradient algorithm (Lin, 1997).
197!-----------------------------------------------------------------------
198!
199! Compute pressure and its vertical integral. Initialize pressure at
200! the free-surface as zero.
201!
202#ifdef ATM_PRESS
203 oneatm=1013.25_r8 ! 1 atm = 1013.25 mb
204 fac=100.0_r8/g
205#endif
206#ifdef WEC_VF
207 fac1=rho0/g
208#endif
209 j_loop : DO j=jstrv-1,jend
210 DO i=istru-1,iend
211 p(i,j,n(ng))=0.0_r8
212#ifdef WEC_VF
213 p(i,j,n(ng))=p(i,j,n(ng))+fac1*zetat(i,j)
214#endif
215#ifdef ATM_PRESS
216 p(i,j,n(ng))=p(i,j,n(ng))+fac*(pair(i,j)-oneatm)
217#endif
218#ifdef TIDE_GENERATING_FORCES
219 p(i,j,n(ng))=p(i,j,n(ng))-g*eq_tide(i,j)
220#endif
221 END DO
222 DO k=n(ng),1,-1
223 DO i=istru-1,iend
224 p(i,j,k-1)=p(i,j,k)+ &
225 & hz(i,j,k)*rho(i,j,k)
226 fx(i,j,k)=0.5_r8*hz(i,j,k)*(p(i,j,k)+p(i,j,k-1))
227 END DO
228 END DO
229!
230! Calculate pressure gradient in the XI-direction (m4/s2).
231!
232 IF (j.ge.jstr) THEN
233 DO i=istru,iend
234 fc(i,n(ng))=0.0_r8
235 END DO
236 cff=0.5_r8*g
237 cff1=g/rho0
238 DO k=n(ng),1,-1
239 DO i=istru,iend
240 dh=z_w(i,j,k-1)-z_w(i-1,j,k-1)
241 fc(i,k-1)=0.5_r8*dh*(p(i,j,k-1)+p(i-1,j,k-1))
242 ru(i,j,k,nrhs)=(cff*(hz(i-1,j,k)+ &
243 & hz(i ,j,k))* &
244 & (z_w(i-1,j,n(ng))- &
245 & z_w(i ,j,n(ng)))+ &
246 & cff1*(fx(i-1,j,k)- &
247 & fx(i ,j,k)+ &
248 & fc(i,k )- &
249 & fc(i,k-1)))*on_u(i,j)
250#ifdef WET_DRY
251 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
252#endif
253#ifdef DIAGNOSTICS_UV
254 diaru(i,j,k,nrhs,m3pgrd)=ru(i,j,k,nrhs)
255#endif
256 END DO
257 END DO
258 END IF
259!
260! Calculate pressure gradient in the ETA-direction (m4/s2).
261!
262 IF (j.ge.jstrv) THEN
263 DO i=istr,iend
264 fc(i,n(ng))=0.0_r8
265 END DO
266 cff=0.5_r8*g
267 cff1=g/rho0
268 DO k=n(ng),1,-1
269 DO i=istr,iend
270 dh=z_w(i,j,k-1)-z_w(i,j-1,k-1)
271 fc(i,k-1)=0.5_r8*dh*(p(i,j,k-1)+p(i,j-1,k-1))
272 rv(i,j,k,nrhs)=(cff*(hz(i,j-1,k)+ &
273 & hz(i,j ,k))* &
274 & (z_w(i,j-1,n(ng))- &
275 & z_w(i,j ,n(ng)))+ &
276 & cff1*(fx(i,j-1,k)- &
277 & fx(i,j ,k)+ &
278 & fc(i,k )- &
279 & fc(i,k-1)))*om_v(i,j)
280#ifdef WET_DRY
281 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
282#endif
283#ifdef DIAGNOSTICS_UV
284 diarv(i,j,k,nrhs,m3pgrd)=rv(i,j,k,nrhs)
285#endif
286 END DO
287 END DO
288 END IF
289 END DO j_loop
290!
291 RETURN

References mod_scalars::g, mod_scalars::m3pgrd, and mod_scalars::rho0.

◆ prsgrd42_tile()

subroutine prsgrd_mod::prsgrd42_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) eq_tide,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zetat,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pair,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diaru,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diarv,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) ru,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) rv )
private

Definition at line 107 of file prsgrd42.h.

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

References mod_scalars::g, mod_scalars::m3pgrd, and mod_scalars::rho0.

◆ prsgrd44_tile()

subroutine prsgrd_mod::prsgrd44_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) eq_tide,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zetat,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pair,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diaru,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2,ndrhs), intent(inout) diarv,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) ru,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng),2), intent(inout) rv )
private

Definition at line 106 of file prsgrd44.h.

128!***********************************************************************
129!
130 USE mod_param
131 USE mod_scalars
132!
133! Imported variable declarations.
134!
135 integer, intent(in) :: ng, tile
136 integer, intent(in) :: LBi, UBi, LBj, UBj
137 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
138 integer, intent(in) :: nrhs
139!
140#ifdef ASSUMED_SHAPE
141# ifdef WET_DRY
142 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
143 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
144# endif
145 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
146 real(r8), intent(in) :: om_v(LBi:,LBj:)
147 real(r8), intent(in) :: on_u(LBi:,LBj:)
148 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
149 real(r8), intent(in) :: rho(LBi:,LBj:,:)
150# ifdef TIDE_GENERATING_FORCES
151 real(r8), intent(in) :: eq_tide(LBi:,LBj:)
152# endif
153# ifdef WEC_VF
154 real(r8), intent(in) :: zetat(LBi:,LBj:)
155# endif
156# ifdef ATM_PRESS
157 real(r8), intent(in) :: Pair(LBi:,LBj:)
158# endif
159# ifdef DIAGNOSTICS_UV
160 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
161 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
162# endif
163 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
164 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
165#else
166# ifdef WET_DRY
167 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
168 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
169# endif
170 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
171 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
172 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
173 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
174 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
175# ifdef TIDE_GENERATING_FORCES
176 real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
177# endif
178# ifdef WEC_VF
179 real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
180# endif
181# ifdef ATM_PRESS
182 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
183# endif
184# ifdef DIAGNOSTICS_UV
185 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
186 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
187# endif
188 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
189 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
190#endif
191!
192! Local variable declarations.
193!
194 integer :: i, j, k
195
196 real(r8), parameter :: eps = 1.0e-8_r8
197
198 real(r8) :: Ampl, Hdd, cff, cff1, cff2, cff3, cffL, cffR
199 real(r8) :: deltaL, deltaR, dh, delP, limtr, rr
200#ifdef ATM_PRESS
201 real(r8) :: OneAtm, fac
202#endif
203 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: P
204 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: r
205 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: d
206
207 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: FX
208
209 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
210 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aL
211 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: aR
212 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dL
213 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dR
214 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: d1
215 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: r1
216
217#include "set_bounds.h"
218!
219!---------------------------------------------------------------------
220! Finite-volume pressure gradient force algorithm.
221!---------------------------------------------------------------------
222!
223#ifdef ATM_PRESS
224 oneatm=1013.25_r8 ! 1 atm = 1013.25 mb
225 fac=100.0_r8/g
226#endif
227 DO j=jstrv-1,jend
228 DO k=n(ng)-1,1,-1
229 DO i=istru-1,iend
230 fc(i,k)=1.0_r8/(hz(i,j,k+1)+hz(i,j,k))
231 r(i,j,k)=fc(i,k)*(rho(i,j,k+1)*hz(i,j,k )+ &
232 & rho(i,j,k )*hz(i,j,k+1))
233 d(i,j,k)=fc(i,k)*(rho(i,j,k+1)-rho(i,j,k))
234 END DO
235 END DO
236!
237! Parabolic WENO reconstruction of density field. Compute left and
238! right side limits aL and aR for the density assuming monotonized
239! quartic polynomial distributions within each grid box. Also
240! compute dL and dR which are used as a measure of quadratic
241! variation during subsquent WENO reconciliation of side limits.
242!
243 DO k=2,n(ng)-1
244 DO i=istru-1,iend
245 deltar=hz(i,j,k)*d(i,j,k )
246 deltal=hz(i,j,k)*d(i,j,k-1)
247 IF ((deltar*deltal).lt.0.0_r8) THEN
248 deltar=0.0_r8
249 deltal=0.0_r8
250 END IF
251 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
252 cffr=cff*d(i,j,k )
253 cffl=cff*d(i,j,k-1)
254 IF (abs(deltar).gt.abs(cffl)) deltar=cffl
255 IF (abs(deltal).gt.abs(cffr)) deltal=cffr
256 cff=(deltar-deltal)/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
257 deltar=deltar-cff*hz(i,j,k+1)
258 deltal=deltal+cff*hz(i,j,k-1)
259 ar(i,k)=rho(i,j,k)+deltar
260 al(i,k)=rho(i,j,k)-deltal
261 dr(i,k)=(2.0_r8*deltar-deltal)**2
262 dl(i,k)=(2.0_r8*deltal-deltar)**2
263 END DO
264 END DO
265!
266 DO i=istru-1,iend
267 al(i,n(ng))=ar(i,n(ng)-1)
268 ar(i,n(ng))=2.0_r8*rho(i,j,n(ng))-al(i,n(ng))
269 dr(i,n(ng))=(2.0_r8*ar(i,n(ng))+al(i,n(ng))- &
270 & 3.0_r8*rho(i,j,n(ng)))**2
271 dl(i,n(ng))=(3.0_r8*rho(i,j,n(ng))- &
272 & 2.0_r8*al(i,n(ng))-ar(i,n(ng)))**2
273 ar(i,1)=al(i,2)
274 al(i,1)=2.0_r8*rho(i,j,1)-ar(i,1)
275 dr(i,1)=(2.0_r8*ar(i,1)+al(i,1)-3.0_r8*rho(i,j,1))**2
276 dl(i,1)=(3.0_r8*rho(i,j,1)-2.0_r8*al(i,1)-ar(i,1))**2
277 END DO
278!
279 DO k=1,n(ng)-1
280 DO i=istru-1,iend
281 deltal=max(dl(i,k ),eps)
282 deltar=max(dr(i,k+1),eps)
283 r1(i,k)=(deltar*ar(i,k)+deltal*al(i,k+1))/(deltar+deltal)
284 END DO
285 END DO
286!
287 DO i=istru-1,iend
288#ifdef NEUMANN
289 r1(i,n(ng))=1.5_r8*rho(i,j,n(ng))-0.5_r8*r1(i,n(ng)-1)
290 r1(i,0)=1.5_r8*rho(i,j,1)-0.5_r8*r1(i,1)
291#else
292 r1(i,n(ng))=2.0_r8*rho(i,j,n(ng))-r1(i,n(ng)-1)
293 r1(i,0)=2.0_r8*rho(i,j,1)-r1(i,1)
294#endif
295 END DO
296!
297! Power-law reconciliation step. It starts with the computation of
298! side limits dR and dL of the first derivative assuming parabolic
299! distributions within each grid box. In this version of the code,
300! before doing so (see "else" branch of 3-way switch below), in the
301! situation when interfacial deviations deltaR and deltaL differ by
302! more than a factor of two (hence monotonic parabolic fit becomes
303! impossible), the parabolic assumption is switched to power-law
304! function, such that its derivative is zero at one end and,
305! consequently, larger than that of (would be) limited parabolic
306! on the other end. The basic parabolic version of the code is
307! commented out, but left here for reference.
308!
309 DO k=1,n(ng)
310 DO i=istru-1,iend
311!! cff=2.0_r8/Hz(i,j,k)
312!! dR(i,k)=cff*(2.0_r8*r1(i,k)+r1(i,k-1)-3.0_r8*rho(i,j,k))
313!! dL(i,k)=cff*(3.0_r8*rho(i,j,k)-2.0_r8*r1(i,k-1)-r1(i,k))
314!! cff=r(i,j,k)-r(i,j,k-1)
315!! if (cff*dR(i,k).lt.0.0_r8) dR(i,k)=0.0_r8
316!! if (cff*dL(i,k).lt.0.0_r8) dL(i,k)=0.0_r8
317 deltar=r1(i,k)-rho(i,j,k)
318 deltal=rho(i,j,k)-r1(i,k-1)
319 cff=deltar*deltal
320 IF (cff.gt.eps) THEN
321 cff=(deltar+deltal)/cff
322 ELSE
323 cff=0.0_r8
324 END IF
325 cffl=cff*deltal
326 cffr=cff*deltar
327 IF (cffl.gt.3.0_r8) THEN
328 cffl=cffl*deltal
329 cffr=0.0_r8
330 ELSE IF (cffr.gt.3.0_r8) THEN
331 cffl=0.0_r8
332 cffr=cffr*deltar
333 ELSE
334 cffl=4.0_r8*deltal-2.0_r8*deltar
335 cffr=4.0_r8*deltar-2.0_r8*deltal
336 END IF
337 cff=1.0_r8/hz(i,j,k)
338 dr(i,k)=cff*cffr
339 dl(i,k)=cff*cffl
340 END DO
341 END DO
342!
343! Compute final value of derivative at each interface by reconciling
344! two side limits dR(k) and dL(k+1) coming from adjacent grid boxes.
345! The difference between these two also causes change of interfacial
346! value r(k) by Ampl. The commented code (left here for reference)
347! computes the exact value of Ampl assuming power law reconciliation
348! and solving associated quadratic equation. The code segment below
349! corresponds to Pade fit to exact solution, which avoids computation
350! of SQRT for the sake of computational efficiency.
351!
352 DO k=n(ng)-1,1,-1
353 DO i=istru-1,iend
354 d(i,j,k)=fc(i,k)*(hz(i,j,k+1)*dl(i,k+1)+hz(i,j,k)*dr(i,k))
355 cffr=8.0_r8*(dr(i,k )+2.0_r8*dl(i,k ))
356 cffl=8.0_r8*(dl(i,k+1)+2.0_r8*dr(i,k+1))
357 IF (abs(d(i,j,k)).gt.abs(cffr)) d(i,j,k)=cffr
358 IF (abs(d(i,j,k)).gt.abs(cffl)) d(i,j,k)=cffl
359 IF ((dl(i,k+1)-dr(i,k))* &
360 & (rho(i,j,k+1)-rho(i,j,k)).gt.0.0_r8) THEN
361 hdd=hz(i,j,k)*(d(i,j,k)-dr(i,k))
362 rr=rho(i,j,k)-r1(i,k-1)
363 ELSE
364 hdd=hz(i,j,k+1)*(dl(i,k+1)-d(i,j,k))
365 rr=r1(i,k+1)-rho(i,j,k+1)
366 END IF
367 rr=abs(rr)
368!! Ampl=0.4_r8*Hdd*rr
369!! Hdd=ABS(Hdd)
370!! cff=rr*(rr+0.16_r8*Hdd)
371!! if (cff.gt.eps) Ampl=Ampl/(rr+sqrt(cff))
372 ampl=0.2_r8*hdd*rr
373 hdd=abs(hdd)
374 cff=rr*rr+0.0763636363636363636_r8*hdd* &
375 & (rr+0.004329004329004329_r8*hdd)
376 IF (cff.gt.eps) THEN
377 ampl=ampl*(rr+0.0363636363636363636_r8*hdd)/cff
378 ELSE
379 ampl=0.0_r8
380 END IF
381 r(i,j,k)=r1(i,k)+ampl
382 END DO
383 END DO
384 DO i=istru-1,iend
385#ifdef NEUMANN
386 r(i,j,0)=1.5_r8*rho(i,j,1)-0.5_r8*r(i,j,1)
387 r(i,j,n(ng))=1.5_r8*rho(i,j,n(ng))-0.5_r8*r(i,j,n(ng)-1)
388 d(i,j,0)=0.0_r8
389 d(i,j,n(ng))=0.0_r8
390#else
391 r(i,j,0)=2.0_r8*rho(i,j,1)-r(i,j,1)
392 r(i,j,n(ng))=2.0_r8*rho(i,j,n(ng))-r(i,j,n(ng)-1)
393 d(i,j,0)=d(i,j,1)
394 d(i,j,n(ng))=d(i,j,n(ng)-1)
395#endif
396 END DO
397!
398! Compute pressure (P) and lateral pressure force (FX). Initialize
399! pressure at the free-surface as zero
400!
401 DO i=istru-1,iend
402 p(i,j,n(ng))=0.0_r8
403#ifdef WEC_VF
404 p(i,j,n(ng))=p(i,j,n(ng))+zetat(i,j)
405#endif
406#ifdef ATM_PRESS
407 p(i,j,n(ng))=p(i,j,n(ng))+fac*(pair(i,j)-oneatm)
408#endif
409#ifdef TIDE_GENERATING_FORCES
410 p(i,j,n(ng))=p(i,j,n(ng))-g*eq_tide(i,j)
411#endif
412 END DO
413 cff3=1.0_r8/12.0_r8
414 DO k=n(ng),1,-1
415 DO i=istru-1,iend
416 p(i,j,k-1)=p(i,j,k)+hz(i,j,k)*rho(i,j,k)
417 fx(i,j,k)=0.5_r8*hz(i,j,k)* &
418 & (p(i,j,k)+p(i,j,k-1)+ &
419 & 0.2_r8*hz(i,j,k)* &
420 & (r(i,j,k)-r(i,j,k-1)- &
421 & cff3*hz(i,j,k)*(d(i,j,k)+d(i,j,k-1))))
422 END DO
423 END DO
424!
425! Compute net pressure gradient forces in the XI- and ETA-directions.
426!
427 IF (j.ge.jstr) THEN
428 DO i=istru,iend
429 fc(i,n(ng))=0.0_r8
430 END DO
431 cff=0.5_r8*g
432 cff1=g/rho0
433 cff2=1.0_r8/6.0_r8
434 cff3=1.0_r8/12.0_r8
435 DO k=n(ng),1,-1
436 DO i=istru,iend
437 dh=z_w(i,j,k-1)-z_w(i-1,j,k-1)
438 delp=p(i-1,j,k-1)-p(i,j,k-1)
439 rr=0.5_r8*dh*(r(i,j,k-1)+r(i-1,j,k-1)- &
440 & cff2*dh*(d(i,j,k-1)-d(i-1,j,k-1)))
441 limtr=2.0_r8*delp*rr
442 rr=rr*rr+delp*delp
443 IF (limtr.gt.eps*rr) THEN
444 limtr=limtr/rr
445 ELSE
446 limtr=0.0_r8
447 END IF
448 fc(i,k-1)=0.5_r8*dh* &
449 & (p(i,j,k-1)+p(i-1,j,k-1)+ &
450 & limtr*0.2_r8*dh* &
451 & (r(i,j,k-1)-r(i-1,j,k-1)- &
452 & cff3*dh*(d(i,j,k-1)+d(i-1,j,k-1))))
453 ru(i,j,k,nrhs)=(cff*(hz(i-1,j,k)+hz(i,j,k))* &
454 & (z_w(i-1,j,n(ng))-z_w(i,j,n(ng)))+ &
455 & cff1*(fx(i-1,j,k)-fx(i,j,k)+ &
456 & fc(i,k)-fc(i,k-1)))*on_u(i,j)
457#ifdef WET_DRY
458 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
459#endif
460#ifdef DIAGNOSTICS_UV
461 diaru(i,j,k,nrhs,m3pgrd)=ru(i,j,k,nrhs)
462#endif
463 END DO
464 END DO
465 END IF
466!
467 IF (j.ge.jstrv) THEN
468 DO i=istr,iend
469 fc(i,n(ng))=0.0_r8
470 END DO
471 cff=0.5_r8*g
472 cff1=g/rho0
473 cff2=1.0_r8/6.0_r8
474 cff3=1.0_r8/12.0_r8
475 DO k=n(ng),1,-1
476 DO i=istr,iend
477 dh=z_w(i,j,k-1)-z_w(i,j-1,k-1)
478 delp=p(i,j-1,k-1)-p(i,j,k-1)
479 rr=0.5_r8*dh*(r(i,j,k-1)+r(i,j-1,k-1)- &
480 & cff2*dh*(d(i,j,k-1)-d(i,j-1,k-1)))
481 limtr=2.0_r8*delp*rr
482 rr=rr*rr+delp*delp
483 IF (limtr.gt.eps*rr) THEN
484 limtr=limtr/rr
485 ELSE
486 limtr=0.0_r8
487 END IF
488 fc(i,k-1)=0.5_r8*dh* &
489 & (p(i,j,k-1)+p(i,j-1,k-1)+ &
490 & limtr*0.2_r8*dh* &
491 & (r(i,j,k-1)-r(i,j-1,k-1)- &
492 & cff3*dh*(d(i,j,k-1)+d(i,j-1,k-1))))
493 rv(i,j,k,nrhs)=(cff*(hz(i,j-1,k)+hz(i,j,k))* &
494 & (z_w(i,j-1,n(ng))-z_w(i,j,n(ng)))+ &
495 & cff1*(fx(i,j-1,k)-fx(i,j,k)+ &
496 & fc(i,k)-fc(i,k-1)))*om_v(i,j)
497#ifdef WET_DRY
498 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
499#endif
500#ifdef DIAGNOSTICS_UV
501 diarv(i,j,k,nrhs,m3pgrd)=rv(i,j,k,nrhs)
502#endif
503 END DO
504 END DO
505 END IF
506 END DO
507!
508 RETURN

References mod_scalars::g, mod_scalars::m3pgrd, and mod_scalars::rho0.