source: LMDZ6/trunk/libf/phylmd/Dust/dustemission_mod.f90 @ 5821

Last change on this file since 5821 was 5636, checked in by fhourdin, 7 months ago

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