LMDZ
cltracrn.F90
Go to the documentation of this file.
1 !$Id $
2 
3 SUBROUTINE 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 
288 END SUBROUTINE cltracrn
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
traclmdz_mod::id_rn
integer, save id_rn
Definition: traclmdz_mod.F90:41
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
traclmdz_mod::id_pb
integer, save id_pb
Definition: traclmdz_mod.F90:41
pplay
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL pplay
Definition: calcul_STDlev.h:26
cltracrn
subroutine cltracrn(itr, dtime, u1lay, v1lay, cdrag, coef, t, ftsol, pctsrf, tr, trs, paprs, pplay, delp, masktr, fshtr, hsoltr, tautr, vdeptr, lat, d_tr, d_trs)
Definition: cltracrn.F90:8
paprs
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL &zphi geo500!IM on interpole a chaque pas de temps le paprs
Definition: calcul_STDlev.h:144
indice_sol_mod::nbsrf
integer, parameter nbsrf
Definition: indice_sol_mod.F90:2
cdrag
subroutine cdrag(knon, nsrf, speed, t1, q1, zgeop1, psol, tsurf, qsurf, z0m, z0h, pcfm, pcfh, zri, pref)
Definition: cdrag.F90:8
traclmdz_mod
Definition: traclmdz_mod.F90:3
dimphy
Definition: dimphy.F90:1
indice_sol_mod
Definition: indice_sol_mod.F90:1
rg
real rg
Definition: comcstphy.h:1