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

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

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