source: LMDZ.3.3/trunk/libf/phylmd/cltracrn.F @ 1176

Last change on this file since 1176 was 2, checked in by lmdz, 25 years ago

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