source: LMDZ6/branches/Amaury_dev/libf/phylmd/cltracrn.F90 @ 5409

Last change on this file since 5409 was 5144, checked in by abarral, 5 months ago

Put YOMCST.h into modules

  • 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.0 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 lmdz_yomcst
13
14  IMPLICIT NONE
15  !======================================================================
16  ! Auteur(s): Alex/LMD) date:  fev 99
17  !            inspire de clqh + clvent
18  ! Objet: diffusion verticale de traceurs avec quantite de traceur ds
19  !        le sol ( reservoir de sol de radon )
20
21  ! note : pour l'instant le traceur dans le sol et le flux sont
22  !        calcules mais ils ne servent que de diagnostiques
23  !        seule la tendance sur le traceur est sortie (d_tr)
24  !---------------------------------------------------------------------
25  ! Arguments:
26  ! itr......input-R-  le type de traceur : id_rn(radon), id_pb(plomb)
27  ! dtime....input-R-  intervalle du temps (en secondes) ~ pdtphys
28  ! u1lay....input-R-  vent u de la premiere couche (m/s)
29  ! v1lay....input-R-  vent v de la premiere couche (m/s)
30  ! cdrag....input-R-  cdrag
31  ! coef.....input-R-  le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
32  ! t........input-R-  temperature (K)
33  ! paprs....input-R-  pression a inter-couche (Pa)
34  ! pplay....input-R-  pression au milieu de couche (Pa)
35  ! delp.....input-R-  epaisseur de couche (Pa)
36  ! ftsol....input-R-  temperature du sol (en Kelvin)
37  ! tr.......input-R-  traceurs
38  ! trs......input-R-  traceurs dans le sol
39  ! masktr...input-R-  Masque reservoir de sol traceur (1 = reservoir)
40  ! fshtr....input-R-  Flux surfacique de production dans le sol
41  ! tautr....input-R-  Constante de decroissance du traceur
42  ! vdeptr...input-R-  Vitesse de depot sec dans la couche brownienne
43  ! hsoltr...input-R-  Epaisseur equivalente du reservoir de sol
44  ! lat......input-R-  latitude en degree
45  ! d_tr.....output-R- le changement de "tr"
46  ! d_trs....output-R- le changement de "trs"
47  !======================================================================
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)) == 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==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=70.).OR.      &
200            (itr==id_pb.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=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==id_rn.AND.NINT(masktr(i))==0).OR.       &
226            (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)<-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==id_rn.AND.NINT(masktr(i))==0).OR.    &
235            (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>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==id_rn.AND.NINT(masktr(i))==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.