source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/cltracrn.F @ 5427

Last change on this file since 5427 was 524, checked in by lmdzadmin, 21 years ago

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