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

Last change on this file since 5208 was 5182, checked in by abarral, 10 days ago

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