source: trunk/LMDZ.MARS/libf/phymars/initracer.F90 @ 4063

Last change on this file since 4063 was 4052, checked in by emillour, 3 weeks ago

Mars PCM:
Add posibilty for user to input lifted dust effective radius in callphys.def
as "reff_lift_dust" (default is the previously hard-coded value of 3.0e6 m)
EM
File size: 34.1 KB
Line 
1MODULE initracer_mod
2
3implicit none
4
5contains
6
7SUBROUTINE initracer(ngrid,nq)
8
9use tracer_mod
10use comcstfi_h,     only: pi
11use dust_param_mod, only: doubleq, submicron, dustbin
12use callkeys_mod,   only: water, microphys, scavenging, rdstorm, &
13                          topflows, photochem, callthermos, hdo, &
14                          callnlte, nltemodel, dustinjection,    &
15                          co2clouds, co2useh2o, meteo_flux
16use ioipsl_getin_p_mod, only : getin_p
17
18implicit none
19
20! == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
21!   subject:
22!   --------
23!   Initialization related to tracer
24!   (transported dust, water, chemical species, ice...)
25!
26!   Name of the tracer
27!
28!   Test of dimension :
29!   Initialize tracer related data in tracer_mod, using tracer names provided
30!   by the dynamics in "infotrac"
31!
32!
33!   author: F.Forget
34!   ------
35!    Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003)
36!            Ehouarn Millour (oct. 2008) identify tracers by their names
37!            JB Clement (juillet 2025), upgrading to F90 with module
38! == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
39
40
41integer, intent(in) :: ngrid ! number of atmospheric columns
42integer, intent(in) :: nq    ! number of tracers
43
44integer :: iq,ig,count
45real    :: r0_lift, reff_lift, nueff_lift
46real    :: r0_storm, reff_storm
47! Ratio of small over large dust particles (used when both
48!   doubleq and the submicron mode are active); In Montmessin
49!   et al. (2002), a value of 25 has been deduced;
50real, parameter                  :: popratio = 25.
51character(len=30)                :: txt ! to store some text
52character(len=30), dimension(nq) :: idtracers
53
54!-----------------------------------------------------------------------
55!  radius(nq)      ! aerosol particle radius (m)
56!  rho_q(nq)       ! tracer densities (kg.m-3)
57!  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
58!  alpha_devil(nq) ! lifting coeeficient by dust devil
59!  rho_dust          ! Mars dust density
60!  rho_ice           ! Water ice density
61!  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
62!  varian            ! Characteristic variance of log-normal distribution
63!-----------------------------------------------------------------------
64
65! Sanity check: check that we are running with at least 1 tracer:
66if (nq<1) then
67    write(*,*) "initracer error, nq=",nq," must be >=1!"
68    call abort_physic("initracer","nq<1",1)
69endif
70!------------------------------------------------------------
71!         NAME and molar mass of the tracer
72!------------------------------------------------------------
73
74! Identify tracers by their names: (and set corresponding values of mmol)
75! 0. initialize tracer indexes to zero:
76igcm_dustbin(1:nq)=0
77igcm_co2_ice=0
78igcm_ccnco2_mass=0
79igcm_ccnco2_number=0
80igcm_ccnco2_meteor_mass=0
81igcm_ccnco2_meteor_number=0
82igcm_ccnco2_h2o_mass_ice=0
83igcm_ccnco2_h2o_mass_ccn=0
84igcm_ccnco2_h2o_number=0
85igcm_dust_mass=0
86igcm_dust_number=0
87igcm_ccn_mass=0
88igcm_ccn_number=0
89igcm_dust_submicron=0
90igcm_h2o_vap=0
91igcm_h2o_ice=0
92igcm_hdo_vap=0
93igcm_hdo_ice=0
94igcm_stormdust_mass=0
95igcm_stormdust_number=0
96igcm_topdust_mass=0
97igcm_topdust_number=0
98igcm_co2=0
99igcm_co=0
100igcm_o=0
101igcm_o1d=0
102igcm_o2=0
103igcm_o3=0
104igcm_h=0
105igcm_d=0
106igcm_hd=0
107igcm_h2=0
108igcm_od=0
109igcm_do2=0
110igcm_hdo2=0
111igcm_oh=0
112igcm_ho2=0
113igcm_h2o2=0
114igcm_ch4=0
115igcm_n2=0
116igcm_ar=0
117igcm_ar_n2=0
118igcm_n=0
119igcm_no=0
120igcm_no2=0
121igcm_n2d=0
122igcm_he=0
123igcm_co2plus=0
124igcm_oplus=0
125igcm_o2plus=0
126igcm_coplus=0
127igcm_cplus=0
128igcm_nplus=0
129igcm_noplus=0
130igcm_n2plus=0
131igcm_hplus=0
132igcm_hco2plus=0
133igcm_hcoplus=0
134igcm_h2oplus=0
135igcm_h3oplus=0
136igcm_ohplus=0
137igcm_elec=0
138
139! 1. find dust tracers
140count=0
141if (dustbin>0) then
142    do iq=1,nq
143        txt=" "
144        write(txt,'(a4,i2.2)')'dust',count+1
145        if (noms(iq) == txt) then
146        count = count + 1
147        idtracers(count) = txt
148        igcm_dustbin(count)=iq
149        mmol(iq)=100.
150        endif
151    enddo !do iq=1,nq
152endif ! of if (dustbin>0)
153if (doubleq) then
154    do iq=1,nq
155        if (noms(iq) == "dust_mass") then
156            igcm_dust_mass=iq
157            count = count + 1
158            idtracers(count) = "dust_mass"
159        endif
160        if (noms(iq) == "dust_number") then
161            igcm_dust_number=iq
162            count = count + 1
163            idtracers(count) = "dust_number"
164        endif
165    enddo
166endif ! of if (doubleq)
167if (microphys) then
168    do iq=1,nq
169        if (noms(iq) == "ccn_mass") then
170            igcm_ccn_mass=iq
171            count = count + 1
172            idtracers(count) = "ccn_mass"
173        endif
174        if (noms(iq) == "ccn_number") then
175            igcm_ccn_number=iq
176            count = count + 1
177            idtracers(count) = "ccn_number"
178        endif
179    enddo
180endif ! of if (microphys)
181if (submicron) then
182    do iq=1,nq
183        if (noms(iq) == "dust_submicron") then
184            igcm_dust_submicron=iq
185            mmol(iq)=100.
186            count = count + 1
187            idtracers(count) = "dust_submicron"
188        endif
189    enddo
190endif ! of if (submicron)
191if (rdstorm) then
192    do iq=1,nq
193        if (noms(iq) == "stormdust_mass") then
194            igcm_stormdust_mass=iq
195            count = count + 1
196            idtracers(count) = "stormdust_mass"
197        endif
198        if (noms(iq) == "stormdust_number") then
199            igcm_stormdust_number=iq
200            count = count + 1
201            idtracers(count) = "stormdust_number"
202        endif
203    enddo
204endif ! of if (rdstorm)
205if (topflows) then
206    do iq=1,nq
207        if (noms(iq) == "topdust_mass") then
208            igcm_topdust_mass=iq
209            count = count + 1
210            idtracers(count) = "topdust_mass"
211        endif
212        if (noms(iq) == "topdust_number") then
213            igcm_topdust_number=iq
214            count = count + 1
215            idtracers(count) = "topdust_number"
216        endif
217    enddo
218endif ! of if (topflows)
219! 2. find chemistry and water tracers
220do iq=1,nq
221    if (noms(iq) == "co2") then
222        igcm_co2=iq
223        mmol(igcm_co2)=44.
224        count = count + 1
225        idtracers(count) = "co2"
226    endif
227    if (noms(iq) == "co") then
228        igcm_co=iq
229        mmol(igcm_co)=28.
230        count = count + 1
231        idtracers(count) = "co"
232    endif
233    if (noms(iq) == "o") then
234        igcm_o=iq
235        mmol(igcm_o)=16.
236        count = count + 1
237        idtracers(count) = "o"
238    endif
239    if (noms(iq) == "o1d") then
240        igcm_o1d=iq
241        mmol(igcm_o1d)=16.
242        count = count + 1
243        idtracers(count) = "o1d"
244    endif
245    if (noms(iq) == "o2") then
246        igcm_o2=iq
247        mmol(igcm_o2)=32.
248        count = count + 1
249        idtracers(count) = "o2"
250    endif
251    if (noms(iq) == "o3") then
252        igcm_o3=iq
253        mmol(igcm_o3)=48.
254        count = count + 1
255        idtracers(count) = "o3"
256    endif
257    if (noms(iq) == "h") then
258        igcm_h=iq
259        mmol(igcm_h)=1.
260        count = count + 1
261        idtracers(count) = "h"
262    endif
263    if (noms(iq) == "h2") then
264        igcm_h2=iq
265        mmol(igcm_h2)=2.
266        count = count + 1
267        idtracers(count) = "h2"
268    endif
269    if (noms(iq) == "oh") then
270        igcm_oh=iq
271        mmol(igcm_oh)=17.
272        count = count + 1
273        idtracers(count) = "oh"
274    endif
275    if (noms(iq) == "ho2") then
276        igcm_ho2=iq
277        mmol(igcm_ho2)=33.
278        count = count + 1
279        idtracers(count) = "ho2"
280    endif
281    if (noms(iq) == "h2o2") then
282        igcm_h2o2=iq
283        mmol(igcm_h2o2)=34.
284        count = count + 1
285        idtracers(count) = "h2o2"
286    endif
287    if (noms(iq) == "n2") then
288        igcm_n2=iq
289        mmol(igcm_n2)=28.
290        count = count + 1
291        idtracers(count) = "n2"
292    endif
293    if (noms(iq) == "ch4") then
294        igcm_ch4=iq
295        mmol(igcm_ch4)=16.
296        count = count + 1
297        idtracers(count) = "ch4"
298    endif
299    if (noms(iq) == "ar") then
300        igcm_ar=iq
301        mmol(igcm_ar)=40.
302        count = count + 1
303        idtracers(count) = "ar"
304    endif
305    if (noms(iq) == "n") then
306        igcm_n=iq
307        mmol(igcm_n)=14.
308        count = count + 1
309        idtracers(count) = "n"
310    endif
311    if (noms(iq) == "no") then
312        igcm_no=iq
313        mmol(igcm_no)=30.
314        count = count + 1
315        idtracers(count) = "no"
316    endif
317    if (noms(iq) == "no2") then
318        igcm_no2=iq
319        mmol(igcm_no2)=46.
320        count = count + 1
321        idtracers(count) = "no2"
322    endif
323    if (noms(iq) == "n2d") then
324        igcm_n2d=iq
325        mmol(igcm_n2d)=28.
326        count = count + 1
327        idtracers(count) = "n2d"
328    endif
329    if (noms(iq) == "he") then
330        igcm_he=iq
331        mmol(igcm_he)=4.
332        count = count + 1
333        idtracers(count) = "he"
334    endif
335    if (noms(iq) == "co2plus") then
336        igcm_co2plus=iq
337        mmol(igcm_co2plus)=44.
338        count = count + 1
339        idtracers(count) = "co2plus"
340    endif
341    if (noms(iq) == "oplus") then
342        igcm_oplus=iq
343        mmol(igcm_oplus)=16.
344        count = count + 1
345        idtracers(count) = "oplus"
346    endif
347    if (noms(iq) == "o2plus") then
348        igcm_o2plus=iq
349        mmol(igcm_o2plus)=32.
350        count = count + 1
351        idtracers(count) = "o2plus"
352    endif
353    if (noms(iq) == "coplus") then
354        igcm_coplus=iq
355        mmol(igcm_coplus)=28.
356        count = count + 1
357        idtracers(count) = "coplus"
358    endif
359    if (noms(iq) == "cplus") then
360        igcm_cplus=iq
361        mmol(igcm_cplus)=12.
362        count = count + 1
363        idtracers(count) = "cplus"
364    endif
365    if (noms(iq) == "nplus") then
366        igcm_nplus=iq
367        mmol(igcm_nplus)=14.
368        count = count + 1
369        idtracers(count) = "nplus"
370    endif
371    if (noms(iq) == "noplus") then
372        igcm_noplus=iq
373        mmol(igcm_noplus)=30.
374        count = count + 1
375        idtracers(count) = "noplus"
376    endif
377    if (noms(iq) == "n2plus") then
378        igcm_n2plus=iq
379        mmol(igcm_n2plus)=28.
380        count = count + 1
381        idtracers(count) = "n2plus"
382    endif
383    if (noms(iq) == "hplus") then
384        igcm_hplus=iq
385        mmol(igcm_hplus)=1.
386        count = count + 1
387        idtracers(count) = "hplus"
388    endif
389    if (noms(iq) == "hco2plus") then
390        igcm_hco2plus=iq
391        mmol(igcm_hco2plus)=45.
392        count = count + 1
393        idtracers(count) = "hco2plus"
394    endif
395    if (noms(iq) == "hcoplus") then
396        igcm_hcoplus=iq
397        mmol(igcm_hcoplus)=29.
398        count = count + 1
399        idtracers(count) = "hcoplus"
400    endif
401    if (noms(iq) == "h2oplus") then
402        igcm_h2oplus=iq
403        mmol(igcm_h2oplus)=18.
404        count = count + 1
405        idtracers(count) = "h2oplus"
406    endif
407    if (noms(iq) == "h3oplus") then
408        igcm_h3oplus=iq
409        mmol(igcm_h3oplus)=19.
410        count = count + 1
411        idtracers(count) = "h3oplus"
412    endif
413    if (noms(iq) == "ohplus") then
414        igcm_ohplus=iq
415        mmol(igcm_ohplus)=17.
416        count = count + 1
417        idtracers(count) = "ohplus"
418    endif
419    if (noms(iq) == "elec") then
420        igcm_elec=iq
421        mmol(igcm_elec)=1./1822.89
422        count = count + 1
423        idtracers(count) = "elec"
424    endif
425    if (noms(iq) == "h2o_vap") then
426        igcm_h2o_vap=iq
427        mmol(igcm_h2o_vap)=18.
428        count = count + 1
429        idtracers(count) = "h2o_vap"
430    endif
431    if (noms(iq) == "hdo_vap") then
432        igcm_hdo_vap=iq
433        mmol(igcm_hdo_vap)=19.
434        count = count + 1
435        idtracers(count) = "hdo_vap"
436    endif
437    if (noms(iq) == "od") then
438        igcm_od=iq
439        mmol(igcm_od)=18.
440        count = count + 1
441        idtracers(count) = "od"
442    endif
443    if (noms(iq) == "d") then
444        igcm_d=iq
445        mmol(igcm_d)=2.
446        count = count + 1
447        idtracers(count) = "d"
448    endif
449    if (noms(iq) == "hd") then
450        igcm_hd=iq
451        mmol(igcm_hd)=3.
452        count = count + 1
453        idtracers(count) = "hd"
454    endif
455    if (noms(iq) == "do2") then
456        igcm_do2=iq
457        mmol(igcm_do2)=34.
458        count = count + 1
459        idtracers(count) = "do2"
460    endif
461    if (noms(iq) == "hdo2") then
462        igcm_hdo2=iq
463        mmol(igcm_hdo2)=35.
464        count = count + 1
465        idtracers(count) = "hdo2"
466    endif
467    if (noms(iq) == "co2_ice") then
468        igcm_co2_ice=iq
469        mmol(igcm_co2_ice)=44.
470        count = count + 1
471        idtracers(count) = "co2_ice"
472    endif
473    if (noms(iq) == "h2o_ice") then
474        igcm_h2o_ice=iq
475        mmol(igcm_h2o_ice)=18.
476        count = count + 1
477        idtracers(count) = "h2o_ice"
478    endif
479    if (noms(iq) == "hdo_ice") then
480        igcm_hdo_ice=iq
481        mmol(igcm_hdo_ice)=19.
482        count = count + 1
483        idtracers(count) = "hdo_ice"
484    endif
485    ! Other stuff: e.g. for simulations using co2 + neutral gaz
486    if (noms(iq) == "Ar_N2") then
487        igcm_ar_n2=iq
488        mmol(igcm_ar_n2)=30.
489        count = count + 1
490        idtracers(count) = "Ar_N2"
491    endif
492    if (co2clouds) then
493        if (noms(iq) == "ccnco2_mass") then
494            igcm_ccnco2_mass=iq
495            count = count + 1
496            idtracers(count) = "ccnco2_mass"
497        endif
498        if (noms(iq) == "ccnco2_number") then
499            igcm_ccnco2_number=iq
500            count = count + 1
501            idtracers(count) = "ccnco2_number"
502        endif
503        if (meteo_flux) then
504            if (noms(iq) == "ccnco2_meteor_mass") then
505                igcm_ccnco2_meteor_mass=iq
506                count = count + 1
507                idtracers(count) = "ccnco2_meteor_mass"
508            endif
509            if (noms(iq) == "ccnco2_meteor_number") then
510                igcm_ccnco2_meteor_number=iq
511                count = count + 1
512                idtracers(count) = "ccnco2_meteor_number"
513            endif
514        end if
515        if (co2useh2o) then
516            if (noms(iq) == "ccnco2_h2o_number") then
517                igcm_ccnco2_h2o_number=iq
518                count = count + 1
519                idtracers(count) = "ccnco2_h2o_number"
520            endif
521            if (noms(iq) == "ccnco2_h2o_mass_ice") then
522                igcm_ccnco2_h2o_mass_ice=iq
523                count = count + 1
524                idtracers(count) = "ccnco2_h2o_mass_ice"
525            endif
526            if (noms(iq) == "ccnco2_h2o_mass_ccn") then
527                igcm_ccnco2_h2o_mass_ccn=iq
528                count = count + 1
529                idtracers(count) = "ccnco2_h2o_mass_ccn"
530            endif
531        end if
532    endif
533enddo ! of do iq=1,nq
534
535! check that we identified all tracers:
536if (count /= nq) then
537    write(*,*) "initracer: found only",count,"tracers but",nq,"are expected according to traceur.def!"
538    write(*,*) "initracer: the only identified tracers are:"
539    do iq = 1,count
540        write(*,*) '      ',iq,trim(idtracers(iq))
541    enddo
542    call abort_physic("initracer","tracer mismatch",1)
543else
544    write(*,*) "initracer: found all expected tracers, namely:"
545    do iq = 1,nq
546        write(*,*) '      ',iq,' ',trim(noms(iq))
547    enddo
548endif
549
550! Additional test required for HDO
551! We need to compute some things for H2O before HDO
552if (hdo) then
553    if (igcm_h2o_vap>igcm_hdo_vap) then
554        call abort_physic("initracer","Tracer H2O must be initialized before HDO",1)
555    else if ((nqfils(igcm_h2o_ice)<1) .or. (nqfils(igcm_h2o_vap)<1)) then
556        write(*,*) "HDO must be transported as a son of H2O: complete traceur.def"
557        call abort_physic("initracer","adapt your tracer.def",1)
558    else if ((igcm_hdo_ice < nq-2) .or. (igcm_hdo_vap < nq-2)) then
559        write(*,*) "The isotopes (HDO) must be placed at the end of the list in traceur.def"
560        call abort_physic("initracer","adapt your tracer.def",1)
561    endif
562endif
563
564!------------------------------------------------------------
565!     Initialize tracers .... (in tracer_mod)
566!------------------------------------------------------------
567! start by setting everything to (default) zero
568rho_q(1:nq)=0     ! tracer density (kg.m-3)
569radius(1:nq)=0.   ! tracer particle radius (m)
570alpha_lift(1:nq) =0.  ! tracer saltation vertical flux/horiz flux ratio (m-1)
571alpha_devil(1:nq)=0.  ! tracer lifting coefficient by dust devils
572
573
574! some reference values
575rho_dust=2500.  ! Mars dust density (kg.m-3)
576rho_ice=920.    ! Water ice density (kg.m-3)
577rho_ice_co2=1650.
578
579if (doubleq) then
580    ! verify that we indeed have dust_mass and dust_number tracers
581    if (igcm_dust_mass == 0) then
582           write(*,*) "initracer: error!"
583            write(*,*) "  Cannot use doubleq option without a dust_mass tracer!"
584           call abort_physic("initracer","doubleq issue",1)
585    endif
586    if (igcm_dust_number == 0) then
587        write(*,*) "initracer: error!"
588        write(*,*) "  Cannot use doubleq option without a dust_number tracer!"
589        call abort_physic("initracer","doubleq issue",1)
590    endif
591
592    ! "doubleq" technique
593    ! -------------------
594    ! transport of mass and number mixing ratio)
595    ! iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
596
597    if( (nq<2).or.(water.and.(nq<4)) .or.(hdo.and.(nq<6) )) then
598        write(*,*)'initracer: nq is too low : nq=', nq
599        write(*,*)'water= ',water,' doubleq= ',doubleq
600    end if
601
602    nueff_lift = 0.5
603    varian=sqrt(log(1.+nueff_lift))
604
605    rho_q(igcm_dust_mass)=rho_dust
606    rho_q(igcm_dust_number)=rho_dust
607
608    ! Intermediate calcul for computing geometric mean radius r0
609    ! as a function of mass and number mixing ratio Q and N
610    ! (r0 = (r3n_q * Q/ N)^(1/3))
611    r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
612
613    ! Intermediate calcul for computing effective radius reff
614    ! from geometric mean radius r0
615    ! (reff = ref_r0 * r0)
616    ref_r0 = exp(2.5*varian**2)
617
618    ! lifted dust :
619    ! '''''''''''
620    reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
621    call getin_p("reff_lift_dust",reff_lift)
622   
623    alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
624    ! alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff
625
626    ! default lifting settings
627    ! -- GCM: alpha_lift not zero because large-scale lifting by default
628    ! -- MESOSCALE: alpha_lift zero because no lifting at all in mesoscale by default
629#ifdef MESOSCALE
630    alpha_lift(igcm_dust_mass)=0.0
631#else
632    alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
633    if (dustinjection >= 1) then
634        !reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)! already set
635        alpha_lift(igcm_dust_mass)=(4/3.)*reff_lift*rho_dust/2.4
636    endif
637#endif
638
639    r0_lift = reff_lift/ref_r0
640    alpha_devil(igcm_dust_number)=r3n_q*alpha_devil(igcm_dust_mass)/r0_lift**3
641    alpha_lift(igcm_dust_number)=r3n_q*alpha_lift(igcm_dust_mass)/r0_lift**3
642
643    radius(igcm_dust_mass) = reff_lift
644    radius(igcm_dust_number) = reff_lift
645
646    write(*,*) "initracer: doubleq_param reff_lift:", reff_lift
647    write(*,*) "initracer: doubleq_param nueff_lift:", nueff_lift
648    write(*,*) "initracer: doubleq_param alpha_lift:",alpha_lift(igcm_dust_mass)
649    ! ----------------------------------------------------------------------
650    ! rocket dust storm scheme
651    ! lifting tracer stormdust using same distribution than
652    ! normal dust
653    if (rdstorm) then
654        reff_storm=reff_lift !3.e-6
655        r0_storm=reff_storm/ref_r0
656        rho_q(igcm_stormdust_mass)=rho_dust
657        rho_q(igcm_stormdust_number)=rho_dust
658
659        alpha_devil(igcm_stormdust_mass)=9.e-9   ! dust devil lift mass coeff
660        alpha_lift(igcm_stormdust_mass)=4./3./2.4*reff_storm*rho_dust
661
662        write(*,*) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)
663
664        alpha_devil(igcm_stormdust_number)=r3n_q*alpha_devil(igcm_stormdust_mass)/r0_storm**3
665        alpha_lift(igcm_stormdust_number)=r3n_q*alpha_lift(igcm_stormdust_mass)/r0_storm**3
666
667        radius(igcm_stormdust_mass) = reff_storm
668        radius(igcm_stormdust_number) = reff_storm
669    end if !(rdstorm)
670    ! ----------------------------------------------------------------------
671    ! mountain top dust flows scheme
672    ! you need a radius value for topdust to active its sedimentation
673    ! we take the same value as for the normal dust
674    if (topflows) then
675        rho_q(igcm_topdust_mass)=rho_dust
676        rho_q(igcm_topdust_number)=rho_dust
677        radius(igcm_topdust_mass) = reff_lift !3.e-6
678        radius(igcm_topdust_number) = reff_lift !3.e-6
679    end if !(topflows)
680    ! ----------------------------------------------------------------------
681
682else
683
684    ! initialize varian, which may be used (e.g. by surfacearea)
685    ! even with conrath dust
686    nueff_lift = 0.5
687    varian=sqrt(log(1.+nueff_lift))
688
689    if (dustbin>1) then
690        print*,'initracer: STOP! properties of dust need to be set in initracer!!!'
691        call abort_physic("initracer","dustbin properties issue",1)
692
693    else if (dustbin == 1) then
694
695!       This will be used for 1 dust particle size:
696!       ------------------------------------------
697        radius(igcm_dustbin(1))=3.e-6
698        alpha_lift(igcm_dustbin(1))=0.0e-6
699        alpha_devil(igcm_dustbin(1))=7.65e-9
700        rho_q(igcm_dustbin(1))=rho_dust
701
702    endif
703end if    ! (doubleq)
704
705
706! Scavenging of dust particles by H2O clouds:
707! ------------------------------------------
708! Initialize the two tracers used for the CCNs
709if (water.and.doubleq.and.scavenging) then
710    radius(igcm_ccn_mass) = radius(igcm_dust_mass)
711    alpha_lift(igcm_ccn_mass) = 1e-30
712    alpha_devil(igcm_ccn_mass) = 1e-30
713    rho_q(igcm_ccn_mass) = rho_dust
714
715    radius(igcm_ccn_number) = radius(igcm_ccn_mass)
716    alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
717    alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
718    rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
719endif ! of if (water.AND.doubleq.AND.scavenging)
720
721! Submicron dust mode:
722! --------------------
723if (submicron) then
724    radius(igcm_dust_submicron)=0.1e-6
725    rho_q(igcm_dust_submicron)=rho_dust
726    if (doubleq) then
727!        If doubleq is also active, we use the population ratio:
728        alpha_lift(igcm_dust_submicron) = alpha_lift(igcm_dust_number)*popratio*rho_q(igcm_dust_submicron)*4./3.*pi*radius(igcm_dust_submicron)**3.
729        alpha_devil(igcm_dust_submicron)=1.e-30
730    else
731        alpha_lift(igcm_dust_submicron)=1e-6
732        alpha_devil(igcm_dust_submicron)=1.e-30
733    endif ! (doubleq)
734end if  ! (submicron)
735
736! Initialization for water vapor
737! ------------------------------
738if (water) then
739    radius(igcm_h2o_vap)=0.
740    alpha_lift(igcm_h2o_vap) =0.
741    alpha_devil(igcm_h2o_vap)=0.
742    if (water.and.(nq>=2)) then
743        radius(igcm_h2o_ice)=3.e-6
744        rho_q(igcm_h2o_ice)=rho_ice
745        alpha_lift(igcm_h2o_ice) =0.
746        alpha_devil(igcm_h2o_ice)=0.
747    else if(water.and.(nq<2)) then
748        write(*,*) 'nq is too low : nq=', nq
749        write(*,*) 'water= ',water
750    endif
751end if  ! (water)
752
753! Initialization for hdo vapor
754! ------------------------------
755if (hdo) then
756    radius(igcm_hdo_vap)=0.
757    alpha_lift(igcm_hdo_vap) =0.
758    alpha_devil(igcm_hdo_vap)=0.
759    if (water.and.(nq>=2)) then
760        radius(igcm_hdo_ice)=3.e-6
761        rho_q(igcm_hdo_ice)=rho_ice
762        alpha_lift(igcm_hdo_ice) =0.
763        alpha_devil(igcm_hdo_ice)=0.
764    else if (hdo.and.(nq<6)) then
765        write(*,*) 'nq is too low : nq=', nq
766        write(*,*) 'hdo= ',hdo
767    endif
768end if  ! (hdo)
769
770! Output for records:
771! ~~~~~~~~~~~~~~~~~~
772write(*,*)
773Write(*,*) '******** initracer : dust transport parameters :'
774write(*,*) 'alpha_lift = ', alpha_lift
775write(*,*) 'alpha_devil = ', alpha_devil
776write(*,*) 'radius  = ', radius
777if(doubleq) then
778    write(*,*) 'reff_lift (um) =  ', reff_lift
779    write(*,*) 'size distribution variance  = ', varian
780    write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
781end if
782
783!
784! some extra (possibly redundant) sanity checks for tracers:
785! ---------------------------------------------------------
786if ((.not.doubleq).and.(dustbin>0)) then
787    ! verify that we indeed have 'dustbin' dust tracers
788    count=0
789    do iq=1,dustbin
790        if (igcm_dustbin(iq)/=0) then
791            count = count + 1
792        endif
793    enddo
794    if (count/=dustbin) then
795        write(*,*) "initracer: error!"
796        write(*,*) "  dustbin is set to ",dustbin," but we only have the following dust tracers:"
797        do iq=1,count
798            write(*,*)"   ",trim(noms(igcm_dustbin(iq)))
799        enddo
800        call abort_physic("initracer","dustbin issue",1)
801    endif
802endif
803
804if (water) then
805    ! verify that we indeed have h2o_vap and h2o_ice tracers
806    if (igcm_h2o_vap == 0) then
807        write(*,*) "initracer: error!"
808        write(*,*) "  Cannot use water option without an h2o_vap tracer!"
809        call abort_physic("initracer","water cycle issue",1)
810    endif
811    if (igcm_h2o_ice == 0) then
812        write(*,*) "initracer: error!"
813        write(*,*) "  Cannot use water option without an h2o_ice tracer!"
814        call abort_physic("initracer","water cycle issue",1)
815    endif
816endif
817
818if (hdo) then
819    ! verify that we indeed have hdo_vap and hdo_ice tracers
820    if (igcm_hdo_vap == 0) then
821        write(*,*) "initracer: error!"
822        write(*,*) "  Cannot use hdo option without an hdo_vap tracer!"
823        call abort_physic("initracer","hdo cycle issue",1)
824    endif
825    if (igcm_hdo_ice == 0) then
826        write(*,*) "initracer: error!"
827        write(*,*) "  Cannot use hdo option without an hdo_ice tracer!"
828        call abort_physic("initracer","hdo cycle issue",1)
829    endif
830endif
831
832if (microphys) then
833    ! verify that we indeed have ccn_number and ccn_mass tracers
834    if (igcm_ccn_number == 0) then
835        write(*,*) "initracer: error!"
836        write(*,*) "  Cannot use microphys option without a ccn_number tracer!"
837        call abort_physic("initracer","microphys issue",1)
838    endif
839    if (igcm_ccn_mass == 0) then
840        write(*,*) "initracer: error!"
841        write(*,*) "  Cannot use microphys option without a ccn_mass tracer!"
842        call abort_physic("initracer","microphys issue",1)
843    endif
844endif
845
846if (co2clouds) then
847    !verify that we have co2_ice and co2 tracers
848    if (igcm_co2 == 0) then
849        write(*,*) "initracer: error!"
850        write(*,*) "  Cannot use co2 clouds option without a co2 tracer!"
851        call abort_physic("initracer","co2 clouds issue",1)
852    end if
853    if (igcm_co2_ice == 0) then
854        write(*,*) "initracer: error!"
855        write(*,*) "  Cannot use co2 clouds option without a co2_ice tracer!"
856        call abort_physic("initracer","co2 clouds issue",1)
857    end if
858    if (igcm_ccnco2_number == 0) then
859        write(*,*) "initracer: error!"
860        write(*,*) "  Cannot use co2 clouds option without a ccnco2_number tracer!"
861        call abort_physic("initracer","co2 clouds issue",1)
862    end if
863    if (igcm_ccnco2_mass == 0) then
864        write(*,*) "initracer: error!"
865        write(*,*) "  Cannot use co2 clouds option without a ccnco2_mass tracer!"
866        call abort_physic("initracer","co2 clouds issue",1)
867    end if
868    if (co2useh2o) then
869        if (igcm_ccnco2_h2o_number == 0) then
870            write(*,*) "initracer: error!"
871            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_number tracer!"
872            call abort_physic("initracer","co2 clouds issue",1)
873        end if
874        if (igcm_ccnco2_h2o_mass_ice == 0) then
875            write(*,*) "initracer: error!"
876            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ice tracer!"
877            call abort_physic("initracer","co2 clouds issue",1)
878        end if
879        if (igcm_ccnco2_h2o_mass_ccn == 0) then
880            write(*,*) "initracer: error!"
881            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ccn tracer!"
882            call abort_physic("initracer","co2 clouds issue",1)
883        end if
884    end if
885    if (meteo_flux) then
886        if (igcm_ccnco2_meteor_number == 0) then
887            write(*,*) "initracer: error!"
888            write(*,*) "  Cannot use co2 clouds option without a ccnco2_meteor_number tracer!"
889            call abort_physic("initracer","co2 clouds issue",1)
890        end if
891        if (igcm_ccnco2_meteor_mass == 0) then
892            write(*,*) "initracer: error!"
893            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ice tracer!"
894            call abort_physic("initracer","co2 clouds issue",1)
895        end if
896        if (igcm_ccnco2_meteor_mass == 0) then
897            write(*,*) "initracer: error!"
898            write(*,*) "  Cannot use co2 clouds option without a ccnco2_meteor_mass tracer!"
899            call abort_physic("initracer","co2 clouds issue",1)
900        end if
901    end if
902endif
903
904if (rdstorm) then
905    ! verify that we indeed have stormdust_mass and stormdust_number tracers
906    if (igcm_stormdust_mass == 0) then
907        write(*,*) "initracer: error!"
908        write(*,*) "  Cannot use rdstorm option without a stormdust_mass tracer!"
909        call abort_physic("initracer","rdstorm issue",1)
910    endif
911    if (igcm_stormdust_number == 0) then
912        write(*,*) "initracer: error!"
913        write(*,*) "  Cannot use rdstorm option without a stormdust_number tracer!"
914        call abort_physic("initracer","rdstorm issue",1)
915    endif
916endif
917
918if (topflows) then
919    ! verify that we indeed have topdust_mass and topdust_number tracers
920    if (igcm_topdust_mass == 0) then
921        write(*,*) "initracer: error!"
922        write(*,*) "  Cannot use topflows option without a topdust_mass tracer!"
923        call abort_physic("initracer","topflows issue",1)
924    endif
925    if (igcm_topdust_number == 0) then
926        write(*,*) "initracer: error!"
927        write(*,*) "  Cannot use topflows option without a topdust_number tracer!"
928        call abort_physic("initracer","topflows issue",1)
929    endif
930endif
931
932if (callnlte) then ! NLTE requirements
933    if (nltemodel>=1) then
934        ! check that co2, co, o and n2 tracers are available
935        if (igcm_co2 == 0) then
936            write(*,*) "initracer: error!"
937            write(*,*) "  with nltemodel>0, we need the co2 tracer!"
938            call abort_physic("initracer","missing co2 tracer",1)
939        endif
940        if (igcm_co == 0) then
941            write(*,*) "initracer: error!"
942            write(*,*) "  with nltemodel>0, we need the co tracer!"
943            call abort_physic("initracer","missing co tracer",1)
944        endif
945        if (igcm_o == 0) then
946            write(*,*) "initracer: error!"
947            write(*,*) "  with nltemodel>0, we need the o tracer!"
948            call abort_physic("initracer","missing o tracer",1)
949        endif
950        if (igcm_n2 == 0) then
951            write(*,*) "initracer: error!"
952            write(*,*) "  with nltemodel>0, we need the n2 tracer!"
953            call abort_physic("initracer","missing n2 tracer",1)
954        endif
955    endif
956endif
957
958if (scavenging) then
959    ! verify that we indeed have ccn_mass and ccn_number tracers
960    if (igcm_ccn_mass == 0 .and. igcm_ccnco2_mass == 0) then
961        write(*,*) "initracer: error!"
962        write(*,*) "  Cannot use scavenging option without a ccn_mass or ccnco2_mass tracer!"
963        call abort_physic("initracer","scavenging issue",1)
964    endif
965    if (igcm_ccn_number == 0 .and. igcm_ccnco2_number == 0 ) then
966        write(*,*) "initracer: error!"
967        write(*,*) "  Cannot use scavenging option without a ccn_number or ccnco2_number tracer!"
968        call abort_physic("initracer","scavenging issue",1)
969    endif
970endif ! of if (scavenging)
971
972if (photochem .or. callthermos) then
973    ! verify that we indeed have the chemistry tracers
974    if (igcm_co2 == 0) then
975        write(*,*) "initracer: error!"
976        write(*,*) "  Cannot use chemistry option without a co2 tracer!"
977        call abort_physic("initracer","missing co2 tracer",1)
978    endif
979    if (igcm_co == 0) then
980        write(*,*) "initracer: error!"
981        write(*,*) "  Cannot use chemistry option without a co tracer!"
982        call abort_physic("initracer","missing co tracer",1)
983    endif
984    if (igcm_o == 0) then
985        write(*,*) "initracer: error!"
986        write(*,*) "  Cannot use chemistry option without a o tracer!"
987        call abort_physic("initracer","missing o tracer",1)
988    endif
989    if (igcm_o1d == 0) then
990        write(*,*) "initracer: error!"
991        write(*,*) "  Cannot use chemistry option without a o1d tracer!"
992        call abort_physic("initracer","missing o1d tracer",1)
993    endif
994    if (igcm_o2 == 0) then
995        write(*,*) "initracer: error!"
996        write(*,*) "  Cannot use chemistry option without an o2 tracer!"
997        call abort_physic("initracer","missing o2 tracer",1)
998    endif
999    if (igcm_o3 == 0) then
1000        write(*,*) "initracer: error!"
1001        write(*,*) "  Cannot use chemistry option without an o3 tracer!"
1002        call abort_physic("initracer","missing o3 tracer",1)
1003    endif
1004    if (igcm_h == 0) then
1005        write(*,*) "initracer: error!"
1006        write(*,*) "  Cannot use chemistry option without a h tracer!"
1007        call abort_physic("initracer","missing h tracer",1)
1008    endif
1009    if (igcm_h2 == 0) then
1010        write(*,*) "initracer: error!"
1011        write(*,*) "  Cannot use chemistry option without a h2 tracer!"
1012        call abort_physic("initracer","missing h2 tracer",1)
1013    endif
1014    if (igcm_oh == 0) then
1015        write(*,*) "initracer: error!"
1016        write(*,*) "  Cannot use chemistry option without an oh tracer!"
1017        call abort_physic("initracer","missing oh tracer",1)
1018    endif
1019    if (igcm_ho2 == 0) then
1020        write(*,*) "initracer: error!"
1021        write(*,*) "  Cannot use chemistry option without a ho2 tracer!"
1022        call abort_physic("initracer","missing ho2 tracer",1)
1023    endif
1024    if (igcm_h2o2 == 0) then
1025        write(*,*) "initracer: error!"
1026        write(*,*) "  Cannot use chemistry option without a h2o2 tracer!"
1027        call abort_physic("initracer","missing h2o2 tracer",1)
1028    endif
1029    if (igcm_n2 == 0) then
1030        write(*,*) "initracer: error!"
1031        write(*,*) "  Cannot use chemistry option without a n2 tracer!"
1032        call abort_physic("initracer","missing n2 tracer",1)
1033    endif
1034    if (igcm_ar == 0) then
1035        write(*,*) "initracer: error!"
1036        write(*,*) "  Cannot use chemistry option without an ar tracer!"
1037        call abort_physic("initracer","missing ar tracer",1)
1038    endif
1039endif ! of if (photochem .or. callthermos)
1040
1041! Initialisation for CO2 clouds
1042if (co2clouds) then
1043    radius(igcm_ccnco2_mass) = radius(igcm_dust_mass)
1044    alpha_lift(igcm_ccnco2_mass) = 1e-30
1045    alpha_devil(igcm_ccnco2_mass) = 1e-30
1046    rho_q(igcm_ccnco2_mass) = rho_dust
1047    radius(igcm_ccnco2_number) = radius(igcm_ccnco2_mass)
1048    alpha_lift(igcm_ccnco2_number) = alpha_lift(igcm_ccnco2_mass)
1049    alpha_devil(igcm_ccnco2_number) = alpha_devil(igcm_ccnco2_mass)
1050    rho_q(igcm_ccnco2_number) = rho_q(igcm_ccnco2_mass)
1051
1052    radius(igcm_co2)=0.
1053    alpha_lift(igcm_co2) =0.
1054    alpha_devil(igcm_co2)=0.
1055    radius(igcm_co2_ice)=1.e-8
1056    rho_q(igcm_co2_ice)=rho_ice_co2
1057    alpha_lift(igcm_co2_ice) =0.
1058    alpha_devil(igcm_co2_ice)=0.
1059endif
1060
1061END SUBROUTINE initracer
1062
1063END MODULE initracer_mod
1064
Note: See TracBrowser for help on using the repository browser.