source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/dustemission_mod.F90 @ 2241

Last change on this file since 2241 was 2217, checked in by jescribano, 10 years ago

Bugs corrections. Included a correction/tunning factor for the Chimere-dust emissions, Constant of MB95 equal to 2.61 as in MB95. No spurious increase of u* before horizontal flux calculations in the dust emission scheme. Values of AG00 binding energies fixed as the original AG00 divided by 3 as is Sow et al 2011 ACPD.
File size: 47.2 KB
Line 
1MODULE dustemission_mod
2
3  IMPLICIT NONE
4  !Parameters   
5!  INTEGER, PARAMETER     :: nbins=12  ! number of aerosol bins: original
6!  INTEGER, PARAMETER     :: nbins=800  ! number of aerosol bins: for spla
7!  INTEGER, PARAMETER     :: nbins=8000  ! number of aerosol bins: for spla
8 
9  INTEGER, PARAMETER     :: flag_feff=1 ! 0: deactivate feff (drag partition scheme)
10  INTEGER, PARAMETER     :: nbins=800  ! number of aerosol bins: for spla
11  INTEGER, PARAMETER     :: nmode=3   ! number of soil-dust modes
12  INTEGER, PARAMETER     :: ntyp=5   ! number of soil types
13  INTEGER, PARAMETER     :: nwb=12   ! number of points for the 10m wind
14! speed weibull distribution (>=2)
15  real   ,parameter     :: z10m=1000. !10m in cm
16  REAL,PARAMETER         :: kref=3. !weibull parameter
17  INTEGER, PARAMETER    :: nats=14 !number of mineral types (14 here for sand,
18                                   ! silt, clay etc.)
19  integer, parameter :: nclass=200000
20
21
22  real   , parameter :: dmin=0.0001
23  real   , parameter :: dmax=0.2
24  integer, parameter :: nspe=nmode*3+1
25  real   ,parameter     :: vkarm=0.41
26!JE20150202 : updating scheme to chimere13b <<<
27! original values
28!  integer, parameter :: div1=3.
29!  integer, parameter :: div2=3.
30!  integer, parameter :: div3=3.
31!  real   , parameter :: e1=3.61/div1
32!  real   , parameter :: e2=3.52/div2
33!  real   , parameter :: e3=3.46/div3
34!  real   , parameter :: rop=2.65 ! particle density g/m3
35!  real   , parameter :: roa=0.001227  ! air density g/m3
36!  real   , parameter :: pi=3.14159  !!
37!  real   , parameter :: gravity=981. !! cm!!
38!  real   , parameter :: cd=1.*roa/gravity
39! new values
40!  logical, parameter :: ok_splatunning=.true.
41! Div=3 from S. Alfaro (Sow et al ACPD 2011)
42!JE 20150206
43  integer, parameter :: div1=3.
44  integer, parameter :: div2=3.
45  integer, parameter :: div3=3.
46  real   , parameter :: e1=3.61/div1
47  real   , parameter :: e2=3.52/div2
48  real   , parameter :: e3=3.46/div3
49  real   , parameter :: rop=2.65 ! particle density g/m3
50  real   , parameter :: roa=0.001227  ! air density g/m3
51  real   , parameter :: pi=3.14159  !!
52  real   , parameter :: gravity=981. !! cm!!
53! C=2.61 from Marticorena and Bergametti 1995 instead of Gillete and Chen 2001
54! (recommended C=1.1  in supply-limited dust source area.. )
55   real   , parameter :: cd=2.61*roa/gravity
56!  real   , parameter :: cd=1.0*roa/gravity
57!JE20150202>>>>
58  real,parameter     :: beta=16300.
59  real, parameter, dimension(3) :: diam=(/1.5,6.7,14.2/)
60  INTEGER, PARAMETER     :: ndistb=3
61  real, parameter, dimension(3) :: sig=(/1.7,1.6,1.5/)
62
63!   INTEGER, PARAMETER     :: nbinsHR=3000 !original
64   INTEGER, PARAMETER     :: nbinsHR=30000
65!min and max dust size in um
66      REAL, PARAMETER     :: sizedustmin=0.0599 ! for spla
67      REAL, PARAMETER     :: sizedustmax=63.
68!      REAL, PARAMETER     :: sizedustmin=0.09  ! original
69!      REAL, PARAMETER     :: sizedustmax=63.
70
71
72  ! Calc variables
73  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: massfrac
74  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: binsHR
75  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: binsHRcm
76  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: itv
77  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: sizedust !size dust bin (in um)
78  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: szdcm !the same but (in cm)
79
80  !soil inputs from file donnees_lisa.nc . Should be in the same grid as the
81  !model (regridded with nearest neighborhood from surfnew.nc)
82  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: sol
83  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: P
84  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: zos
85  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: z01
86  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: z02
87  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: D
88  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: A
89  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: solspe
90  INTEGER,DIMENSION(:,:),ALLOCATABLE,SAVE :: masklisa
91!!!  INTEGER,DIMENSION(:),ALLOCATABLE,SAVE :: maskdust
92  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: feff
93  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: feffdbg
94
95  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: sizeclass
96  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: sizeclass2
97  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: uth
98  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: uth2
99  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srel
100  REAL,DIMENSION(:,:),  ALLOCATABLE,SAVE :: srel2
101
102  INTEGER :: nat ! SOL data inside the loop (use as soil type index?)
103  REAL :: ustarsalt
104  REAL :: var3a,var3b
105  INTEGER  :: ns,nd,nsi,npi,ni ! counters
106  INTEGER :: ncl
107
108! outputs
109  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: m1dflux !fluxes for each soil mode
110  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: m2dflux
111  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: m3dflux
112
113
114
115!$OMP THREADPRIVATE(m1dflux)
116!$OMP THREADPRIVATE(m2dflux)
117!$OMP THREADPRIVATE(m3dflux)
118!$OMP THREADPRIVATE(massfrac)
119!$OMP THREADPRIVATE(binsHR)
120!$OMP THREADPRIVATE(binsHRcm)
121!$OMP THREADPRIVATE(itv)
122!$OMP THREADPRIVATE(sizedust)
123!$OMP THREADPRIVATE(szdcm)
124!$OMP THREADPRIVATE(sol)
125!$OMP THREADPRIVATE(P)
126!$OMP THREADPRIVATE(zos)
127!$OMP THREADPRIVATE(z01)
128!$OMP THREADPRIVATE(z02)
129!$OMP THREADPRIVATE(D)
130!$OMP THREADPRIVATE(A)
131!$OMP THREADPRIVATE(solspe)
132!$OMP THREADPRIVATE(masklisa)
133!!!!$OMP THREADPRIVATE(maskdust)
134!$OMP THREADPRIVATE(feff)
135!$OMP THREADPRIVATE(feffdbg)
136!$OMP THREADPRIVATE(sizeclass)
137!$OMP THREADPRIVATE(sizeclass2)
138!$OMP THREADPRIVATE(uth)
139!$OMP THREADPRIVATE(uth2)
140!$OMP THREADPRIVATE(srel)
141!$OMP THREADPRIVATE(srel2)
142
143
144
145  CONTAINS
146
147!--------------------------------------------------------------------------------------
148!======================================================================================
149!**************************************************************************************
150!======================================================================================
151!--------------------------------------------------------------------------------------
152
153  SUBROUTINE dustemission( debutphy, xlat, xlon, &    !Input
154                          pctsrf,zu10m,zv10m,wstar, & !Input
155                          ale_bl,ale_wake, &          !Input
156                          param_wstarBL, param_wstarWAKE, &  !Input
157                          emdustacc,emdustcoa,emdustsco,maskdust)    !Output
158  USE dimphy
159  USE infotrac
160  USE write_field_phy
161  USE mod_grid_phy_lmdz
162  USE mod_phys_lmdz_para
163
164  IMPLICIT NONE
165  ! input :
166  ! output: flux_sparam_ddfine,flux_sparam_ddcoa,
167  ! first:
168  ! Model grid parameters
169  REAL,DIMENSION(klon),     INTENT(IN)     :: xlat
170  REAL,DIMENSION(klon),     INTENT(IN)     :: xlon
171  REAL,DIMENSION(klon),     INTENT(IN)     :: pctsrf
172  REAL,DIMENSION(klon),INTENT(IN)          :: zu10m   ! 10m zonal wind
173  REAL,DIMENSION(klon),INTENT(IN)          :: zv10m   ! meridional 10m wind
174  REAL,DIMENSION(klon),INTENT(IN)          :: wstar
175  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
176  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
177  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarWAKE
178  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
179 
180 
181  LOGICAL  :: debutphy ! First physiqs run or not
182  ! Intermediate variable: 12 bins emissions
183  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbinloc ! vertical emission fluxes
184
185  !OUT variables
186  REAL,DIMENSION(klon)  :: emdustacc,emdustcoa,emdustsco ! emission in spla dust bins:
187! old radio : acc=0.03-0.5 micrometers, ccoa:0.5-10 micrometers
188! new acc=0.03-0.5 micrometers, coa:0.5-3 micrometers ,sco:3-15 um
189  INTEGER,DIMENSION(klon) :: maskdust ! where the emissions were calculated
190!  INTEGER,DIMENSION(klon_glo) :: maskdust_glo ! auxiliar
191!  REAL,DIMENSION(klon_glo) :: raux_klon_glo ! auxiliar
192
193!$OMP THREADPRIVATE(emisbinloc)
194!!!!!!$OMP THREADPRIVATE(maskdust)
195  IF (debutphy) THEN
196     ALLOCATE( emisbinloc(klon,nbins) )
197  ENDIF
198
199  IF( debutphy ) THEN 
200     CALL initdust(xlat,xlon,pctsrf)
201  ENDIF
202
203!JE20141124  CALL  calcdustemission(debutphy,zu10m,zv10m,wstar,ale_bl,ale_wake,emisbinloc)
204  CALL  calcdustemission(debutphy,zu10m,zv10m,wstar,ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & !I
205                         emisbinloc)   !O
206
207  CALL makemask(maskdust)
208
209  IF( debutphy ) THEN 
210!       call gather(maskdust,maskdust_glo)
211!     !$OMP MASTER
212!     IF (is_mpi_root .AND. is_omp_root) THEN
213!       CALL writefield_phy("maskdust",float(maskdust_glo),1)
214       CALL writefield_phy("maskdust",float(maskdust),1)
215!     ENDIF
216!     !$OMP END MASTER
217!     !$OMP BARRIER
218  ENDIF
219
220  CALL adaptdustemission(debutphy,emisbinloc,emdustacc,emdustcoa,emdustsco)
221  ! output in kg/m2/s
222
223
224  END SUBROUTINE dustemission
225
226!--------------------------------------------------------------------------------------
227!======================================================================================
228!**************************************************************************************
229!======================================================================================
230!--------------------------------------------------------------------------------------
231
232 SUBROUTINE makemask(maskdustloc)
233  USE dimphy
234  USE infotrac
235  IMPLICIT NONE
236  !Input
237  INTEGER,DIMENSION(klon) :: maskdustloc
238  INTEGER :: i,j,k
239  integer :: iaux
240 
241
242  do k=1,klon
243      maskdustloc(k)=0
244      do i=1,ntyp
245         if (masklisa(k,i)>0) then
246             maskdustloc(k)=1
247         endif
248      enddo
249  enddo
250
251  END SUBROUTINE makemask
252
253!--------------------------------------------------------------------------------------
254!======================================================================================
255!**************************************************************************************
256!======================================================================================
257!--------------------------------------------------------------------------------------
258
259 SUBROUTINE adaptdustemission(debutphy,emisbinlocal, &
260                               emdustacc,emdustcoa,emdustsco)
261
262  USE dimphy
263  USE infotrac
264  USE write_field_phy
265  USE mod_grid_phy_lmdz
266  USE mod_phys_lmdz_para
267
268  IMPLICIT NONE
269  !Input
270  REAL,DIMENSION(klon) :: emdustacc,emdustcoa,emdustsco
271  REAL, DIMENSION(klon,nbins) :: emisbinlocal
272  ! Local
273  INTEGER  :: i,j,k
274!!!  INTEGER,DIMENSION(2) :: iminacclow,iminacchigh,imincoalow,imincoahigh ! in
275!case of small nbins.... not ready
276  INTEGER,SAVE ::iminacclow,iminacchigh,imincoalow
277  INTEGER,SAVE ::imincoahigh,iminscohigh,iminscolow
278!  real,parameter :: sizeacclow=0.03
279!  real,parameter :: sizeacchigh=0.5
280!  real,parameter :: sizecoalow=0.5
281!  real,parameter :: sizecoahigh=10.  ! in micrometers
282  real,parameter :: sizeacclow=0.06
283  real,parameter :: sizeacchigh=1.0
284  real,parameter :: sizecoalow=1.0
285  real,parameter :: sizecoahigh=6.  !20 ! diameter in micrometers
286  real,parameter :: sizescolow=6.
287  real,parameter :: sizescohigh=30.  ! in micrometers
288
289  real,parameter :: tunningfactor=4.5  ! factor for fine bins!!! important!!
290!JE20140915  real,parameter :: sizeacclow=0.06
291!JE20140915  real,parameter :: sizeacchigh=1.0
292!JE20140915  real,parameter :: sizecoalow=1.0
293!JE20140915  real,parameter :: sizecoahigh=10.  !20 ! diameter in micrometers
294!JE20140915  real,parameter :: sizescolow=10.
295!JE20140915  real,parameter :: sizescohigh=30.  ! in micrometers
296
297
298
299  logical ::  debutphy
300  real :: diff, auxr1,auxr2,auxr3,auxr4
301  real,dimension(klon,nbins) :: itvmean
302  real,dimension(klon,nbins+1) :: itv2
303!  real,dimension(klon_glo,nbins) :: itvmean_glo
304!  real,dimension(:,:) , allocatable  :: itvmean_glo
305!  real,dimension(:,:), allocatable :: itv2_glo
306 
307  integer, save :: counter,counter1 !dbg
308  REAL, DIMENSION(:,:),ALLOCATABLE,SAVE :: emisbinlocalmean,emisbinlocalmean2 !dbg
309  REAL, DIMENSION(:,:),ALLOCATABLE :: emisbinlocalmean2_glo
310  logical :: writeaerosoldistrib
311!$OMP THREADPRIVATE(iminacclow,iminacchigh,imincoalow,imincoahigh)
312
313writeaerosoldistrib=.false.
314if (debutphy) then
315
316  if (sizedustmin>sizeacclow .or. sizedustmax<sizescohigh) then
317   call abort_gcm('adaptdustemission', 'Dust range problem',1)
318  endif
319  print *,'ALL DUST BIN: tunning EMISSION factor= ',tunningfactor
320  auxr1=9999.
321  auxr2=9999.
322  auxr3=9999.
323  auxr4=9999.
324  do i=1,nbins+1
325   if (abs(sizeacclow-itv(i))<auxr1) then
326          auxr1=abs( sizeacclow-itv(i))
327          iminacclow=i
328   endif
329   if (abs(sizeacchigh-itv(i))<auxr2) then
330          auxr2=abs( sizeacchigh-itv(i))
331          iminacchigh=i
332          imincoalow=i
333   endif
334   if (abs(sizecoahigh-itv(i))<auxr3) then
335          auxr3=abs( sizecoahigh-itv(i))
336          imincoahigh=i
337          iminscolow=i
338   endif
339   if (abs(sizescohigh-itv(i))<auxr4) then
340          auxr4=abs( sizescohigh-itv(i))
341          iminscohigh=i
342   endif
343  enddo
344if (writeaerosoldistrib) then
345!JEdbg<<
346  do j=1,klon
347  do i=1,nbins
348    itvmean(j,i)=(itv(i)+itv(i+1))/2.
349    itv2(j,i)=itv(i)
350    !print*, itv(i),itvmean(i),itv(i+1)
351    !print*, sizedust(i)
352  enddo
353  itv2(j,nbins+1)=itv(nbins+1)
354  enddo
355!  allocate(itv2_glo(klon_glo,nbins+1))
356!  allocate(itvmean_glo(klon_glo,nbins))
357!  ALLOCATE(emisbinlocalmean2_glo(klon_glo,nbins))
358!  call gather(itv2,itv2_glo)
359!  call gather(itvmean,itvmean_glo)
360!!$OMP MASTER
361!  IF (is_mpi_root .AND. is_omp_root) THEN
362!  CALL writefield_phy("itv2",itv2_glo,nbins+1)
363!  CALL writefield_phy("itvmean",itvmean_glo,nbins)
364  CALL writefield_phy("itv2",itv2,nbins+1)
365  CALL writefield_phy("itvmean",itvmean,nbins)
366!  ENDIF
367!!$OMP END MASTER
368!!$OMP BARRIER
369  ALLOCATE(emisbinlocalmean(klon,nbins))
370  ALLOCATE(emisbinlocalmean2(klon,nbins))
371  do i=1,nbins
372    do j=1,klon
373       emisbinlocalmean(j,i)=0.0
374       emisbinlocalmean2(j,i)=0.0
375    enddo
376  enddo
377  counter=1
378  counter1=0
379!JEdbg>>
380endif
381endif
382
383
384!print *,'JE'
385!print *,iminacclow,iminacchigh,imincoalow,imincoahigh
386
387! estimate and integrate bins into only accumulation and coarse
388
389
390do k=1,klon
391auxr1=0.0
392auxr2=0.0
393auxr3=0.0
394  do i=iminacclow,iminacchigh-1
395   auxr1=auxr1+emisbinlocal(k,i)
396  enddo
397  emdustacc(k)= auxr1*tunningfactor
398  do i=imincoalow,imincoahigh-1
399    auxr2=auxr2+emisbinlocal(k,i)
400  enddo
401  emdustcoa(k)=auxr2*tunningfactor
402  do i=iminscolow,iminscohigh-1
403    auxr3=auxr3+emisbinlocal(k,i)
404  enddo
405  emdustsco(k)=auxr3*tunningfactor
406enddo
407
408!JEdbg<<
409if (writeaerosoldistrib) then
410  do i=1,nbins
411    do j=1,klon
412        emisbinlocalmean(j,i)=emisbinlocalmean(j,i)+emisbinlocal(j,i)
413    enddo
414  enddo
415counter1=counter1+1
416! 4 timesteps of physics=4*15min=1hour..
417! 1440 = 15 days
418! 480 = 5 days
419if (MOD(counter,1440).eq. 0) THEN
420   !if (MOD(counter,480).eq. 0) THEN
421   do k = 1,klon
422   do i=1,nbins
423     emisbinlocalmean2(k,i)=emisbinlocalmean(k,i)/float(counter1)
424   enddo
425   enddo
426   counter1=0
427!   call gather(emisbinlocalmean2,emisbinlocalmean2_glo)
428!!$OMP MASTER
429!   IF (is_mpi_root .AND. is_omp_root) THEN
430!   CALL writefield_phy("emisbinlocalmean2",emisbinlocalmean2_glo,nbins)
431   CALL writefield_phy("emisbinlocalmean2",emisbinlocalmean2,nbins)
432!   ENDIF
433!!$OMP END MASTER
434!!$OMP BARRIER
435   do i=1,nbins
436     do j=1,klon
437        emisbinlocalmean(j,i)=0.0
438     enddo
439   enddo
440endif
441counter=counter+1
442endif
443!JEdbg>>
444
445!CALL abort_gcm('calcdustemission', 'OK1',1)
446 
447END SUBROUTINE adaptdustemission
448
449
450!--------------------------------------------------------------------------------------
451!======================================================================================
452!**************************************************************************************
453!======================================================================================
454!--------------------------------------------------------------------------------------
455
456
457!--------------------------------------------------------
458  SUBROUTINE initdust(xlat,xlon,pctsrf)
459  USE dimphy
460  USE infotrac
461  USE write_field_phy
462  USE mod_grid_phy_lmdz
463  USE mod_phys_lmdz_para
464
465  IMPLICIT NONE
466  !Inputs
467  REAL,DIMENSION(klon),     INTENT(IN)     :: xlat
468  REAL,DIMENSION(klon),     INTENT(IN)     :: xlon
469  REAL,DIMENSION(klon),     INTENT(IN)     :: pctsrf
470
471  !Local
472  REAL,DIMENSION(klon,ntyp)                :: solini
473  REAL,DIMENSION(klon,ntyp)                :: Pini
474  REAL,DIMENSION(klon,ntyp)                :: zosini
475  REAL,DIMENSION(klon,ntyp)                :: z01ini
476  REAL,DIMENSION(klon,ntyp)                :: z02ini
477  REAL,DIMENSION(klon,ntyp)                :: Dini
478  REAL,DIMENSION(klon,ntyp)                :: Aini
479
480  REAL,DIMENSION(nclass)                   :: newstep
481
482 !TEMPORAL/ MAKE MASK!!!
483  REAL, PARAMETER        :: longmin=-20.
484  REAL, PARAMETER        :: longmax=77.          !35. 
485  REAL, PARAMETER        :: latmin=10.
486  REAL, PARAMETER        :: latmax=40.
487 !TEMPORAL/ MAKE MASK!!!
488  !Local
489  REAL, PARAMETER        :: eps=1.e-5
490  REAL, PARAMETER        :: aeff=0.35
491  REAL, PARAMETER        :: xeff=10.
492  REAL, PARAMETER        :: trctunt=0.310723
493 
494  INTEGER     :: i,j,k,nts
495  REAL :: dp,dstep
496! Parametres pour le calcul de la repartition de l energie feff(klon,ntyp)
497  REAL :: aa,bb,cc,dd,ee,ff
498  REAL,DIMENSION(klon,ntyp) :: tmp1
499  REAL :: p3t,var1,var2,et
500
501! Parametres pour le calcul de la surface couverte par les particule de diametre
502! dp srel(nats,nclass)
503  REAL :: stotale,su,su_loc,xk,xl,xm,xn
504  REAL,DIMENSION(nclass)                   :: subsoildist
505
506! isolog and massfrac calcs
507  INTEGER :: nbinsout,nb,miniso,nbins1,nbins2,istart,no
508  REAL :: b1,b2,stepbin,minisograd,exden,d2min,d2max,numin,numax
509  REAL                   :: delta1,delta2,ldmd
510!  REAL, PARAMETER     :: sizedustmin=0.09
511  REAL,DIMENSION(nbinsHR):: binsISOGRAD
512  REAL,DIMENSION(nbinsHR):: vdISOGRAD
513  REAL,DIMENSION(nbinsHR):: logvdISOGRAD
514  REAL                   :: curvd
515  REAL                   :: avdhr
516  REAL                   :: bvdhr
517  REAL,DIMENSION(nbins)  :: dmn1
518  REAL,DIMENSION(nbins)  :: dmn2
519  REAL,DIMENSION(nbins)  :: dmn3
520  REAL,DIMENSION(nbinsHR) :: vdHR
521  REAL,DIMENSION(nbinsHR)  :: vdout
522
523
524
525! Parametres pour le calcul de la vitesse de friction seuil uth(nclass)
526! calcul de la vitesse seuil selon la parametrisation de Shao and Lu
527! 2000(icuth=1).
528!  INTEGER                :: ich1
529  REAL                   :: an
530  REAL                   :: gam
531  REAL                   :: sigshao
532  REAL                   :: x1
533  REAL                   :: x2
534! Cas de Iversen and White 1982 (icuth=0) si necessaire.
535  REAL, PARAMETER        :: adust=1331.647
536  REAL, PARAMETER        :: bdust=0.38194
537  REAL, PARAMETER        :: xdust=1.561228
538  REAL, PARAMETER        :: ddust=0.006/(rop*gravity)
539
540  CHARACTER(LEN=10) :: varname
541  CHARACTER(LEN=80) :: fnsolspe
542  CHARACTER(LEN=200) :: line
543 
544
545! nats: number of mineral types (14 here for sand, silt, clay etc.)
546    ALLOCATE( binsHR(nbinsHR) )
547    ALLOCATE( binsHRcm(nbinsHR) )
548    ALLOCATE( itv(nbins+1) )
549    ALLOCATE( sizedust(nbins) )
550    ALLOCATE( szdcm(nbins) )
551    ALLOCATE( massfrac(ndistb,nbins) )
552    ALLOCATE( sol(klon,ntyp) )
553    ALLOCATE( P(klon,ntyp) )
554    ALLOCATE( zos(klon,ntyp) )
555    ALLOCATE( z01(klon,ntyp) )
556    ALLOCATE( z02(klon,ntyp) )
557    ALLOCATE( D(klon,ntyp) )
558    ALLOCATE( A(klon,ntyp) )
559    ALLOCATE( masklisa(klon,ntyp) )
560    ALLOCATE( solspe(nats,nspe) )
561    ALLOCATE( feff(klon,ntyp) )
562    ALLOCATE( feffdbg(klon,ntyp) )
563    ALLOCATE( sizeclass(nclass) )
564    ALLOCATE( sizeclass2(nclass) )
565    ALLOCATE( uth(nclass) )
566    ALLOCATE( uth2(nclass) )
567    ALLOCATE( srel(nats,nclass) )
568    ALLOCATE( srel2(nats,nclass) )
569    ALLOCATE( m1dflux(klon) )
570    ALLOCATE( m2dflux(klon) )
571    ALLOCATE( m3dflux(klon) )
572
573
574
575  ! read input data from donnees_lisa.nc
576  varname='SOL'
577  CALL read_surface(varname,solini)
578  varname='P'
579  CALL read_surface(varname,Pini)
580  varname='ZOS_'
581  CALL read_surface(varname,zosini)
582  varname='ZO1_'
583  CALL read_surface(varname,z01ini)
584  varname='ZO2_'
585  CALL read_surface(varname,z02ini)
586  varname='D'
587  CALL read_surface(varname,Dini)
588  varname='A'
589  CALL read_surface(varname,Aini)
590print *,'beforewritephy',mpi_rank,omp_rank
591  CALL writefield_phy("SOLinit",solini,5)
592  CALL writefield_phy("Pinit",Pini,5)
593  CALL writefield_phy("ZOSinit",zosini,5)
594  CALL writefield_phy("ZO1init",z01ini,5)
595  CALL writefield_phy("ZO2init",z02ini,5)
596  CALL writefield_phy("Dinit",Dini,5)
597  CALL writefield_phy("Ainit",Aini,5)
598print *,'afterwritephy',mpi_rank,omp_rank
599
600  DO i=1,klon
601    DO nts=1,ntyp
602        masklisa(i,nts)=0
603        IF(Pini(i,nts)>=0.) THEN
604               masklisa(i,nts)=1
605        ENDIF
606    ENDDO
607  ENDDO
608!print *,'JEOK1',mpi_rank,omp_rank
609  DO i=1,klon
610    !print *,Pini(i,1),Pini(i,2),Pini(i,3),Pini(i,4),Pini(i,5)
611    DO nts=1,ntyp
612      !IF(xlon(i).ge.longmin.and.xlon(i).le.longmax.and. &
613      !&      xlat(i).ge.latmin.and.xlat(i).le.latmax    &
614      !&      .and.pctsrf(i)>0.5.and.Pini(i,nts)>0.)THEN
615      IF(pctsrf(i)>0.5.and.Pini(i,nts)>0.)THEN
616         
617           sol(i,nts) = solini(i,nts)
618             P(i,nts) = Pini(i,nts)
619           zos(i,nts) = zosini(i,nts)
620           z01(i,nts) = z01ini(i,nts)
621           z02(i,nts) = z02ini(i,nts)
622             D(i,nts) = Dini(i,nts)
623             A(i,nts) = Aini(i,nts)
624!             masklisa(i,nts)=1
625      ELSE
626              sol(i,nts) =0.0
627               P(i,nts)  =0.0
628              zos(i,nts) =0.0
629              z01(i,nts) =0.0
630              z02(i,nts) =0.0
631               D(i,nts)  =0.0
632               A(i,nts)  =0.0
633!            masklisa(i,nts)=0
634      ENDIF
635    ENDDO
636  ENDDO
637
638! print *,'JEOK2',mpi_rank,omp_rank
639if ( 1==1 ) then
640
641! print *,'JEOK4',mpi_rank,omp_rank
642   CALL writefield_phy("SOL",sol,5)
643   CALL writefield_phy("P"  ,P  ,5)
644   CALL writefield_phy("ZOS",zos,5)
645   CALL writefield_phy("ZO1",z01,5)
646   CALL writefield_phy("ZO2",z02,5)
647   CALL writefield_phy("D"  ,D  ,5)
648   CALL writefield_phy("A"  ,A  ,5)
649   CALL writefield_phy("masklisa",float(masklisa),5)
650   CALL writefield_phy("pctsrf",pctsrf,1)
651   CALL writefield_phy("xlon",xlon,1)
652   CALL writefield_phy("xlat",xlat,1)
653!print *,'JEOK5',mpi_rank,omp_rank
654!print *,'JEOK6',mpi_rank,omp_rank
655
656endif
657
658  !CALL abort_gcm('initdustemission', 'OK1',1)
659
660! Lecture des donnees du LISA indiquant le type de sol utilise.
661       ! in lines: landuse
662! nats: number of mineral types (14 here for sand, silt, clay etc.)
663! data format in columns
664! mmd1  sigma1  p1  mmd2  sigma2  p2 mmd3 ... alpha
665!print *,'JEOK7',mpi_rank,omp_rank
666!$OMP MASTER
667IF (is_mpi_root .AND. is_omp_root) THEN
668!print *,'JEOK9',mpi_rank,omp_rank
669 fnsolspe='SOILSPEC.data'
670  PRINT*,'  o Reading ',fnsolspe(1:40)
671  OPEN(10,file=fnsolspe)
672  READ(10,*)line
673  DO i=1,nats
674     READ(10,*)line
675     READ(10,*)(solspe(i,j),j=1,nspe)
676!JE     alfa(i)=solspe(i,10)
677!     PRINT*,'i,alfa(i)',i,alfa(i)
678!     WRITE(18,*)i,alfa(i)
679  END DO
680!     print*,'solspe(14,10)= ',solspe(14,10)
681  CLOSE(10)
682ENDIF
683!$OMP END MASTER
684!$OMP BARRIER
685!print *,'JEOK10',mpi_rank,omp_rank
686  call bcast(solspe)
687! Calcul de la distribution en taille des particules de Dust
688! Elle depend du nombre de classe des particules nclass.
689!c making full resolution soil diameter table
690      dstep=0.0005
691      dp=dmin
692      do i=1,nclass
693         dp=dp*exp(dstep)
694         sizeclass(i)=dp
695         if(dp.ge.dmax+eps)goto 30
696         newstep(i)=dstep
697        ! WRITE(18,*)i,sizeclass(i)
698      enddo
69930   continue
700      print*,'IK5'
701      ncl=i-1
702      print*,'   soil size classes used   ',ncl,' / ',nclass
703      print*,'   soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl)
704      if(ncl.gt.nclass)stop
705
706! Threshold velocity:
707if (.false.) then
708!if (.true.) then
709!c 0: Iversen and White 1982
710       print *,'Using  Iversen and White 1982 Uth'
711         do i=1,ncl
712            bb=adust*(sizeclass(i)**xdust)+bdust
713            cc=sqrt(1+ddust*(sizeclass(i)**(-2.5)))
714            xk=sqrt(abs(rop*gravity*sizeclass(i)/roa))
715            if (bb.lt.10.) then
716               dd=sqrt(1.928*(bb**0.092)-1.)
717               uth(i)=0.129*xk*cc/dd
718            else
719               ee=-0.0617*(bb-10.)
720               ff=1.-0.0858*exp(ee)
721               uth(i)=0.12*xk*cc*ff
722            endif
723         enddo
724endif
725if(.true.) then
726! 1: Shao and Lu 2000
727       print *,'Using  Shao and Lu 2000 Uth'
728            an=0.0123
729            gam=0.3
730            do i=1,ncl
731               sigshao=rop/roa
732               x1=sigshao*gravity*sizeclass(i)
733               x2=gam/(roa*sizeclass(i))
734               uth(i)=sqrt(an*(x1+x2))
735            enddo
736endif
737
738
739
740!Calcul de la surface couverte par les particules de diametre dp
741!srel(nats,nclass)
742
743!  OPEN(18,file='srel.dat',form='formatted',access='sequential')
744  DO ns=1,nats
745     stotale=0.
746     DO i=1,ncl
747        su=0.
748        DO j=1,nmode
749           nd=((j-1)*3)+1
750           nsi=((j-1)*3)+2
751           npi=((j-1)*3)+3
752           IF (solspe(ns,nd).EQ.0.)THEN
753              su_loc=0.
754           ELSE
755              xk=solspe(ns,npi)/(sqrt(2.*pi)*log(solspe(ns,nsi)))
756              xl=((log(sizeclass(i))-log(solspe(ns,nd)))**2) &
757     &              /(2.*(log(solspe(ns,nsi)))**2)
758              xm=xk*exp(-xl)
759              xn=rop*(2./3.)*(sizeclass(i)/2.)
760              su_loc=xm*newstep(i)/xn
761           END IF
762           su=su+su_loc
763        END DO
764        subsoildist(i)=su
765        stotale=stotale+su
766     END DO
767     DO i=1,ncl
768        IF (subsoildist(i).gt.0..and.stotale.gt.0.)THEN
769            srel(ns,i)=subsoildist(i)/stotale
770
771        ELSE
772            srel(ns,i)=0.
773        END IF
774!        WRITE(18,*)i,srel(1,i)
775     END DO
776!     PRINT*,'ns , srel(10000) ',ns,srel(ns,10000)
777  END DO
778
779
780! Calcul de la repartition de l energie feff(klon,ntyp)
781
782  !efficient fraction for the friction velocities as a function
783  !of z0s, Zo1 et Zo2 to retrieve the apparent threshold
784  !wind friction velocity.
785      !    feff(:,:)=0.
786       do i=1,klon
787          do k=1,ntyp
788     !     print*,'IKKK ',i,klon,k,ntyp
789             if (zos(i,k).eq.0..or.z01(i,k).eq.0.) then
790     !       if (zos(i,k)<=0..or.z01(i,k)<=0.) then
791!              if (zos(i,k)<0..or.z01(i,k)<0.) then
792     !            print*,'INI DUST WARNING zos ou z01<0',zos(i,k),z01(i,k)
793!              endif
794              feff(i,k)=0.
795              feffdbg(i,k)=0.
796 !         print*,'IKKK A ',i,klon,k,ntyp
797            else
798! drag partition betzeen the erodable surface and zo1
799     !     print*,'IKKK B0 ',i,klon,k,ntyp,z01(i,k),zos(i,k),xeff,aeff
800              aa=log(z01(i,k)/zos(i,k))
801              tmp1(i,k)=aa
802              bb=log(aeff*(xeff/zos(i,k))**0.8)
803              cc=1.-aa/bb
804              feffdbg(i,k)=cc
805       !   print*,'IKKK B1 ',i,klon,k,ntyp
806! drag partition between zo1 and zo2
807! feff: total efficient fraction
808              if(D(i,k).eq.0.)then
809                 feff(i,k)=cc
810   !       print*,'IKKK C ',i,klon,k,ntyp
811              else
812                 dd=log(z02(i,k)/z01(i,k))
813                 ee=log(aeff*(D(i,k)/z01(i,k))**0.8)
814                 feff(i,k)=(1.-dd/ee)*cc
815   !       print*,'IKKK D ',i,klon,k,ntyp
816              endif
817              if (feff(i,k).lt.0.)feff(i,k)=0.
818              if (feffdbg(i,k).lt.0.)feffdbg(i,k)=0.
819              if (feff(i,k).gt.1.)feff(i,k)=1.
820              if (feffdbg(i,k).gt.1.)feffdbg(i,k)=1.
821    !      print*,'IKKK E ',i,klon,k,ntyp
822            endif
823          enddo
824        enddo
825! JE20150120<<
826  if (flag_feff .eq. 0) then
827    print *,'JE_dbg FORCED deactivated feff'
828    do i=1,klon
829      do k=1,ntyp
830        feff(i,k)=1.
831      enddo
832    enddo
833  endif
834! JE20150120>>
835
836
837if (1==1) then
838!  !  CALL writefield_phy("AA",tmp1(1:klon,1:5),5)
839!
840    CALL writefield_phy("REPART5",feff(1:klon,1:5),5)
841    CALL writefield_phy("REPART5dbg",feffdbg(1:klon,1:5),5)
842endif
843
844
845!c---------------FOR def_prepag01modif(var3a,var3b)-----------------------
846        p3t=0.0001
847        var1=e3**(1./3.)
848        var2=(rop*pi)**(1./3.)
849        var3a=trctunt*var1/var2
850        et=e1+(sqrt(p3t*(e3*e3+e1*e2-e2*e3-e1*e3))/p3t)
851        var1=et**(1./3.)
852        var2=(rop*pi)**(1./3.)
853        var3b=trctunt*var1/var2
854
855! JE20140902: comment all isograd distributions and make own massfrac:
856
857
858!  if (.false.) then
859!!**************L718
860!
861!!c------------------------------------------------------------------------
862!! isolog distrib and massfrac calculations.
863!
864!      nbinsout=nbins+1
865!      b1=log(sizedustmin)
866!      b2=log(sizedustmax)
867!! restricted ISOLOG bins distributions
868!
869!!      step=(b2-b1)/(nbinsout-1)
870!!      DO ni=1,nbinsout
871!!         itv(ni)=exp(b1+(ni-1)*step)
872!!      ENDDO
873!!OPEN(18,file='vdepothrsizbin.dat',form='formatted',access='sequential')
874!! Restricted ISOGRADIENT bins distributions
875!!!!!!!Making HIGH RESOLUTION dust size distribution (in um)!!!!!!
876!  stepbin=(b2-b1)/(nbinsHR-1)
877!  DO nb=1,nbinsHR
878!     binsHR(nb)=exp(b1+(nb-1)*stepbin)
879!  END DO
880!
881!  DO nb=1,nbinsHR
882!     binsHRcm(nb)=binsHR(nb)*1.e-4
883!  END DO
884!! Making HIGH RESOLUTION dry deposition velocity
885!     CALL calcvd(vdout)
886!
887!
888! DO nb=1,nbinsHR
889!     vdHR(nb)=vdout(nb)
890!!  WRITE(18,*),binsHR(nb),vdHR(nb)
891!  END DO
892!
893!   !searching for minimum value of dry deposition velocity
894!  minisograd=1.e20
895!  DO nb=1,nbinsHR
896!     IF(vdHR(nb).le.minisograd)THEN
897!       minisograd=vdHR(nb)
898!       miniso=nb
899!     END IF
900!  END DO
901!
902!! searching for optimal number of bins in positive slope Vd part
903!
904!  nbins1=1
905!  nbins2=nbinsout-1
906!     DO k=1,1000
907!        delta1=abs(log(vdHR(1))-log(vdHR(miniso)) )/nbins1
908!        delta2=abs(log(vdHR(nbinsHR))-log(vdHR(miniso)))/(nbins2-1)
909!        IF(delta2.GE.delta1)THEN
910!        GOTO 50
911!
912!        ELSE
913!           nbins2=nbins2-1
914!           nbins1=nbins1+1
915!        END IF
916!        IF(nbins2.EQ.1)THEN
917!           PRINT*,'!! warning: lower limit was found: STOP'
918!           STOP
919!        END IF
920!     END DO
921! 50  CONTINUE
922!print*,'IK10'
923!! building the optimized distribution
924!  logvdISOGRAD(1)=log(vdHR(1))
925!  binsISOGRAD(1)=binsHR(1)
926!  DO k=2,nbins1
927!     logvdISOGRAD(k)=logvdISOGRAD(1)-(k-1)*delta1
928!  END DO
929!
930!  logvdISOGRAD(nbins1+1)=log(minisograd)
931!
932!  DO k=1,nbins2
933!     logvdISOGRAD(nbins1+1+k)=logvdISOGRAD(nbins1+1)+k*delta2
934!  END DO
935!
936!  DO k=1,nbinsout
937!     vdISOGRAD(k)=exp(logvdISOGRAD(k))
938!  END DO
939!! sequential search of the correspondig particle diameter
940!  istart=1
941!  DO k=2,nbinsout-1
942!     curvd=vdISOGRAD(k)
943!     DO nb=istart,nbinsHR
944!        avdhr=vdHR(nb)
945!        bvdhr=vdHR(nb+1)
946!        IF(nb.LT.(miniso-1).AND.curvd.LT.avdhr.AND. &
947!           curvd.GT.bvdhr)THEN
948!           binsISOGRAD(k)=binsHR(nb)
949!           istart=nb
950!           GOTO 60
951!        END IF
952!        IF(nb.eq.miniso)THEN
953!           binsISOGRAD(k)=binsHR(nb)
954!           istart=nb+1
955!           GOTO 60
956!        END IF
957!        IF(nb.GT.miniso.AND.curvd.GT.avdhr.AND. &
958!                              curvd.LT.bvdhr)THEN
959!           binsISOGRAD(k)=binsHR(nb)
960!           istart=nb
961!           GOTO 60
962!        END IF
963!     END DO
964! 60  CONTINUE
965!  END DO
966!print*,'IK11'
967!  binsISOGRAD(nbinsout)=binsHR(nbinsHR)
968!  vdISOGRAD(nbinsout)=vdHR(nbinsHR)
969!  DO nb=1,nbinsout
970!     itv(nb)=binsISOGRAD(nb)
971!  END DO
972!  end ! JE20140902 if false
973
974! Making dust size distribution (in um)
975!
976      nbinsout=nbins+1
977      b1=log(sizedustmin)
978      b2=log(sizedustmax)
979      stepbin=(b2-b1)/(nbinsout-1)
980     DO ni=1,nbinsout
981        itv(ni)=exp(b1+(ni-1)*stepbin)
982     ENDDO
983
984!  stepbin=(b2-b1)/(nbinsHR-1)
985!  DO nb=1,nbinsHR
986!     binsHR(nb)=exp(b1+(nb-1)*stepbin)
987!  END DO
988!
989!  DO nb=1,nbinsHR
990!     binsHRcm(nb)=binsHR(nb)*1.e-4
991!  END DO
992
993
994
995      DO nb=1,nbins
996         ldmd=(log(itv(nb+1))-log(itv(nb)))/2.
997         sizedust(nb)=exp(log(itv(nb))+ldmd)
998        PRINT*,nb, itv(nb),' < ',sizedust(nb),' < ',itv(nb+1)
999      ENDDO
1000! Making dust size distribution (in cm) ???
1001      DO nb=1,nbins
1002         szdcm(nb)=sizedust(nb)*1.e-4
1003      ENDDO
1004!c preparation of emissions reaffectation.
1005      DO k=1,ndistb
1006         exden=sqrt(2.)*log(sig(k))
1007         DO nb=1,nbins
1008            d2min=itv(nb)
1009            d2max=itv(nb+1)
1010            numin=log(d2min/diam(k))
1011            numax=log(d2max/diam(k))
1012            massfrac(k,nb)=0.5*(erf(numax/exden)-erf(numin/exden))
1013            !print *,k,nb,massfrac(k,nb)
1014         ENDDO
1015      ENDDO
1016
1017!$OMP MASTER
1018IF (is_mpi_root .AND. is_omp_root) THEN
1019      OPEN (unit=15001, file='massfrac')
1020      DO k=1,ndistb
1021        DO nb=1,nbins
1022            write(15001,*),k,nb,massfrac(k,nb)
1023        ENDDO
1024      ENDDO
1025ENDIF
1026!$OMP END MASTER
1027!$OMP BARRIER
1028
1029 !CALL abort_gcm('calcdustemission', 'OK1',1)
1030
1031  !------------
1032
1033
1034  END SUBROUTINE initdust
1035
1036!--------------------------------------------------------------------------------------
1037!======================================================================================
1038!**************************************************************************************
1039!======================================================================================
1040!--------------------------------------------------------------------------------------
1041
1042  SUBROUTINE calcdustemission(debutphy,zu10m,zv10m,wstar, &
1043                              ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, &
1044                              emisbin)
1045  ! emisions over 12 dust bin
1046  USE dimphy
1047  USE infotrac
1048
1049  IMPLICIT NONE
1050  ! Input
1051  LOGICAL, INTENT(IN)                   :: debutphy ! First physiqs run or not
1052  REAL,DIMENSION(klon),INTENT(IN)          :: zu10m   ! 10m zonal wind
1053  REAL,DIMENSION(klon),INTENT(IN)          :: zv10m   ! meridional 10m wind
1054  REAL,DIMENSION(klon),INTENT(IN)          :: wstar   
1055  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
1056  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
1057 
1058  ! Local variables
1059!  INTEGER, PARAMETER :: flag_wstar=150
1060!  INTEGER, PARAMETER :: flag_wstar=120
1061!  INTEGER, PARAMETER :: flag_wstar=125
1062  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarWAKE
1063  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
1064  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE :: fluxdust ! horizonal emission fluxes in UNITS for the nmod soil aerosol modes
1065  REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10ms   ! 10m wind distribution in m/s
1066  REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10cm   ! 10m wind distribution in cm/s
1067  REAL,DIMENSION(klon)                  :: zwstar   
1068  REAL,DIMENSION(nwb)                :: probu
1069!  REAL, DIMENSION(nmode) :: fluxN,ftN,adN,fdpN,pN,eN ! in the original code N=1,2,3
1070  REAL :: flux1,flux2,flux3,ft1,ft2,ft3
1071
1072
1073  REAL, PARAMETER        :: umin=21.
1074  REAL, PARAMETER        :: woff=0.5  ! min value of 10m wind speed accepted for emissions
1075  REAL :: pdfcum,U10mMOD,pdfu,weilambda
1076  REAL :: z0salt,ceff,cerod,cpcent
1077  REAL :: cdnms,ustarns,modwm,utmin
1078!JE20150202   REAL :: cdnms,ustarns,modwm
1079  REAL :: fdp1,fdp2,ad1,ad2,ad3,flux_diam
1080  REAL :: dfec1,dfec2,dfec3,t1,t2,t3,p1,p2,p3,dec,ec
1081  ! auxiliar counters
1082  INTEGER                               :: kwb
1083  INTEGER                               :: i,j,k,l,n
1084  INTEGER  :: kfin,ideb,ifin,kfin2,istep
1085  REAL :: auxreal
1086
1087  ! Output
1088  !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbin ! vertical emission fluxes in UNITS for the 12 bins
1089  REAL,DIMENSION(klon,nbins)  :: emisbin ! vertical emission fluxes in UNITS for the 12 bins
1090!$OMP THREADPRIVATE(fluxdust)
1091!$OMP THREADPRIVATE(wind10ms)
1092!$OMP THREADPRIVATE(wind10cm)
1093
1094  !----------------------------------------------------
1095  ! initialization
1096  !----------------------------------------------------
1097  IF (debutphy) THEN
1098!   ALLOCATE( emisbin(klon,nbins) )
1099   ALLOCATE( fluxdust(klon,nmode) )
1100   ALLOCATE( wind10ms(nwb) )
1101   ALLOCATE( wind10cm(nwb) )
1102  ENDIF !debutphy
1103
1104 
1105  DO i=1,klon
1106   DO j=1,nbins
1107    emisbin(i,j)  = 0.0
1108   ENDDO !j,nbind
1109   DO j=1,nmode
1110    fluxdust(i,j) = 0.0
1111   ENDDO !j,nmode
1112  zwstar(i) = 0.0
1113  ENDDO !i,klon
1114  !----------------------------------------------------
1115  ! calculations
1116  !----------------------------------------------------
1117
1118  ! including BL processes..
1119!JE201041120  zwstar=sqrt(2.*(ale_bl+0.01*(flag_wstar-100)*ale_wake))
1120!JE20141124  zwstar=sqrt(2.*(flag_wstarBL*ale_bl+0.01*(flag_wstar-100)*ale_wake))
1121!  print
1122!*,'zwstar=sqrt(2.*(',flag_wstarBL,'ale_bl+0.01*(',flag_wstar,'-100)*ale_wake))'
1123  !
1124    DO i=1,klon  ! main loop
1125     zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)+param_wstarWAKE(i)*ale_wake(i)))
1126     U10mMOD=MAX(woff,sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i)))
1127     pdfcum=0.
1128     ! Wind weibull distribution:
1129
1130           DO kwb=1,nwb
1131                flux1=0.
1132                flux2=0.
1133                flux3=0.
1134!JE20141205<< differences in weibull distribution: expectance of the distribution is
1135!0.89*U10mMODD instead of U10mMOD
1136! now: lambda parameter of weibull ditribution is estimated as
1137! lambda=U10mMOD/gamma(1+1/kref)
1138! gamma function estimated with stirling formula
1139                auxreal=1.+1./kref
1140                weilambda = U10mMOD/exp(auxreal*log(auxreal)-auxreal &
1141                         - 0.5*log(auxreal/(2.*pi))+1./(12.*auxreal) &
1142                         -1./(360.*(auxreal**3.))+1./(1260.*(auxreal**5.)))
1143                IF(nwb.gt.1)THEN
1144                   wind10ms(kwb)=kwb*2.*U10mMOD/nwb
1145!original
1146!                   pdfu=(kref/U10mMOD)*(wind10ms(kwb)/U10mMOD)**(kref-1) &
1147!                      *exp(-(wind10ms(kwb)/U10mMOD)**kref)
1148                   pdfu=(kref/weilambda)*(wind10ms(kwb)/weilambda)**(kref-1) &
1149                      *exp(-(wind10ms(kwb)/weilambda)**kref)
1150!                   !print *,'JEdbg  U10mMOD weilambda  ',U10mMOD,weilambda
1151!JE20141205>>
1152
1153                   probu(kwb)=pdfu*2.*U10mMOD/nwb
1154                   pdfcum=pdfcum+probu(kwb)
1155                      IF(probu(kwb).le.1.e-2)GOTO 70
1156                ELSE
1157                   wind10ms(kwb)=U10mMOD
1158                   probu(kwb)=1.
1159                ENDIF
1160             wind10cm(kwb)=wind10ms(kwb)*100.
1161             DO n=1,ntyp
1162                   ft1=0.
1163                   ft2=0.
1164                   ft3=0.
1165!JE20150129<<<<
1166
1167             IF(.FALSE.) THEN
1168!                  nat=int(sol(i,n))
1169!                    print *,i,n
1170                    IF(sol(i,n).gt.1..and.sol(i,n).lt.15.) nat=int(sol(i,n))
1171!JE20140526<<
1172!                    print *,'JE: WARNING: nat=0 forced to nat=99!! and doing nothing'
1173                   IF(sol(i,n).lt.0.5) THEN
1174                      nat=99
1175                      GOTO 80
1176                    ENDIF
1177!JE20140526>>
1178 
1179
1180                 !IF(n.eq.1.and.nat.eq.99)GOTO 80
1181             !      if(n.eq.1) print*,'nat1=',nat,'sol1=',sol(i,n)
1182                   IF(n.eq.1.and.nat.eq.99)GOTO 80
1183
1184             ENDIF
1185             IF(.TRUE.) THEN
1186                nat=int(sol(i,n))
1187                if(n == 1 .and. nat >= 14 .or. nat < 1 .or. nat > 19) GOTO 80
1188             ENDIF
1189!JE20150129>>>>
1190
1191                      z0salt=z02(i,n)
1192                      ceff=feff(i,n)
1193                      cerod=A(i,n)
1194                      cpcent=P(i,n)
1195                      ustarsalt=0.
1196                   IF(ceff.le.0..or.z0salt.eq.0.)GOTO 80
1197                   IF(cerod.eq.0.or.cpcent.eq.0.)GOTO 80
1198! in cm: utmin, umin, z10m, z0salt, ustarns
1199! in meters: modwm
1200! no dimension for: cdnms
1201! Cas ou wsta=0.
1202                      cdnms=vkarm/(log(z10m/z0salt))
1203                      modwm=sqrt((wind10ms(kwb)**2)+(1.2*zwstar(i))**2)
1204                      ustarns=cdnms*modwm*100.
1205!JE20150202 <<
1206! Do not have too much sense.. and is not anymore in the chimere14b version.
1207!
1208!                      utmin=umin/(cdnms*ceff)
1209!                   IF(wind10cm(kwb).ge.utmin)THEN
1210!                      ustarsalt=ustarns+  &
1211!                    (0.3*(wind10cm(kwb)/100.-utmin/100.)**2.)
1212!                   ELSE
1213!                      ustarsalt=ustarns
1214!                   ENDIF
1215! ustarsalt should be :
1216                    ustarsalt=ustarns
1217!JE20150202 >>
1218
1219
1220                   IF(ustarsalt.lt.umin/ceff)GOTO 80
1221!                      print*,'ustarsalt = ',ustarsalt
1222!----------------------------------------
1223                    CALL def_copyncl(kfin)
1224
1225!                  CALL def_ag01(kfin,ft1,ft2,ft3)
1226      do ni=1,kfin
1227         fdp1=1.-(uth2(ni)/(ceff*ustarsalt))
1228         if (fdp1.le.0..or.srel2(nat,ni).eq.0.) then
1229            ad1=0.
1230            ad2=0.
1231            ad3=0.
1232         else
1233       
1234            fdp2=(1.+(uth2(ni)/(ceff*ustarsalt)))**2.
1235            flux_diam=cd*srel2(nat,ni)*fdp1*fdp2*(ustarsalt**3)
1236            dec=flux_diam*beta
1237            ec=(pi/3.)*1.e+2*rop* (sizeclass2(ni)**3) *(ustarsalt**2)
1238            dfec1=(ec-e1)
1239            dfec2=(ec-e2)
1240            dfec3=(ec-e3)
1241            t1=0.
1242            t2=0.
1243            t3=0.
1244            if(ec.ge.e1)t1=1.
1245            if(ec.ge.e2)t2=1.
1246            if(ec.ge.e3)t3=1.
1247            if(dfec3.ne.0.)then
1248               p1=t1*dfec1/dfec3
1249               p2=t2*(1.-p1)*dfec2/dfec3
1250            else
1251               p1=0.
1252               p2=0.
1253            endif
1254            p3=t3*(1.-p1-p2)
1255            ad1=(pi/6.)*dec*rop*p1*((diam(1)*1.e-4)**3.)/e1
1256            ad2=(pi/6.)*dec*rop*p2*((diam(2)*1.e-4)**3.)/e2
1257            ad3=(pi/6.)*dec*rop*p3*((diam(3)*1.e-4)**3.)/e3
1258
1259         endif
1260         ft1=ft1+ad1
1261         ft2=ft2+ad2
1262         ft3=ft3+ad3
1263      enddo ! of loop over nc
1264
1265!....
1266
1267        flux1=flux1+ft1*cpcent*cerod
1268        flux2=flux2+ft2*cpcent*cerod
1269        flux3=flux3+ft3*cpcent*cerod
1270!       print *,'JEflux :',kwb,n,flux1,flux2,flux3
127180 CONTINUE
1272             ENDDO !n=1,ntyp
127370 CONTINUE
1274        fluxdust(i,1)=fluxdust(i,1)+flux1*probu(kwb)
1275        fluxdust(i,2)=fluxdust(i,2)+flux2*probu(kwb)
1276        fluxdust(i,3)=fluxdust(i,3)+flux3*probu(kwb)
1277   ENDDO !kwb=1,nwb
1278      m1dflux(i)=10.*fluxdust(i,1)
1279      m2dflux(i)=10.*fluxdust(i,2)          ! tous en Kg/m2/s
1280      m3dflux(i)=10.*fluxdust(i,3)
1281
1282
1283
1284    ENDDO  ! i, klon
1285
1286
1287
1288   
1289  !from horizontal to  vertical fluxes for each bin
1290  DO k=1,klon
1291   DO i=1,ndistb
1292    DO j=1,nbins
1293!JE20150202 <<
1294    emisbin(k,j) = emisbin(k,j)+10*fluxdust(k,i)*massfrac(i,j) 
1295!JE20150202 >>
1296    ENDDO !j, nbind
1297   ENDDO  !i, nmode
1298  ENDDO !k,klon
1299
1300
1301!print *,' JE fluxdust in calcdust'
1302! DO k=1,klon
1303!   DO i=1,ndistb
1304!!print *,k,i,fluxdust(k,i)
1305!enddo
1306!enddo
1307!print *,' JE emisbin in calcdust'
1308!do k=1,klon
1309!do j=1,nbins
1310!!print *,k,j,emisbin(k,j)
1311!enddo
1312!enddo
1313!  CALL abort_gcm('calcdustemission', 'OK1',1)
1314
1315  END SUBROUTINE calcdustemission
1316!--------------------------------------------------------------------------------------
1317!======================================================================================
1318!**************************************************************************************
1319!======================================================================================
1320!--------------------------------------------------------------------------------------
1321
1322
1323
1324SUBROUTINE def_copyncl(kfin)
1325      implicit none
1326
1327 integer i,n,kfin,ideb,ifin,istep,kfin2
1328    real dsmin,dsmax
1329
1330! estimation of the reduced soil size distribution
1331         dsmin=var3a*(ustarsalt**(-2./3.))
1332         dsmax=var3b*(ustarsalt**(-2./3.))
1333!      print*,'ustarsalt = ',ustarsalt,'dsmin=',dsmin,'dsmax=',dsmax
1334! dichotomy
1335         call def_dichotomy(sizeclass,nclass,1,ncl,dsmin,ideb)
1336   !      print*,'ideb = ',ideb
1337         call def_dichotomy(sizeclass,nclass,ideb,ncl,dsmax,ifin)
1338   !      print*,'ifin = ',ifin
1339! readaptation of large sizes particles
1340         kfin=0
1341         do i=ideb,ifin
1342            kfin=kfin+1
1343            sizeclass2(kfin)=sizeclass(i)
1344            uth2(kfin)=uth(i)
1345            srel2(nat,kfin)=srel(nat,i)
1346         enddo
1347  !          print*,'je suis la'
1348         kfin2=kfin
1349         istep=50
1350         do i=ifin,ncl,istep
1351            kfin=kfin+1
1352            sizeclass2(kfin)=sizeclass(i)
1353            uth2(kfin)=uth(i)
1354            srel2(nat,kfin)=srel(nat,i)*istep
1355         enddo
1356         if(kfin.ge.nclass)then
1357            print*,'$$$$ Tables dimension problem:',kfin,'>',nclass
1358         endif
1359!---------------       
1360
1361      return
1362END SUBROUTINE def_copyncl
1363
1364!--------------------------------------------------------------------------------------
1365!======================================================================================
1366!**************************************************************************************
1367!======================================================================================
1368!--------------------------------------------------------------------------------------
1369
1370subroutine def_dichotomy(siz,nclass,i1,i2,ds,iout)
1371!c---------------------------------------------------------------
1372!c 'size' is the table to scan
1373!c      of dimension 'nclass', and reduced size of interest [i1:i2]
1374!c 'ds' is the value to find in 'size'
1375!c 'iout' is the number in the table 'size' correspondig to 'ds'
1376!c---------------------------------------------------------------
1377
1378      implicit none
1379      integer i1,i2,nclass,iout,ismin,ismax,k2,ihalf,idiff
1380      real siz(nclass),ds
1381!c-----------------------------
1382      iout=0
1383      ismin=i1
1384      ismax=i2
1385      ihalf=int((ismax+ismin)/2.)
1386      do k2=1,1000000
1387          if(ds.gt.siz(ihalf))then
1388             ismin=ihalf
1389          else
1390             ismax=ihalf
1391          endif
1392          ihalf=int((ismax+ismin)/2.)
1393          idiff=ismax-ismin
1394          if(idiff.le.1)then
1395             iout=ismin
1396             goto 52
1397          endif
1398      enddo
1399 52   continue
1400      if(iout.eq.0)then
1401        print*,'$$$$ Tables dimension problem: ',iout
1402      endif
1403
1404      end subroutine def_dichotomy
1405
1406!--------------------------------------------------------------------------------------
1407!======================================================================================
1408!**************************************************************************************
1409!======================================================================================
1410!--------------------------------------------------------------------------------------
1411
1412
1413  SUBROUTINE calcvd(vdout)
1414  IMPLICIT NONE
1415  INTEGER                     :: nb
1416!JE already def in module  INTEGER,PARAMETER           :: nbinsHR=3000
1417  INTEGER,PARAMETER           :: idiffusi=1
1418  INTEGER,PARAMETER           :: idrydep=2
1419  REAL                        :: dmn1,dmn2,dmn3,ra,rb
1420  REAL                        :: rhod,muair,nuair,R,airmolw
1421  REAL                        :: bolz,rhoair,gravity,karman,zed,z0m
1422  REAL                        :: ustarbin,temp,pres,rac,lambda,ccexp
1423  REAL                        :: Cc,rexp,pi,St
1424  REAL,DIMENSION(nbinsHR)     :: setvel
1425  REAL,DIMENSION(nbinsHR)     :: diffmol1
1426  REAL,DIMENSION(nbinsHR)     :: diffmol2
1427  REAL,DIMENSION(nbinsHR)     :: schmidtnumb
1428  REAL,DIMENSION(nbinsHR)     :: diffmole
1429
1430  REAL,DIMENSION(nbinsHR),   INTENT(OUT)              :: vdout
1431! users physical parameters
1432       temp=273.15+25. ! in Kelvin degrees
1433       pres=1013. ! in hPa
1434! calculation of Ra
1435       zed=10.
1436       z0m=0.05
1437       karman=0.41
1438       ustarbin=30. ! cm/s
1439       ra=log(zed/z0m)/(karman*ustarbin)
1440!c physical constants (all must in units g / cm / s)
1441       rhod=2.65                 ! g.cm-3
1442       muair=1.72e-4             ! g.cm-1.s-1
1443       nuair=0.1461              ! air kinematic viscosity [cm2.s-1]
1444       R=8.314                   ! molar gas constant J.mol-1.K-1
1445       airmolw=28.8              ! molar weight of air
1446       bolz=1.38e-16                    ! g.cm2.s-2
1447       rhoair=1.225e-3                  ! g.cm-3
1448       gravity=981.                ! gravity [cm.s-2]
1449       pi=3.14159
1450       rac=sqrt(8.0*airmolw/(pi*R*temp))
1451       lambda=(2.*muair)/(pres*rac)
1452! c Cc and Vs
1453!       binsHRcm=binsize*1.e-4
1454       DO nb=1,nbinsHR
1455       ccexp=exp((-0.55*binsHRcm(nb))/lambda)
1456       Cc=1.+(2.*lambda)/binsHRcm(nb)*(1.257+0.4*ccexp)
1457       setvel(nb)=(1./18.)*((binsHRcm(nb))**2*rhod*gravity*Cc)/muair
1458!c Molecular diffusivity (s/cm2) and Schmidt number
1459!c Note that the molecular diffusivity uses diameter in micro-m
1460!c to give a result directly in s/cm2
1461       dmn1=2.38e-7/binsHR(nb)
1462       dmn2=0.163/binsHR(nb)
1463       dmn3=0.0548*exp(-6.66*binsHR(nb))/binsHR(nb)
1464       diffmol1(nb)=dmn1*(1.+dmn2+dmn3)
1465       diffmol2(nb)=bolz*temp*Cc/(3.*pi*muair*binsHRcm(nb))
1466       IF(idiffusi.EQ.1)diffmole(nb)=diffmol1(nb)
1467       IF(idiffusi.EQ.2)diffmole(nb)=diffmol2(nb)
1468       schmidtnumb(nb)=nuair/diffmole(nb)
1469       St=setvel(nb)*ustarbin*ustarbin/(gravity*nuair)
1470       rb=1./(ustarbin*((schmidtnumb(nb))**(-2./3.)+10.**(-3./St)))
1471!c wesely (primarily designed for gases)
1472       IF(idrydep.EQ.1)THEN
1473          vdout(nb)=1./(ra+rb+ra*rb*setvel(nb))+setvel(nb)
1474       END IF
1475!c venkatram and pleim (more adaptated to particles but numerically unstable)
1476       IF(idrydep.EQ.2)THEN
1477        rexp=exp(-(ra+rb)*setvel(nb))
1478        vdout(nb)=setvel(nb)/(1.-rexp)
1479       END IF
1480       END DO
1481!c-----------------
1482       RETURN
1483      END SUBROUTINE calcvd
1484
1485
1486
1487END MODULE dustemission_mod
1488
Note: See TracBrowser for help on using the repository browser.