source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/cltracrn.F90 @ 5441

Last change on this file since 5441 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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