ROMS
Loading...
Searching...
No Matches
hsimt_tvd.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3 MODULE hsimt_tvd_mod
4#if defined NONLINEAR && defined SOLVE3D
5!
6!git $Id$
7!=======================================================================
8! Copyright (c) 2002-2025 The ROMS Group !
9! Licensed under a MIT/X style license Hui Wu !
10! See License_ROMS.md Tarandeep Kalra !
11!==================================================== John C. Warner ===
12! !
13! This routine computes horizontal tracer advection fluxes FX and FE !
14! using tthe third High-order Spatial Interpolation at the Middle !
15! Temporal level (HSIMT; Wu and Zhu, 2010) with a Total Variation !
16! Diminishing limiter scheme. !
17! !
18! Reference: !
19! !
20! Hui Wu and Jianrong Zhu, 2010: Advection scheme with 3rd !
21! high-order spatial interpolation at the middle temporal !
22! level and its application to saltwater intrusion in the !
23! Changjiang Estuary, Ocean Modelling 33, 33-51, !
24! doi:10.1016/j.ocemod.2009.12.001 !
25! !
26!=======================================================================
27!
28 implicit none
29
30 PUBLIC :: hsimt_tvd_tile
31
32 CONTAINS
33!
34!***********************************************************************
35 SUBROUTINE hsimt_tvd_tile (ng, tile, &
36 & LBi, UBi, LBj, UBj, &
37 & IminS, ImaxS, JminS, JmaxS, &
38# ifdef MASKING
39 & rmask, umask, vmask, &
40# endif
41# ifdef WET_DRY
42 & rmask_wet, umask_wet, vmask_wet, &
43# endif
44 & pm, pn, &
45 & Huon, Hvom, oHz, t, &
46 & FX, FE)
47!***********************************************************************
48!
49 USE mod_param
50 USE mod_ncparam
51 USE mod_scalars
52!
53! Imported variable declarations.
54!
55 integer, intent(in) :: ng, tile
56 integer, intent(in) :: lbi, ubi, lbj, ubj
57 integer, intent(in) :: imins, imaxs, jmins, jmaxs
58!
59# ifdef ASSUMED_SHAPE
60# ifdef MASKING
61 real(r8), intent(in) :: rmask(lbi:,lbj:)
62 real(r8), intent(in) :: umask(lbi:,lbj:)
63 real(r8), intent(in) :: vmask(lbi:,lbj:)
64# endif
65# ifdef WET_DRY
66 real(r8), intent(in) :: rmask_wet(lbi:,lbj:)
67 real(r8), intent(in) :: umask_wet(lbi:,lbj:)
68 real(r8), intent(in) :: vmask_wet(lbi:,lbj:)
69# endif
70 real(r8), intent(in) :: pm(lbi:,lbj:)
71 real(r8), intent(in) :: pn(lbi:,lbj:)
72 real(r8), intent(in) :: huon(lbi:,lbj:)
73 real(r8), intent(in) :: hvom(lbi:,lbj:)
74 real(r8), intent(in) :: ohz(imins:,jmins:)
75 real(r8), intent(in) :: t(lbi:,lbj:)
76 real(r8), intent(out) :: fx(imins:,jmins:)
77 real(r8), intent(out) :: fe(imins:,jmins:)
78# else
79# ifdef MASKING
80 real(r8), intent(in) :: rmask(lbi:ubi,lbj:ubj)
81 real(r8), intent(in) :: umask(lbi:ubi,lbj:ubj)
82 real(r8), intent(in) :: vmask(lbi:ubi,lbj:ubj)
83# endif
84# ifdef WET_DRY
85 real(r8), intent(in) :: rmask_wet(lbi:ubi,lbj:ubj)
86 real(r8), intent(in) :: umask_wet(lbi:ubi,lbj:ubj)
87 real(r8), intent(in) :: vmask_wet(lbi:ubi,lbj:ubj)
88# endif
89 real(r8), intent(in) :: pm(lbi:ubi,lbj:ubj)
90 real(r8), intent(in) :: pn(lbi:ubi,lbj:ubj)
91 real(r8), intent(in) :: huon(lbi:ubi,lbj:ubj)
92 real(r8), intent(in) :: hvom(lbi:ubi,lbj:ubj)
93 real(r8), intent(in) :: ohz(imins:imaxs,jmins:jmaxs)
94 real(r8), intent(in) :: t(lbi:ubi,lbj:ubj)
95 real(r8), intent(out) :: fx(imins:imaxs,jmins:jmaxs)
96 real(r8), intent(out) :: fe(imins:imaxs,jmins:jmaxs)
97# endif
98!
99! Local variable declarations.
100!
101 integer :: i, is, j, k, ii, jj
102
103 real(r8) :: cc1 = 0.25_r8
104 real(r8) :: cc2 = 0.5_r8
105 real(r8) :: cc3 = 1.0_r8/12.0_r8
106 real(r8) :: eps1 = 1.0e-12_r8
107
108 real(r8) :: cff, cff1
109 real(r8) :: betal, betar, betad, betau
110 real(r8) :: rl, rr, rd, ru, rkal, rkar, rkad, rkau
111 real(r8) :: a1, b1, sw_eta, sw_xi
112
113 real(r8), dimension(IminS:ImaxS) :: gradx, kax, okax
114 real(r8), dimension(JminS:JmaxS) :: grade, kae, okae
115
116# include "set_bounds.h"
117!
118!-----------------------------------------------------------------------
119! Compute tracer horizontal aadvective fluxes using the HSIMT scheme
120! (Wu and Zhu, 2010).
121!-----------------------------------------------------------------------
122!
123 DO j=jstr,jend
124 DO i=istru-1,iendp2
125 cff=0.125_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))*dt(ng)
126 cff1=cff*(ohz(i-1,j)+ohz(i,j))
127 gradx(i)=t(i,j)-t(i-1,j)
128 kax(i)=1.0_r8-abs(huon(i,j)*cff1)
129# ifdef MASKING
130 gradx(i)=gradx(i)*umask(i,j)
131 kax(i)=kax(i)*umask(i,j)
132# endif
133 END DO
134 IF (.not.ewperiodic(ng)) THEN
135 IF (domain(ng)%Western_Edge(tile)) THEN
136 IF (huon(istr,j).ge.0.0_r8) THEN
137 gradx(istr-1)=0.0_r8
138 kax(istr-1)=0.0_r8
139 END IF
140 END IF
141 IF (domain(ng)%Eastern_Edge(tile)) THEN
142 IF (huon(iend+1,j).lt.0.0_r8) THEN
143 gradx(iend+2)=0.0_r8
144 kax(iend+2)=0.0_r8
145 END IF
146 END IF
147 END IF
148 DO i=istr,iend+1
149 IF (kax(i).le.eps1) THEN
150 okax(i)=0.0_r8
151 ELSE
152 okax(i)=1.0_r8/max(kax(i),eps1)
153 END IF
154 IF (huon(i,j).ge.0.0_r8) THEN
155 IF (abs(gradx(i)).le.eps1) THEN
156 rl=0.0_r8
157 rkal=0.0_r8
158 ELSE
159 rl=gradx(i-1)/gradx(i)
160 rkal=kax(i-1)*okax(i)
161 END IF
162 a1= cc1*kax(i)+cc2-cc3*okax(i)
163 b1=-cc1*kax(i)+cc2+cc3*okax(i)
164 betal=a1+b1*rl
165 cff=0.5_r8*max(0.0_r8, &
166 & min(2.0_r8, 2.0_r8*rl*rkal, betal))* &
167 & gradx(i)*kax(i)
168# ifdef MASKING
169 ii=max(i-2,0)
170 cff=cff*rmask(ii,j)
171# endif
172 sw_xi=t(i-1,j)+cff
173 ELSE
174 IF (abs(gradx(i)).le.eps1) THEN
175 rr=0.0_r8
176 rkar=0.0_r8
177 ELSE
178 rr=gradx(i+1)/gradx(i)
179 rkar=kax(i+1)*okax(i)
180 END IF
181 a1= cc1*kax(i)+cc2-cc3*okax(i)
182 b1=-cc1*kax(i)+cc2+cc3*okax(i)
183 betar=a1+b1*rr
184 cff=0.5_r8*max(0.0_r8, &
185 & min(2.0_r8, 2.0_r8*rr*rkar, betar))* &
186 & gradx(i)*kax(i)
187# ifdef MASKING
188 ii=min(i+1,lm(ng)+1)
189 cff=cff*rmask(ii,j)
190# endif
191 sw_xi=t(i,j)-cff
192 END IF
193 fx(i,j)=sw_xi*huon(i,j)
194 END DO
195 END DO
196!
197 DO i=istr,iend
198 DO j=jstrv-1,jendp2
199 cff=0.125_r8*(pn(i,j)+pn(i,j-1))*(pm(i,j)+pm(i,j-1))*dt(ng)
200 cff1=cff*(ohz(i,j)+ohz(i,j-1))
201 grade(j)=t(i,j)-t(i,j-1)
202 kae(j)=1.0_r8-abs(hvom(i,j)*cff1)
203# ifdef MASKING
204 grade(j)=grade(j)*vmask(i,j)
205 kae(j)=kae(j)*vmask(i,j)
206# endif
207 END DO
208 IF (.not.nsperiodic(ng)) THEN
209 IF (domain(ng)%Southern_Edge(tile)) THEN
210 IF (hvom(i,jstr).ge.0.0_r8) THEN
211 grade(jstr-1)=0.0_r8
212 kae(jstr-1)=0.0_r8
213 END IF
214 END IF
215 IF (domain(ng)%Northern_Edge(tile)) THEN
216 IF (hvom(i,jend+1).lt.0.0_r8) THEN
217 grade(jend+2)=0.0_r8
218 kae(jend+2)=0.0_r8
219 END IF
220 END IF
221 END IF
222 DO j=jstr,jend+1
223 IF (kae(j).le.eps1) THEN
224 okae(j)=0.0_r8
225 ELSE
226 okae(j)=1.0_r8/max(kae(j),eps1)
227 END IF
228 IF (hvom(i,j).ge.0.0_r8) THEN
229 IF (abs(grade(j)).le.eps1) THEN
230 rd=0.0_r8
231 rkad=0.0_r8
232 ELSE
233 rd=grade(j-1)/grade(j)
234 rkad=kae(j-1)*okae(j)
235 END IF
236 a1= cc1*kae(j)+cc2-cc3*okae(j)
237 b1=-cc1*kae(j)+cc2+cc3*okae(j)
238 betad=a1+b1*rd
239 cff=0.5_r8*max(0.0_r8, &
240 & min(2.0_r8, 2.0_r8*rd*rkad, betad))* &
241 & grade(j)*kae(j)
242# ifdef MASKING
243 jj=max(j-2,0)
244 cff=cff*rmask(i,jj)
245# endif
246 sw_eta=t(i,j-1)+cff
247 ELSE
248 IF (abs(grade(j)).le.eps1) THEN
249 ru=0.0_r8
250 rkau=0.0_r8
251 ELSE
252 ru=grade(j+1)/grade(j)
253 rkau=kae(j+1)*okae(j)
254 END IF
255 a1= cc1*kae(j)+cc2-cc3*okae(j)
256 b1=-cc1*kae(j)+cc2+cc3*okae(j)
257 betau=a1+b1*ru
258 cff=0.5*max(0.0_r8, &
259 & min(2.0_r8, 2.0_r8*ru*rkau, betau))* &
260 & grade(j)*kae(j)
261# ifdef MASKING
262 jj=min(j+1,mm(ng)+1)
263 cff=cff*rmask(i,jj)
264# endif
265 sw_eta=t(i,j)-cff
266 END IF
267 fe(i,j)=sw_eta*hvom(i,j)
268 END DO
269 END DO
270!
271 RETURN
272 END SUBROUTINE hsimt_tvd_tile
273#endif
274 END MODULE hsimt_tvd_mod
hsimt_tvd_mod
Definition hsimt_tvd.F:3
hsimt_tvd_mod::hsimt_tvd_tile
subroutine, public hsimt_tvd_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, pm, pn, huon, hvom, ohz, t, fx, fe)
Definition hsimt_tvd.F:47
mod_ncparam
Definition mod_ncparam.F:2
mod_param
Definition mod_param.F:2
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::cc3
real(r8) cc3
Definition mod_scalars.F:395
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::cc1
real(r8) cc1
Definition mod_scalars.F:393
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::cc2
real(r8) cc2
Definition mod_scalars.F:394