source: LMDZ6/trunk/libf/phylmd/Dust/dustemission_mod.F90 @ 3786

Last change on this file since 3786 was 3786, checked in by asima, 4 years ago

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