ROMS
Loading...
Searching...
No Matches
ana_psource.h
Go to the documentation of this file.
1!!
2 SUBROUTINE ana_psource (ng, tile, model)
3!
4!! git $Id$
5!!======================================================================
6!! Copyright (c) 2002-2025 The ROMS Group !
7!! Licensed under a MIT/X style license !
8!! See License_ROMS.md !
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_stepping
21!
22! Imported variable declarations
23!
24 integer, intent(in) :: ng, tile, model
25!
26! Local variable declarations.
27!
28 character (len=*), parameter :: MyFile = &
29 & __FILE__
30!
31#include "tile.h"
32!
33 CALL ana_psource_tile (ng, tile, model, &
34 & lbi, ubi, lbj, ubj, &
35 & imins, imaxs, jmins, jmaxs, &
36 & nnew(ng), knew(ng), &
37 & ocean(ng) % zeta, &
38 & ocean(ng) % ubar, &
39 & ocean(ng) % vbar, &
40#ifdef SOLVE3D
41 & ocean(ng) % u, &
42 & ocean(ng) % v, &
43 & grid(ng) % z_w, &
44#endif
45 & grid(ng) % h, &
46 & grid(ng) % on_u, &
47 & grid(ng) % om_v)
48!
49! Set analytical header file name used.
50!
51#ifdef DISTRIBUTE
52 IF (lanafile) THEN
53#else
54 IF (lanafile.and.(tile.eq.0)) THEN
55#endif
56 ananame(20)=myfile
57 END IF
58!
59 RETURN
60 END SUBROUTINE ana_psource
61!
62!***********************************************************************
63 SUBROUTINE ana_psource_tile (ng, tile, model, &
64 & LBi, UBi, LBj, UBj, &
65 & IminS, ImaxS, JminS, JmaxS, &
66 & nnew, knew, &
67 & zeta, ubar, vbar, &
68#ifdef SOLVE3D
69 & u, v, z_w, &
70#endif
71 & h, on_u, om_v)
72!***********************************************************************
73!
74 USE mod_param
75 USE mod_parallel
76 USE mod_scalars
77#ifdef SEDIMENT
78 USE mod_sediment
79#endif
80 USE mod_sources
81
82#ifdef DISTRIBUTE
83!
86#endif
87!
88! Imported variable declarations.
89!
90 integer, intent(in) :: ng, tile, model
91 integer, intent(in) :: LBi, UBi, LBj, UBj
92 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
93 integer, intent(in) :: nnew, knew
94!
95#ifdef ASSUMED_SHAPE
96 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
97 real(r8), intent(in) :: ubar(LBi:,LBj:,:)
98 real(r8), intent(in) :: vbar(LBi:,LBj:,:)
99# ifdef SOLVE3D
100 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
101 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
102 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
103# endif
104 real(r8), intent(in) :: h(LBi:,LBj:)
105 real(r8), intent(in) :: on_u(LBi:,LBj:)
106 real(r8), intent(in) :: om_v(LBi:,LBj:)
107#else
108 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
109 real(r8), intent(in) :: ubar(LBi:UBi,LBj:UBj,:)
110 real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,:)
111# ifdef SOLVE3D
112 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
113 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
114 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
115# endif
116 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
117 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
118 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
119#endif
120!
121! Local variable declarations.
122!
123 integer :: Npts, NSUB, is, i, j, k, ised
124!
125 real(r8) :: Pspv = 0.0_r8
126 real(r8), save :: area_east, area_west
127 real(r8) :: cff, fac, my_area_east, my_area_west
128
129#if defined DISTRIBUTE && defined SOLVE3D
130 real(r8), dimension(Msrc(ng)*N(ng)) :: Pwrk
131#endif
132#if defined DISTRIBUTE
133 real(r8), dimension(2) :: rbuffer
134!
135 character (len=3), dimension(2) :: io_handle
136#endif
137
138#include "set_bounds.h"
139!
140!-----------------------------------------------------------------------
141! If initialization, set point Sources and/or Sinks locations.
142!-----------------------------------------------------------------------
143!
144 IF ((iic(ng).eq.ntstart(ng)).or.(iic(ng).eq.0)) THEN
145!
146! Set-up point Sources/Sink number (Nsrc), direction (Dsrc), I- and
147! J-grid locations (Isrc,Jsrc). Currently, the direction can be along
148! XI-direction (Dsrc=0) or along ETA-direction (Dsrc=1). The
149! mass sources are located at U- or V-points so the grid locations
150! should range from 1 =< Isrc =< L and 1 =< Jsrc =< M.
151!
152! Vertical mass sources can be added my setting a W-direction (Dsrc=2)
153! and mass sources are located at Rho-points so the grid locations
154! should range from 0 =< Isrc =< L and 0 =< Jsrc =< M.
155!
156#if defined RIVERPLUME1
157 IF (master.and.domain(ng)%SouthWest_Test(tile)) THEN
158 nsrc(ng)=1
159 sources(ng)%Dsrc(nsrc(ng))=0.0_r8 ! horizontal, LuvSrc=T
160 sources(ng)%Isrc(nsrc(ng))=2 ! i = 2
161 sources(ng)%Jsrc(nsrc(ng))=50 ! j = 50
162!
163 nsrc(ng)=nsrc(ng)+10 ! Add rainfall location
164 DO is=2,6
165 sources(ng)%Dsrc(is)=2.0_r8 ! vertical influx, LwSrc=T
166 sources(ng)%Isrc(is)=6 ! i = 6
167 sources(ng)%Jsrc(is)=is+25 ! j = 27 to 31
168 END DO
169 DO is=7,11
170 sources(ng)%Dsrc(is)=2.0_r8 ! vertical inclux, LwSrc=T
171 sources(ng)%Isrc(is)=7 ! i = 7
172 sources(ng)%Jsrc(is)=is+20 ! j = 27 to 31
173 END DO
174 END IF
175#elif defined RIVERPLUME2
176 IF (master.and.domain(ng)%SouthWest_Test(tile)) THEN
177 nsrc(ng)=1+lm(ng)*2
178 DO is=1,(nsrc(ng)-1)/2
179 sources(ng)%Dsrc(is)=1.0_r8
180 sources(ng)%Isrc(is)=is
181 sources(ng)%Jsrc(is)=1
182 END DO
183 DO is=(nsrc(ng)-1)/2+1,nsrc(ng)-1
184 sources(ng)%Dsrc(is)=1.0_r8
185 sources(ng)%Isrc(is)=is-lm(ng)
186 sources(ng)%Jsrc(is)=mm(ng)+1
187 END DO
188 sources(ng)%Dsrc(nsrc(ng))=0.0_r8
189 sources(ng)%Isrc(nsrc(ng))=1
190 sources(ng)%Jsrc(nsrc(ng))=60
191 END IF
192#elif defined SED_TEST1
193 IF (master.and.domain(ng)%SouthWest_Test(tile)) THEN
194 nsrc(ng)=mm(ng)*2
195 DO is=1,nsrc(ng)/2
196 sources(ng)%Dsrc(is)=0.0_r8
197 sources(ng)%Isrc(is)=1
198 sources(ng)%Jsrc(is)=is
199 END DO
200 DO is=nsrc(ng)/2+1,nsrc(ng)
201 sources(ng)%Dsrc(is)=0.0_r8
202 sources(ng)%Isrc(is)=lm(ng)+1
203 sources(ng)%Jsrc(is)=is-mm(ng)
204 END DO
205 END IF
206#else
207 ana_psource.h: no values provided for nsrc, dsrc, isrc, jsrc.
208#endif
209
210#ifdef DISTRIBUTE
211!
212! Broadcast point sources/sinks information to all nodes.
213!
214 CALL mp_bcasti (ng, inlm, nsrc(ng))
215 CALL mp_bcasti (ng, inlm, sources(ng)%Isrc)
216 CALL mp_bcasti (ng, inlm, sources(ng)%Jsrc)
217 CALL mp_bcastf (ng, inlm, sources(ng)%Dsrc)
218#endif
219 END IF
220!
221!-----------------------------------------------------------------------
222! Set momentum point Sources and/or Sinks.
223!-----------------------------------------------------------------------
224!
225 momentum : IF (luvsrc(ng).or.lwsrc(ng)) THEN
226
227#ifdef SOLVE3D
228!
229! If appropriate, set-up nondimensional shape function to distribute
230! mass point sources/sinks vertically. It must add to unity!!.
231!
232# ifdef DISTRIBUTE
233 sources(ng)%Qshape=pspv
234# endif
235 npts=msrc(ng)*n(ng)
236
237!$OMP BARRIER
238
239# if defined SED_TEST1
240 DO k=1,n(ng)
241 DO is=1,nsrc(ng)
242 i=sources(ng)%Isrc(is)
243 j=sources(ng)%Jsrc(is)
244 IF (((istrt.le.i).and.(i.le.iendt)).and. &
245 & ((jstrt.le.j).and.(j.le.jendt))) THEN
246 IF (ubar(i,j,knew).ne.0.0_r8) THEN
247 cff=abs(u(i,j,k,nnew)/ubar(i,j,knew))
248 ELSE
249 cff=1.0_r8
250 END IF
251 sources(ng)%Qshape(is,k)=cff* &
252 & (z_w(i-1,j,k )- &
253 & z_w(i-1,j,k-1 )+ &
254 & z_w(i ,j,k )- &
255 & z_w(i ,j,k-1 ))/ &
256 & (z_w(i-1,j,n(ng))- &
257 & z_w(i-1,j,0 )+ &
258 & z_w(i ,j,n(ng))- &
259 & z_w(i ,j,0 ))
260 END IF
261 END DO
262 END DO
263# ifdef DISTRIBUTE
264 pwrk=reshape(sources(ng)%Qshape,(/npts/))
265 CALL mp_collect (ng, inlm, npts, pspv, pwrk)
266 sources(ng)%Qshape=reshape(pwrk,(/msrc(ng),n(ng)/))
267# endif
268
269# elif defined RIVERPLUME1
270
271 IF (domain(ng)%NorthEast_Test(tile)) THEN
272 DO k=1,n(ng)
273 DO is=1,nsrc(ng)
274 sources(ng)%Qshape(is,k)=1.0_r8/real(n(ng),r8)
275 END DO
276 END DO
277 END IF
278
279# elif defined RIVERPLUME2
280 DO k=1,n(ng)
281 DO is=1,nsrc(ng)-1
282 i=sources(ng)%Isrc(is)
283 j=sources(ng)%Jsrc(is)
284 IF (((istrt.le.i).and.(i.le.iendt)).and. &
285 & ((jstrt.le.j).and.(j.le.jendt))) THEN
286 IF (vbar(i,j,knew).ne.0.0_r8) THEN
287 cff=abs(v(i,j,k,nnew)/vbar(i,j,knew))
288 ELSE
289 cff=1.0_r8
290 END IF
291 sources(ng)%Qshape(is,k)=cff* &
292 & (z_w(i,j-1,k )- &
293 & z_w(i,j-1,k-1 )+ &
294 & z_w(i,j ,k )- &
295 & z_w(i,j ,k-1 ))/ &
296 & (z_w(i,j-1,n(ng))- &
297 & z_w(i,j-1,0 )+ &
298 & z_w(i,j ,n(ng))- &
299 & z_w(i,j ,0 ))
300 END IF
301 END DO
302 END DO
303 IF (master.and.domain(ng)%SouthWest_Test(tile)) THEN
304 DO k=1,n(ng)
305 sources(ng)%Qshape(nsrc(ng),k)=1.0_r8/real(n(ng),r8)
306 END DO
307 END IF
308# ifdef DISTRIBUTE
309 pwrk=reshape(sources(ng)%Qshape,(/npts/))
310 CALL mp_collect (ng, inlm, npts, pspv, pwrk)
311 sources(ng)%Qshape=reshape(pwrk,(/msrc(ng),n(ng)/))
312# endif
313
314# else
315!!
316!! Notice that there is not need for distributed-memory communications
317!! here since the computation below does not have a parallel tile
318!! dependency. All the nodes are computing this simple statement.
319!!
320 IF (domain(ng)%NorthEast_Test(tile)) THEN
321 DO k=1,n(ng)
322 DO is=1,nsrc(ng)
323 sources(ng)%Qshape(is,k)=1.0_r8/real(n(ng),r8)
324 END DO
325 END DO
326 END IF
327# endif
328#endif
329!
330! Set-up vertically integrated mass transport (m3/s) of point
331! Sources/Sinks (positive in the positive U- or V-direction and
332! viceversa).
333!
334#ifdef DISTRIBUTE
335 sources(ng)%Qbar=pspv
336#endif
337
338!$OMP BARRIER
339
340#if defined RIVERPLUME1
341 IF ((tdays(ng)-dstart).lt.0.5_r8) THEN
342 fac=1.0_r8+tanh((time(ng)-43200.0_r8)/43200.0_r8)
343 ELSE
344 fac=1.0_r8
345 END IF
346 DO is=1,1 ! horizontal influx only
347 sources(ng)%Qbar(is)=fac*1500.0_r8
348 END DO
349!
350! Rainfall is 10 cm/hour.
351!
352! Qbar = 3000m * 1500m * 0.10 m/hr * 1hr/3600s = 125 m3/s
353!
354 DO is=2,nsrc(ng)
355 sources(ng)%Qbar(is)=125.0_r8
356 END DO
357
358#elif defined RIVERPLUME2
359 DO is=1,(nsrc(ng)-1)/2 ! North end
360 i=sources(ng)%Isrc(is)
361 j=sources(ng)%Jsrc(is)
362 IF (((istrt.le.i).and.(i.le.iendt)).and. &
363 & ((jstrt.le.j).and.(j.le.jendt))) THEN
364 sources(ng)%Qbar(is)=-0.05_r8*om_v(i,j)* &
365 & (0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
366 & zeta(i ,j,knew)+h(i ,j)))
367 END IF
368 END DO
369 DO is=(nsrc(ng)-1)/2+1,nsrc(ng)-1 ! South end
370 i=sources(ng)%Isrc(is)
371 j=sources(ng)%Jsrc(is)
372 IF (((istrt.le.i).and.(i.le.iendt)).and. &
373 & ((jstrt.le.j).and.(j.le.jendt))) THEN
374 sources(ng)%Qbar(is)=-0.05_r8*om_v(i,j)* &
375 & (0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
376 & zeta(i ,j,knew)+h(i ,j)))
377 END IF
378 END DO
379 IF (master.and.domain(ng)%SouthWest_Test(tile)) THEN
380 sources(ng)%Qbar(nsrc(ng))=1500.0_r8 ! West wall
381 END IF
382# ifdef DISTRIBUTE
383 CALL mp_collect (ng, inlm, msrc(ng), pspv, sources(ng)%Qbar)
384# endif
385
386#elif defined SED_TEST1
387 my_area_west=0.0_r8 ! West end
388 fac=-36.0_r8*10.0_r8*1.0_r8
389 DO is=1,nsrc(ng)/2
390 i=sources(ng)%Isrc(is)
391 j=sources(ng)%Jsrc(is)
392 IF (((istrt.le.i).and.(i.le.iendt)).and. &
393 & ((jstrt.le.j).and.(j.le.jendt))) THEN
394 cff=0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
395 & zeta(i ,j,knew)+h(i ,j))*on_u(i,j)
396 sources(ng)%Qbar(is)=fac*cff
397 my_area_west=my_area_west+cff
398 END IF
399 END DO
400!
401 my_area_east=0.0_r8 ! East end
402 fac=-36.0_r8*10.0_r8*1.0_r8
403 DO is=nsrc(ng)/2+1,nsrc(ng)
404 i=sources(ng)%Isrc(is)
405 j=sources(ng)%Jsrc(is)
406 IF (((istrt.le.i).and.(i.le.iendt)).and. &
407 & ((jstrt.le.j).and.(j.le.jendt))) THEN
408 cff=0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
409 & zeta(i ,j,knew)+h(i ,j))*on_u(i,j)
410 sources(ng)%Qbar(is)=fac*cff
411 my_area_east=my_area_east+cff
412 END IF
413 END DO
414!
415# ifdef DISTRIBUTE
416 nsub=1 ! distributed-memory
417# else
418 IF (domain(ng)%SouthWest_Corner(tile).and. &
419 & domain(ng)%NorthEast_Corner(tile)) THEN
420 nsub=1 ! non-tiled application
421 ELSE
422 nsub=ntilex(ng)*ntilee(ng) ! tiled application
423 END IF
424# endif
425!$OMP CRITICAL (PSOURCE)
426 IF (tile_count.eq.0) THEN
427 area_west=0.0_r8
428 area_east=0.0_r8
429 END IF
430 area_west=area_west+my_area_west
431 area_east=area_east+my_area_east
433 IF (tile_count.eq.nsub) THEN
434 tile_count=0
435# ifdef DISTRIBUTE
436 rbuffer(1)=area_west
437 rbuffer(2)=area_east
438 io_handle(1)='SUM'
439 io_handle(2)='SUM'
440 CALL mp_reduce (ng, inlm, 2, rbuffer, io_handle)
441 area_west=rbuffer(1)
442 area_east=rbuffer(1)
443# endif
444 DO is=1,nsrc(ng)/2
445 sources(ng)%Qbar(is)=sources(ng)%Qbar(is)/area_west
446 END DO
447 DO is=nsrc(ng)/2+1,nsrc(ng)
448 sources(ng)%Qbar(is)=sources(ng)%Qbar(is)/area_east
449 END DO
450 END IF
451!$OMP END CRITICAL (PSOURCE)
452
453# ifdef DISTRIBUTE
454 CALL mp_collect (ng, inlm, msrc(ng), pspv, sources(ng)%Qbar)
455# endif
456#else
457 ana_psource.h: no values provided for qbar.
458#endif
459
460#ifdef SOLVE3D
461!
462! Set-up mass transport profile (m3/s) of point Sources/Sinks.
463!
464!$OMP BARRIER
465
466 IF (domain(ng)%NorthEast_Test(tile)) THEN
467 DO k=1,n(ng)
468 DO is=1,nsrc(ng)
469 sources(ng)%Qsrc(is,k)=sources(ng)%Qbar(is)* &
470 & sources(ng)%Qshape(is,k)
471 END DO
472 END DO
473 END IF
474#endif
475 END IF momentum
476
477#ifdef SOLVE3D
478!
479!-----------------------------------------------------------------------
480! Set tracers point Sources and/or Sinks.
481!-----------------------------------------------------------------------
482!
483 tracers : IF (any(ltracersrc(:,ng))) THEN
484 sources(ng)%Tsrc=0.0_r8 ! initialize
485!
486! Set-up tracer (tracer units) point Sources/Sinks.
487!
488# if defined RIVERPLUME1
489 IF (domain(ng)%NorthEast_Test(tile)) THEN
490 DO k=1,n(ng)
491 DO is=1,nsrc(ng)
492 IF (is.eq.1) THEN
493 sources(ng)%Tsrc(is,k,itemp)=10.0_r8
494 ELSE
495 sources(ng)%Tsrc(is,k,itemp)=t0(ng)
496 END IF
497# ifdef SALINITY
498 IF (is.eq.1) THEN
499 sources(ng)%Tsrc(is,k,isalt)=0.0_r8
500 ELSE
501 sources(ng)%Tsrc(is,k,isalt)=s0(ng)
502 END IF
503# endif
504 END DO
505 END DO
506 END IF
507
508# elif defined RIVERPLUME2
509 IF (domain(ng)%NorthEast_Test(tile)) THEN
510 DO k=1,n(ng)
511 DO is=1,nsrc(ng)-1
512 sources(ng)%Tsrc(is,k,itemp)=t0(ng)
513# ifdef SALINITY
514 sources(ng)%Tsrc(is,k,isalt)=s0(ng)
515# endif
516 END DO
517 sources(ng)%Tsrc(nsrc(ng),k,itemp)=t0(ng)
518# ifdef SALINITY
519 sources(ng)%Tsrc(nsrc(ng),k,isalt)=s0(ng)
520# endif
521 END DO
522 END IF
523
524# elif defined SED_TEST1
525!
526! No tracers point sources.
527!
528# else
529 ana_psource.h: no values provided for tsrc.
530# endif
531 END IF tracers
532#endif
533!
534 RETURN
535 END SUBROUTINE ana_psource_tile
subroutine ana_psource_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nnew, knew, zeta, ubar, vbar, u, v, z_w, h, on_u, om_v)
Definition ana_psource.h:72
subroutine ana_psource(ng, tile, model)
Definition ana_psource.h:3
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
logical lanafile
character(len=256), dimension(39) ananame
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
logical master
integer tile_count
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer, dimension(:), allocatable lm
Definition mod_param.F:455
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
integer, dimension(:), allocatable mm
Definition mod_param.F:456
logical, dimension(:), allocatable luvsrc
logical, dimension(:,:), allocatable ltracersrc
integer, dimension(:), allocatable iic
real(r8), dimension(:), allocatable t0
real(r8), dimension(:), allocatable s0
real(dp), dimension(:), allocatable tdays
real(dp) dstart
logical, dimension(:), allocatable lwsrc
integer isalt
integer itemp
real(dp), dimension(:), allocatable time
integer, dimension(:), allocatable ntstart
type(t_sources), dimension(:), allocatable sources
Definition mod_sources.F:90
integer, dimension(:), allocatable nsrc
Definition mod_sources.F:97
integer, dimension(:), allocatable msrc
Definition mod_sources.F:96
integer, dimension(:), allocatable knew
integer, dimension(:), allocatable nnew