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

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

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