source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/dustemission_mod.F90 @ 5112

Last change on this file since 5112 was 5112, checked in by abarral, 2 months ago

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