source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cltracrn.F90 @ 3269

Last change on this file since 3269 was 1279, checked in by Laurent Fairhead, 15 years ago

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