1 | SUBROUTINE ana_psource (ng, tile, model)
|
---|
2 | !
|
---|
3 | !! svn $Id: ana_psource.h 48 2007-05-09 16:15:44Z arango $
|
---|
4 | !!======================================================================
|
---|
5 | !! Copyright (c) 2002-2007 The ROMS/TOMS Group !
|
---|
6 | !! Licensed under a MIT/X style license !
|
---|
7 | !! See License_ROMS.txt !
|
---|
8 | !! !
|
---|
9 | !=======================================================================
|
---|
10 | ! !
|
---|
11 | ! This subroutine sets analytical tracer and mass point Sources !
|
---|
12 | ! and/or Sinks. River runoff can be consider as a point source. !
|
---|
13 | ! !
|
---|
14 | !=======================================================================
|
---|
15 | !
|
---|
16 | USE mod_param
|
---|
17 | USE mod_grid
|
---|
18 | USE mod_ncparam
|
---|
19 | USE mod_ocean
|
---|
20 | USE mod_sources
|
---|
21 | USE mod_stepping
|
---|
22 | !
|
---|
23 | integer, intent(in) :: ng, tile, model
|
---|
24 |
|
---|
25 | integer :: LBi, UBi, LBj, UBj
|
---|
26 | !
|
---|
27 | LBi=LBOUND(GRID(ng)%h,DIM=1)
|
---|
28 | UBi=UBOUND(GRID(ng)%h,DIM=1)
|
---|
29 | LBj=LBOUND(GRID(ng)%h,DIM=2)
|
---|
30 | UBj=UBOUND(GRID(ng)%h,DIM=2)
|
---|
31 | !
|
---|
32 | CALL ana_psource_grid (ng, model, LBi, UBi, LBj, UBj, &
|
---|
33 | & nnew(ng), knew(ng), Nsrc(ng), &
|
---|
34 | & OCEAN(ng) % zeta, &
|
---|
35 | & OCEAN(ng) % ubar, &
|
---|
36 | & OCEAN(ng) % vbar, &
|
---|
37 | #ifdef SOLVE3D
|
---|
38 | & OCEAN(ng) % u, &
|
---|
39 | & OCEAN(ng) % v, &
|
---|
40 | & GRID(ng) % z_w, &
|
---|
41 | #endif
|
---|
42 | & GRID(ng) % h, &
|
---|
43 | & GRID(ng) % on_u, &
|
---|
44 | & GRID(ng) % om_v, &
|
---|
45 | & SOURCES(ng) % Isrc, &
|
---|
46 | & SOURCES(ng) % Jsrc, &
|
---|
47 | & SOURCES(ng) % Lsrc, &
|
---|
48 | & SOURCES(ng) % Dsrc, &
|
---|
49 | #ifdef SOLVE3D
|
---|
50 | # if defined UV_PSOURCE || defined Q_PSOURCE
|
---|
51 | & SOURCES(ng) % Qshape, &
|
---|
52 | & SOURCES(ng) % Qsrc, &
|
---|
53 | # endif
|
---|
54 | # ifdef TS_PSOURCE
|
---|
55 | & SOURCES(ng) % Tsrc, &
|
---|
56 | # endif
|
---|
57 | #endif
|
---|
58 | & SOURCES(ng) % Qbar)
|
---|
59 | !
|
---|
60 | ! Set analytical header file name used.
|
---|
61 | !
|
---|
62 | IF (Lanafile) THEN
|
---|
63 | ANANAME(20)='ROMS/Functionals/ana_psource.h'
|
---|
64 | END IF
|
---|
65 |
|
---|
66 | RETURN
|
---|
67 | END SUBROUTINE ana_psource
|
---|
68 | !
|
---|
69 | !***********************************************************************
|
---|
70 | SUBROUTINE ana_psource_grid (ng, model, LBi, UBi, LBj, UBj, &
|
---|
71 | & nnew, knew, Nsrc, &
|
---|
72 | & zeta, ubar, vbar, &
|
---|
73 | #ifdef SOLVE3D
|
---|
74 | & u, v, z_w, &
|
---|
75 | #endif
|
---|
76 | & h, on_u, om_v, &
|
---|
77 | & Isrc, Jsrc, Lsrc, Dsrc, &
|
---|
78 | #ifdef SOLVE3D
|
---|
79 | # if defined UV_PSOURCE || defined Q_PSOURCE
|
---|
80 | & Qshape, Qsrc, &
|
---|
81 | # endif
|
---|
82 | # ifdef TS_PSOURCE
|
---|
83 | & Tsrc, &
|
---|
84 | # endif
|
---|
85 | #endif
|
---|
86 | & Qbar)
|
---|
87 | !***********************************************************************
|
---|
88 | !
|
---|
89 | USE mod_param
|
---|
90 | USE mod_scalars
|
---|
91 | #ifdef SEDIMENT
|
---|
92 | USE mod_sediment
|
---|
93 | #endif
|
---|
94 | #ifndef ASSUMED_SHAPE
|
---|
95 | USE mod_sources, ONLY : Msrc
|
---|
96 | #endif
|
---|
97 | !
|
---|
98 | ! Imported variable declarations.
|
---|
99 | !
|
---|
100 | integer, intent(in) :: ng, model, LBi, UBi, LBj, UBj
|
---|
101 | integer, intent(in) :: nnew, knew
|
---|
102 |
|
---|
103 | integer, intent(out) :: Nsrc
|
---|
104 | !
|
---|
105 | #ifdef ASSUMED_SHAPE
|
---|
106 | logical, intent(out) :: Lsrc(:,:)
|
---|
107 |
|
---|
108 | integer, intent(out) :: Isrc(:)
|
---|
109 | integer, intent(out) :: Jsrc(:)
|
---|
110 |
|
---|
111 | real(r8), intent(in) :: zeta(LBi:,LBj:,:)
|
---|
112 | real(r8), intent(in) :: ubar(LBi:,LBj:,:)
|
---|
113 | real(r8), intent(in) :: vbar(LBi:,LBj:,:)
|
---|
114 | # ifdef SOLVE3D
|
---|
115 | real(r8), intent(in) :: u(LBi:,LBj:,:,:)
|
---|
116 | real(r8), intent(in) :: v(LBi:,LBj:,:,:)
|
---|
117 | real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
|
---|
118 | # endif
|
---|
119 | real(r8), intent(in) :: h(LBi:,LBj:)
|
---|
120 | real(r8), intent(in) :: on_u(LBi:,LBj:)
|
---|
121 | real(r8), intent(in) :: om_v(LBi:,LBj:)
|
---|
122 |
|
---|
123 | real(r8), intent(out) :: Dsrc(:)
|
---|
124 | real(r8), intent(out) :: Qbar(:)
|
---|
125 | # ifdef SOLVE3D
|
---|
126 | # if defined UV_PSOURCE || defined Q_PSOURCE
|
---|
127 | real(r8), intent(out) :: Qshape(:,:)
|
---|
128 | real(r8), intent(out) :: Qsrc(:,:)
|
---|
129 | # endif
|
---|
130 | # ifdef TS_PSOURCE
|
---|
131 | real(r8), intent(out) :: Tsrc(:,:,:)
|
---|
132 | # endif
|
---|
133 | # endif
|
---|
134 | #else
|
---|
135 | logical, intent(out) :: Lsrc(Msrc,NT(ng))
|
---|
136 |
|
---|
137 | integer, intent(out) :: Isrc(Msrc)
|
---|
138 | integer, intent(out) :: Jsrc(Msrc)
|
---|
139 |
|
---|
140 | real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
|
---|
141 | real(r8), intent(in) :: ubar(LBi:UBi,LBj:UBj,3)
|
---|
142 | real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,3)
|
---|
143 | # ifdef SOLVE3D
|
---|
144 | real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
|
---|
145 | real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
|
---|
146 | real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
|
---|
147 | # endif
|
---|
148 | real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
|
---|
149 | real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
|
---|
150 | real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
|
---|
151 |
|
---|
152 | real(r8), intent(out) :: Dsrc(Msrc)
|
---|
153 | real(r8), intent(out) :: Qbar(Msrc)
|
---|
154 | # ifdef SOLVE3D
|
---|
155 | # if defined UV_PSOURCE || defined Q_PSOURCE
|
---|
156 | real(r8), intent(out) :: Qshape(Msrc,N(ng))
|
---|
157 | real(r8), intent(out) :: Qsrc(Msrc,N(ng))
|
---|
158 | # endif
|
---|
159 | # ifdef TS_PSOURCE
|
---|
160 | real(r8), intent(out) :: Tsrc(Msrc,N(ng),NT(ng))
|
---|
161 | # endif
|
---|
162 | # endif
|
---|
163 | #endif
|
---|
164 | !
|
---|
165 | ! Local variable declarations.
|
---|
166 | !
|
---|
167 | integer :: is, i, j, k, ised
|
---|
168 | real(r8) :: fac, my_area, ramp
|
---|
169 | !
|
---|
170 | !-----------------------------------------------------------------------
|
---|
171 | ! Set tracer and/or mass point sources and/or sink.
|
---|
172 | !-----------------------------------------------------------------------
|
---|
173 | !
|
---|
174 | IF (iic(ng).eq.ntstart(ng)) THEN
|
---|
175 | !
|
---|
176 | ! Set-up point Sources/Sink number (Nsrc), direction (Dsrc), I- and
|
---|
177 | ! J-grid locations (Isrc,Jsrc), and logical switch for type of tracer
|
---|
178 | ! to apply (Lsrc). Currently, the direction can be along XI-direction
|
---|
179 | ! (Dsrc = 0) or along ETA-direction (Dsrc > 0). The mass sources are
|
---|
180 | ! located at U- or V-points so the grid locations should range from
|
---|
181 | ! 1 =< Isrc =< L and 1 =< Jsrc =< M.
|
---|
182 | !
|
---|
183 | Lsrc(:,:)=.FALSE.
|
---|
184 | #if defined RIVERPLUME1
|
---|
185 | Nsrc=1
|
---|
186 | Dsrc(Nsrc)=0.0_r8
|
---|
187 | Isrc(Nsrc)=1
|
---|
188 | Jsrc(Nsrc)=50
|
---|
189 | Lsrc(Nsrc,itemp)=.TRUE.
|
---|
190 | Lsrc(Nsrc,isalt)=.TRUE.
|
---|
191 | #elif defined RIVERPLUME2
|
---|
192 | Nsrc=1+Lm(ng)*2
|
---|
193 | DO is=1,(Nsrc-1)/2
|
---|
194 | Dsrc(is)=1.0_r8
|
---|
195 | Isrc(is)=is
|
---|
196 | Jsrc(is)=1
|
---|
197 | Lsrc(is,itemp)=.TRUE.
|
---|
198 | Lsrc(is,isalt)=.TRUE.
|
---|
199 | END DO
|
---|
200 | DO is=(Nsrc-1)/2+1,Nsrc-1
|
---|
201 | Dsrc(is)=1.0_r8
|
---|
202 | Isrc(is)=is-Lm(ng)
|
---|
203 | Jsrc(is)=Mm(ng)+1
|
---|
204 | Lsrc(is,itemp)=.TRUE.
|
---|
205 | Lsrc(is,isalt)=.TRUE.
|
---|
206 | END DO
|
---|
207 | Dsrc(Nsrc)=0.0_r8
|
---|
208 | Isrc(Nsrc)=1
|
---|
209 | Jsrc(Nsrc)=60
|
---|
210 | Lsrc(Nsrc,itemp)=.TRUE.
|
---|
211 | Lsrc(Nsrc,isalt)=.TRUE.
|
---|
212 | #elif defined SED_TEST1
|
---|
213 | Nsrc=Mm(ng)*2
|
---|
214 | DO is=1,Nsrc/2
|
---|
215 | Dsrc(is)=0.0_r8
|
---|
216 | Isrc(is)=1
|
---|
217 | Jsrc(is)=is
|
---|
218 | END DO
|
---|
219 | DO is=Nsrc/2+1,Nsrc
|
---|
220 | Dsrc(is)=0.0_r8
|
---|
221 | Isrc(is)=Lm(ng)+1
|
---|
222 | Jsrc(is)=is-Mm(ng)
|
---|
223 | END DO
|
---|
224 | #else
|
---|
225 | ana_psource.h: No values provided for Dsrc, Isrc, Jsrc, Lsrc.
|
---|
226 | #endif
|
---|
227 | END IF
|
---|
228 | #if defined UV_PSOURCE || defined Q_PSOURCE
|
---|
229 | # ifdef SOLVE3D
|
---|
230 | !
|
---|
231 | ! If appropriate, set-up nondimensional shape function to distribute
|
---|
232 | ! mass point sources/sinks vertically. It must add to unity!!.
|
---|
233 | !
|
---|
234 | # if defined SED_TEST1
|
---|
235 | DO k=1,N(ng)
|
---|
236 | DO is=1,Nsrc
|
---|
237 | i=Isrc(is)
|
---|
238 | j=Jsrc(is)
|
---|
239 | Qshape(is,k)=ABS(u(i,j,k,nnew)/ubar(i,j,knew))* &
|
---|
240 | & (z_w(i-1,Mm(ng)/2,k)-z_w(i-1,Mm(ng)/2,k-1)+ &
|
---|
241 | & z_w(i ,Mm(ng)/2,k)-z_w(i ,Mm(ng)/2,k-1))/ &
|
---|
242 | & (z_w(i-1,Mm(ng)/2,N(ng))-z_w(i-1,Mm(ng)/2,0)+ &
|
---|
243 | & z_w(i ,Mm(ng)/2,N(ng))-z_w(i ,Mm(ng)/2,0))
|
---|
244 | END DO
|
---|
245 | END DO
|
---|
246 | # elif defined RIVERPLUME2
|
---|
247 | DO k=1,N(ng)
|
---|
248 | DO is=1,Nsrc-1
|
---|
249 | i=Isrc(is)
|
---|
250 | j=Jsrc(is)
|
---|
251 | Qshape(is,k)=ABS(v(i,j,k,nnew)/vbar(i,j,knew))* &
|
---|
252 | & (z_w(i,j-1,k)-z_w(i,j-1,k-1)+ &
|
---|
253 | & z_w(i ,j,k)-z_w(i ,j,k-1))/ &
|
---|
254 | & (z_w(i,j-1,N(ng))-z_w(i,j-1,0)+ &
|
---|
255 | & z_w(i ,j,N(ng))-z_w(i ,j,0))
|
---|
256 | END DO
|
---|
257 | Qshape(Nsrc,k)=1.0_r8/REAL(N(ng),r8)
|
---|
258 | END DO
|
---|
259 | # else
|
---|
260 | DO k=1,N(ng)
|
---|
261 | DO is=1,Nsrc
|
---|
262 | Qshape(is,k)=1.0_r8/REAL(N(ng),r8)
|
---|
263 | END DO
|
---|
264 | END DO
|
---|
265 | # endif
|
---|
266 | # endif
|
---|
267 | !
|
---|
268 | ! Set-up vertically integrated mass transport (m3/s) of point
|
---|
269 | ! Sources/Sinks (positive in the positive U- or V-direction and
|
---|
270 | ! viceversa).
|
---|
271 | !
|
---|
272 | # if defined RIVERPLUME1
|
---|
273 | IF ((tdays(ng)-dstart).lt.0.5_r8) THEN
|
---|
274 | fac=1.0_r8+TANH((time(ng)-43200.0_r8)/43200.0_r8)
|
---|
275 | ELSE
|
---|
276 | fac=1.0_r8
|
---|
277 | END IF
|
---|
278 | DO is=1,Nsrc
|
---|
279 | Qbar(is)=fac*1500.0_r8
|
---|
280 | END DO
|
---|
281 | # elif defined RIVERPLUME2
|
---|
282 | DO is=1,(Nsrc-1)/2 ! North end
|
---|
283 | i=Isrc(is)
|
---|
284 | j=Jsrc(is)
|
---|
285 | Qbar(is)=-0.05_r8*om_v(i,j)*(0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
|
---|
286 | & zeta(i ,j,knew)+h(i ,j)))
|
---|
287 | END DO
|
---|
288 | ! South End
|
---|
289 | DO is=(Nsrc-1)/2+1,Nsrc-1 ! South end
|
---|
290 | i=Isrc(is)
|
---|
291 | j=Jsrc(is)
|
---|
292 | Qbar(is)=-0.05_r8*om_v(i,j)*(0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
|
---|
293 | & zeta(i ,j,knew)+h(i ,j)))
|
---|
294 | END DO
|
---|
295 | Qbar(Nsrc)=1500.0_r8 ! West wall
|
---|
296 | # elif defined SED_TEST1
|
---|
297 | my_area=0.0_r8 ! West end
|
---|
298 | DO is=1,Nsrc/2
|
---|
299 | i=Isrc(is)
|
---|
300 | j=Jsrc(is)
|
---|
301 | my_area=my_area+0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
|
---|
302 | & zeta(i ,j,knew)+h(i ,j))*on_u(i,j)
|
---|
303 | END DO
|
---|
304 | fac=-36.0_r8*10.0_r8*1.0_r8
|
---|
305 | DO is=1,Nsrc/2
|
---|
306 | i=Isrc(is)
|
---|
307 | j=Jsrc(is)
|
---|
308 | Qbar(is)=fac*(0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
|
---|
309 | & zeta(i ,j,knew)+h(i ,j)))* &
|
---|
310 | & on_u(i,j)/my_area
|
---|
311 | END DO
|
---|
312 | my_area=0.0_r8 ! East end
|
---|
313 | DO is=Nsrc/2+1,Nsrc
|
---|
314 | i=Isrc(is)
|
---|
315 | j=Jsrc(is)
|
---|
316 | my_area=my_area+0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
|
---|
317 | & zeta(i ,j,knew)+h(i ,j))*on_u(i,j)
|
---|
318 | END DO
|
---|
319 | fac=-36.0_r8*10.0_r8*1.0_r8
|
---|
320 | DO is=Nsrc/2+1,Nsrc
|
---|
321 | i=Isrc(is)
|
---|
322 | j=Jsrc(is)
|
---|
323 | Qbar(is)=fac*(0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
|
---|
324 | & zeta(i ,j,knew)+h(i ,j)))* &
|
---|
325 | & on_u(i,j)/my_area
|
---|
326 | END DO
|
---|
327 | # else
|
---|
328 | ana_psource.h: No values provided for Qbar.
|
---|
329 | # endif
|
---|
330 | # ifdef SOLVE3D
|
---|
331 | !
|
---|
332 | ! Set-up mass transport profile (m3/s) of point Sources/Sinks.
|
---|
333 | !
|
---|
334 | DO k=1,N(ng)
|
---|
335 | DO is=1,Nsrc
|
---|
336 | Qsrc(is,k)=Qbar(is)*Qshape(is,k)
|
---|
337 | END DO
|
---|
338 | END DO
|
---|
339 | # endif
|
---|
340 | #endif
|
---|
341 | #if defined TS_PSOURCE && defined SOLVE3D
|
---|
342 | !
|
---|
343 | ! Set-up tracer (tracer units) point Sources/Sinks.
|
---|
344 | !
|
---|
345 | # if defined RIVERPLUME1
|
---|
346 | DO k=1,N(ng)
|
---|
347 | DO is=1,Nsrc
|
---|
348 | Tsrc(is,k,itemp)=T0(ng)
|
---|
349 | Tsrc(is,k,isalt)=0.0_r8
|
---|
350 | # ifdef SEDIMENT
|
---|
351 | DO ised=1,NST
|
---|
352 | Tsrc(is,k,ised+2)=0.0_r8
|
---|
353 | END DO
|
---|
354 | # endif
|
---|
355 | END DO
|
---|
356 | END DO
|
---|
357 | # elif defined RIVERPLUME2
|
---|
358 | DO k=1,N(ng)
|
---|
359 | DO is=1,Nsrc-1
|
---|
360 | Tsrc(is,k,itemp)=T0(ng)
|
---|
361 | Tsrc(is,k,isalt)=S0(ng)
|
---|
362 | END DO
|
---|
363 | Tsrc(Nsrc,k,itemp)=T0(ng)
|
---|
364 | Tsrc(Nsrc,k,isalt)=0.0_r8
|
---|
365 | END DO
|
---|
366 | # else
|
---|
367 | ana_psource.h: No values provided for Tsrc.
|
---|
368 | # endif
|
---|
369 | #endif
|
---|
370 | RETURN
|
---|
371 | END SUBROUTINE ana_psource_grid
|
---|