Context Navigation

analytical.F @ 5

Last change on this file since 5 was 5, checked in by cholod, 13 years ago
patch_romsagrif_12_04_2011.tar : Bug fixes
File size: 53.3 KB

Line
1	! $Id: analytical.F 703 2011-04-11 15:57:49Z gcambon $
2	!
3	!======================================================================
4	! ROMS_AGRIF is a branch of ROMS developped at IRD and INRIA, in France
5	! The two other branches from UCLA (Shchepetkin et al)
6	! and Rutgers University (Arango et al) are under MIT/X style license.
7	! ROMS_AGRIF specific routines (nesting) are under CeCILL-C license.
8	!
9	! ROMS_AGRIF website : http://roms.mpl.ird.fr
10	!======================================================================
11	!
12	#include "cppdefs.h"
13	!
14	! ANALYTICAL PACKAGE:
15	!--------------------------------------------------------------------
16	!
17	! This package is used to provide various analytical fields to the
18	! model when appropriate.
19	!
20	! Routines:
21	!
22	! ana_bmflux_tile Analytical kinematic bottom momentum flux.
23	! ana_btflux_tile Analytical kinematic bottom flux of tracer
24	! type variables.
25	! ana_bsedim_tile Analytical bottom sediment grain size
26	! and density.
27	! ana_smflux_tile Analytical kinematic surface momentum flux
28	! (wind stress).
29	! ana_srflux_tile Analytical kinematic surface shortwave
30	! radiation.
31	! ana_ssh_tile Analytical sea surface height climatology.
32	! ana_sst_tile Analytical sea surface temperature and dQdSST
33	! which are used during heat flux correction.
34	! ana_sss_tile Analytical sea surface salinity which is used
35	! during salt flux correction.
36	! ana_stflux_tile Analytical kinematic surface flux of tracer type
37	! variables.
38	! ana_tclima_tile Analytical tracer climatology fields.
39	! ana_uclima_tile Analytical tracer climatology fields.
40	! ana_wwave_tile Analytical wind induced wave amplitude,
41	! direction and period.
42	! ana_sediment_tile Analytical sediment
43	! ana_psource_tile Analytical point source
44	! ana_bry_tile Analytical boundary forcing.
45	!
46	!-------------------------------------------------------------------
47	!
48
49	# if !defined OPENMP
50	integer function omp_get_thread_num()
51	omp_get_thread_num=0
52	return
53	end
54	integer function omp_get_num_threads()
55	omp_get_num_threads=1
56	return
57	end
58	# endif
59
60	#ifdef ANA_BMFLUX
61	subroutine ana_bmflux_tile (Istr,Iend,Jstr,Jend)
62	!
63	!---------------------------------------------------------------------
64	! This routine sets kinematic bottom momentum flux (bottom stress)
65	! "bustr" and "bvstr" [m^2/s^2] using an analytical expression.
66	!---------------------------------------------------------------------
67	!
68	implicit none
69	# include "param.h"
70	# include "grid.h"
71	# include "forces.h"
72	# include "scalars.h"
73	integer Istr,Iend,Jstr,Jend, i,j
74	!
75	# include "compute_auxiliary_bounds.h"
76	!
77	do j=JstrR,JendR
78	do i=Istr,IendR
79	bustr(i,j)=???
80	enddo
81	enddo
82	do j=Jstr,JendR
83	do i=IstrR,IendR
84	bvstr(i,j)=???
85	enddo
86	enddo
87	return
88	end
89	#endif /* ANA_BMFLUX */
90	#ifdef SOLVE3D
91	# if defined ANA_BTFLUX \|\| defined ANA_BSFLUX \|\| defined ANA_BPFLUX
92	subroutine ana_btflux_tile (Istr,Iend,Jstr,Jend, itrc)
93	!
94	!---------------------------------------------------------------------
95	! This routine sets kinematic bottom flux of tracer type variables
96	! [tracer units m/s].
97	!
98	! On Input:
99	! itrc Tracer type array index.
100	!---------------------------------------------------------------------
101	!
102	implicit none
103	# include "param.h"
104	# include "grid.h"
105	# include "forces.h"
106	# include "scalars.h"
107	integer itrc, Istr,Iend,Jstr,Jend, i,j
108	!
109	# include "compute_auxiliary_bounds.h"
110	!
111	if (itrc.eq.itemp) then
112	!
113	! Set kinematic bottom heat flux [degC m/s] at horizontal RHO-points.
114	!--------------------------------------------------------------------
115	!
116	# if defined BASIN \|\| defined CANYON_A \|\| defined CANYON_B \
117	\|\| defined EQUATOR \|\| defined GRAV_ADJ \|\| defined INNERSHELF \
118	\|\| defined OVERFLOW \|\| defined REGIONAL \|\| defined RIVER \
119	\|\| defined SEAMOUNT \|\| defined SED_TEST2 \|\| defined SHELFRONT \
120	\|\| defined UPWELLING \|\| defined VORTEX \|\| defined INTERNAL
121	do j=JstrR,JendR
122	do i=IstrR,IendR
123	btflx(i,j,itemp)=0.
124	enddo
125	enddo
126	# else
127	do j=JstrR,JendR
128	do i=IstrR,IendR
129	btflx(i,j,itemp)=???
130	enddo
131	enddo
132	# endif
133	# ifdef SALINITY
134	elseif (itrc.eq.isalt) then
135	!
136	! Set kinematic bottom salt flux (m/s) at horizontal RHO-points,
137	! scaling by bottom salinity is done in STEP3D.
138	!---------------------------------------------------------------------
139	!
140	# if defined EQUATOR \|\| defined INNERSHELF \|\| defined REGIONAL \
141	\|\| defined RIVER \|\| defined SED_TEST2 \|\| defined SHELFRONT \
142	\|\| defined UPWELLING \|\| defined SEAMOUNT \|\| defined VORTEX
143	do j=JstrR,JendR
144	do i=IstrR,IendR
145	btflx(i,j,isalt)=0.
146	enddo
147	enddo
148	# else
149	do j=JstrR,JendR
150	do i=IstrR,IendR
151	btflx(i,j,isalt)=???
152	enddo
153	enddo
154	# endif
155	# endif /* SALINITY */
156	else
157	!
158	!---------------------------------------------------------------------
159	! Set kinematic surface flux of additional tracers,
160	! for example sediments, bio..., to zero
161	!---------------------------------------------------------------------
162	!
163	do j=JstrR,JendR
164	do i=IstrR,IendR
165	btflx(i,j,itrc)=0.
166	enddo
167	enddo
168	endif
169	return
170	end
171	# endif /* ANA_BTFLUX */
172	#endif /* SOLVE3D */
173	!
174	#if defined ANA_BSEDIM && defined BBL
175	subroutine ana_bsedim (tile)
176	implicit none
177	# include "param.h"
178	integer tile
179	# include "compute_tile_bounds.h"
180	call ana_bsedim_tile (Istr,Iend,Jstr,Jend)
181	return
182	end
183	subroutine ana_bsedim_tile (Istr,Iend,Jstr,Jend)
184	!
185	!---------------------------------------------------------------------
186	! This routine sets initial bottom sediment grain diameter size [m]
187	! and density used in the bottom boundary formulation [kg/m^3].
188	!---------------------------------------------------------------------
189	!
190	implicit none
191	# include "param.h"
192	# include "bbl.h"
193	# include "grid.h"
194	# include "scalars.h"
195	integer Istr,Iend,Jstr,Jend, i,j
196	!
197	# include "compute_auxiliary_bounds.h"
198	!
199	# if defined SED_TEST2 \|\| defined REGIONAL
200	!
201	! taucb=critical threshold stress for initiation of motion
202	! (=bedload for coarse grains).
203	! critical suspension stress: ustar_crit=0.8*w_set
204	!
205	! determine taucb from Shields curve, fit provided by
206	! Soulsby & Whitehouse 1997, Threshold of sediment motion
207	! in coastal environments, Proc. Pacific Coasts and Ports
208	! '97 Conf., pp 149--154, Univ Canterbury, Nw Zealand.
209	!
210	! visk=1.3e-3/rhow; (kinem. visc., nu=mu/rhow)
211	! D=d50(g(Sdens/rhow-1)/(visk^2))^0.33333
212	! thetcr=0.3./(1+1.2D) + 0.055(1-exp(-0.02*D))
213	! taucb=thetcr.(g(sdens-rhow).*d50);
214	!
215	! Souslby's (1997) estimate of settling velocity
216	! w_set = visk(sqrt(10.36^2+1.049D^3)-10.36)/d50 [m/s]
217	! with D as above
218	!
219	do j=JstrR,JendR
220	do i=IstrR,IendR
221	Ssize(i,j)=1.5e-4 ! d50 [m]
222	Sdens(i,j)=2650.0 ! rho sediment [kg/m^3]
223	taucb(i,j)=0.16 ! critical bedload stress [N/m^2]
224	w_set(i,j)=0.013 ! analytical settling velocity [m/s]
225	Hripple(i,j)=0.01 ! analytical initial ripple height [m]
226	Lripple(i,j)=0.10 ! analytical initial ripple length [m]
227	enddo
228	enddo
229	# else
230	ANA_BSEDIM: no values provided for SSIZE and SDENS.
231	# endif
232	return
233	end
234	#endif /* ANA_BSEDIM && BBL */
235
236	#ifdef ANA_SMFLUX
237	subroutine ana_smflux_tile (Istr,Iend,Jstr,Jend)
238	!
239	! Sets kinematic surface momentum flux (wind stress) "sustr" and "svstr"
240	! [m^2/s^2] using an analytical expression.
241	!
242	# ifdef AGRIF
243	use Agrif_UTIL
244	# endif
245	implicit none
246	# include "param.h"
247	# include "grid.h"
248	# include "forces.h"
249	# include "scalars.h"
250	integer Istr,Iend,Jstr,Jend, i,j
251	real Ewind, Nwind, dircoef, windamp
252	real cff1, cff2
253	! data windamp /0./
254	data Ewind, Nwind, dircoef /0., 0., 0./
255	! save windamp
256	!
257	#include "compute_extended_bounds.h"
258	!
259	! Set kinematic surface momentum flux (wind stress) component in the
260	! XI-direction (m^2/s^2) at horizontal U-points.
261	!
262
263	windamp = 0.
264
265	# ifdef BASIN
266	cff1=0.0001 * 0.5(1.+tanh((time-6.86400.)/(3.*86400.)))
267	cff2=2.*pi/el
268	do j=JstrR,JendR
269	do i=IstrR,IendR
270	sustr(i,j)=-cff1cos(cff2yr(i,j))
271	enddo
272	enddo
273	# elif defined CANYON_A \|\| defined CANYON_B
274	do j=JstrR,JendR
275	do i=IstrR,IendR
276	sustr(i,j)=0.00010.5sin(2.pitdays/10.)*
277	& (1.-tanh((yr(i,j)-0.5*el)/10000.))
278	enddo
279	enddo
280	# elif defined EQUATOR
281	do j=JstrR,JendR
282	do i=IstrR,IendR
283	sustr(i,j)=-0.05/rho0
284	enddo
285	enddo
286	# elif defined SED_TEST2
287	do j=JstrR,JendR
288	do i=IstrR,IendR
289	windamp=0.5+
290	& 0.5*TANH((time-user(9))/user(10))
291	sustr(i,j)=windamp*user(1)
292	enddo
293	enddo
294	# elif defined UPWELLING
295	if (tdays.le.2.) then
296	windamp=-0.1sin(pitdays/4.)/rho0
297	else
298	windamp=-0.1/rho0
299	endif
300	do j=JstrR,JendR
301	do i=IstrR,IendR
302	sustr(i,j)=windamp
303	enddo
304	enddo
305	# elif defined GRAV_ADJ \|\| defined OVERFLOW \|\| defined SEAMOUNT \
306	\|\| defined SHELFRONT \|\| defined SOLITON \|\| defined INNERSHELF \
307	\|\| defined RIVER \|\| defined VORTEX \|\| defined REGIONAL \
308	\|\| defined INTERNAL
309	do j=JstrR,JendR
310	do i=IstrR,IendR
311	sustr(i,j)=0.
312	enddo
313	enddo
314	# else
315	do j=JstrR,JendR
316	do i=IstrR,IendR
317	sustr(i,j)=???
318	enddo
319	enddo
320	# endif
321	!
322	! Set kinematic surface momentum flux (wind stress) component in the
323	! ETA-direction (m^2/s^2) at horizontal V-points.
324	!
325	# if defined BASIN \|\| defined CANYON_A \|\| defined CANYON_B \
326	\|\| defined EQUATOR \|\| defined GRAV_ADJ \|\| defined OVERFLOW \
327	\|\| defined REGIONAL \|\| defined RIVER \|\| defined SEAMOUNT \
328	\|\| defined SHELFRONT \|\| defined SOLITON \|\| defined UPWELLING \
329	\|\| defined VORTEX \|\| defined INTERNAL
330	do j=JstrR,JendR
331	do i=IstrR,IendR
332	svstr(i,j)=0.
333	enddo
334	enddo
335	# elif defined SED_TEST2
336	do j=JstrR,JendR
337	do i=IstrR,IendR
338	windamp=0.5+
339	& 0.5*TANH((time-user(9))/user(10))
340	svstr(i,j)= windamp*user(2)
341	enddo
342	enddo
343	# elif defined INNERSHELF
344	do j=JstrR,JendR
345	do i=IstrR,IendR
346	svstr(i,j)=-1.e-4
347	enddo
348	enddo
349	# else
350	do j=JstrR,JendR
351	do i=IstrR,IendR
352	svstr(i,j)=???
353	enddo
354	enddo
355	# endif
356	return
357	end
358	#endif /* ANA_SMFLUX */
359	#ifdef ANA_SRFLUX
360	subroutine ana_srflux_tile (Istr,Iend,Jstr,Jend)
361	!
362	!---------------------------------------------------------------------
363	! This subroutine sets kinematic surface solar shortwave radiation
364	! flux "srflx" (degC m/s) using an analytical expression.
365	!---------------------------------------------------------------------
366	!
367	implicit none
368	# include "param.h"
369	# include "grid.h"
370	# include "forces.h"
371	# include "scalars.h"
372	integer Istr,Iend,Jstr,Jend, i,j
373	!
374	# include "compute_auxiliary_bounds.h"
375	!
376	! Set kinematic surface solar shortwave radiation [degC m/s] at
377	! horizontal RHO-points.
378	!
379	# if defined EQUATOR \|\| defined INNERSHELF \|\| defined OVERFLOW \
380	\|\| defined REGIONAL \|\| defined RIVER \|\| defined SEAMOUNT \
381	\|\| defined SHELFRONT \|\| defined UPWELLING \|\| defined VORTEX \
382	\|\| defined INTERNAL
383	do j=JstrR,JendR
384	do i=IstrR,IendR
385	srflx(i,j)=0.
386	enddo
387	enddo
388	# else
389	do j=JstrR,JendR
390	do i=IstrR,IendR
391	srflx(i,j)=???
392	enddo
393	enddo
394	# endif
395	return
396	end
397	#endif /* ANA_SRFLUX */
398
399	#if defined ANA_SSH && defined ZCLIMATOLOGY
400	subroutine ana_ssh (tile)
401	implicit none
402	# include "param.h"
403	integer tile
404	# include "compute_tile_bounds.h"
405	call ana_ssh_tile (Istr,Iend,Jstr,Jend)
406	return
407	end
408	!
409	subroutine ana_ssh_tile (Istr,Iend,Jstr,Jend)
410	!
411	!---------------------------------------------------------------------
412	! This routine sets analytical sea surface height climatology [m].
413	!---------------------------------------------------------------------
414	!
415	implicit none
416	# include "param.h"
417	# include "grid.h"
418	# include "climat.h"
419	# include "scalars.h"
420	integer Istr,Iend,Jstr,Jend, i,j
421	# ifdef INTERNAL
422	real U0,omega,kwave,ETA0
423	# endif
424	!
425	# include "compute_auxiliary_bounds.h"
426	!
427	! Set sea surface height (meters).
428	!
429	# if defined REGIONAL
430	do j=JstrR,JendR
431	do i=IstrR,IendR
432	ssh(i,j)=0.
433	enddo
434	enddo
435	# elif defined INTERNAL
436	U0=0.02
437	omega=2.pi/(12.43600)
438	kwave=((omegaomega)-(f(1,1)f(1,1)))/(g*h(1,1))
439	ETA0=kwaveh(1,1)U0/omega
440	do j=JstrR,JendR
441	do i=IstrR,IendR
442	ssh(i,j)=ETA0sin(omegatime-kwave*(xr(i,j)))
443	enddo
444	enddo
445	# else
446	do j=JstrR,JendR
447	do i=IstrR,IendR
448	ssh(i,j)=???
449	enddo
450	enddo
451	# endif
452	return
453	end
454	#endif /* ANA_SSH && ZCLIMATOLOGY */
455
456	#if defined ANA_SST && defined QCORRECTION
457	subroutine ana_sst_tile (Istr,Iend,Jstr,Jend)
458	!
459	!--------------------------------------------------------------------
460	! This routine sets sea surface temperature SST[Celsius] and surface
461	! net heat flux sensitivity dQdSTT to sea surface temperature using
462	! analytical expressions. dQdSTT is usually computed in units of
463	! [Watts/m^2/degC]. It needs to be scaled to [m/s] by dividing by
464	! rho0*Cp. These forcing fields are used when the heat flux
465	! correction is activated:
466	!
467	! Q_model ~ Q + dQdSST * (T_model - SST)
468	!--------------------------------------------------------------------
469	!
470	implicit none
471	# include "param.h"
472	# include "grid.h"
473	# include "forces.h"
474	# include "scalars.h"
475	integer Istr,Iend,Jstr,Jend, i,j
476	# if defined EQUATOR
477	real y1,y2,sst1,sst2
478	# endif
479	!
480	# include "compute_auxiliary_bounds.h"
481	!
482	# if defined EQUATOR
483	! SST = 25C and lineraly decreases to 10C 1200km from the Equator.
484	y1=1200.E+3
485	y2=1500.E+3
486	sst1=10.
487	sst2=25.
488	do j=JstrR,JendR
489	do i=IstrR,IendR
490	sst(i,j)=sst1
491	if ((yr(i,j).gt.-y1).and.(yr(i,j).lt.y1)) then
492	sst(i,j)=sst2
493	else
494	if ((yr(i,j).gt.-y2).and.(yr(i,j).lt.-y1)) then
495	sst(i,j)=((sst2-sst1)yr(i,j)-sst1y1+y2*sst2)/(y2-y1)
496	endif
497	if ((yr(i,j).gt.y1).and.(yr(i,j).lt.y2)) then
498	sst(i,j)=((sst2-sst1)yr(i,j)+sst1y1-y2*sst2)/(y1-y2)
499	endif
500	endif
501	dqdt(i,j)=-50.0/(rho0*Cp)
502	enddo
503	enddo
504	# else
505	do j=JstrR,JendR
506	do i=IstrR,IendR
507	sst(i,j)=???
508	dqdt(i,j)=???
509	enddo
510	enddo
511	# endif
512	return
513	end
514	#endif /* ANA_SST && QCORRECTION */
515	#if defined SALINITY && defined SFLX_CORR && defined ANA_SSS
516	subroutine ana_sss_tile (Istr,Iend,Jstr,Jend)
517	!
518	!--------------------------------------------------------------------
519	! This routine sets sea surface salinity SSS[PSU] using
520	! analytical expressions. This forcing field is used when the
521	! salt flux correction is activated:
522	!
523	! SSSFLX_model ~ SSS(E-P) + CST (SSS_model - SSS)
524	!
525	! we use DQDSST for CST....
526	!
527	!--------------------------------------------------------------------
528	!
529	implicit none
530	# include "param.h"
531	# include "grid.h"
532	# include "forces.h"
533	# include "scalars.h"
534	integer Istr,Iend,Jstr,Jend, i,j
535	!
536	# include "compute_auxiliary_bounds.h"
537	!
538	do j=JstrR,JendR
539	do i=IstrR,IendR
540	sss(i,j)=???
541	enddo
542	enddo
543	return
544	end
545	#endif /* SALINITY && SFLX_CORR && ANA_SSS */
546	#if defined ANA_STFLUX \|\| defined ANA_SSFLUX
547	subroutine ana_stflux_tile (Istr,Iend,Jstr,Jend, itrc)
548	!
549	!--------------------------------------------------------------------
550	! This routine sets kinematic surface flux of tracer type variables
551	! "stflx" (tracer units m/s) using analytical expressions.
552	!
553	! On Input:
554	! itrc Tracer type array index.
555	!--------------------------------------------------------------------
556	!
557	implicit none
558	# include "param.h"
559	# include "grid.h"
560	# include "forces.h"
561	# include "scalars.h"
562	integer itrc, Istr,Iend,Jstr,Jend, i,j
563	c
564	# include "compute_auxiliary_bounds.h"
565	c
566	if (itrc.eq.itemp) then
567	!
568	! Set kinematic surface heat flux [degC m/s] at horizontal
569	! RHO-points.
570	!
571	#if defined BASIN \|\| defined CANYON_A \|\| defined CANYON_B \
572	\|\| defined EQUATOR \|\| defined GRAV_ADJ \|\| defined INNERSHELF \
573	\|\| defined OVERFLOW \|\| defined REGIONAL \|\| defined RIVER \
574	\|\| defined SEAMOUNT \|\| defined SHELFRONT \|\| defined UPWELLING \
575	\|\| defined VORTEX \|\| defined INTERNAL
576
577	do j=JstrR,JendR
578	do i=IstrR,IendR
579	stflx(i,j,itemp)=0.
580	enddo
581	enddo
582	#else
583	do j=JstrR,JendR
584	do i=IstrR,IendR
585	stflx(i,j,itemp)=???
586	enddo
587	enddo
588	#endif
589	#ifdef SALINITY
590	elseif (itrc.eq.isalt) then
591	!
592	! Set kinematic surface freshwater flux (m/s) at horizontal
593	! RHO-points, scaling by surface salinity is done in STEP3D.
594	!
595	# if defined EQUATOR \|\| defined INNERSHELF \|\| defined REGIONAL \
596	\|\| defined RIVER \|\| defined SEAMOUNT \|\| defined SHELFRONT \
597	\|\| defined UPWELLING
598
599	do j=JstrR,JendR
600	do i=IstrR,IendR
601	stflx(i,j,isalt)=0.
602	enddo
603	enddo
604	# else
605	do j=JstrR,JendR
606	do i=IstrR,IendR
607	stflx(i,j,isalt)=???
608	enddo
609	enddo
610	# endif
611	#endif /* SALINITY */
612	else
613	!
614	! Set kinematic surface flux of additional tracers, if any.
615	!
616	do j=JstrR,JendR
617	do i=IstrR,IendR
618	stflx(i,j,itrc)=0.
619	enddo
620	enddo
621	endif
622	return
623	end
624	#endif /* ANA_STFLUX \|\| ANA_SSFLUX */
625	!---------------------------------------------------------------------
626	#if defined TCLIMATOLOGY
627	subroutine ana_tclima (tile)
628	implicit none
629	# include"param.h"
630	integer tile
631	#include "compute_tile_bounds.h"
632	call ana_tclima_tile (Istr,Iend,Jstr,Jend)
633	return
634	end
635	subroutine ana_tclima_tile (Istr,Iend,Jstr,Jend)
636	!
637	!---------------------------------------------------------------------
638	! This routine sets analytical ACTIVE (T&S) tracer climatology fields.
639	!---------------------------------------------------------------------
640	!
641	implicit none
642	# include "param.h"
643	# include "grid.h"
644	# include "climat.h"
645	# include "ocean3d.h"
646	# include "scalars.h"
647	# include "sediment.h"
648	integer Istr,Iend,Jstr,Jend, i,j,k, itrc
649	real cff,cff1
650	!
651	# include "compute_auxiliary_bounds.h"
652	!
653	! Set climatology fields for tracer type variables.
654	!---------------------------------------------------------------------
655	!
656	# ifdef ANA_TCLIMA
657	do k=1,N
658	do j=JstrR,JendR
659	do i=IstrR,IendR
660	tclm(i,j,k,itemp)=???
661	# ifdef SALINITY
662	tclm(i,j,k,isalt)=???
663	# endif /* SALINITY */
664	enddo
665	enddo
666	enddo
667	# endif
668
669	# ifdef BIOLOGY
670	# define temp cff
671	# define SiO4 cff1
672	do k=1,N
673	do j=JstrR,JendR
674	do i=IstrR,IendR
675	# ifdef ANA_TCLIMA
676	temp=t(i,j,k,1,itemp)
677	if (temp.lt.8.) then
678	SiO4=30.
679	elseif (temp.ge.8. .and. temp.le.11.) then
680	SiO4=30.-((temp-8.)*(20./3.))
681	elseif (temp.gt.11. .and. temp.le.13.) then
682	SiO4=10.-((temp-11.)*(8./2.))
683	elseif (temp.gt.13. .and. temp.le.16.) then
684	SiO4=2.-((temp-13.)*(2./3.))
685	elseif (temp.gt.16.) then
686	SiO4=0.
687	endif
688	tclm(i,j,k,iNO3_)=1.67+0.5873SiO4+0.0144SiO4**2
689	& +0.0003099SiO4*3
690	# ifdef PISCES
691	tclm(i,j,k,iDIC_)=2150.
692	tclm(i,j,k,iTAL_)=2350.
693	tclm(i,j,k,iOXY_)=200.
694	tclm(i,j,k,iCAL_)=0.01
695	tclm(i,j,k,iPO4_)=tclm(i,j,k,iNO3_)/16.
696	tclm(i,j,k,iPOC_)=0.01
697	tclm(i,j,k,iSIL_)=91.51
698	tclm(i,j,k,iPHY_)=0.01
699	tclm(i,j,k,iZOO_)=0.01
700	tclm(i,j,k,iDOC_)=5.
701	tclm(i,j,k,iDIA_)=0.01
702	tclm(i,j,k,iMES_)=0.01
703	tclm(i,j,k,iBSI_)=1.5e-3
704	tclm(i,j,k,iFER_)=6.e-4
705	tclm(i,j,k,iBFE_)=1.e-2*5.e-6
706	tclm(i,j,k,iGOC_)=0.01
707	tclm(i,j,k,iSFE_)=1.e-2*5.e-6
708	tclm(i,j,k,iDFE_)=1.e-2*5.e-6
709	tclm(i,j,k,iDSI_)=1.e-2*0.15
710	tclm(i,j,k,iNFE_)=1.e-2*5.e-6
711	tclm(i,j,k,iNCH_)=1.e-2*12./55.
712	tclm(i,j,k,iDCH_)=1.e-2*12./55.
713	tclm(i,j,k,iNH4_)=1.e-2
714	# elif defined BIO_NChlPZD
715	tclm(i,j,k,iChla)=0.08
716	tclm(i,j,k,iPhy1)=0.1
717	tclm(i,j,k,iZoo1)=0.06
718	tclm(i,j,k,iDet1)=0.02
719	# elif defined BIO_N2ChlPZD2
720	tclm(i,j,k,iNH4_)=0.1
721	tclm(i,j,k,iChla)=0.08
722	tclm(i,j,k,iPhy1)=0.06
723	tclm(i,j,k,iZoo1)=0.04
724	tclm(i,j,k,iDet1)=0.02
725	tclm(i,j,k,iDet2)=0.02
726	# elif defined BIO_N2P2Z2D2
727	tclm(i,j,k,iNH4_)=0.1
728	tclm(i,j,k,iPhy1)=0.04
729	tclm(i,j,k,iPhy2)=0.06
730	tclm(i,j,k,iZoo1)=0.04
731	tclm(i,j,k,iZoo2)=0.04
732	tclm(i,j,k,iDet1)=0.02
733	tclm(i,j,k,iDet2)=0.02
734	# endif
735	# else
736	if (.not.got_tclm(iNO3_)) then
737	temp=t(i,j,k,1,itemp)
738	if (temp.lt.8.) then
739	SiO4=30.
740	elseif (temp.ge.8. .and. temp.le.11.) then
741	SiO4=30.-((temp-8.)*(20./3.))
742	elseif (temp.gt.11. .and. temp.le.13.) then
743	SiO4=10.-((temp-11.)*(8./2.))
744	elseif (temp.gt.13. .and. temp.le.16.) then
745	SiO4=2.-((temp-13.)*(2./3.))
746	elseif (temp.gt.16.) then
747	SiO4=0.
748	endif
749	tclm(i,j,k,iNO3_)=1.67+0.5873SiO4+0.0144SiO4**2
750	& +0.0003099SiO4*3
751	endif
752	# ifdef PISCES
753	if (.not.got_tclm(iDIC_)) tclm(i,j,k,iDIC_)=2150.
754	if (.not.got_tclm(iTAL_)) tclm(i,j,k,iTAL_)=2350.
755	if (.not.got_tclm(iOXY_)) tclm(i,j,k,iOXY_)=200.
756	if (.not.got_tclm(iCAL_)) tclm(i,j,k,iCAL_)=0.01
757	if (.not.got_tclm(iPO4_)) then
758	temp=t(i,j,k,1,itemp)
759	if (temp.lt.8.) then
760	SiO4=30.
761	elseif (temp.ge.8. .and. temp.le.11.) then
762	SiO4=30.-((temp-8.)*(20./3.))
763	elseif (temp.gt.11. .and. temp.le.13.) then
764	SiO4=10.-((temp-11.)*(8./2.))
765	elseif (temp.gt.13. .and. temp.le.16.) then
766	SiO4=2.-((temp-13.)*(2./3.))
767	elseif (temp.gt.16.) then
768	SiO4=0.
769	endif
770	tclm(i,j,k,iPO4_)=(1.67+0.5873SiO4+0.0144SiO4**2
771	& +0.0003099SiO4*3)/16.
772	endif
773	if (.not.got_tclm(iPOC_)) tclm(i,j,k,iPOC_)=0.01
774	if (.not.got_tclm(iSIL_)) tclm(i,j,k,iSIL_)=91.51
775	if (.not.got_tclm(iPHY_)) tclm(i,j,k,iPHY_)=0.01
776	if (.not.got_tclm(iZOO_)) tclm(i,j,k,iZOO_)=0.01
777	if (.not.got_tclm(iDOC_)) tclm(i,j,k,iDOC_)=5.
778	if (.not.got_tclm(iDIA_)) tclm(i,j,k,iDIA_)=0.01
779	if (.not.got_tclm(iMES_)) tclm(i,j,k,iMES_)=0.01
780	if (.not.got_tclm(iBSI_)) tclm(i,j,k,iBSI_)=0.0015
781	if (.not.got_tclm(iFER_)) tclm(i,j,k,iFER_)=6.e-4
782	if (.not.got_tclm(iBFE_)) tclm(i,j,k,iBFE_)=5.e-8
783	if (.not.got_tclm(iGOC_)) tclm(i,j,k,iGOC_)=0.01
784	if (.not.got_tclm(iSFE_)) tclm(i,j,k,iSFE_)=5.e-8
785	if (.not.got_tclm(iDFE_)) tclm(i,j,k,iDFE_)=5.e-8
786	if (.not.got_tclm(iDSI_)) tclm(i,j,k,iDSI_)=0.0015
787	if (.not.got_tclm(iNFE_)) tclm(i,j,k,iNFE_)=5.e-8
788	if (.not.got_tclm(iNCH_)) tclm(i,j,k,iNCH_)=1.e-2*12./55.
789	if (.not.got_tclm(iDCH_)) tclm(i,j,k,iDCH_)=1.e-2*12./55.
790	if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.01
791	# elif defined BIO_NChlPZD
792	if (.not.got_tclm(iChla)) tclm(i,j,k,iChla)=0.08
793	if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.1
794	if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.06
795	if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
796	# elif defined BIO_N2ChlPZD2
797	if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.1
798	if (.not.got_tclm(iChla)) tclm(i,j,k,iChla)=0.08
799	if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.06
800	if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.04
801	if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
802	if (.not.got_tclm(iDet2)) tclm(i,j,k,iDet2)=0.02
803	# elif defined BIO_N2P2Z2D2
804	if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.1
805	if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.04
806	if (.not.got_tclm(iPhy2)) tclm(i,j,k,iPhy2)=0.06
807	if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.04
808	if (.not.got_tclm(iZoo2)) tclm(i,j,k,iZoo2)=0.04
809	if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
810	if (.not.got_tclm(iDet2)) tclm(i,j,k,iDet2)=0.02
811	# endif
812	# endif /* ANA_TCLIMA */
813	enddo
814	enddo
815	enddo
816	# undef SiO4
817	# undef temp
818	# endif /* BIOLOGY */
819
820	# ifdef SEDIMENT
821	do k=1,N
822	do j=JstrR,JendR
823	do i=IstrR,IendR
824	# ifdef ANA_TCLIMA
825	tclm(i,j,k,isand)=Csed(1)
826	tclm(i,j,k,isilt)=Csed(2)
827	# else
828	if (.not.got_tclm(isand)) then
829	tclm(i,j,k,isand)=0. !Csed(1)
830	endif
831	if (.not.got_tclm(isilt)) then
832	tclm(i,j,k,isilt)=0. !Csed(2)
833	endif
834	# endif
835	enddo
836	enddo
837	enddo
838	# endif /* SEDIMENT */
839
840	# ifdef PASSIVE_TRACER
841	do k=1,N
842	do j=JstrR,JendR
843	do i=IstrR,IendR
844	# ifdef ANA_TCLIMA
845	tclm(i,j,k,itpas)=0.0
846	# else
847	if (.not.got_tclm(itpas)) then
848	tclm(i,j,k,itpas)=0.0
849	endif
850	# endif
851	enddo
852	enddo
853	enddo
854	# endif
855
856	# if defined EW_PERIODIC \|\| defined NS_PERIODIC \|\| defined MPI
857	do itrc=1,NT
858	# ifndef ANA_TCLIMA
859	if (.not.got_tclm(itrc)) then
860	# endif
861	call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
862	& tclm(START_2D_ARRAY,1,itrc))
863	# ifndef ANA_TCLIMA
864	endif
865	# endif
866	enddo
867	# endif
868
869	return
870	end
871	#endif /* TCLIMATOLOGY */
872	!
873	!====================================================================
874	! subroutine ana_uclima
875	!====================================================================
876	!
877	#if defined ANA_M2CLIMA && defined M2CLIMATOLOGY \|\|\
878	(defined ANA_M3CLIMA && defined M3CLIMATOLOGY)
879	subroutine ana_uclima (tile)
880	implicit none
881	# include "param.h"
882	integer tile
883	# include "compute_tile_bounds.h"
884	call ana_uclima_tile (Istr,Iend,Jstr,Jend)
885	return
886	end
887	!---------------------------------------------------------------------
888	!
889	subroutine ana_uclima_tile (Istr,Iend,Jstr,Jend)
890	!
891	!---------------------------------------------------------------------
892	! This routine sets analytical momentum climatology fields.
893	!---------------------------------------------------------------------
894	!
895	implicit none
896	# include "param.h"
897	# include "grid.h"
898	# include "climat.h"
899	# include "scalars.h"
900	integer Istr,Iend,Jstr,Jend, i,j,k
901	# ifdef INTERNAL
902	real U0,omega,kwave,V0
903	# endif
904	!
905	# include "compute_auxiliary_bounds.h"
906	!
907	# ifdef EW_PERIODIC
908	# define IU_RANGE Istr,Iend
909	# define IV_RANGE Istr,Iend
910	# else
911	# define IU_RANGE Istr,IendR
912	# define IV_RANGE IstrR,IendR
913	# endif
914
915	# ifdef NS_PERIODIC
916	# define JU_RANGE Jstr,Jend
917	# define JV_RANGE Jstr,Jend
918	# else
919	# define JU_RANGE JstrR,JendR
920	# define JV_RANGE Jstr,JendR
921	# endif
922	!
923	# if defined ANA_M2CLIMA && defined M2CLIMATOLOGY
924	# if defined REGIONAL
925	do j=JstrR,JendR
926	do i=IstrR,IendR
927	ubclm(i,j)=0.
928	vbclm(i,j)=0.
929	enddo
930	enddo
931	# elif defined INTERNAL
932	U0=0.02
933	omega=2.pi/(12.43600)
934	kwave=sqrt(((omegaomega)-(f(1,1)f(1,1)))/(g*h(1,1)))
935	V0=f(1,1)*U0/omega
936	do j=JU_RANGE
937	do i=IU_RANGE
938	ubclm(i,j)=U0sin(omegatime-kwave0.5(xr(i,j)+xr(i-1,j)))
939	enddo
940	enddo
941	do j=JV_RANGE
942	do i=IV_RANGE
943	vbclm(i,j)=V0cos(omegatime-kwave0.5(xr(i,j)+xr(i,j-1)))
944	enddo
945	enddo
946	# else
947	do j=JstrR,JendR
948	do i=IstrR,IendR
949	ubclm(i,j)=???
950	vbclm(i,j)=???
951	enddo
952	enddo
953	# endif
954	# if defined EW_PERIODIC \|\| defined NS_PERIODIC \|\| defined MPI
955	call exchange_u2d_tile (Istr,Iend,Jstr,Jend, ubclm)
956	call exchange_v2d_tile (Istr,Iend,Jstr,Jend, vbclm)
957	# endif
958	# endif
959	# if defined ANA_M3CLIMA && defined M3CLIMATOLOGY && defined SOLVE3D
960	# if defined REGIONAL
961	do k=1,N
962	do j=JstrR,JendR
963	do i=IstrR,IendR
964	uclm(i,j,k)=0.
965	vclm(i,j,k)=0.
966	enddo
967	enddo
968	enddo
969	# else
970	do k=1,N
971	do j=JstrR,JendR
972	do i=IstrR,IendR
973	uclm(i,j,k)=???
974	vclm(i,j,k)=???
975	enddo
976	enddo
977	enddo
978	# endif
979	# endif
980	# undef IU_RANGE
981	# undef JU_RANGE
982	# undef IV_RANGE
983	# undef JV_RANGE
984	return
985	end
986	#endif /* ANA_M2CLIMA && M2CLIMATOLOGY \|\| (ANA_M3CLIMA && M3CLIMATOLOGY) */
987	!
988	!====================================================================
989	! subroutine ana_wwave
990	!====================================================================
991	!
992	#if defined ANA_WWAVE && defined BBL
993	subroutine ana_wwave (tile)
994	implicit none
995	# include "param.h"
996	integer tile
997	# include "compute_tile_bounds.h"
998	call ana_wwave_tile (Istr,Iend,Jstr,Jend)
999	return
1000	end
1001	!
1002	subroutine ana_wwave_tile (Istr,Iend,Jstr,Jend)
1003	!
1004	!---------------------------------------------------------------------
1005	! This routine sets wind induced wave amplitude, direction
1006	! and period used in the bottom boundary layer formulation.
1007	!---------------------------------------------------------------------
1008	!
1009	implicit none
1010	# include "param.h"
1011	# include "grid.h"
1012	# include "forces.h"
1013	# include "scalars.h"
1014	integer Istr,Iend,Jstr,Jend, i,j
1015	!
1016	# include "compute_auxiliary_bounds.h"
1017	!
1018	! Set wind induced wave amplitude (m), direction (radians) and
1019	! period (s) at RHO-points.
1020	!
1021	# if defined SED_TEST2 \|\| defined REGIONAL
1022	do j=JstrR,JendR
1023	do i=IstrR,IendR
1024	Awave(i,j)=1.0
1025	Dwave(i,j)=270.*deg2rad
1026	Pwave(i,j)=10.
1027	enddo
1028	enddo
1029	# else
1030	ANA_WWAVE: no values provided for AWAVE, DWAVE, and PWAVE.
1031	do j=JstrR,JendR
1032	do i=IstrR,IendR
1033	Awave(i,j)=???
1034	Dwave(i,j)=???
1035	Pwave(i,j)=???
1036	enddo
1037	enddo
1038	# endif
1039	return
1040	end
1041	#endif /* ANA_WWAVE && BBL */
1042
1043	#if defined SEDIMENT
1044	subroutine ana_sediment (tile)
1045	implicit none
1046	# include "param.h"
1047	integer tile
1048	# include "compute_tile_bounds.h"
1049	call ana_sediment_tile (Istr,Iend,Jstr,Jend)
1050	return
1051	end
1052	!
1053	subroutine ana_sediment_tile (Istr,Iend,Jstr,Jend)
1054	!
1055	!---------------------------------------------------------------------
1056	! This routine sets sediment ripple and bed parameters
1057	!---------------------------------------------------------------------
1058	!
1059	implicit none
1060	# include "param.h"
1061	# include "grid.h"
1062	# include "scalars.h"
1063	# include "sediment.h"
1064	# include "bbl.h"
1065	integer Istr,Iend,Jstr,Jend, i,j, ilay, itrc
1066	!
1067	#undef DEBUG
1068	#ifdef DEBUG
1069	integer ick,jck
1070	parameter (ick=60, jck=35)
1071	#endif
1072	# include "compute_auxiliary_bounds.h"
1073	# if defined SED_TEST2 \|\| defined REGIONAL
1074	do j=JstrR,JendR
1075	do i=IstrR,IendR
1076	#ifdef BBL
1077	# ifndef ANA_SEDIMENT
1078	if(.not.got_inibed(1)) then
1079	# endif
1080	Hripple(i,j)=Hrip ! initial ripple height [m] from .in file
1081	# ifndef ANA_SEDIMENT
1082	endif
1083	if(.not.got_inibed(2)) then
1084	# endif
1085	Lripple(i,j)=Lrip ! initial ripple length [m] from .in file
1086	# ifndef ANA_SEDIMENT
1087	endif
1088	# endif
1089	#endif
1090	do ilay=1,NLAY
1091	#ifndef ANA_SEDIMENT
1092	if(.not.got_inised(1)) then
1093	#endif
1094	bed_thick(i,j,ilay)=Bthk(ilay)
1095	#ifndef ANA_SEDIMENT
1096	endif
1097	if(.not.got_inised(2)) then
1098	#endif
1099	bed_poros(i,j,ilay)=Bpor(ilay)
1100	#ifndef ANA_SEDIMENT
1101	endif
1102	if(.not.got_inised(3)) then
1103	#endif
1104	do itrc=1,NST
1105	bed_frac(i,j,ilay,itrc)=Bfr(ilay,itrc)
1106	enddo
1107	#ifndef ANA_SEDIMENT
1108	endif
1109	#endif
1110	enddo
1111
1112	#ifdef DEBUG
1113	if(j.eq.jck.and.i.eq.ick) then
1114	write(6,) '******* ANA_SEDIMENT *********'
1115	do itrc=1,NST
1116	write(6,*) 'Sd(itrc)',Sd(itrc)
1117	do ilay=1,NLAY
1118	write(6,*) 'i,j,ilay,itrc',i,j,ilay,itrc
1119	write(6,*) 'bed_frac',bed_frac(i,j,ilay,itrc)
1120	enddo
1121	enddo
1122	do ilay=1,NLAY
1123	write(6,*) 'bed_thick',bed_thick(i,j,ilay)
1124	write(6,*) 'bed_por',bed_poros(i,j,ilay)
1125	enddo
1126	endif
1127	#endif
1128	enddo
1129	enddo
1130	# else
1131	ana_sediment: no values provided
1132	do j=JstrR,JendR
1133	do i=IstrR,IendR
1134	do ilay=1,NLAY
1135	bed_thick(i,j,ilay)=???
1136	bed_poros(i,j,ilay)=???
1137	do itrc=1,NST
1138	bed_frac(i,j,ilay,itrc)=???
1139	enddo
1140	enddo
1141	enddo
1142	enddo
1143	# endif
1144	return
1145	end
1146	#endif /* defined SEDIMENT */
1147
1148	#if defined PSOURCE && defined ANA_PSOURCE && defined SOLVE3D
1149	!
1150	!-----------------------------------------------------------
1151	! Set analytical tracer and mass point sources and sinks
1152	!-----------------------------------------------------------
1153	!
1154	subroutine ana_psource_tile (Istr,Iend,Jstr,Jend)
1155	implicit none
1156	# include "param.h"
1157	# include "scalars.h"
1158	# include "sources.h"
1159	# include "ocean3d.h"
1160	# include "grid.h"
1161	!
1162	integer is, k, Istr,Iend,Jstr,Jend, i,j
1163	real cff, cff1, cff2, ramp, Hs
1164	# include "compute_auxiliary_bounds.h"
1165
1166	if (iic.eq.ntstart) then
1167
1168	!
1169	! Set-up nondimensional shape Qshape, must add to unity
1170	!
1171	# if defined RIVER
1172	# define EXP_SHAPE
1173	# ifdef CST_SHAPE
1174	cff=1./float(N)
1175	do k=1,N ! Uniform vertical
1176	do is=1,Nsrc ! distribution
1177	Qshape(is,k)=cff
1178	enddo
1179	enddo
1180	# elif defined EXP_SHAPE
1181	do is=1,Nsrc ! Exponential vertical
1182	Hs=h(Isrc(is),Jsrc(is)) ! distribution
1183	cff=5. ! Hs/z0 (z0 surface layer depth)
1184	cff1=cff/(1-exp(-cff))
1185	cff2=0.
1186	do k=1,N
1187	Qshape(is,k)=cff1exp(z_r(Isrc(is),Jsrc(is),k)cff/Hs)*
1188	& (z_w(Isrc(is),Jsrc(is),k)-z_w(Isrc(is),Jsrc(is),k-1))/Hs
1189	cff2=cff2+Qshape(is,k)
1190	enddo
1191	do k=1,N
1192	Qshape(is,k)=Qshape(is,k)/cff2
1193	enddo
1194	enddo
1195	# elif defined AL_SHAPE
1196	do is=1,Nsrc ! Set-up nondimensional shape
1197	do k=1,10
1198	Qshape(is,k)=0.0
1199	enddo
1200	do k=11,14
1201	Qshape(is,k)=0.05
1202	enddo
1203	Qshape(is,15)=0.1 ! These most add to unity!
1204	Qshape(is,16)=0.1 ! These most add to unity!
1205	Qshape(is,17)=0.1 ! These most add to unity!
1206	Qshape(is,18)=0.1 ! These most add to unity!
1207	Qshape(is,19)=0.2
1208	Qshape(is,20)=0.2
1209	enddo
1210	# endif
1211
1212	# elif defined REGIONAL
1213	do is=1,Nsrc ! Exponential vertical
1214	# ifdef MPI
1215	i=Isrc_mpi(is,mynode)
1216	j=Jsrc_mpi(is,mynode)
1217	# else
1218	i=Isrc(is)
1219	j=Jsrc(is)
1220	# endif
1221	Hs=h(i,j) ! distribution
1222	cff=5. ! Hs/z0 (z0 surface layer depth)
1223	cff1=cff/(1-exp(-cff))
1224	cff2=0.
1225	do k=1,N
1226	Qshape(is,k)=cff1exp(z_r(Isrc(is),Jsrc(is),k)cff/Hs)*
1227	& (z_w(Isrc(is),Jsrc(is),k)-z_w(Isrc(is),Jsrc(is),k-1))/Hs
1228	cff2=cff2+Qshape(is,k)
1229	enddo
1230	do k=1,N
1231	Qshape(is,k)=Qshape(is,k)/cff2
1232	enddo
1233	enddo
1234
1235	# else
1236	ERROR ### CPP-key 'ANA_PSOURCE' is defined, but no code
1237	ERROR ### is provided to set up Dsrc, Isrc, Jsrc, Lsrc.
1238	# endif /* REGIONAL */
1239
1240	endif ! iic.eq.ntstart
1241
1242	!
1243	! Set-up vertically integrated mass transport [m3/s] of point
1244	! sources (these may be time-dependent; positive in the positive U-
1245	! or V-direction and vice-versa) and vertically distribute them
1246	! according to mass transport profile chosen above.
1247	!
1248	# if defined RIVER
1249	ramp=1 !TANH(dtsec2dayfloat(iic-ntstart))
1250	do is=1,Nsrc
1251	Qbar(is)=ramp*Qbar(is)
1252	enddo
1253	# elif defined REGIONAL
1254	ramp=1 !TANH(dtsec2dayfloat(iic-ntstart))
1255	do is=1,Nsrc
1256	Qbar(is)=ramp*Qbar(is)
1257	enddo
1258	# else
1259	ERROR ### CPP-key 'ANA_PSOURCE' is defined, but no code
1260	ERROR ### is provided to set up Qbar(is) analytically.
1261	# endif
1262	do is=1,Nsrc
1263	do k=1,N
1264	Qsrc(is,k)=Qbar(is)*Qshape(is,k)
1265	enddo
1266	enddo
1267	!
1268	! Set-up tracer (tracer units) point Sources/Sinks.
1269	!
1270	# if defined RIVER
1271	do k=1,N
1272	do is=1,Nsrc
1273	Tsrc(is,k,itemp)=Tsrc0(is,itemp)
1274	! Tsrc(is,k,itemp)=4.+10.*exp(z_r(Isrc(is),Jsrc(is),k)/50.)
1275	Tsrc(is,k,isalt)=Tsrc0(is,isalt)
1276	# if defined PASSIVE_TRACER
1277	Tsrc(is,k,ipas)=Tsrc0(is,itpas)
1278	# endif
1279	enddo
1280	enddo
1281	# elif defined REGIONAL
1282	do k=1,N
1283	do is=1,Nsrc
1284	Tsrc(is,k,itemp)=Tsrc0(is,itemp)
1285	Tsrc(is,k,isalt)=Tsrc0(is,isalt)
1286	# if defined PASSIVE_TRACER
1287	Tsrc(is,k,itpas)=Tsrc0(is,itpas)
1288	# endif
1289	# if defined BIOLOGY
1290	Tsrc(is,k,iNO3_)=Tsrc0(is,iNO3_)
1291	# endif
1292	enddo
1293	enddo
1294	# else
1295	ERROR ### CPP-key 'ANA_PSOURCE' is defined, but no code
1296	ERROR ### is provided to set up Tsrc(is) analytically.
1297	# endif
1298	return
1299	end
1300	#endif /* PSOURCE && ANA_PSOURCE SOLVE3D*/
1301
1302	#ifdef ANA_BRY
1303	!
1304	!---------------------------------------------------------------------
1305	! Set analytical boundary forcing
1306	!---------------------------------------------------------------------
1307	!
1308	subroutine ana_bry (tile)
1309	implicit none
1310	# include "param.h"
1311	integer tile
1312	# include "compute_tile_bounds.h"
1313	call ana_bry_tile (Istr,Iend,Jstr,Jend)
1314	return
1315	end
1316	!
1317	subroutine ana_bry_tile (Istr,Iend,Jstr,Jend)
1318	implicit none
1319	integer Istr,Iend,Jstr,Jend, i,j,k, itrc
1320	# include "param.h"
1321	# include "boundary.h"
1322	!
1323	# include "compute_auxiliary_bounds.h"
1324	!
1325	# ifdef OBC_WEST
1326	if (WESTERN_EDGE) then
1327	# ifdef Z_FRC_BRY
1328	do j=JstrR,JendR
1329	zetabry_west(j)=0.
1330	enddo
1331	# endif
1332	# ifdef M2_FRC_BRY
1333	do j=JstrR,JendR
1334	ubarbry_west(j)=0.
1335	vbarbry_west(j)=0.
1336	enddo
1337	# endif
1338	# ifdef SOLVE3D && (defined M3_FRC_BRY \|\| defined T_FRC_BRY)
1339	do k=1,N
1340	do j=JstrR,JendR
1341	# ifdef M3_FRC_BRY
1342	ubry_west(j,k)=0.
1343	vbry_west(j,k)=0.
1344	# endif
1345	# ifdef T_FRC_BRY
1346	do itrc=1,NT
1347	tbry_west(j,k,itrc)=0.
1348	enddo
1349	# endif
1350	enddo
1351	enddo
1352	# endif /* SOLVE3D && (M3_FRC_BRY \|\| T_FRC_BRY)*/
1353	endif
1354	# endif /* OBC_WEST */
1355	!
1356	# ifdef OBC_EAST
1357	if (EASTERN_EDGE) then
1358	# ifdef Z_FRC_BRY
1359	do j=JstrR,JendR
1360	zetabry_east(j)=0.
1361	enddo
1362	# endif
1363	# ifdef M2_FRC_BRY
1364	do j=JstrR,JendR
1365	ubarbry_east(j)=0.
1366	vbarbry_east(j)=0.
1367	enddo
1368	# endif
1369	# ifdef SOLVE3D && (defined M3_FRC_BRY \|\| defined T_FRC_BRY)
1370	do k=1,N
1371	do j=JstrR,JendR
1372	# ifdef M3_FRC_BRY
1373	ubry_east(j,k)=0.
1374	vbry_east(j,k)=0.
1375	# endif
1376	# ifdef T_FRC_BRY
1377	do itrc=1,NT
1378	tbry_east(j,k,itrc)=0.
1379	enddo
1380	# endif
1381	enddo
1382	enddo
1383	# endif /* SOLVE3D && (M3_FRC_BRY \|\| T_FRC_BRY)*/
1384	endif
1385	# endif /* OBC_EAST */
1386	!
1387	# ifdef OBC_SOUTH
1388	if (SOUTHERN_EDGE) then
1389	# ifdef Z_FRC_BRY
1390	do i=IstrR,IendR
1391	zetabry_south(i)=0.
1392	enddo
1393	# endif
1394	# ifdef M2_FRC_BRY
1395	do i=IstrR,IendR
1396	ubarbry_south(i)=0.
1397	vbarbry_south(i)=0.
1398	enddo
1399	# endif
1400	# ifdef SOLVE3D && (defined M3_FRC_BRY \|\| defined T_FRC_BRY)
1401	do k=1,N
1402	do i=IstrR,IendR
1403	# ifdef M3_FRC_BRY
1404	ubry_south(i,k)=0.
1405	vbry_south(i,k)=0.
1406	# endif
1407	# ifdef T_FRC_BRY
1408	do itrc=1,NT
1409	tbry_south(i,k,itrc)=0.
1410	enddo
1411	# endif
1412	enddo
1413	enddo
1414	# endif /* SOLVE3D && (M3_FRC_BRY \|\| T_FRC_BRY)*/
1415	endif
1416	# endif /* OBC_SOUTH */
1417	!
1418	# ifdef OBC_NORTH
1419	if (NORTHERN_EDGE) then
1420	# ifdef Z_FRC_BRY
1421	do i=IstrR,IendR
1422	zetabry_north(i)=0.
1423	enddo
1424	# endif
1425	# ifdef M2_FRC_BRY
1426	do i=IstrR,IendR
1427	ubarbry_north(i)=0.
1428	vbarbry_north(i)=0.
1429	enddo
1430	# endif
1431	# ifdef SOLVE3D && (defined M3_FRC_BRY \|\| defined T_FRC_BRY)
1432	do k=1,N
1433	do i=IstrR,IendR
1434	# ifdef M3_FRC_BRY
1435	ubry_north(i,k)=0.
1436	vbry_north(i,k)=0.
1437	# endif
1438	# ifdef T_FRC_BRY
1439	do itrc=1,NT
1440	tbry_north(i,k,itrc)=0.
1441	enddo
1442	# endif
1443	enddo
1444	enddo
1445	# endif /* SOLVE3D && (M3_FRC_BRY \|\| T_FRC_BRY)*/
1446	endif
1447	# endif /* OBC_NORTH */
1448	!
1449	return
1450	end
1451	#endif /* defined ANA_BRY */
1452
1453	#if defined BIOLOGY && defined T_FRC_BRY
1454	!
1455	!---------------------------------------------------------------------
1456	! Set analytical boundary forcing for biological tracers
1457	!---------------------------------------------------------------------
1458	!
1459	subroutine ana_bry_bio (tile)
1460	implicit none
1461	integer tile
1462	# include "param.h"
1463	# include "compute_tile_bounds.h"
1464	call ana_bry_bio_tile (Istr,Iend,Jstr,Jend)
1465	return
1466	end
1467	!
1468	subroutine ana_bry_bio_tile (Istr,Iend,Jstr,Jend)
1469	implicit none
1470	integer Istr,Iend,Jstr,Jend, i,j,k, itrc
1471	real xno3,temp,SiO4
1472	# include "param.h"
1473	# include "boundary.h"
1474	# include "scalars.h"
1475	!
1476	# include "compute_auxiliary_bounds.h"
1477	!
1478	# ifdef PISCES
1479	# ifdef OBC_WEST
1480	if (WESTERN_EDGE) then
1481	do k=1,N
1482	do j=JstrR,JendR
1483	temp=tbry_west(j,k,itemp)
1484	if (temp.lt.8.) then
1485	SiO4=30.
1486	elseif (temp.ge.8. .and. temp.le.11.) then
1487	SiO4=30.-((temp-8.)*(20./3.))
1488	elseif (temp.gt.11. .and. temp.le.13.) then
1489	SiO4=10.-((temp-11.)*(8./2.))
1490	elseif (temp.gt.13. .and. temp.le.16.) then
1491	SiO4=2.-((temp-13.)*(2./3.))
1492	elseif (temp.gt.16.) then
1493	SiO4=0.
1494	endif
1495	xno3=1.67+0.5873SiO4+0.0144SiO4**2
1496	& +0.0003099SiO4*3
1497	itrc=iCAL_
1498	if(.not.got_tbry(itrc)) then
1499	tbry_west(j,k,itrc)=0.01
1500	endif
1501	itrc=iPOC_
1502	if(.not.got_tbry(itrc)) then
1503	tbry_west(j,k,itrc)=0.01
1504	endif
1505	itrc=iPHY_
1506	if(.not.got_tbry(itrc)) then
1507	tbry_west(j,k,itrc)=0.01
1508	endif
1509	itrc=iZOO_
1510	if(.not.got_tbry(itrc)) then
1511	tbry_west(j,k,itrc)=0.01
1512	endif
1513	itrc=iDIA_
1514	if(.not.got_tbry(itrc)) then
1515	tbry_west(j,k,itrc)=0.01
1516	endif
1517	itrc=iBSI_
1518	if(.not.got_tbry(itrc)) then
1519	tbry_west(j,k,itrc)=1.5e-3
1520	endif
1521	itrc=iBFE_
1522	if(.not.got_tbry(itrc)) then
1523	tbry_west(j,k,itrc)=1.e-2*5.e-6
1524	endif
1525	itrc=iGOC_
1526	if(.not.got_tbry(itrc)) then
1527	tbry_west(j,k,itrc)=0.01
1528	endif
1529	itrc=iSFE_
1530	if(.not.got_tbry(itrc)) then
1531	tbry_west(j,k,itrc)=1.e-2*5.e-6
1532	endif
1533	itrc=iDFE_
1534	if(.not.got_tbry(itrc)) then
1535	tbry_west(j,k,itrc)=1.e-2*5.e-6
1536	endif
1537	itrc=iDSI_
1538	if(.not.got_tbry(itrc)) then
1539	tbry_west(j,k,itrc)=1.e-2*0.15
1540	endif
1541	itrc=iNFE_
1542	if(.not.got_tbry(itrc)) then
1543	tbry_west(j,k,itrc)=1.e-2*5e-6
1544	endif
1545	itrc=iNCH_
1546	if(.not.got_tbry(itrc)) then
1547	tbry_west(j,k,itrc)=1.e-2*12./55.
1548	endif
1549	itrc=iDCH_
1550	if(.not.got_tbry(itrc)) then
1551	tbry_west(j,k,itrc)=1.e-2*12./55.
1552	endif
1553	itrc=iNH4_
1554	if(.not.got_tbry(itrc)) then
1555	tbry_west(j,k,itrc)=1.e-2
1556	endif
1557	enddo
1558	enddo
1559	endif
1560	# endif /* OBC_WEST */
1561	!
1562	# ifdef OBC_EAST
1563	if (EASTERN_EDGE) then
1564	do k=1,N
1565	do j=JstrR,JendR
1566	temp=tbry_east(j,k,itemp)
1567	if (temp.lt.8.) then
1568	SiO4=30.
1569	elseif (temp.ge.8. .and. temp.le.11.) then
1570	SiO4=30.-((temp-8.)*(20./3.))
1571	elseif (temp.gt.11. .and. temp.le.13.) then
1572	SiO4=10.-((temp-11.)*(8./2.))
1573	elseif (temp.gt.13. .and. temp.le.16.) then
1574	SiO4=2.-((temp-13.)*(2./3.))
1575	elseif (temp.gt.16.) then
1576	SiO4=0.
1577	endif
1578	xno3=1.67+0.5873SiO4+0.0144SiO4**2
1579	& +0.0003099SiO4*3
1580	itrc=iCAL_
1581	if(.not.got_tbry(itrc)) then
1582	tbry_east(j,k,itrc)=0.01
1583	endif
1584	itrc=iPOC_
1585	if(.not.got_tbry(itrc)) then
1586	tbry_east(j,k,itrc)=0.01
1587	endif
1588	itrc=iPHY_
1589	if(.not.got_tbry(itrc)) then
1590	tbry_east(j,k,itrc)=0.01
1591	endif
1592	itrc=iZOO_
1593	if(.not.got_tbry(itrc)) then
1594	tbry_east(j,k,itrc)=0.01
1595	endif
1596	itrc=iDIA_
1597	if(.not.got_tbry(itrc)) then
1598	tbry_east(j,k,itrc)=0.01
1599	endif
1600	itrc=iBSI_
1601	if(.not.got_tbry(itrc)) then
1602	tbry_east(j,k,itrc)=1.5e-3
1603	endif
1604	itrc=iBFE_
1605	if(.not.got_tbry(itrc)) then
1606	tbry_east(j,k,itrc)=1.e-2*5.e-6
1607	endif
1608	itrc=iGOC_
1609	if(.not.got_tbry(itrc)) then
1610	tbry_east(j,k,itrc)=0.01
1611	endif
1612	itrc=iSFE_
1613	if(.not.got_tbry(itrc)) then
1614	tbry_east(j,k,itrc)=1.e-2*5.e-6
1615	endif
1616	itrc=iDFE_
1617	if(.not.got_tbry(itrc)) then
1618	tbry_east(j,k,itrc)=1.e-2*5.e-6
1619	endif
1620	itrc=iDSI_
1621	if(.not.got_tbry(itrc)) then
1622	tbry_east(j,k,itrc)=1.e-2*0.15
1623	endif
1624	itrc=iNFE_
1625	if(.not.got_tbry(itrc)) then
1626	tbry_east(j,k,itrc)=1.e-2*5e-6
1627	endif
1628	itrc=iNCH_
1629	if(.not.got_tbry(itrc)) then
1630	tbry_east(j,k,itrc)=1.e-2*12./55.
1631	endif
1632	itrc=iDCH_
1633	if(.not.got_tbry(itrc)) then
1634	tbry_east(j,k,itrc)=1.e-2*12./55.
1635	endif
1636	itrc=iNH4_
1637	if(.not.got_tbry(itrc)) then
1638	tbry_east(j,k,itrc)=1.e-2
1639	endif
1640	enddo
1641	enddo
1642	endif
1643	# endif /* OBC_EAST */
1644	!
1645	# ifdef OBC_NORTH
1646	if (NORTHERN_EDGE) then
1647	do k=1,N
1648	do j=IstrR,IendR
1649	temp=tbry_north(j,k,itemp)
1650	if (temp.lt.8.) then
1651	SiO4=30.
1652	elseif (temp.ge.8. .and. temp.le.11.) then
1653	SiO4=30.-((temp-8.)*(20./3.))
1654	elseif (temp.gt.11. .and. temp.le.13.) then
1655	SiO4=10.-((temp-11.)*(8./2.))
1656	elseif (temp.gt.13. .and. temp.le.16.) then
1657	SiO4=2.-((temp-13.)*(2./3.))
1658	elseif (temp.gt.16.) then
1659	SiO4=0.
1660	endif
1661	xno3=1.67+0.5873SiO4+0.0144SiO4**2
1662	& +0.0003099SiO4*3
1663	itrc=iCAL_
1664	if(.not.got_tbry(itrc)) then
1665	tbry_north(j,k,itrc)=0.01
1666	endif
1667	itrc=iPOC_
1668	if(.not.got_tbry(itrc)) then
1669	tbry_north(j,k,itrc)=0.01
1670	endif
1671	itrc=iPHY_
1672	if(.not.got_tbry(itrc)) then
1673	tbry_north(j,k,itrc)=0.01
1674	endif
1675	itrc=iZOO_
1676	if(.not.got_tbry(itrc)) then
1677	tbry_north(j,k,itrc)=0.01
1678	endif
1679	itrc=iDIA_
1680	if(.not.got_tbry(itrc)) then
1681	tbry_north(j,k,itrc)=0.01
1682	endif
1683	itrc=iBSI_
1684	if(.not.got_tbry(itrc)) then
1685	tbry_north(j,k,itrc)=1.5e-3
1686	endif
1687	itrc=iBFE_
1688	if(.not.got_tbry(itrc)) then
1689	tbry_north(j,k,itrc)=1.e-2*5.e-6
1690	endif
1691	itrc=iGOC_
1692	if(.not.got_tbry(itrc)) then
1693	tbry_north(j,k,itrc)=0.01
1694	endif
1695	itrc=iSFE_
1696	if(.not.got_tbry(itrc)) then
1697	tbry_north(j,k,itrc)=1.e-2*5.e-6
1698	endif
1699	itrc=iDFE_
1700	if(.not.got_tbry(itrc)) then
1701	tbry_north(j,k,itrc)=1.e-2*5.e-6
1702	endif
1703	itrc=iDSI_
1704	if(.not.got_tbry(itrc)) then
1705	tbry_north(j,k,itrc)=1.e-2*0.15
1706	endif
1707	itrc=iNFE_
1708	if(.not.got_tbry(itrc)) then
1709	tbry_north(j,k,itrc)=1.e-2*5e-6
1710	endif
1711	itrc=iNCH_
1712	if(.not.got_tbry(itrc)) then
1713	tbry_north(j,k,itrc)=1.e-2*12./55.
1714	endif
1715	itrc=iDCH_
1716	if(.not.got_tbry(itrc)) then
1717	tbry_north(j,k,itrc)=1.e-2*12./55.
1718	endif
1719	itrc=iNH4_
1720	if(.not.got_tbry(itrc)) then
1721	tbry_north(j,k,itrc)=1.e-2
1722	endif
1723	enddo
1724	enddo
1725	endif
1726	# endif /* OBC_NORTH */
1727
1728	# ifdef OBC_SOUTH
1729	if (SOUTHERN_EDGE) then
1730	do k=1,N
1731	do j=IstrR,IendR
1732	temp=tbry_south(j,k,itemp)
1733	if (temp.lt.8.) then
1734	SiO4=30.
1735	elseif (temp.ge.8. .and. temp.le.11.) then
1736	SiO4=30.-((temp-8.)*(20./3.))
1737	elseif (temp.gt.11. .and. temp.le.13.) then
1738	SiO4=10.-((temp-11.)*(8./2.))
1739	elseif (temp.gt.13. .and. temp.le.16.) then
1740	SiO4=2.-((temp-13.)*(2./3.))
1741	elseif (temp.gt.16.) then
1742	SiO4=0.
1743	endif
1744	xno3=1.67+0.5873SiO4+0.0144SiO4**2
1745	& +0.0003099SiO4*3
1746	itrc=iCAL_
1747	if(.not.got_tbry(itrc)) then
1748	tbry_south(j,k,itrc)=0.01
1749	endif
1750	itrc=iPOC_
1751	if(.not.got_tbry(itrc)) then
1752	tbry_south(j,k,itrc)=0.01
1753	endif
1754	itrc=iPHY_
1755	if(.not.got_tbry(itrc)) then
1756	tbry_south(j,k,itrc)=0.01
1757	endif
1758	itrc=iZOO_
1759	if(.not.got_tbry(itrc)) then
1760	tbry_south(j,k,itrc)=0.01
1761	endif
1762	itrc=iDIA_
1763	if(.not.got_tbry(itrc)) then
1764	tbry_south(j,k,itrc)=0.01
1765	endif
1766	itrc=iBSI_
1767	if(.not.got_tbry(itrc)) then
1768	tbry_south(j,k,itrc)=1.5e-3
1769	endif
1770	itrc=iBFE_
1771	if(.not.got_tbry(itrc)) then
1772	tbry_south(j,k,itrc)=1.e-2*5.e-6
1773	endif
1774	itrc=iGOC_
1775	if(.not.got_tbry(itrc)) then
1776	tbry_south(j,k,itrc)=0.01
1777	endif
1778	itrc=iSFE_
1779	if(.not.got_tbry(itrc)) then
1780	tbry_south(j,k,itrc)=1.e-2*5.e-6
1781	endif
1782	itrc=iDFE_
1783	if(.not.got_tbry(itrc)) then
1784	tbry_south(j,k,itrc)=1.e-2*5.e-6
1785	endif
1786	itrc=iDSI_
1787	if(.not.got_tbry(itrc)) then
1788	tbry_south(j,k,itrc)=1.e-2*0.15
1789	endif
1790	itrc=iNFE_
1791	if(.not.got_tbry(itrc)) then
1792	tbry_south(j,k,itrc)=1.e-2*5e-6
1793	endif
1794	itrc=iNCH_
1795	if(.not.got_tbry(itrc)) then
1796	tbry_south(j,k,itrc)=1.e-2*12./55.
1797	endif
1798	itrc=iDCH_
1799	if(.not.got_tbry(itrc)) then
1800	tbry_south(j,k,itrc)=1.e-2*12./55.
1801	endif
1802	itrc=iNH4_
1803	if(.not.got_tbry(itrc)) then
1804	tbry_south(j,k,itrc)=1.e-2
1805	endif
1806	enddo
1807	enddo
1808	endif
1809	# endif /* OBC_SOUTH */
1810	# endif /* PISCES */
1811	return
1812	end
1813	#else
1814	subroutine empty_analytical
1815	return
1816	end
1817	#endif /* BIOLOGY && defined T_FRC_BRY */
1818
1819

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: Roms_tools/Roms_Agrif/analytical.F @ 5

Download in other formats: