source: LMDZ4/branches/V3_test/libf/phylmd/cltracrn.F @ 3947

Last change on this file since 3947 was 704, checked in by Laurent Fairhead, 18 years ago

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