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

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

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