ROMS
Loading...
Searching...
No Matches
ana_initial.h
Go to the documentation of this file.
1!!
2 SUBROUTINE ana_initial (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 initial conditions for momentum and tracer !
12! type variables using analytical expressions. !
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 IF (model.eq.inlm) THEN
34 CALL ana_nlminitial_tile (ng, tile, model, &
35 & lbi, ubi, lbj, ubj, &
36 & imins, imaxs, jmins, jmaxs, &
37 & grid(ng) % h, &
38#ifdef SPHERICAL
39 & grid(ng) % lonr, &
40 & grid(ng) % latr, &
41#else
42 & grid(ng) % xr, &
43 & grid(ng) % yr, &
44#endif
45#ifdef SOLVE3D
46 & grid(ng) % z_r, &
47 & ocean(ng) % u, &
48 & ocean(ng) % v, &
49 & ocean(ng) % t, &
50#endif
51 & ocean(ng) % ubar, &
52 & ocean(ng) % vbar, &
53 & ocean(ng) % zeta)
54#ifdef TANGENT
55 ELSE IF ((model.eq.itlm).or.(model.eq.irpm)) THEN
56 CALL ana_tlminitial_tile (ng, tile, model, &
57 & lbi, ubi, lbj, ubj, &
58 & imins, imaxs, jmins, jmaxs, &
59 & kstp(ng), &
60# ifdef SOLVE3D
61 & nstp(ng), &
62 & ocean(ng) % tl_u, &
63 & ocean(ng) % tl_v, &
64 & ocean(ng) % tl_t, &
65# endif
66 & ocean(ng) % tl_ubar, &
67 & ocean(ng) % tl_vbar, &
68 & ocean(ng) % tl_zeta)
69#endif
70#ifdef ADJOINT
71 ELSE IF (model.eq.iadm) THEN
72 CALL ana_adminitial_tile (ng, tile, model, &
73 & lbi, ubi, lbj, ubj, &
74 & imins, imaxs, jmins, jmaxs, &
75 & knew(ng), &
76# ifdef SOLVE3D
77 & nstp(ng), &
78 & ocean(ng) % ad_u, &
79 & ocean(ng) % ad_v, &
80 & ocean(ng) % ad_t, &
81# endif
82 & ocean(ng) % ad_ubar, &
83 & ocean(ng) % ad_vbar, &
84 & ocean(ng) % ad_zeta)
85#endif
86 END IF
87!
88! Set analytical header file name used.
89!
90#ifdef DISTRIBUTE
91 IF (lanafile) THEN
92#else
93 IF (lanafile.and.(tile.eq.0)) THEN
94#endif
95 ananame(10)=myfile
96 END IF
97!
98 RETURN
99 END SUBROUTINE ana_initial
100!
101!***********************************************************************
102 SUBROUTINE ana_nlminitial_tile (ng, tile, model, &
103 & LBi, UBi, LBj, UBj, &
104 & IminS, ImaxS, JminS, JmaxS, &
105 & h, &
106#ifdef SPHERICAL
107 & lonr, latr, &
108#else
109 & xr, yr, &
110#endif
111#ifdef SOLVE3D
112 & z_r, &
113 & u, v, t, &
114#endif
115 & ubar, vbar, zeta)
116!***********************************************************************
117!
118 USE mod_param
119 USE mod_parallel
120 USE mod_grid
121 USE mod_ncparam
122 USE mod_iounits
123 USE mod_scalars
124!
125#ifdef CHANNEL
126# ifdef DISTRIBUTE
127 USE distribute_mod, ONLY : mp_bcasti
128# endif
129 USE erf_mod, ONLY : erf
130#endif
131 USE stats_mod, ONLY : stats_2dfld
132#ifdef SOLVE3D
133 USE stats_mod, ONLY : stats_3dfld
134#endif
135!
136! Imported variable declarations.
137!
138 integer, intent(in) :: ng, tile, model
139 integer, intent(in) :: LBi, UBi, LBj, UBj
140 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
141!
142#ifdef ASSUMED_SHAPE
143 real(r8), intent(in) :: h(LBi:,LBj:)
144# ifdef SPHERICAL
145 real(r8), intent(in) :: lonr(LBi:,LBj:)
146 real(r8), intent(in) :: latr(LBi:,LBj:)
147# else
148 real(r8), intent(in) :: xr(LBi:,LBj:)
149 real(r8), intent(in) :: yr(LBi:,LBj:)
150# endif
151# ifdef SOLVE3D
152 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
153
154 real(r8), intent(out) :: u(LBi:,LBj:,:,:)
155 real(r8), intent(out) :: v(LBi:,LBj:,:,:)
156 real(r8), intent(out) :: t(LBi:,LBj:,:,:,:)
157# endif
158 real(r8), intent(out) :: ubar(LBi:,LBj:,:)
159 real(r8), intent(out) :: vbar(LBi:,LBj:,:)
160 real(r8), intent(out) :: zeta(LBi:,LBj:,:)
161#else
162# ifdef SPHERICAL
163 real(r8), intent(in) :: lonr(LBi:UBi,LBj:UBj)
164 real(r8), intent(in) :: latr(LBi:UBi,LBj:UBj)
165# else
166 real(r8), intent(in) :: xr(LBi:UBi,LBj:UBj)
167 real(r8), intent(in) :: yr(LBi:UBi,LBj:UBj)
168# endif
169 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
170# ifdef SOLVE3D
171 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
172
173 real(r8), intent(out) :: u(LBi:UBi,LBj:UBj,N(ng),2)
174 real(r8), intent(out) :: v(LBi:UBi,LBj:UBj,N(ng),2)
175 real(r8), intent(out) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
176# endif
177 real(r8), intent(out) :: ubar(LBi:UBi,LBj:UBj,:)
178 real(r8), intent(out) :: vbar(LBi:UBi,LBj:UBj,:)
179 real(r8), intent(out) :: zeta(LBi:UBi,LBj:UBj,:)
180#endif
181!
182! Local variable declarations.
183!
184 logical, save :: first = .true.
185!
186 integer :: Iless, Iplus, i, itrc, j, k
187!
188#ifdef CHANNEL
189 real(r8), parameter :: guscale = 40.0e+03_r8
190 real(r8), parameter :: u0 = 1.6_r8
191#endif
192 real(r8) :: depth, dx, val1, val2, val3, val4, x, x0, y, y0
193!
194 TYPE (T_STATS), save :: Stats(7) ! ubar, vbar, zeta, u, v, t, s
195
196#include "set_bounds.h"
197!
198!-----------------------------------------------------------------------
199! Initialize field statistics structure.
200!-----------------------------------------------------------------------
201!
202 IF (first) THEN
203 first=.false.
204 DO i=1,SIZE(stats,1)
205 stats(i) % checksum=0_i8b
206 stats(i) % count=0
207 stats(i) % min=large
208 stats(i) % max=-large
209 stats(i) % avg=0.0_r8
210 stats(i) % rms=0.0_r8
211 END DO
212 END IF
213!
214!-----------------------------------------------------------------------
215! Initial conditions for 2D momentum (m/s) components.
216!-----------------------------------------------------------------------
217!
218#if defined CHANNEL && !defined ONLY_TS_IC
219 y0=0.5_r8*el(ng)
220 DO j=jstrt,jendt
221 DO i=istrp,iendt
222 val1=(grid(ng)%yu(i,j)-y0)/guscale
223 ubar(i,j,1)=u0*exp(-val1*val1)/6.0_r8
224 END DO
225 END DO
226 DO j=jstrp,jendt
227 DO i=istrt,iendt
228 vbar(i,j,1)=0.0_r8
229 END DO
230 END DO
231#elif defined SOLITON
232 x0=2.0_r8*xl(ng)/3.0_r8
233 y0=0.5_r8*el(ng)
234 val1=0.395_r8
235 val2=0.771_r8*(val1*val1)
236 IF (ng.eq.1) THEN
237 DO j=jstrt,jendt
238 DO i=istrp,iendt
239 x=0.5_r8*(xr(i-1,j)+xr(i,j))-x0
240 y=0.5_r8*(yr(i-1,j)+yr(i,j))-y0
241 val3=exp(-val1*x)
242 val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2)
243 ubar(i,j,1)=0.25_r8*val4*(6.0_r8*y*y-9.0_r8)* &
244 & exp(-0.5_r8*y*y)
245 END DO
246 END DO
247 DO j=jstrp,jendt
248 DO i=istrt,iendt
249 x=0.5_r8*(xr(i,j-1)+xr(i,j))-x0
250 y=0.5_r8*(yr(i,j-1)+yr(i,j))-y0
251 val3=exp(-val1*x)
252 val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2)
253 vbar(i,j,1)=2.0_r8*val4*y*(-2.0_r8*val1*tanh(val1*x))* &
254 & exp(-0.5_r8*y*y)
255 END DO
256 END DO
257 ELSE
258 DO j=jstrt,jendt
259 DO i=istrp,iendt
260 ubar(i,j,1)=0.0_r8
261 END DO
262 END DO
263 DO j=jstrp,jendt
264 DO i=istrt,iendt
265 vbar(i,j,1)=0.0_r8
266 END DO
267 END DO
268 END IF
269#elif defined RIVERPLUME2
270 DO j=jstrt,jendt
271 DO i=istrp,iendt
272 ubar(i,j,1)=0.0_r8
273 END DO
274 END DO
275 DO j=jstrp,jendt
276 DO i=istrt,iendt
277 vbar(i,j,1)=-0.05_r8
278 END DO
279 END DO
280#elif defined SED_TEST1
281 val1=100.0_r8
282 DO j=jstrt,jendt
283 DO i=istrp,iendt
284 ubar(i,j,1)=-10.0_r8/(10.0_r8+9.0e-06_r8*real(i,r8)*val1)
285 END DO
286 END DO
287 DO j=jstrp,jendt
288 DO i=istrt,iendt
289 vbar(i,j,1)=0.0_r8
290 END DO
291 END DO
292#elif defined SED_TOY
293 DO j=jstrt,jendt
294 DO i=istrp,iendt
295 ubar(i,j,1)=0.0_r8
296 END DO
297 END DO
298 DO j=jstrp,jendt
299 DO i=istrt,iendt
300 vbar(i,j,1)=0.0_r8
301 END DO
302 END DO
303#elif defined TEST_CHAN
304 val1=100.0_r8
305 DO j=jstrt,jendt
306 DO i=istrp,iendt
307 ubar(i,j,1)=0.0_r8
308 END DO
309 END DO
310 DO j=jstrp,jendt
311 DO i=istrt,iendt
312 vbar(i,j,1)=0.0_r8
313 END DO
314 END DO
315#else
316 DO j=jstrt,jendt
317 DO i=istrp,iendt
318 ubar(i,j,1)=0.0_r8
319 END DO
320 END DO
321 DO j=jstrp,jendt
322 DO i=istrt,iendt
323 vbar(i,j,1)=0.0_r8
324 END DO
325 END DO
326#endif
327!
328! Report statistics.
329!
330 CALL stats_2dfld (ng, tile, inlm, u2dvar, stats(1), 0, &
331 & lbi, ubi, lbj, ubj, ubar(:,:,1))
332 IF (domain(ng)%NorthEast_Corner(tile)) THEN
333 WRITE (stdout,10) trim(vname(2,idubar))//': '// &
334 & trim(vname(1,idubar)), &
335 & ng, stats(1)%min, stats(1)%max
336 END IF
337 CALL stats_2dfld (ng, tile, inlm, v2dvar, stats(2), 0, &
338 & lbi, ubi, lbj, ubj, vbar(:,:,1))
339 IF (domain(ng)%NorthEast_Corner(tile)) THEN
340 WRITE (stdout,10) trim(vname(2,idvbar))//': '// &
341 & trim(vname(1,idvbar)), &
342 & ng, stats(2)%min, stats(2)%max
343 END IF
344!
345!-----------------------------------------------------------------------
346! Initial conditions for free-surface (m).
347!-----------------------------------------------------------------------
348!
349#if defined CHANNEL && !defined ONLY_TS_IC
350 y0=0.5_r8*el(ng)
351# ifdef SOLVE3D
352 DO j=jstrt,jendt
353 DO i=istrt,iendt
354 val1=(yr(i,j)-y0)/guscale
355 val2=-u0*guscale*grid(ng)%f(i,j)*sqrt(pi)/(12.0_r8*g)
356 zeta(i,j,1)=val2*erf(val1)
357 END DO
358 END DO
359# else
360 DO j=jstrt,jendt
361 DO i=istrt,iendt
362 val1=(yr(i,j)-y0)/guscale
363 val2=-0.5_r8*u0*guscale*grid(ng)%f(i,j)*sqrt(pi)/g
364 zeta(i,j,1)=val2*erf(val1)
365 END DO
366 END DO
367# endif
368# ifdef DISTRIBUTE
369 CALL mp_bcasti (ng, model, exit_flag) ! in case of error in ERF
370# endif
371#elif defined KELVIN
372!! val1=1.0_r8 ! zeta0
373!! val2=2.0_r8*pi/(12.42_r8*3600.0_r8) ! M2 Tide period
374 DO j=jstrt,jendt
375 DO i=istrt,iendt
376!! zeta(i,j,1)=val1* &
377!! & EXP(-GRID(ng)%f(i,j)*GRID(ng)%yp(i,j)/ &
378!! & SQRT(g*GRID(ng)%h(i,j)))* &
379!! & COS(val2*GRID(ng)%xp(i,j)/ &
380!! & SQRT(g*GRID(ng)%h(i,j)))
381 zeta(i,j,1)=0.0_r8
382 END DO
383 END DO
384#elif defined SOLITON
385 IF (ng.eq.1) THEN
386 x0=2.0_r8*xl(ng)/3.0_r8
387 y0=0.5_r8*el(ng)
388 val1=0.395_r8
389 val2=0.771_r8*(val1*val1)
390 DO j=jstrt,jendt
391 DO i=istrt,iendt
392 x=xr(i,j)-x0
393 y=yr(i,j)-y0
394 val3=exp(-val1*x)
395 val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2)
396 zeta(i,j,1)=0.25_r8*val4*(6.0_r8*y*y+3.0_r8)* &
397 & exp(-0.5_r8*y*y)
398 END DO
399 END DO
400 ELSE
401 DO j=jstrt,jendt
402 DO i=istrt,iendt
403 zeta(i,j,1)=0.0_r8
404 END DO
405 END DO
406 END IF
407#elif defined SED_TEST1
408 val1=100.0_r8
409 DO j=jstrt,jendt
410 DO i=istrt,iendt
411 zeta(i,j,1)=9.0e-06_r8*real(i,r8)*val1
412 END DO
413 END DO
414#elif defined TEST_CHAN
415 DO j=jstrt,jendt
416 DO i=istrt,iendt
417 zeta(i,j,1)=-0.4040_r8*real(i,r8)/real(lm(ng)+1,r8)
418 END DO
419 END DO
420#else
421 DO j=jstrt,jendt
422 DO i=istrt,iendt
423 zeta(i,j,1)=0.0_r8
424 END DO
425 END DO
426#endif
427!
428! Report statistics.
429!
430 CALL stats_2dfld (ng, tile, inlm, r2dvar, stats(3), 0, &
431 & lbi, ubi, lbj, ubj, zeta(:,:,1))
432 IF (domain(ng)%NorthEast_Corner(tile)) THEN
433 WRITE (stdout,10) trim(vname(2,idfsur))//': '// &
434 & trim(vname(1,idfsur)), &
435 & ng, stats(3)%min, stats(3)%max
436 END IF
437
438#ifdef SOLVE3D
439!
440!-----------------------------------------------------------------------
441! Initial conditions for 3D momentum components (m/s).
442!-----------------------------------------------------------------------
443!
444# if defined CHANNEL && !defined ONLY_TS_IC
445 y0=0.5_r8*el(ng)
446 DO k=1,n(ng)
447 DO j=jstrt,jendt
448 DO i=istrp,iendt
449 val1=(grid(ng)%yu(i,j)-y0)/guscale
450 val2=(z_r(i,j,k)+z_r(i-1,j,k))/(h(i,j)+h(i-1,j))
451 val3=u0*(0.5_r8+val2+(0.5_r8*val2*val2))*exp(-val1*val1)
452 u(i,j,k,1)=val3
453 END DO
454 END DO
455 END DO
456 DO k=1,n(ng)
457 DO j=jstrp,jendt
458 DO i=istrt,iendt
459 v(i,j,k,1)=0.0_r8
460 END DO
461 END DO
462 END DO
463# elif defined RIVERPLUME2
464 DO k=1,n(ng)
465 DO j=jstrt,jendt
466 DO i=istrp,iendt
467 u(i,j,k,1)=0.0_r8
468 END DO
469 END DO
470 DO j=jstrp,jendt
471 DO i=istrt,iendt
472 v(i,j,k,1)=-0.05_r8*log((h(i,j)+z_r(i,j,k))/zob(ng))/ &
473 & (log(h(i,j)/zob(ng))-1.0_r8+zob(ng)/h(i,j))
474 END DO
475 END DO
476 END DO
477# elif defined SED_TEST1
478 DO k=1,n(ng)
479 DO j=jstrt,jendt
480 DO i=istrp,iendt
481 u(i,j,k,1)=-1.0_r8*log((h(i,j)+z_r(i,j,k))/zob(ng))/ &
482 & (log(h(i,j)/zob(ng))-1.0_r8+zob(ng)/h(i,j))
483 END DO
484 END DO
485 DO j=jstrp,jendt
486 DO i=istrt,iendt
487 v(i,j,k,1)=0.0_r8
488 END DO
489 END DO
490 END DO
491# elif defined SED_TOY
492 DO k=1,n(ng)
493 DO j=jstrt,jendt
494 DO i=istrp,iendt
495 u(i,j,k,1)=1.0_r8
496!! u(i,j,k,1)=-1.0_r8
497!! u(i,j,k,1)=0.0_r8
498 END DO
499 END DO
500 DO j=jstrp,jendt
501 DO i=istrt,iendt
502 v(i,j,k,1)=0.0_r8
503 END DO
504 END DO
505 END DO
506# else
507 DO k=1,n(ng)
508 DO j=jstrt,jendt
509 DO i=istrp,iendt
510 u(i,j,k,1)=0.0_r8
511 END DO
512 END DO
513 DO j=jstrp,jendt
514 DO i=istrt,iendt
515 v(i,j,k,1)=0.0_r8
516 END DO
517 END DO
518 END DO
519# endif
520!
521! Report statistics.
522!
523 CALL stats_3dfld (ng, tile, inlm, u3dvar, stats(4), 0, &
524 & lbi, ubi, lbj, ubj, 1, n(ng), u(:,:,:,1))
525 IF (domain(ng)%NorthEast_Corner(tile)) THEN
526 WRITE (stdout,10) trim(vname(2,iduvel))//': '// &
527 & trim(vname(1,iduvel)), &
528 & ng, stats(4)%min, stats(4)%max
529 END IF
530 CALL stats_3dfld (ng, tile, inlm, v3dvar, stats(5), 0, &
531 & lbi, ubi, lbj, ubj, 1, n(ng), v(:,:,:,1))
532 IF (domain(ng)%NorthEast_Corner(tile)) THEN
533 WRITE (stdout,10) trim(vname(2,idvvel))//': '// &
534 & trim(vname(1,idvvel)), &
535 & ng, stats(5)%min, stats(5)%max
536 END IF
537!
538!-----------------------------------------------------------------------
539! Initial conditions for tracer type variables.
540!-----------------------------------------------------------------------
541!
542! Set initial conditions for potential temperature (Celsius) and
543! salinity (PSU).
544!
545# if defined BENCHMARK
546 val1=(44.69_r8/39.382_r8)**2
547 val2=val1*(rho0*800.0_r8/g)*(5.0e-05_r8/((42.689_r8/44.69_r8)**2))
548 DO k=1,n(ng)
549 DO j=jstrt,jendt
550 DO i=istrt,iendt
551 t(i,j,k,1,itemp)=val2*exp(z_r(i,j,k)/800.0_r8)* &
552 & (0.6_r8-0.4_r8*tanh(z_r(i,j,k)/800.0_r8))
553# ifdef SALINITY
554 t(i,j,k,1,isalt)=35.0_r8
555# endif
556 END DO
557 END DO
558 END DO
559# elif defined BASIN
560 val1=(44.69_r8/39.382_r8)**2
561 val2=val1*(rho0*800.0_r8/g)*(5.0e-05_r8/((42.689_r8/44.69_r8)**2))
562 DO k=1,n(ng)
563 DO j=jstrt,jendt
564 DO i=istrt,iendt
565 t(i,j,k,1,itemp)=val2*exp(z_r(i,j,k)/800.0_r8)* &
566 & (0.6_r8-0.4_r8*tanh(z_r(i,j,k)/800.0_r8))
567 END DO
568 END DO
569 END DO
570# elif defined BL_TEST
571 DO k=1,n(ng)
572 DO j=jstrt,jendt
573 DO i=istrt,iendt
574 val1=tanh(1.1_r8*z_r(i,j,k)+11.0_r8)
575 t(i,j,k,1,itemp)=t0(ng)+6.25_r8*val1
576# ifdef SALINITY
577 t(i,j,k,1,isalt)=s0(ng)-0.75_r8*val1
578# endif
579 END DO
580 END DO
581 END DO
582# elif defined CHANNEL
583 y0=0.5_r8*el(ng)
584 DO k=1,n(ng)
585 DO j=jstrt,jendt
586 DO i=istrt,iendt
587 val1=(yr(i,j)-y0)/guscale
588 val2=-0.5_r8*u0*guscale*grid(ng)%f(i,j)*sqrt(pi)/ &
589 & (tcoef(ng)*g*h(i,j))
590 val3=(val2*erf(val1)+t0(ng))*(1.0_r8+z_r(i,j,k)/h(i,j))
591 t(i,j,k,1,itemp)=val3
592 END DO
593 END DO
594 END DO
595# ifdef DISTRIBUTE
596 CALL mp_bcasti (ng, model, exit_flag) ! in case of error in ERF
597# endif
598# elif defined CANYON
599 DO k=1,n(ng)
600 DO j=jstrt,jendt
601 DO i=istrt,iendt
602 t(i,j,k,1,itemp)=3.488_r8*exp(z_r(i,j,k)/800.0_r8)* &
603 & (1.0_r8-(2.0_r8/3.0_r8)* &
604 & tanh(z_r(i,j,k)/800.0_r8))
605 END DO
606 END DO
607 END DO
608# elif defined CHANNEL_NECK
609 DO k=1,n(ng)
610 DO j=jstrt,jendt
611 DO i=istrt,iendt
612 t(i,j,k,1,itemp)=20.0_r8
613!! t(i,j,k,1,itemp)=14.0_r8+8.0_r8*EXP(z_r(i,j,k)/50.0_r8)
614 END DO
615 END DO
616 END DO
617# elif defined COUPLING_TEST
618 val1=40.0_r8
619 DO k=1,n(ng)
620 DO j=jstrt,jendt
621 DO i=istrt,iendt
622 t(i,j,k,1,itemp)=val1*exp(z_r(i,j,k)/800.0_r8)* &
623 & (0.6_r8-0.4_r8*tanh(z_r(i,j,k)/800.0_r8))+ &
624 & 1.5_r8
625# ifdef SALINITY
626 t(i,j,k,1,isalt)=35.0_r8
627# endif
628 END DO
629 END DO
630 END DO
631# elif defined DOUBLE_GYRE
632 val1=(44.69_r8/39.382_r8)**2
633 val2=val1*(rho0*100.0_r8/g)*(5.0e-5_r8/((42.689_r8/44.69_r8)**2))
634 DO k=1,n(ng)
635 DO j=jstrt,jendt
636 DO i=istrt,iendt
637 val3=t0(ng)+val2*exp(z_r(i,j,k)/100.0_r8)* &
638 & (10.0_r8-0.4_r8*tanh(z_r(i,j,k)/100.0_r8))
639 val4=yr(i,j)/el(ng)
640 t(i,j,k,1,itemp)=val3-3.0_r8*val4
641# ifdef SALINITY
642 t(i,j,k,1,isalt)=34.5_r8-0.001_r8*z_r(i,j,k)-val4
643# endif
644 END DO
645 END DO
646 END DO
647# elif defined ESTUARY_TEST
648 DO k=1,n(ng)
649 DO j=jstrt,jendt
650 DO i=istrt,iendt
651 t(i,j,k,1,itemp)=10.0_r8
652# ifdef SALINITY
653 IF (xr(i,j).le.30000.0_r8) then
654 t(i,j,k,1,isalt)=30.0_r8
655 ELSEIF (xr(i,j).le.80000.0_r8) then
656 t(i,j,k,1,isalt)=(80000.0_r8-xr(i,j))/50000.0_r8*30.0_r8
657 ELSE
658 t(i,j,k,1,isalt)=0.0_r8
659 END IF
660# endif
661 END DO
662 END DO
663 END DO
664# elif defined FLT_TEST
665 DO k=1,n(ng)
666 DO j=jstrt,jendt
667 DO i=istrt,iendt
668 t(i,j,k,1,itemp)=t0(ng)
669 END DO
670 END DO
671 END DO
672# elif defined GRAV_ADJ
673 DO k=1,n(ng)
674 DO j=jstrt,jendt
675 DO i=istrt,min((lm(ng)+1)/2,iendt)
676 t(i,j,k,1,itemp)=t0(ng)+5.0_r8
677# ifdef SALINITY
678 t(i,j,k,1,isalt)=0.0_r8
679# endif
680 END DO
681 DO i=max((lm(ng)+1)/2+1,istrt),iendt
682 t(i,j,k,1,itemp)=t0(ng)
683# ifdef SALINITY
684 t(i,j,k,1,isalt)=s0(ng)
685# endif
686 END DO
687!! DO i=IstrT,IendT
688!! IF (i.lt.Lm(ng)/2) THEN
689!! t(i,j,k,1,itemp)=T0(ng)+5.0_r8
690!! ELSE IF (i.eq.Lm(ng)/2) THEN
691!! t(i,j,k,1,itemp)=T0(ng)+4.0_r8
692!! ELSE IF (i.eq.Lm(ng)/2+1) THEN
693!! t(i,j,k,1,itemp)=T0(ng)+1.0_r8
694!! ELSE
695!! t(i,j,k,1,itemp)=T0(ng)
696!! END IF
697!! END DO
698 END DO
699 END DO
700# elif defined LAB_CANYON
701 DO k=1,n(ng)
702 DO j=jstrt,jendt
703 DO i=istrt,iendt
704 t(i,j,k,1,itemp)=-659.34183_r8*z_r(i,j,k)
705 END DO
706 END DO
707 END DO
708# elif defined LAKE_SIGNELL
709 DO k=1,n(ng)
710 DO j=jstrt,jendt
711 DO i=istrt,iendt
712 t(i,j,k,1,itemp)=10.0_r8
713# ifdef SALINITY
714 t(i,j,k,1,isalt)=30.0_r8
715# endif
716 END DO
717 END DO
718 END DO
719# elif defined LMD_TEST
720 DO k=1,n(ng)
721 DO j=jstrt,jendt
722 DO i=istrt,iendt
723 t(i,j,k,1,itemp)=min(13.0_r8, &
724 & 7.0_r8+0.2_r8*(z_r(i,j,k)+50.0_r8))
725# ifdef SALINITY
726 t(i,j,k,1,isalt)=35.0_r8
727# endif
728 END DO
729 END DO
730 END DO
731# elif defined MIXED_LAYER
732 DO k=1,n(ng)
733 DO j=jstrt,jendt
734 DO i=istrt,iendt
735 t(i,j,k,1,itemp)=10.0_r8+3.0_r8*(z_r(i,j,k)+h(i,j))/ &
736 & h(i,j)
737# ifdef SALINITY
738 t(i,j,k,1,isalt)=s0(ng)
739# endif
740 END DO
741 END DO
742 END DO
743# elif defined NJ_BIGHT
744 DO k=1,n(ng)
745 DO j=jstrt,jendt
746 DO i=istrt,iendt
747 depth=z_r(i,j,k)
748 IF (depth.ge.-15.0_r8) THEN
749 t(i,j,k,1,itemp)= 2.049264257728403e+01_r8-depth* &
750 & (2.640850848793918e-01_r8+depth* &
751 & (2.751125328535212e-01_r8+depth* &
752 & (9.207489761648872e-02_r8+depth* &
753 & (1.449075725742839e-02_r8+depth* &
754 & (1.078215685912076e-03_r8+depth* &
755 & (3.240318053903974e-05_r8+ &
756 & 1.262826857690271e-07_r8*depth))))))
757# ifdef SALINITY
758 t(i,j,k,1,isalt)= 3.066489149193135e+01_r8-depth* &
759 & (1.476725262946735e-01_r8+depth* &
760 & (1.126455760313399e-01_r8+depth* &
761 & (3.900923281871022e-02_r8+depth* &
762 & (6.939014937447098e-03_r8+depth* &
763 & (6.604436696792939e-04_r8+depth* &
764 & (3.191792361954220e-05_r8+ &
765 & 6.177352634409320e-07_r8*depth))))))
766# endif
767 ELSE
768 t(i,j,k,1,itemp)=14.6_r8+ &
769 & 6.70_r8*tanh(1.1_r8*depth+15.9_r8)
770# ifdef SALINITY
771 t(i,j,k,1,isalt)=31.3_r8- &
772 & 0.55_r8*tanh(1.1_r8*depth+15.9_r8)
773# endif
774 END IF
775 END DO
776 END DO
777 END DO
778# elif defined OVERFLOW
779 DO k=1,n(ng)
780 DO j=jstrt,jendt
781 DO i=istrt,iendt
782 t(i,j,k,1,itemp)=t0(ng)-0.5_r8*t0(ng)*(1.0_r8+ &
783 & tanh((yr(i,j)-60000.0_r8)/2000.0_r8))
784 END DO
785 END DO
786 END DO
787# elif defined RIVERPLUME1
788 DO k=1,n(ng)
789 DO j=jstrt,jendt
790 DO i=istrt,iendt
791 t(i,j,k,1,itemp)=t0(ng)+0.01_r8*real(k,r8)
792# ifdef SALINITY
793 t(i,j,k,1,isalt)=s0(ng)
794# endif
795 END DO
796 END DO
797 END DO
798# elif defined RIVERPLUME2
799 DO k=1,n(ng)
800 DO j=jstrt,jendt
801 DO i=istrt,iendt
802 t(i,j,k,1,itemp)=t0(ng)
803# ifdef SALINITY
804 t(i,j,k,1,isalt)=s0(ng)
805# endif
806 END DO
807 END DO
808 END DO
809# elif defined SEAMOUNT
810 DO k=1,n(ng)
811 DO j=jstrt,jendt
812 DO i=istrt,iendt
813 t(i,j,k,1,itemp)=t0(ng)+7.5_r8*exp(z_r(i,j,k)/1000.0_r8)
814 END DO
815 END DO
816 END DO
817# elif defined SED_TEST1
818 DO k=1,n(ng)
819 DO j=jstrt,jendt
820 DO i=istrt,iendt
821 t(i,j,k,1,itemp)=20.0_r8
822# ifdef SALINITY
823 t(i,j,k,1,isalt)=0.0_r8
824# endif
825 END DO
826 END DO
827 END DO
828# elif defined UPWELLING
829 DO k=1,n(ng)
830 DO j=jstrt,jendt
831 DO i=istrt,iendt
832 t(i,j,k,1,itemp)=t0(ng)+8.0_r8*exp(z_r(i,j,k)/50.0_r8)
833!! t(i,j,k,1,itemp)=T0(ng)+(z_r(i,j,k)+75.0_r8)/150.0_r8+
834!! & 4.0_r8*(1.0_r8+TANH((z_r(i,j,k)+35.0_r8)/
835!! & 6.5_r8))
836# ifdef SALINITY
837!! t(i,j,k,1,isalt)=1.0E-04_r8*yr(i,j)-S0(ng)
838 t(i,j,k,1,isalt)=s0(ng)
839!! IF (j.lt.Mm(ng)/2) THEN
840!! t(i,j,k,1,isalt)=0.0_r8
841!! ELSE IF (j.eq.Mm(ng)/2) THEN
842!! t(i,j,k,1,isalt)=0.5_r8
843!! ELSE IF (j.gt.Mm(ng)/2) THEN
844!! t(i,j,k,1,isalt)=1.0_r8
845!! END IF
846# endif
847 END DO
848 END DO
849 END DO
850# elif defined WINDBASIN
851 DO k=1,n(ng)
852 DO j=jstrt,jendt
853 DO i=istrt,iendt
854 t(i,j,k,1,itemp)=20.0_r8 ! homogeneous
855!! t(i,j,k,1,itemp)=14.0_r8+8.0_r8*EXP(z_r(i,j,k)/50.0_r8)- &
856!! & T0(ng) ! stratified
857 END DO
858 END DO
859 END DO
860# else
861 DO k=1,n(ng)
862 DO j=jstrt,jendt
863 DO i=istrt,iendt
864 t(i,j,k,1,itemp)=t0(ng)
865# ifdef SALINITY
866 t(i,j,k,1,isalt)=s0(ng)
867# endif
868 END DO
869 END DO
870 END DO
871# endif
872!
873! Report statistics.
874!
875 DO itrc=1,nat
876 CALL stats_3dfld (ng, tile, inlm, r3dvar, stats(itrc+5), 0, &
877 & lbi, ubi, lbj, ubj, 1, n(ng), t(:,:,:,1,itrc))
878 IF (domain(ng)%NorthEast_Corner(tile)) THEN
879 WRITE (stdout,10) trim(vname(2,idtvar(itrc)))//': '// &
880 & trim(vname(1,idtvar(itrc))), &
881 & ng, stats(itrc+5)%min, stats(itrc+5)%max
882 END IF
883 END DO
884#endif
885!
886 10 FORMAT (3x,' ANA_INITIAL - ',a,/,19x, &
887 & '(Grid = ',i2.2,', Min = ',1p,e15.8,0p, &
888 & ' Max = ',1p,e15.8,0p,')')
889!
890 RETURN
891 END SUBROUTINE ana_nlminitial_tile
892
893#ifdef TANGENT
894!
895!***********************************************************************
896 SUBROUTINE ana_tlminitial_tile (ng, tile, model, &
897 & LBi, UBi, LBj, UBj, &
898 & IminS, ImaxS, JminS, JmaxS, &
899 & kstp, &
900# ifdef SOLVE3D
901 & nstp, &
902 & tl_u, tl_v, tl_t, &
903# endif
904 & tl_ubar, tl_vbar, tl_zeta)
905!***********************************************************************
906!
907 USE mod_param
908 USE mod_scalars
909!
910! Imported variable declarations.
911!
912 integer, intent(in) :: ng, tile, model
913 integer, intent(in) :: LBi, UBi, LBj, UBj
914 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
915 integer, intent(in) :: kstp
916# ifdef SOLVE3D
917 integer, intent(in) :: nstp
918# endif
919!
920# ifdef ASSUMED_SHAPE
921# ifdef SOLVE3D
922 real(r8), intent(out) :: tl_u(LBi:,LBj:,:,:)
923 real(r8), intent(out) :: tl_v(LBi:,LBj:,:,:)
924 real(r8), intent(out) :: tl_t(LBi:,LBj:,:,:,:)
925# endif
926 real(r8), intent(out) :: tl_ubar(LBi:,LBj:,:)
927 real(r8), intent(out) :: tl_vbar(LBi:,LBj:,:)
928 real(r8), intent(out) :: tl_zeta(LBi:,LBj:,:)
929# else
930# ifdef SOLVE3D
931 real(r8), intent(out) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
932 real(r8), intent(out) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
933 real(r8), intent(out) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
934# endif
935 real(r8), intent(out) :: tl_ubar(LBi:UBi,LBj:UBj,:)
936 real(r8), intent(out) :: tl_vbar(LBi:UBi,LBj:UBj,:)
937 real(r8), intent(out) :: tl_zeta(LBi:UBi,LBj:UBj,:)
938# endif
939!
940! Local variable declarations.
941!
942 integer :: i, itrc, j, k
943
944# include "set_bounds.h"
945!
946!-----------------------------------------------------------------------
947! Initial conditions for tangent linear 2D momentum (s/m) components.
948!-----------------------------------------------------------------------
949!
950 DO j=jstrt,jendt
951 DO i=istrp,iendt
952 tl_ubar(i,j,kstp)=0.0_r8
953 END DO
954 END DO
955 DO j=jstrp,jendt
956 DO i=istrt,iendt
957 tl_vbar(i,j,kstp)=0.0_r8
958 END DO
959 END DO
960!
961!-----------------------------------------------------------------------
962! Initial conditions for tangent linear free-surface (1/m).
963!-----------------------------------------------------------------------
964!
965 DO j=jstrt,jendt
966 DO i=istrt,iendt
967 tl_zeta(i,j,kstp)=0.0_r8
968 END DO
969 END DO
970# ifdef SOLVE3D
971!
972!-----------------------------------------------------------------------
973! Initial conditions for tangent linear 3D momentum components (s/m).
974!-----------------------------------------------------------------------
975!
976 DO k=1,n(ng)
977 DO j=jstrt,jendt
978 DO i=istrp,iendt
979 tl_u(i,j,k,nstp)=0.0_r8
980 END DO
981 END DO
982 DO j=jstrp,jendt
983 DO i=istrt,iendt
984 tl_v(i,j,k,nstp)=0.0_r8
985 END DO
986 END DO
987 END DO
988!
989!-----------------------------------------------------------------------
990! Initial conditions for tangent linear active tracers (1/Tunits).
991!-----------------------------------------------------------------------
992!
993 DO itrc=1,nat
994 DO k=1,n(ng)
995 DO j=jstrt,jendt
996 DO i=istrt,iendt
997 tl_t(i,j,k,nstp,itrc)=0.0_r8
998 END DO
999 END DO
1000 END DO
1001 END DO
1002# endif
1003!
1004 RETURN
1005 END SUBROUTINE ana_tlminitial_tile
1006#endif
1007
1008#ifdef ADJOINT
1009!
1010!***********************************************************************
1011 SUBROUTINE ana_adminitial_tile (ng, tile, model, &
1012 & LBi, UBi, LBj, UBj, &
1013 & IminS, ImaxS, JminS, JmaxS, &
1014 & knew, &
1015# ifdef SOLVE3D
1016 & nstp, &
1017 & ad_u, ad_v, ad_t, &
1018# endif
1019 & ad_ubar, ad_vbar, ad_zeta)
1020!***********************************************************************
1021!
1022 USE mod_param
1023 USE mod_scalars
1024!
1025! Imported variable declarations.
1026!
1027 integer, intent(in) :: ng, tile, model
1028 integer, intent(in) :: LBi, UBi, LBj, UBj
1029 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1030 integer, intent(in) :: knew
1031# ifdef SOLVE3D
1032 integer, intent(in) :: nstp
1033# endif
1034!
1035# ifdef ASSUMED_SHAPE
1036# ifdef SOLVE3D
1037 real(r8), intent(out) :: ad_u(LBi:,LBj:,:,:)
1038 real(r8), intent(out) :: ad_v(LBi:,LBj:,:,:)
1039 real(r8), intent(out) :: ad_t(LBi:,LBj:,:,:,:)
1040# endif
1041 real(r8), intent(out) :: ad_ubar(LBi:,LBj:,:)
1042 real(r8), intent(out) :: ad_vbar(LBi:,LBj:,:)
1043 real(r8), intent(out) :: ad_zeta(LBi:,LBj:,:)
1044# else
1045# ifdef SOLVE3D
1046 real(r8), intent(out) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1047 real(r8), intent(out) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1048 real(r8), intent(out) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1049# endif
1050 real(r8), intent(out) :: ad_ubar(LBi:UBi,LBj:UBj,:)
1051 real(r8), intent(out) :: ad_vbar(LBi:UBi,LBj:UBj,:)
1052 real(r8), intent(out) :: ad_zeta(LBi:UBi,LBj:UBj,:)
1053# endif
1054!
1055! Local variable declarations.
1056!
1057 integer :: i, itrc, j, k
1058
1059# include "set_bounds.h"
1060!
1061!-----------------------------------------------------------------------
1062! Initial conditions for adjoint 2D momentum (s/m) components.
1063!-----------------------------------------------------------------------
1064!
1065 DO j=jstrt,jendt
1066 DO i=istrp,iendt
1067 ad_ubar(i,j,knew)=0.0_r8
1068 END DO
1069 END DO
1070 DO j=jstrp,jendt
1071 DO i=istrt,iendt
1072 ad_vbar(i,j,knew)=0.0_r8
1073 END DO
1074 END DO
1075!
1076!-----------------------------------------------------------------------
1077! Initial conditions for adjoint free-surface (1/m).
1078!-----------------------------------------------------------------------
1079!
1080 DO j=jstrt,jendt
1081 DO i=istrt,iendt
1082 ad_zeta(i,j,knew)=0.0_r8
1083 END DO
1084 END DO
1085# ifdef SOLVE3D
1086!
1087!-----------------------------------------------------------------------
1088! Initial conditions for adjoint 3D momentum components (s/m).
1089!-----------------------------------------------------------------------
1090!
1091 DO k=1,n(ng)
1092 DO j=jstrt,jendt
1093 DO i=istrp,iendt
1094 ad_u(i,j,k,nstp)=0.0_r8
1095 END DO
1096 END DO
1097 DO j=jstrp,jendt
1098 DO i=istrt,iendt
1099 ad_v(i,j,k,nstp)=0.0_r8
1100 END DO
1101 END DO
1102 END DO
1103!
1104!-----------------------------------------------------------------------
1105! Initial conditions for adjoint active tracers (1/Tunits).
1106!-----------------------------------------------------------------------
1107!
1108 DO itrc=1,nat
1109 DO k=1,n(ng)
1110 DO j=jstrt,jendt
1111 DO i=istrt,iendt
1112 ad_t(i,j,k,nstp,itrc)=0.0_r8
1113 END DO
1114 END DO
1115 END DO
1116 END DO
1117# endif
1118!
1119 RETURN
1120 END SUBROUTINE ana_adminitial_tile
1121#endif
subroutine ana_adminitial_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, knew, nstp, ad_u, ad_v, ad_t, ad_ubar, ad_vbar, ad_zeta)
subroutine ana_initial(ng, tile, model)
Definition ana_initial.h:3
subroutine ana_nlminitial_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, h, lonr, latr, xr, yr, z_r, u, v, t, ubar, vbar, zeta)
subroutine ana_tlminitial_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kstp, nstp, tl_u, tl_v, tl_t, tl_ubar, tl_vbar, tl_zeta)
Definition erf.F:2
real(r8) function, public erf(x)
Definition erf.F:43
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
integer stdout
integer idubar
integer idvvel
integer idfsur
integer, dimension(:), allocatable idtvar
logical lanafile
integer iduvel
character(len=maxlen), dimension(6, 0:nv) vname
character(len=256), dimension(39) ananame
integer idvbar
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer nat
Definition mod_param.F:499
integer, parameter inlm
Definition mod_param.F:662
integer, parameter irpm
Definition mod_param.F:664
integer, parameter r3dvar
Definition mod_param.F:721
integer, parameter iadm
Definition mod_param.F:665
integer, parameter u3dvar
Definition mod_param.F:722
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer, dimension(:), allocatable lm
Definition mod_param.F:455
integer, parameter u2dvar
Definition mod_param.F:718
integer, parameter itlm
Definition mod_param.F:663
integer, parameter r2dvar
Definition mod_param.F:717
integer, parameter v2dvar
Definition mod_param.F:719
integer, parameter v3dvar
Definition mod_param.F:723
real(r8), dimension(:), allocatable t0
real(r8), dimension(:), allocatable el
real(dp), parameter large
real(r8), dimension(:), allocatable s0
real(r8), dimension(:), allocatable tcoef
integer exit_flag
integer isalt
integer itemp
real(r8), dimension(:), allocatable zob
real(dp) g
real(r8), dimension(:), allocatable xl
real(dp) rho0
real(dp), parameter pi
integer, dimension(:), allocatable kstp
integer, dimension(:), allocatable knew
integer, dimension(:), allocatable nstp
subroutine, public stats_2dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, f, fmask, debug)
Definition stats.F:47
subroutine, public stats_3dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, lbk, ubk, f, fmask, debug)
Definition stats.F:342