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