source: LMDZ6/trunk/libf/phylmd/Dust/sediment_mod.f90 @ 5248

Last change on this file since 5248 was 5246, checked in by abarral, 7 weeks ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)
File size: 9.5 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 "dimphy.h"
17  INCLUDE "YOMCST.h"
18  INCLUDE "YOECUMF.h"
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  !
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/
67  !
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/
70  !
71  !
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
98  !
99  !--------- Air viscosity (poise=0.1 kg/m-sec)-----------
100  !
101   DO k=1, klev
102   DO i=1, klon
103  !
104   zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
105  !
106   temp=t_seri(i,k)-RTT
107  !
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
113  !
114   ENDDO
115   ENDDO
116  !
117  !--------- for Sea Salt -------------------
118  !
119  !
120  !
121   IF(id_coss>0) THEN
122   DO k=1, klev
123   DO i=1,klon
124  !
125  !---cal. correction factor hygroscopic growth of aerosols
126  !
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))
132  !
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))
138  !
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.)
142  !
143   lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
144  !
145   CC=1.0+1.257*lambda/(mmd_ss*ss_g)/1.e6  ! C-correction factor
146  !
147   v_sed=v_stokes*CC                       ! m/sec !orig
148  !
149  !---------check for v_sed*dt<zdz
150  !
151   IF (v_sed*time_step.GT.zdz(i,k)) THEN
152     v_sed=zdz(i,k)/time_step
153   ENDIF
154  !
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 !!!!!!!
159  !
160   ENDDO          !klon
161   ENDDO          !klev
162  !
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
172  !
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            sed_ss3D(i,k)*time_step
185  ENDDO
186  ENDDO
187
188  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189  !
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
198  !
199  !
200
201  !--------- For dust ------------------
202  !
203  !
204   IF(id_codu>0) THEN
205   DO k=1, klev
206   DO i=1,klon
207  !
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.)
211  !
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
215  !
216  !---------check for v_sed*dt<zdz
217  !
218   IF (v_sed*time_step.GT.zdz(i,k)) THEN
219     v_sed=zdz(i,k)/time_step
220   ENDIF
221
222  !
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 !!!!!!!
226  !
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
240  !
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
249  !
250  DO k = 1, klev
251  DO i = 1, klon
252     tr_seri(i,k,id_codu)=tr_seri(i,k,id_codu)+ &
253           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
266  !
267
268
269  !--------- For scoarse  dust ------------------
270  !
271  !
272   IF(id_scdu>0) THEN
273   DO k=1, klev
274   DO i=1,klon
275  !
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.)
279  !
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
283  !
284  !---------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
291  !
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 !!!!!!!
295  !
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
308  !
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         sed_dustsco3D(i,k)*time_step
320  ENDDO
321  ENDDO
322  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
323
324
325  !
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
334  !
335
336
337
338
339  !
340   RETURN
341END SUBROUTINE sediment_mod
Note: See TracBrowser for help on using the repository browser.