source: trunk/LMDZ.MARS/libf/phymars/initracer.F @ 3026

Last change on this file since 3026 was 3026, checked in by jbclement, 16 months ago

Mars PCM/PEM 1D:
Small fixes to be able to run the Mars PCM 1D without "water" + Improvements/addition of scripts in deftank/pem to run the PEM 1D model according to Laskar orbital parameters.
JBC

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