source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/sediment_mod.F @ 5099

Last change on this file since 5099 was 5099, checked in by abarral, 2 months ago

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line
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"
16      INCLUDE "YOMCST.h"
17      INCLUDE "YOECUMF.h"
18c
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
39c
40c------local variables
41c
42       INTEGER i, k, nbre_RH
43       PARAMETER(nbre_RH=12)
44c
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
56c
57       INTEGER RH_num
58       REAL RH_MAX, DELTA, rh, RH_tab(nbre_RH)
59       PARAMETER (RH_MAX=95.)
60c
61       DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
62c
63c
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/
66c
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/
69c
70c
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
96       rho_dust=2600. !kg/m3
97c
98c--------- Air viscosity (poise=0.1 kg/m-sec)-----------
99c
100       DO k=1, klev
101       DO i=1, klon
102c
103       zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
104c
105       temp=t_seri(i,k)-RTT
106c
107       IF (temp<0.) THEN
108         air_visco(i,k)=(1.718+0.0049*temp-1.2e-5*temp*temp)*1.e-4
109       ELSE
110         air_visco(i,k)=(1.718+0.0049*temp)*1.e-4
111       ENDIF
112c
113       ENDDO
114       ENDDO
115c
116c--------- for Sea Salt -------------------
117c
118c
119c
120       IF(id_coss>0) THEN
121       DO k=1, klev
122       DO i=1,klon
123c
124c---cal. correction factor hygroscopic growth of aerosols
125c
126        rh=MIN(RHcl(i,k)*100.,RH_MAX)
127        RH_num = INT( rh/10. + 1.)
128        IF (rh>85.) RH_num=10
129        IF (rh>90.) RH_num=11
130        DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
131c
132        ss_g=ss_growth_f(rh_num) +
133     .       DELTA*(ss_growth_f(RH_num+1)-ss_growth_f(RH_num))
134
135        rho_ss1=rho_ss(rh_num) +
136     .       DELTA*(rho_ss(RH_num+1)-rho_ss(RH_num))             
137c
138        v_stokes=RG*(rho_ss1-zrho(i,k))*      !m/sec
139     .           (mmd_ss*ss_g)*(mmd_ss*ss_g)*
140     .           1.e-12/(18.0*air_visco(i,k)/10.)
141c
142       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
143c
144       CC=1.0+1.257*lambda/(mmd_ss*ss_g)/1.e6  ! C-correction factor
145c
146       v_sed=v_stokes*CC                       ! m/sec !orig
147c
148c---------check for v_sed*dt<zdz
149c
150       IF (v_sed*time_step>zdz(i,k)) THEN
151         v_sed=zdz(i,k)/time_step     
152       ENDIF
153c
154       v_dep_ss(i,k)= v_sed
155       sed_flux(i,k)= tr_seri(i,k,id_coss)*v_sed !g/cm3*m/sec
156       !sed_ss3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
157      ! conc_sed_ss3D(i,k)=sed_flux(i,k)*1.e6      !g/m3*sec !!!!!!!
158c
159       ENDDO          !klon
160       ENDDO          !klev
161c
162!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
163       sed_ss3D(:,:)=0.0  ! initialisation
164     
165       DO k=1, klev
166       DO i=1, klon
167       sed_ss3D(i,k)=sed_ss3D(i,k)-
168     .        sed_flux(i,k)/zdz(i,k) !!!!!!!!!!!!!!!!!!!!!!
169       ENDDO          !klon
170       ENDDO          !klev
171c
172       DO k=1, klev-1
173       DO i=1, klon
174        sed_ss3D(i,k)=sed_ss3D(i,k)+                   
175     .                  sed_flux(i,k+1)/zdz(i,k) !!!!!!!!
176
177       ENDDO          !klon
178       ENDDO          !klev
179
180      DO k = 1, klev
181      DO i = 1, klon
182          tr_seri(i,k,id_coss)=tr_seri(i,k,id_coss)+
183     s   sed_ss3D(i,k)*time_step
184      ENDDO
185      ENDDO
186
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188c
189       DO i=1, klon
190         sed_ss(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
191       ENDDO          !klon
192       ELSE
193        DO i=1, klon
194          sed_ss(i)=0.
195        ENDDO
196       ENDIF
197c
198c
199
200c--------- For dust ------------------
201c
202c
203       IF(id_codu>0) THEN
204       DO k=1, klev
205       DO i=1,klon
206c
207        v_stokes=RG*(rho_dust-zrho(i,k))*      !m/sec
208     .           mmd_dust*mmd_dust*
209     .           1.e-12/(18.0*air_visco(i,k)/10.)
210c
211       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
212       CC=1.0+1.257*lambda/(mmd_dust)/1.e6        !dimensionless
213       v_sed=v_stokes*CC                       !m/sec
214c
215c---------check for v_sed*dt<zdz
216c
217       IF (v_sed*time_step>zdz(i,k)) THEN
218         v_sed=zdz(i,k)/time_step     
219       ENDIF
220
221c
222       v_dep_dust(i,k)= v_sed
223       sed_flux(i,k)  = tr_seri(i,k,id_codu)*v_sed !g/cm3.m/sec
224       !sed_dust3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
225c
226       ENDDO          !klon
227       ENDDO          !klev
228
229!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
230       sed_dust3D(:,:)=0.0  ! initialisation
231
232       DO k=1, klev
233       DO i=1, klon
234       sed_dust3D(i,k)=sed_dust3D(i,k)-
235     .                  sed_flux(i,k)/zdz(i,k)
236       ENDDO          !klon
237       ENDDO          !klev
238
239c
240!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
241       
242       DO k=1, klev-1
243       DO i=1, klon
244        sed_dust3D(i,k)=sed_dust3D(i,k) +
245     .                  sed_flux(i,k+1)/zdz(i,k)
246       ENDDO          !klon
247       ENDDO          !klev
248c
249      DO k = 1, klev
250      DO i = 1, klon
251         tr_seri(i,k,id_codu)=tr_seri(i,k,id_codu)+
252     s    sed_dust3D(i,k)*time_step
253      ENDDO
254      ENDDO
255
256
257       DO i=1, klon
258         sed_dust(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
259       ENDDO          !klon
260       ELSE
261        DO i=1, klon
262          sed_dust(i)=0.
263        ENDDO
264       ENDIF
265c
266
267
268c--------- For scoarse  dust ------------------
269c
270c
271       IF(id_scdu>0) THEN
272       DO k=1, klev
273       DO i=1,klon
274c
275        v_stokes=RG*(rho_dust-zrho(i,k))*      !m/sec
276     .           mmd_dustsco*mmd_dustsco*
277     .           1.e-12/(18.0*air_visco(i,k)/10.)
278c
279       lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
280       CC=1.0+1.257*lambda/(mmd_dustsco)/1.e6        !dimensionless
281       v_sed=v_stokes*CC                       !m/sec
282c
283c---------check for v_sed*dt<zdz
284
285
286       IF (v_sed*time_step>zdz(i,k)) THEN
287         v_sed=zdz(i,k)/time_step
288       ENDIF
289
290c
291       v_dep_dustsco(i,k)= v_sed
292       sed_flux(i,k)     = tr_seri(i,k,id_scdu)*v_sed !g/cm3.m/sec
293       !sed_dustsco3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
294c
295       ENDDO          !klon
296       ENDDO          !klev
297
298!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
299       sed_dustsco3D(:,:)=0.0  ! initialisation
300
301       DO k=1, klev
302       DO i=1, klon
303       sed_dustsco3D(i,k)=sed_dustsco3D(i,k)-
304     .                  sed_flux(i,k)/zdz(i,k)
305       ENDDO          !klon
306       ENDDO          !klev
307c
308       DO k=1, klev-1
309       DO i=1, klon
310        sed_dustsco3D(i,k)=sed_dustsco3D(i,k) +
311     .                  sed_flux(i,k+1)/zdz(i,k)
312       ENDDO          !klon
313       ENDDO          !klev
314
315      DO k = 1, klev
316      DO i = 1, klon
317       tr_seri(i,k,id_scdu)=tr_seri(i,k,id_scdu)+
318     s  sed_dustsco3D(i,k)*time_step
319      ENDDO
320      ENDDO
321!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
322
323
324c
325       DO i=1, klon
326         sed_dustsco(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
327       ENDDO          !klon
328       ELSE
329        DO i=1, klon
330          sed_dustsco(i)=0.
331        ENDDO
332       ENDIF
333c
334
335
336
337
338c
339       RETURN
340       END
Note: See TracBrowser for help on using the repository browser.