source: LMDZ6/trunk/libf/phylmd/cltracrn.f90 @ 5277

Last change on this file since 5277 was 5274, checked in by abarral, 7 hours ago

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