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

Last change on this file since 2467 was 2463, checked in by emillour, 4 years ago

Mars GCM
Follow-up of previous commit. Also set the molar mass of HDO ice to 19.
EM

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