Context Navigation

cltracrn.F90 @ 1237

Last change on this file since 1237 was 1191, checked in by jghattas, 15 years ago

Reecriture de phytrac et les routines concernes (Anthony Jamelot)

les suffix change de F -> F90 (nflxtr.F90,cltracrn.F90,initrrnpb.F90,cvltr.F90,minmaxqfi.F90,cltrac.F90,phytrac.F90)

Traitement d'un nouveau traceur berelium (optionel, toujours pour des
tests)(Anthony Jamelot)

radiornpb.F change du nom pour radio_decay.F90 car il traite maintenant tout les traceurs radioactives
ajoute init_be.F90

Nouveau interface dans phytrac pour serparer les calculs et appels
specifique a INCA avec les traitements des traceurs specifiques au LMDZ
(JG)

ajoute tracinca_mod.F90 pour les appeles a INCA
ajoute traclmdz_mod.F90 pour les calculs des traceurs specifiques a LMDZ
enleve fichier restartrac et ajoute la variable trs dans restartphy.nc

La convergence numerique a etait rompue uniquement pour les traceurs
LMDZ RN et PB.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 10.7 KB

Line
1	!$Id $
2
3	SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &
4	cdrag,coef,t,ftsol,pctsrf, &
5	tr,trs,paprs,pplay,delp, &
6	masktr,fshtr,hsoltr,tautr,vdeptr, &
7	lat,d_tr,d_trs )
8
9	USE dimphy
10	IMPLICIT NONE
11	!======================================================================
12	! Auteur(s): Alex/LMD) date: fev 99
13	! inspire de clqh + clvent
14	! Objet: diffusion verticale de traceurs avec quantite de traceur ds
15	! le sol ( reservoir de sol de radon )
16	!
17	! note : pour l'instant le traceur dans le sol et le flux sont
18	! calcules mais ils ne servent que de diagnostiques
19	! seule la tendance sur le traceur est sortie (d_tr)
20	!---------------------------------------------------------------------
21	! Arguments:
22	! itr......input-R- le type de traceur 1- Rn 2 - Pb
23	! dtime....input-R- intervalle du temps (en secondes) ~ pdtphys
24	! u1lay....input-R- vent u de la premiere couche (m/s)
25	! v1lay....input-R- vent v de la premiere couche (m/s)
26	! cdrag....input-R- cdrag
27	! coef.....input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
28	! t........input-R- temperature (K)
29	! paprs....input-R- pression a inter-couche (Pa)
30	! pplay....input-R- pression au milieu de couche (Pa)
31	! delp.....input-R- epaisseur de couche (Pa)
32	! ftsol....input-R- temperature du sol (en Kelvin)
33	! tr.......input-R- traceurs
34	! trs......input-R- traceurs dans le sol
35	! masktr...input-R- Masque reservoir de sol traceur (1 = reservoir)
36	! fshtr....input-R- Flux surfacique de production dans le sol
37	! tautr....input-R- Constante de decroissance du traceur
38	! vdeptr...input-R- Vitesse de depot sec dans la couche brownienne
39	! hsoltr...input-R- Epaisseur equivalente du reservoir de sol
40	! lat......input-R- latitude en degree
41	! d_tr.....output-R- le changement de "tr"
42	! d_trs....output-R- le changement de "trs"
43	!======================================================================
44	include "YOMCST.h"
45	include "indicesol.h"
46	!
47	!Entrees
48	INTEGER,INTENT(IN) :: itr
49	REAL,INTENT(IN) :: dtime
50	REAL,DIMENSION(klon),INTENT(IN) :: u1lay, v1lay
51	REAL,DIMENSION(klon),INTENT(IN) :: cdrag
52	REAL,DIMENSION(klon,klev),INTENT(IN) :: coef, t
53	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol, pctsrf
54	REAL,DIMENSION(klon,klev),INTENT(IN) :: tr
55	REAL,DIMENSION(klon),INTENT(IN) :: trs
56	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
57	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay, delp
58	REAL,DIMENSION(klon),INTENT(IN) :: masktr
59	REAL,DIMENSION(klon),INTENT(IN) :: fshtr
60	REAL,INTENT(IN) :: hsoltr
61	REAL,INTENT(IN) :: tautr
62	REAL,INTENT(IN) :: vdeptr
63	REAL,DIMENSION(klon),INTENT(IN) :: lat
64
65	!Sorties
66	REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
67	REAL,DIMENSION(klon),INTENT(OUT) :: d_trs ! (diagnostic) traceur ds le sol
68
69	!Locales
70	REAL,DIMENSION(klon,klev) :: flux_tr ! (diagnostic) flux de traceur
71	INTEGER :: i, k, n, l
72	REAL,DIMENSION(klon) :: rotrhi
73	REAL,DIMENSION(klon,klev) :: zx_coef
74	REAL,DIMENSION(klon) :: zx_buf
75	REAL,DIMENSION(klon,klev) :: zx_ctr
76	REAL,DIMENSION(klon,klev) :: zx_dtr
77	REAL,DIMENSION(klon) :: zx_trs
78	REAL :: zx_a, zx_b
79
80	REAL,DIMENSION(klon,klev) :: local_tr
81	REAL,DIMENSION(klon) :: local_trs
82	REAL,DIMENSION(klon) :: zts ! champ de temperature du sol
83	REAL,DIMENSION(klon) :: zx_alpha1, zx_alpha2
84	!======================================================================
85	!AA Pour l'instant les 4 types de surface ne sont pas pris en compte
86	!AA On fabrique avec zts un champ de temperature de sol
87	!AA que le pondere par la fraction de nature de sol.
88
89	DO i = 1,klon
90	zts(i) = 0.
91	ENDDO
92
93	DO n=1,nbsrf
94	DO i = 1,klon
95	zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
96	ENDDO
97	ENDDO
98
99	DO i = 1,klon
100	rotrhi(i) = RD * zts(i) / hsoltr
101	ENDDO
102
103	DO k = 1, klev
104	DO i = 1, klon
105	local_tr(i,k) = tr(i,k)
106	ENDDO
107	ENDDO
108
109	DO i = 1, klon
110	local_trs(i) = trs(i)
111	ENDDO
112	!======================================================================
113	!AA Attention si dans clmain zx_alf1(i) = 1.0
114	!AA Il doit y avoir coherence (dc la meme chose ici)
115
116	DO i = 1, klon
117	!AA zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
118	zx_alpha1(i) = 1.0
119	zx_alpha2(i) = 1.0 - zx_alpha1(i)
120	ENDDO
121	!======================================================================
122	DO i = 1, klon
123	zx_coef(i,1) = cdrag(i)(1.0+SQRT(u1lay(i)2+v1lay(i)*2)) &
124	pplay(i,1)/(RDt(i,1))
125	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
126	ENDDO
127
128	DO k = 2, klev
129	DO i = 1, klon
130	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
131	(paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
132	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
133	ENDDO
134	ENDDO
135	!======================================================================
136	DO i = 1, klon
137	zx_buf(i) = delp(i,klev) + zx_coef(i,klev)
138	zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
139	zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
140	ENDDO
141
142	DO l = klev-1, 2 , -1
143	DO i = 1, klon
144	zx_buf(i) = delp(i,l)+zx_coef(i,l) &
145	+zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
146
147	zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) &
148	+ zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
149	zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
150	ENDDO
151	ENDDO
152
153	DO i = 1, klon
154	zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2)) &
155	+ masktr(i) * zx_coef(i,1) &
156	( zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2) )
157
158	zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1) &
159	+ zx_ctr(i,2) &
160	*(zx_coef(i,2) &
161	- masktr(i) * zx_coef(i,1) &
162	*zx_alpha2(i) ) ) / zx_buf(i)
163	zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
164	ENDDO
165	!======================================================================
166	! Calculer d'abord local_trs nouvelle quantite dans le reservoir
167	! de sol
168	!=====================================================================
169
170	DO i = 1, klon
171	!-------------------------
172	! Au dessus des continents
173	!--
174	! Le pb peut se deposer partout : vdeptr = 10-3 m/s
175	! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
176	!-------------------------------------------------------------------
177	IF ( NINT(masktr(i)) .EQ. 1 ) THEN
178	zx_trs(i) = local_trs(i)
179	zx_a = zx_trs(i) &
180	+fshtr(i)dtimerotrhi(i) &
181	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
182	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
183	+zx_alpha2(i)*zx_ctr(i,2))
184	! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance
185	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
186	* (1.-zx_dtr(i,1) &
187	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
188	+ dtime / tautr &
189	+ dtime * vdeptr / hsoltr
190	zx_trs(i) = zx_a / zx_b
191	local_trs(i) = zx_trs(i)
192	ENDIF
193	!--------------------------------------------------------
194	! Si on est entre 60N et 70N on divise par 2 l'emanation
195	!--------------------------------------------------------
196
197	IF ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR. &
198	(itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN
199	zx_trs(i) = local_trs(i)
200	zx_a = zx_trs(i) &
201	+(fshtr(i)/2.)dtimerotrhi(i) &
202	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
203	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
204	+zx_alpha2(i)*zx_ctr(i,2))
205	!
206	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
207	* (1.-zx_dtr(i,1) &
208	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
209	+ dtime / tautr &
210	+ dtime * vdeptr / hsoltr
211	!
212	zx_trs(i) = zx_a / zx_b
213	local_trs(i) = zx_trs(i)
214	ENDIF
215
216	!----------------------------------------------
217	! Au dessus des oceans et aux hautes latitudes
218	!--
219	! au dessous de -60S pas d'emission de radon au dessus
220	! des oceans et des continents
221	!---------------------------------------------------------------
222
223	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR. &
224	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
225	zx_trs(i) = 0.
226	local_trs(i) = 0.
227	END IF
228	!--
229	! au dessus de 70 N pas d'emission de radon au dessus
230	! des oceans et des continents
231	!--------------------------------------------------------------
232	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR. &
233	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
234	zx_trs(i) = 0.
235	local_trs(i) = 0.
236	END IF
237	!---------------------------------------------
238	! Au dessus des oceans la source est nulle
239	!--------------------------------------------
240
241	IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
242	zx_trs(i) = 0.
243	local_trs(i) = 0.
244	END IF
245
246	ENDDO ! sur le i=1,klon
247	!
248	!======================================================================
249	! Une fois on a zx_trs, on peut faire l'iteration
250	!======================================================================
251
252	DO i = 1, klon
253	local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
254	ENDDO
255	DO l = 2, klev
256	DO i = 1, klon
257	local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
258	ENDDO
259	ENDDO
260	!======================================================================
261	! Calcul du flux de traceur (flux_tr): UA/(m**2 s)
262	!======================================================================
263	DO i = 1, klon
264	flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG &
265	* (zx_alpha1(i)local_tr(i,1)+zx_alpha2(i)local_tr(i,2) &
266	-zx_trs(i)) / dtime
267	ENDDO
268	DO l = 2, klev
269	DO i = 1, klon
270	flux_tr(i,l) = zx_coef(i,l)/RG &
271	* (local_tr(i,l)-local_tr(i,l-1)) / dtime
272	ENDDO
273	ENDDO
274	!======================================================================
275	! Calcul des tendances du traceur ds le sol et dans l'atmosphere
276	!======================================================================
277	DO l = 1, klev
278	DO i = 1, klon
279	d_tr(i,l) = local_tr(i,l) - tr(i,l)
280	ENDDO
281	ENDDO
282	DO i = 1, klon
283	d_trs(i) = local_trs(i) - trs(i)
284	ENDDO
285
286	END SUBROUTINE cltracrn

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

Download in other formats:

Original Format