source: LMDZ6/branches/LMDZ_DECOUPLE/libf/phylmd/Dust/dustemission_mod.F90 @ 5006

Last change on this file since 5006 was 2630, checked in by fhourdin, 8 years ago

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