source: LMDZ5/branches/testing/libf/phylmd/Dust/sediment_mod.F @ 4538

Last change on this file since 4538 was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch
File size: 10.3 KB
Line 
1c----- This subroutine calculates the sedimentation flux of Tracers
2c
3      SUBROUTINE 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)
8cnhl     .                                       xlon,xlat,
9c
10       USE dimphy
11       USE infotrac
12      IMPLICIT NONE
13c
14#include "dimensions.h"
15#include "chem.h"
16c #include "dimphy.h"
17#include "YOMCST.h"
18#include "YOECUMF.h"
19c
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
40c
41c------local variables
42c
43       INTEGER i, k, nbre_RH
44       PARAMETER(nbre_RH=12)
45c
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
57c
58       INTEGER RH_num
59       REAL RH_MAX, DELTA, rh, RH_tab(nbre_RH)
60       PARAMETER (RH_MAX=95.)
61c
62       DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
63c
64c
65       DATA rho_ss/2160. ,2160. ,2160.,  2160,  1451.6, 1367.9,
66     .             1302.9,1243.2,1182.7, 1149.5,1111.6, 1063.1/
67c
68       DATA ss_growth_f/0.503, 0.503, 0.503, 0.503, 0.724, 0.782,
69     .                  0.838, 0.905, 1.000, 1.072, 1.188, 1.447/
70c
71c
72       mmd_ss=12.7   !dia -um at 80% for bin 0.5-20 um but 90% of real mmd
73!               obsolete      mmd_dust=2.8  !micrometer for bin 0.5-20 and 0.5-10 um
74! 4tracer SPLA:       mmd_dust=11.0  !micrometer for bin 0.5-20 and 0.5-10 um
75!3days       mmd_dust=3.333464  !micrometer for bin 0.5-20 and 0.5-10 um
76!3days       mmd_dustsco=12.91315  !micrometer for bin 0.5-20 and 0.5-10 um
77!JE20140911       mmd_dust=3.002283  !micrometer for bin 0.5-20 and 0.5-10 um
78!JE20140911       mmd_dustsco=13.09771  !micrometer for bin 0.5-20 and 0.5-10 um
79!JE20140911        mmd_dust=5.156346  !micrometer for bin 0.5-20 and 0.5-10 um
80!JE20140911        mmd_dustsco=15.56554  !micrometer for bin 0.5-20 and 0.5-10 um
81        IF (ok_chimeredust) THEN
82!JE20150212<< : changes in ustar in dustmod changes emission distribution
83!        mmd_dust=3.761212  !micrometer for bin 0.5-3 and 0.5-10 um
84!        mmd_dustsco=15.06167  !micrometer for bin 3-20 and 0.5-10 um
85!JE20150212>>
86!JE20150618: Change in div3 of dustmod changes distribution. now is div3=6
87!div=3        mmd_dust=3.983763
88!div=3        mmd_dustsco=15.10854
89        mmd_dust=3.898047
90        mmd_dustsco=15.06167
91        ELSE
92        mmd_dust=11.0  !micrometer for bin 0.5-20 and 0.5-10 um
93        mmd_dustsco=100. ! absurd value, bin not used in this scheme
94        ENDIF
95
96
97       rho_dust=2600. !kg/m3
98c
99c--------- Air viscosity (poise=0.1 kg/m-sec)-----------
100c
101       DO k=1, klev
102       DO i=1, klon
103c
104       zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
105c
106       temp=t_seri(i,k)-RTT
107c
108       IF (temp.LT.0.) THEN
109         air_visco(i,k)=(1.718+0.0049*temp-1.2e-5*temp*temp)*1.e-4
110       ELSE
111         air_visco(i,k)=(1.718+0.0049*temp)*1.e-4
112       ENDIF
113c
114       ENDDO
115       ENDDO
116c
117c--------- for Sea Salt -------------------
118c
119c
120c
121       IF(id_coss>0) THEN
122       DO k=1, klev
123       DO i=1,klon
124c
125c---cal. correction factor hygroscopic growth of aerosols
126c
127        rh=MIN(RHcl(i,k)*100.,RH_MAX)
128        RH_num = INT( rh/10. + 1.)
129        IF (rh.gt.85.) RH_num=10
130        IF (rh.gt.90.) RH_num=11
131        DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
132c
133        ss_g=ss_growth_f(rh_num) +
134     .       DELTA*(ss_growth_f(RH_num+1)-ss_growth_f(RH_num))
135
136        rho_ss1=rho_ss(rh_num) +
137     .       DELTA*(rho_ss(RH_num+1)-rho_ss(RH_num))             
138c
139        v_stokes=RG*(rho_ss1-zrho(i,k))*      !m/sec
140     .           (mmd_ss*ss_g)*(mmd_ss*ss_g)*
141     .           1.e-12/(18.0*air_visco(i,k)/10.)
142c
143       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
144c
145       CC=1.0+1.257*lambda/(mmd_ss*ss_g)/1.e6  ! C-correction factor
146c
147       v_sed=v_stokes*CC                       ! m/sec !orig
148c
149c---------check for v_sed*dt<zdz
150c
151       IF (v_sed*time_step.GT.zdz(i,k)) THEN
152         v_sed=zdz(i,k)/time_step     
153       ENDIF
154c
155       v_dep_ss(i,k)= v_sed
156       sed_flux(i,k)= tr_seri(i,k,id_coss)*v_sed !g/cm3*m/sec
157       !sed_ss3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
158      ! conc_sed_ss3D(i,k)=sed_flux(i,k)*1.e6      !g/m3*sec !!!!!!!
159c
160       ENDDO          !klon
161       ENDDO          !klev
162c
163!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
164       sed_ss3D(:,:)=0.0  ! initialisation
165     
166       DO k=1, klev
167       DO i=1, klon
168       sed_ss3D(i,k)=sed_ss3D(i,k)-
169     .        sed_flux(i,k)/zdz(i,k) !!!!!!!!!!!!!!!!!!!!!!
170       ENDDO          !klon
171       ENDDO          !klev
172c
173       DO k=1, klev-1
174       DO i=1, klon
175        sed_ss3D(i,k)=sed_ss3D(i,k)+                   
176     .                  sed_flux(i,k+1)/zdz(i,k) !!!!!!!!
177
178       ENDDO          !klon
179       ENDDO          !klev
180
181      DO k = 1, klev
182      DO i = 1, klon
183          tr_seri(i,k,id_coss)=tr_seri(i,k,id_coss)+
184     s   sed_ss3D(i,k)*time_step
185      ENDDO
186      ENDDO
187
188!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189c
190       DO i=1, klon
191         sed_ss(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
192       ENDDO          !klon
193       ELSE
194        DO i=1, klon
195          sed_ss(i)=0.
196        ENDDO
197       ENDIF
198c
199c
200
201c--------- For dust ------------------
202c
203c
204       IF(id_codu>0) THEN
205       DO k=1, klev
206       DO i=1,klon
207c
208        v_stokes=RG*(rho_dust-zrho(i,k))*      !m/sec
209     .           mmd_dust*mmd_dust*
210     .           1.e-12/(18.0*air_visco(i,k)/10.)
211c
212       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
213       CC=1.0+1.257*lambda/(mmd_dust)/1.e6        !dimensionless
214       v_sed=v_stokes*CC                       !m/sec
215c
216c---------check for v_sed*dt<zdz
217c
218       IF (v_sed*time_step.GT.zdz(i,k)) THEN
219         v_sed=zdz(i,k)/time_step     
220       ENDIF
221
222c
223       v_dep_dust(i,k)= v_sed
224       sed_flux(i,k)  = tr_seri(i,k,id_codu)*v_sed !g/cm3.m/sec
225       !sed_dust3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
226c
227       ENDDO          !klon
228       ENDDO          !klev
229
230!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
231       sed_dust3D(:,:)=0.0  ! initialisation
232
233       DO k=1, klev
234       DO i=1, klon
235       sed_dust3D(i,k)=sed_dust3D(i,k)-
236     .                  sed_flux(i,k)/zdz(i,k)
237       ENDDO          !klon
238       ENDDO          !klev
239
240c
241!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242       
243       DO k=1, klev-1
244       DO i=1, klon
245        sed_dust3D(i,k)=sed_dust3D(i,k) +
246     .                  sed_flux(i,k+1)/zdz(i,k)
247       ENDDO          !klon
248       ENDDO          !klev
249c
250      DO k = 1, klev
251      DO i = 1, klon
252         tr_seri(i,k,id_codu)=tr_seri(i,k,id_codu)+
253     s    sed_dust3D(i,k)*time_step
254      ENDDO
255      ENDDO
256
257
258       DO i=1, klon
259         sed_dust(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
260       ENDDO          !klon
261       ELSE
262        DO i=1, klon
263          sed_dust(i)=0.
264        ENDDO
265       ENDIF
266c
267
268
269c--------- For scoarse  dust ------------------
270c
271c
272       IF(id_scdu>0) THEN
273       DO k=1, klev
274       DO i=1,klon
275c
276        v_stokes=RG*(rho_dust-zrho(i,k))*      !m/sec
277     .           mmd_dustsco*mmd_dustsco*
278     .           1.e-12/(18.0*air_visco(i,k)/10.)
279c
280       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
281       CC=1.0+1.257*lambda/(mmd_dustsco)/1.e6        !dimensionless
282       v_sed=v_stokes*CC                       !m/sec
283c
284c---------check for v_sed*dt<zdz
285
286
287       IF (v_sed*time_step.GT.zdz(i,k)) THEN
288         v_sed=zdz(i,k)/time_step
289       ENDIF
290
291c
292       v_dep_dustsco(i,k)= v_sed
293       sed_flux(i,k)     = tr_seri(i,k,id_scdu)*v_sed !g/cm3.m/sec
294       !sed_dustsco3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
295c
296       ENDDO          !klon
297       ENDDO          !klev
298
299!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300       sed_dustsco3D(:,:)=0.0  ! initialisation
301
302       DO k=1, klev
303       DO i=1, klon
304       sed_dustsco3D(i,k)=sed_dustsco3D(i,k)-
305     .                  sed_flux(i,k)/zdz(i,k)
306       ENDDO          !klon
307       ENDDO          !klev
308c
309       DO k=1, klev-1
310       DO i=1, klon
311        sed_dustsco3D(i,k)=sed_dustsco3D(i,k) +
312     .                  sed_flux(i,k+1)/zdz(i,k)
313       ENDDO          !klon
314       ENDDO          !klev
315
316      DO k = 1, klev
317      DO i = 1, klon
318       tr_seri(i,k,id_scdu)=tr_seri(i,k,id_scdu)+
319     s  sed_dustsco3D(i,k)*time_step
320      ENDDO
321      ENDDO
322!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
323
324
325c
326       DO i=1, klon
327         sed_dustsco(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
328       ENDDO          !klon
329       ELSE
330        DO i=1, klon
331          sed_dustsco(i)=0.
332        ENDDO
333       ENDIF
334c
335
336
337
338
339c
340       RETURN
341       END
Note: See TracBrowser for help on using the repository browser.