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

Last change on this file since 5104 was 5104, checked in by abarral, 3 months ago

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