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

Last change on this file since 5284 was 5274, checked in by abarral, 3 days ago

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