Ticket #31: ana_psource.h

File ana_psource.h, 12.2 KB (added by m.hadfield, 19 years ago)
Line 
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