source: trunk/LMDZ.EARTH/libf/phylmd/cltracrn.F90 @ 1157

Last change on this file since 1157 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 10.8 KB
Line 
1!$Id $
2
3SUBROUTINE 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
287END SUBROUTINE cltracrn
Note: See TracBrowser for help on using the repository browser.