source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/sediment_mod.f90 @ 5441

Last change on this file since 5441 was 5182, checked in by abarral, 4 months ago

(WIP) Replace REPROBUS CPP KEY by logical
properly name modules
File size: 10.3 KB
Line 
1!----- This SUBROUTINE calculates the sedimentation flux of Tracers
2
3SUBROUTINE sediment_mod(t_seri, pplay, zrho, paprs, time_step, RHcl, &
4        id_coss, id_codu, id_scdu, &
5        ok_chimeredust, &
6        sed_ss, sed_dust, sed_dustsco, &
7        sed_ss3D, sed_dust3D, sed_dustsco3D, tr_seri)
8  !nhl     .                                       xlon,xlat,
9
10  USE dimphy
11  USE lmdz_infotrac
12  USE lmdz_YOECUMF
13  USE lmdz_yomcst
14  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
15  USE lmdz_chem, ONLY: idms, iso2, iso4, ih2s, idmso, imsa, ih2o2, &
16          n_avogadro, masse_s, masse_so4, rho_water, rho_ice
17
18  IMPLICIT NONE
19
20  REAL :: RHcl(klon, klev)     ! humidite relative ciel clair
21  REAL :: tr_seri(klon, klev, nbtr) !conc of tracers
22  REAL :: sed_ss(klon) !sedimentation flux of Sea Salt (g/m2/s)
23  REAL :: sed_dust(klon) !sedimentation flux of dust (g/m2/s)
24  REAL :: sed_dustsco(klon) !sedimentation flux of scoarse  dust (g/m2/s)
25  REAL :: sed_ss3D(klon, klev) !sedimentation flux of Sea Salt (g/m2/s)
26  REAL :: sed_dust3D(klon, klev) !sedimentation flux of dust (g/m2/s)
27  REAL :: sed_dustsco3D(klon, klev) !sedimentation flux of scoarse  dust (g/m2/s)
28  REAL :: t_seri(klon, klev)   !Temperature at mid points of Z (K)
29  REAL :: v_dep_ss(klon, klev)  ! sed. velocity for SS m/s
30  REAL :: v_dep_dust(klon, klev)  ! sed. velocity for dust m/s
31  REAL :: v_dep_dustsco(klon, klev)  ! sed. velocity for dust m/s
32  REAL :: pplay(klon, klev)    !pressure at mid points of Z (Pa)
33  REAL :: zrho(klon, klev)     !Density of air at mid points of Z (kg/m3)
34  REAL :: paprs(klon, klev + 1)    !pressure at interface of layers Z (Pa)
35  REAL :: time_step            !time step (sec)
36  LOGICAL :: ok_chimeredust
37  REAL :: xlat(klon)       ! latitudes pour chaque point
38  REAL :: xlon(klon)       ! longitudes pour chaque point
39  INTEGER :: id_coss, id_codu, id_scdu
40
41  !------local variables
42
43  INTEGER :: i, k, nbre_RH
44  PARAMETER(nbre_RH = 12)
45
46  REAL :: lambda, ss_g
47  REAL :: mmd_ss      !mass median diameter of SS (um)
48  REAL :: mmd_dust          !mass median diameter of dust (um)
49  REAL :: mmd_dustsco          !mass median diameter of scoarse dust (um)
50  REAL :: rho_ss(nbre_RH), rho_ss1 !density of sea salt (kg/m3)
51  REAL :: rho_dust          !density of dust(kg/m3)
52  REAL :: v_stokes, CC, v_sed, ss_growth_f(nbre_RH)
53  REAL :: sed_flux(klon, klev)  ! sedimentation flux g/m2/s
54  REAL :: air_visco(klon, klev)
55  REAL :: zdz(klon, klev)       ! layers height (m)
56  REAL :: temp                 ! temperature in degree Celius
57
58  INTEGER :: RH_num
59  REAL :: RH_MAX, DELTA, rh, RH_tab(nbre_RH)
60  PARAMETER (RH_MAX = 95.)
61
62  DATA RH_tab/0., 10., 20., 30., 40., 50., 60., 70., 80., 85., 90., 95./
63
64  DATA rho_ss/2160., 2160., 2160., 2160, 1451.6, 1367.9, &
65          1302.9, 1243.2, 1182.7, 1149.5, 1111.6, 1063.1/
66
67  DATA ss_growth_f/0.503, 0.503, 0.503, 0.503, 0.724, 0.782, &
68          0.838, 0.905, 1.000, 1.072, 1.188, 1.447/
69
70  mmd_ss = 12.7   !dia -um at 80% for bin 0.5-20 um but 90% of real mmd
71  ! obsolete      mmd_dust=2.8  !micrometer for bin 0.5-20 and 0.5-10 um
72  ! 4tracer SPLA:       mmd_dust=11.0  !micrometer for bin 0.5-20 and 0.5-10 um
73  !3days       mmd_dust=3.333464  !micrometer for bin 0.5-20 and 0.5-10 um
74  !3days       mmd_dustsco=12.91315  !micrometer for bin 0.5-20 and 0.5-10 um
75  !JE20140911       mmd_dust=3.002283  !micrometer for bin 0.5-20 and 0.5-10 um
76  !JE20140911       mmd_dustsco=13.09771  !micrometer for bin 0.5-20 and 0.5-10 um
77  !JE20140911        mmd_dust=5.156346  !micrometer for bin 0.5-20 and 0.5-10 um
78  !JE20140911        mmd_dustsco=15.56554  !micrometer for bin 0.5-20 and 0.5-10 um
79  IF (ok_chimeredust) THEN
80    !JE20150212<< : changes in ustar in dustmod changes emission distribution
81    ! mmd_dust=3.761212  !micrometer for bin 0.5-3 and 0.5-10 um
82    ! mmd_dustsco=15.06167  !micrometer for bin 3-20 and 0.5-10 um
83    !JE20150212>>
84    !JE20150618: Change in div3 of dustmod changes distribution. now is div3=6
85    !div=3        mmd_dust=3.983763
86    !div=3        mmd_dustsco=15.10854
87    mmd_dust = 3.898047
88    mmd_dustsco = 15.06167
89  ELSE
90    mmd_dust = 11.0  !micrometer for bin 0.5-20 and 0.5-10 um
91    mmd_dustsco = 100. ! absurd value, bin not used in this scheme
92  ENDIF
93
94  rho_dust = 2600. !kg/m3
95
96  !--------- Air viscosity (poise=0.1 kg/m-sec)-----------
97
98  DO k = 1, klev
99    DO i = 1, klon
100
101      zdz(i, k) = (paprs(i, k) - paprs(i, k + 1)) / zrho(i, k) / RG
102
103      temp = t_seri(i, k) - RTT
104
105      IF (temp<0.) THEN
106        air_visco(i, k) = (1.718 + 0.0049 * temp - 1.2e-5 * temp * temp) * 1.e-4
107      ELSE
108        air_visco(i, k) = (1.718 + 0.0049 * temp) * 1.e-4
109      ENDIF
110
111    ENDDO
112  ENDDO
113
114  !--------- for Sea Salt -------------------
115
116  IF(id_coss>0) THEN
117    DO k = 1, klev
118      DO i = 1, klon
119
120        !---cal. correction factor hygroscopic growth of aerosols
121
122        rh = MIN(RHcl(i, k) * 100., RH_MAX)
123        RH_num = INT(rh / 10. + 1.)
124        IF (rh>85.) RH_num = 10
125        IF (rh>90.) RH_num = 11
126        DELTA = (rh - RH_tab(RH_num)) / (RH_tab(RH_num + 1) - RH_tab(RH_num))
127
128        ss_g = ss_growth_f(rh_num) + &
129                DELTA * (ss_growth_f(RH_num + 1) - ss_growth_f(RH_num))
130
131        rho_ss1 = rho_ss(rh_num) + &
132                DELTA * (rho_ss(RH_num + 1) - rho_ss(RH_num))
133
134        v_stokes = RG * (rho_ss1 - zrho(i, k)) * & !m/sec
135                (mmd_ss * ss_g) * (mmd_ss * ss_g) * &
136                1.e-12 / (18.0 * air_visco(i, k) / 10.)
137
138        lambda = 6.6 * 1.e-8 * (103125 / pplay(i, k)) * (t_seri(i, k) / 293.15)
139
140        CC = 1.0 + 1.257 * lambda / (mmd_ss * ss_g) / 1.e6  ! C-correction factor
141
142        v_sed = v_stokes * CC                       ! m/sec !orig
143
144        !---------check for v_sed*dt<zdz
145
146        IF (v_sed * time_step>zdz(i, k)) THEN
147          v_sed = zdz(i, k) / time_step
148        ENDIF
149
150        v_dep_ss(i, k) = v_sed
151        sed_flux(i, k) = tr_seri(i, k, id_coss) * v_sed !g/cm3*m/sec
152        !sed_ss3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
153        ! conc_sed_ss3D(i,k)=sed_flux(i,k)*1.e6      !g/m3*sec !!!!!!!
154
155      ENDDO          !klon
156    ENDDO          !klev
157
158    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
159    sed_ss3D(:, :) = 0.0  ! initialisation
160
161    DO k = 1, klev
162      DO i = 1, klon
163        sed_ss3D(i, k) = sed_ss3D(i, k) - &
164                sed_flux(i, k) / zdz(i, k) !!!!!!!!!!!!!!!!!!!!!!
165      ENDDO          !klon
166    ENDDO          !klev
167
168    DO k = 1, klev - 1
169      DO i = 1, klon
170        sed_ss3D(i, k) = sed_ss3D(i, k) + &
171                sed_flux(i, k + 1) / zdz(i, k) !!!!!!!!
172
173      ENDDO          !klon
174    ENDDO          !klev
175
176    DO k = 1, klev
177      DO i = 1, klon
178        tr_seri(i, k, id_coss) = tr_seri(i, k, id_coss) + &
179                sed_ss3D(i, k) * time_step
180      ENDDO
181    ENDDO
182
183    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184
185    DO i = 1, klon
186      sed_ss(i) = sed_flux(i, 1) * 1.e6 * 1.e3    !--unit mg/m2/s
187    ENDDO          !klon
188  ELSE
189    DO i = 1, klon
190      sed_ss(i) = 0.
191    ENDDO
192  ENDIF
193
194  !
195
196  !--------- For dust ------------------
197
198  IF(id_codu>0) THEN
199    DO k = 1, klev
200      DO i = 1, klon
201
202        v_stokes = RG * (rho_dust - zrho(i, k)) * & !m/sec
203                mmd_dust * mmd_dust * &
204                1.e-12 / (18.0 * air_visco(i, k) / 10.)
205
206        lambda = 6.6 * 1.e-8 * (103125 / pplay(i, k)) * (t_seri(i, k) / 293.15)
207        CC = 1.0 + 1.257 * lambda / (mmd_dust) / 1.e6        !dimensionless
208        v_sed = v_stokes * CC                       !m/sec
209
210        !---------check for v_sed*dt<zdz
211
212        IF (v_sed * time_step>zdz(i, k)) THEN
213          v_sed = zdz(i, k) / time_step
214        ENDIF
215
216        v_dep_dust(i, k) = v_sed
217        sed_flux(i, k) = tr_seri(i, k, id_codu) * v_sed !g/cm3.m/sec
218        !sed_dust3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
219
220      ENDDO          !klon
221    ENDDO          !klev
222
223    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224    sed_dust3D(:, :) = 0.0  ! initialisation
225
226    DO k = 1, klev
227      DO i = 1, klon
228        sed_dust3D(i, k) = sed_dust3D(i, k) - &
229                sed_flux(i, k) / zdz(i, k)
230      ENDDO          !klon
231    ENDDO          !klev
232
233
234    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235
236    DO k = 1, klev - 1
237      DO i = 1, klon
238        sed_dust3D(i, k) = sed_dust3D(i, k) + &
239                sed_flux(i, k + 1) / zdz(i, k)
240      ENDDO          !klon
241    ENDDO          !klev
242
243    DO k = 1, klev
244      DO i = 1, klon
245        tr_seri(i, k, id_codu) = tr_seri(i, k, id_codu) + &
246                sed_dust3D(i, k) * time_step
247      ENDDO
248    ENDDO
249
250    DO i = 1, klon
251      sed_dust(i) = sed_flux(i, 1) * 1.e6 * 1.e3    !--unit mg/m2/s
252    ENDDO          !klon
253  ELSE
254    DO i = 1, klon
255      sed_dust(i) = 0.
256    ENDDO
257  ENDIF
258  !
259
260
261  !--------- For scoarse  dust ------------------
262
263  IF(id_scdu>0) THEN
264    DO k = 1, klev
265      DO i = 1, klon
266
267        v_stokes = RG * (rho_dust - zrho(i, k)) * & !m/sec
268                mmd_dustsco * mmd_dustsco * &
269                1.e-12 / (18.0 * air_visco(i, k) / 10.)
270
271        lambda = 6.6 * 1.e-8 * (103125 / pplay(i, k)) * (t_seri(i, k) / 293.15)
272        CC = 1.0 + 1.257 * lambda / (mmd_dustsco) / 1.e6        !dimensionless
273        v_sed = v_stokes * CC                       !m/sec
274
275        !---------check for v_sed*dt<zdz
276
277        IF (v_sed * time_step>zdz(i, k)) THEN
278          v_sed = zdz(i, k) / time_step
279        ENDIF
280
281        v_dep_dustsco(i, k) = v_sed
282        sed_flux(i, k) = tr_seri(i, k, id_scdu) * v_sed !g/cm3.m/sec
283        !sed_dustsco3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
284
285      ENDDO          !klon
286    ENDDO          !klev
287
288    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
289    sed_dustsco3D(:, :) = 0.0  ! initialisation
290
291    DO k = 1, klev
292      DO i = 1, klon
293        sed_dustsco3D(i, k) = sed_dustsco3D(i, k) - &
294                sed_flux(i, k) / zdz(i, k)
295      ENDDO          !klon
296    ENDDO          !klev
297
298    DO k = 1, klev - 1
299      DO i = 1, klon
300        sed_dustsco3D(i, k) = sed_dustsco3D(i, k) + &
301                sed_flux(i, k + 1) / zdz(i, k)
302      ENDDO          !klon
303    ENDDO          !klev
304
305    DO k = 1, klev
306      DO i = 1, klon
307        tr_seri(i, k, id_scdu) = tr_seri(i, k, id_scdu) + &
308                sed_dustsco3D(i, k) * time_step
309      ENDDO
310    ENDDO
311    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312
313    DO i = 1, klon
314      sed_dustsco(i) = sed_flux(i, 1) * 1.e6 * 1.e3    !--unit mg/m2/s
315    ENDDO          !klon
316  ELSE
317    DO i = 1, klon
318      sed_dustsco(i) = 0.
319    ENDDO
320  ENDIF
321  !
322
323
324
325
326  !
327
328END SUBROUTINE sediment_mod
Note: See TracBrowser for help on using the repository browser.