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

Last change on this file since 2494 was 2494, checked in by cmathe, 4 years ago

Mars GCM:
co2_ice as scatterer in radiative transfert. Need co2clouds and

activeco2ice .eqv. True. Files involved:

  • aeropacity_mod.F
  • callradite_mod.F
  • physiq_mod.F
  • updatereffrad_mod.F
  • suaer.F90
  • determine co2_ice density from temperature. Used in riceco2 computation.

Files involved:

  • co2cloud.F90
  • improvedco2clouds_mod.F90
  • updaterad.F90
  • updatereffrad_mod.F
  • co2condens_mod4micro.F: variable initialization
  • initracer.F: add nuiceco2_ref = 0.2
  • phyredem.F: remove co2_ice from qsurf since co2_ice => co2ice
  • watercloud_mod.F: tiny typo

CM

File size: 32.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  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      nuiceco2_ref=0.2    !C.M. Effective variance nueff of the
500                          ! co2-ice size distribution
501      if (doubleq) then
502c       "doubleq" technique
503c       -------------------
504c      (transport of mass and number mixing ratio)
505c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
506
507        if( (nq.lt.2).or.(water.and.(nq.lt.4))
508     *       .or.(hdo.and.(nq.lt.6) )) then
509          write(*,*)'initracer: nq is too low : nq=', nq
510          write(*,*)'water= ',water,' doubleq= ',doubleq   
511        end if
512
513        nueff_lift = 0.5
514        varian=sqrt(log(1.+nueff_lift))
515
516        rho_q(igcm_dust_mass)=rho_dust
517        rho_q(igcm_dust_number)=rho_dust
518
519c       Intermediate calcul for computing geometric mean radius r0
520c       as a function of mass and number mixing ratio Q and N
521c       (r0 = (r3n_q * Q/ N)^(1/3))
522        r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
523
524c       Intermediate calcul for computing effective radius reff
525c       from geometric mean radius r0
526c       (reff = ref_r0 * r0)
527        ref_r0 = exp(2.5*varian**2)
528       
529c       lifted dust :
530c       '''''''''''
531        reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
532        alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
533c       alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff
534
535!! default lifting settings
536!! -- GCM: alpha_lift not zero because large-scale lifting by default
537!! -- MESOSCALE: alpha_lift zero because no lifting at all in mesoscale by default
538#ifdef MESOSCALE
539        alpha_lift(igcm_dust_mass)=0.0
540#else
541        alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
542        IF (dustinjection.ge.1) THEN
543                reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)
544                alpha_lift(igcm_dust_mass)=(4/3.)*reff_lift*rho_dust
545     &                                          /2.4
546        ENDIF
547#endif
548
549        r0_lift = reff_lift/ref_r0
550        alpha_devil(igcm_dust_number)=r3n_q*
551     &                        alpha_devil(igcm_dust_mass)/r0_lift**3
552        alpha_lift(igcm_dust_number)=r3n_q*
553     &                        alpha_lift(igcm_dust_mass)/r0_lift**3
554
555        radius(igcm_dust_mass) = reff_lift
556        radius(igcm_dust_number) = reff_lift
557
558        write(*,*) "initracer: doubleq_param reff_lift:", reff_lift
559        write(*,*) "initracer: doubleq_param nueff_lift:", nueff_lift
560        write(*,*) "initracer: doubleq_param alpha_lift:",
561     &    alpha_lift(igcm_dust_mass)
562!c ----------------------------------------------------------------------
563!c rocket dust storm scheme
564!c lifting tracer stormdust using same distribution than
565!c normal dust
566        if (rdstorm) then
567          reff_storm=3.e-6 ! reff_lift !3.e-6
568          r0_storm=reff_storm/ref_r0
569          rho_q(igcm_stormdust_mass)=rho_dust
570          rho_q(igcm_stormdust_number)=rho_dust
571
572          alpha_devil(igcm_stormdust_mass)=9.e-9   ! dust devil lift mass coeff
573          alpha_lift(igcm_stormdust_mass)=4./3./2.4*reff_storm*rho_dust
574
575          write(*,*) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)
576 
577          alpha_devil(igcm_stormdust_number)=r3n_q*
578     &                      alpha_devil(igcm_stormdust_mass)/r0_storm**3
579          alpha_lift(igcm_stormdust_number)=r3n_q*
580     &                       alpha_lift(igcm_stormdust_mass)/r0_storm**3
581 
582          radius(igcm_stormdust_mass) = reff_storm
583          radius(igcm_stormdust_number) = reff_storm
584        end if !(rdstorm)
585!c ----------------------------------------------------------------------
586!c slope wind scheme
587!c you need a radius value for topdust to active its sedimentation
588!c we take the same value as for the normal dust
589        if (slpwind) then
590          rho_q(igcm_topdust_mass)=rho_dust
591          rho_q(igcm_topdust_number)=rho_dust
592          radius(igcm_topdust_mass) = 3.e-6
593          radius(igcm_topdust_number) = 3.e-6
594        end if !(slpwind)
595!c ----------------------------------------------------------------------
596     
597      else
598
599       ! initialize varian, which may be used (e.g. by surfacearea)
600       ! even with conrath dust
601       nueff_lift = 0.5
602       varian=sqrt(log(1.+nueff_lift))
603
604       if (dustbin.gt.1) then
605        print*,'initracer: STOP!',
606     $   ' properties of dust need to be set in initracer !!!'
607        call abort_physic("initracer","dustbin properties issue",1)
608
609       else if (dustbin.eq.1) then
610
611c       This will be used for 1 dust particle size:
612c       ------------------------------------------
613        radius(igcm_dustbin(1))=3.e-6
614        alpha_lift(igcm_dustbin(1))=0.0e-6
615        alpha_devil(igcm_dustbin(1))=7.65e-9
616        rho_q(igcm_dustbin(1))=rho_dust
617
618       endif
619      end if    ! (doubleq)
620
621
622c     Scavenging of dust particles by H2O clouds:
623c     ------------------------------------------
624c     Initialize the two tracers used for the CCNs
625      if (water.AND.doubleq.AND.scavenging) then
626        radius(igcm_ccn_mass) = radius(igcm_dust_mass)
627        alpha_lift(igcm_ccn_mass) = 1e-30
628        alpha_devil(igcm_ccn_mass) = 1e-30
629        rho_q(igcm_ccn_mass) = rho_dust
630
631        radius(igcm_ccn_number) = radius(igcm_ccn_mass)
632        alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
633        alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
634        rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
635      endif ! of if (water.AND.doubleq.AND.scavenging)
636
637c     Submicron dust mode:
638c     --------------------
639
640      if (submicron) then
641        radius(igcm_dust_submicron)=0.1e-6
642        rho_q(igcm_dust_submicron)=rho_dust
643        if (doubleq) then
644c         If doubleq is also active, we use the population ratio:
645          alpha_lift(igcm_dust_submicron) =
646     &      alpha_lift(igcm_dust_number)*popratio*
647     &      rho_q(igcm_dust_submicron)*4./3.*pi*
648     &      radius(igcm_dust_submicron)**3.
649          alpha_devil(igcm_dust_submicron)=1.e-30
650        else
651          alpha_lift(igcm_dust_submicron)=1e-6
652          alpha_devil(igcm_dust_submicron)=1.e-30
653        endif ! (doubleq)
654      end if  ! (submicron)
655
656c     Initialization for water vapor
657c     ------------------------------
658      if(water) then
659         radius(igcm_h2o_vap)=0.
660         alpha_lift(igcm_h2o_vap) =0.
661         alpha_devil(igcm_h2o_vap)=0.
662         if(water.and.(nq.ge.2)) then
663           radius(igcm_h2o_ice)=3.e-6
664           rho_q(igcm_h2o_ice)=rho_ice
665           alpha_lift(igcm_h2o_ice) =0.
666           alpha_devil(igcm_h2o_ice)=0.
667         elseif(water.and.(nq.lt.2)) then
668            write(*,*) 'nq is too low : nq=', nq
669            write(*,*) 'water= ',water
670         endif
671
672      end if  ! (water)
673
674c     Initialization for hdo vapor
675c     ------------------------------
676      if (hdo) then
677         radius(igcm_hdo_vap)=0.
678         alpha_lift(igcm_hdo_vap) =0.
679         alpha_devil(igcm_hdo_vap)=0.
680         if(water.and.(nq.ge.2)) then
681           radius(igcm_hdo_ice)=3.e-6
682           rho_q(igcm_hdo_ice)=rho_ice
683           alpha_lift(igcm_hdo_ice) =0.
684           alpha_devil(igcm_hdo_ice)=0.
685         elseif(hdo.and.(nq.lt.6)) then
686            write(*,*) 'nq is too low : nq=', nq
687            write(*,*) 'hdo= ',hdo
688         endif
689
690      end if  ! (hdo)
691
692     
693! Initialisation for CO2 clouds
694      if (co2clouds ) then
695        radius(igcm_ccnco2_mass) = radius(igcm_dust_mass)
696        alpha_lift(igcm_ccnco2_mass) = 1e-30
697        alpha_devil(igcm_ccnco2_mass) = 1e-30
698        rho_q(igcm_ccnco2_mass) = rho_dust
699        radius(igcm_ccnco2_number) = radius(igcm_ccnco2_mass)
700        alpha_lift(igcm_ccnco2_number) = alpha_lift(igcm_ccnco2_mass)
701        alpha_devil(igcm_ccnco2_number) = alpha_devil(igcm_ccnco2_mass)
702        rho_q(igcm_ccnco2_number) = rho_q(igcm_ccnco2_mass)
703     
704        radius(igcm_co2)=0.
705        alpha_lift(igcm_co2) =0.
706        alpha_devil(igcm_co2)=0.
707        radius(igcm_co2_ice)=1.e-8
708        rho_q(igcm_co2_ice)=rho_ice_co2
709        alpha_lift(igcm_co2_ice) =0.
710        alpha_devil(igcm_co2_ice)=0.
711
712      endif
713     
714c     Output for records:
715c     ~~~~~~~~~~~~~~~~~~
716      write(*,*)
717      Write(*,*) '******** initracer : dust transport parameters :'
718      write(*,*) 'alpha_lift = ', alpha_lift
719      write(*,*) 'alpha_devil = ', alpha_devil
720      write(*,*) 'radius  = ', radius
721      if(doubleq) then
722        write(*,*) 'reff_lift (um) =  ', reff_lift
723        write(*,*) 'size distribution variance  = ', varian
724        write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
725      end if
726
727!
728!     some extra (possibly redundant) sanity checks for tracers:
729!     ---------------------------------------------------------
730
731       if (doubleq) then
732       ! verify that we indeed have dust_mass and dust_number tracers
733         if (igcm_dust_mass.eq.0) then
734           write(*,*) "initracer: error !!"
735           write(*,*) "  cannot use doubleq option without ",
736     &                "a dust_mass tracer !"
737           call abort_physic("initracer","doubleq issue",1)
738         endif
739         if (igcm_dust_number.eq.0) then
740           write(*,*) "initracer: error !!"
741           write(*,*) "  cannot use doubleq option without ",
742     &                "a dust_number tracer !"
743           call abort_physic("initracer","doubleq issue",1)
744         endif
745       endif
746
747       if ((.not.doubleq).and.(dustbin.gt.0)) then
748       ! verify that we indeed have 'dustbin' dust tracers
749         count=0
750         do iq=1,dustbin
751           if (igcm_dustbin(iq).ne.0) then
752             count=count+1
753           endif
754         enddo
755         if (count.ne.dustbin) then
756           write(*,*) "initracer: error !!"
757           write(*,*) "  dustbin is set to ",dustbin,
758     &                " but we only have the following dust tracers:"
759           do iq=1,count
760             write(*,*)"   ",trim(noms(igcm_dustbin(iq)))
761           enddo
762           call abort_physic("initracer","dustbin issue",1)
763         endif
764       endif
765
766       if (water) then
767       ! verify that we indeed have h2o_vap and h2o_ice tracers
768         if (igcm_h2o_vap.eq.0) then
769           write(*,*) "initracer: error !!"
770           write(*,*) "  cannot use water option without ",
771     &                "an h2o_vap tracer !"
772           call abort_physic("initracer","water cycle issue",1)
773         endif
774         if (igcm_h2o_ice.eq.0) then
775           write(*,*) "initracer: error !!"
776           write(*,*) "  cannot use water option without ",
777     &                "an h2o_ice tracer !"
778           call abort_physic("initracer","water cycle issue",1)
779         endif
780       endif
781
782       if (hdo) then
783       ! verify that we indeed have hdo_vap and hdo_ice tracers
784         if (igcm_hdo_vap.eq.0) then
785           write(*,*) "initracer: error !!"
786           write(*,*) "  cannot use hdo option without ",
787     &                "an hdo_vap tracer !"
788           call abort_physic("initracer","hdo cycle issue",1)
789         endif
790         if (igcm_hdo_ice.eq.0) then
791           write(*,*) "initracer: error !!"
792           write(*,*) "  cannot use hdo option without ",
793     &                "an hdo_ice tracer !"
794           call abort_physic("initracer","hdo cycle issue",1)
795         endif
796       endif
797
798
799       if (co2clouds) then
800          !verify that we have co2_ice and co2 tracers
801          if (igcm_co2 .eq. 0) then
802             write(*,*) "initracer: error !!"
803             write(*,*) "  cannot use co2 clouds option without ",
804     &            "a co2 tracer !"
805             call abort_physic("initracer","co2 clouds issue",1)
806          endif
807          if (igcm_co2_ice .eq. 0) then
808             write(*,*) "initracer: error !!"
809             write(*,*) "  cannot use co2 clouds option without ",
810     &            "a co2_ice tracer !"
811             call abort_physic("initracer","co2 clouds issue",1)
812          endif
813       endif
814 
815       if (rdstorm) then
816       ! verify that we indeed have stormdust_mass and stormdust_number tracers
817         if (igcm_stormdust_mass.eq.0) then
818           write(*,*) "initracer: error !!"
819           write(*,*) "  cannot use rdstorm option without ",
820     &                "a stormdust_mass tracer !"
821           call abort_physic("initracer","rdstorm issue",1)
822         endif
823         if (igcm_stormdust_number.eq.0) then
824           write(*,*) "initracer: error !!"
825           write(*,*) "  cannot use rdstorm option without ",
826     &                "a stormdust_number tracer !"
827           call abort_physic("initracer","rdstorm issue",1)
828         endif
829       endif
830
831       if (slpwind) then
832       ! verify that we indeed have topdust_mass and topdust_number tracers
833         if (igcm_topdust_mass.eq.0) then
834           write(*,*) "initracer: error !!"
835           write(*,*) "  cannot use slpwind option without ",
836     &                "a topdust_mass tracer !"
837           call abort_physic("initracer","slpwind issue",1)
838         endif
839         if (igcm_topdust_number.eq.0) then
840           write(*,*) "initracer: error !!"
841           write(*,*) "  cannot use slpwind option without ",
842     &                "a topdust_number tracer !"
843           call abort_physic("initracer","slpwind issue",1)
844         endif
845       endif
846     
847       if (callnlte) then ! NLTE requirements
848         if (nltemodel.ge.1) then
849           ! check that co2, co, o and n2 tracers are available
850           if (igcm_co2.eq.0) then
851             write(*,*) "initracer: error !!"
852             write(*,*) "  with nltemodel>0, we need the co2 tracer!"
853             call abort_physic("initracer","missing co2 tracer",1)
854           endif
855           if (igcm_co.eq.0) then
856             write(*,*) "initracer: error !!"
857             write(*,*) "  with nltemodel>0, we need the co tracer!"
858             call abort_physic("initracer","missing co tracer",1)
859           endif
860           if (igcm_o.eq.0) then
861             write(*,*) "initracer: error !!"
862             write(*,*) "  with nltemodel>0, we need the o tracer!"
863             call abort_physic("initracer","missing o tracer",1)
864           endif
865           if (igcm_n2.eq.0) then
866             write(*,*) "initracer: error !!"
867             write(*,*) "  with nltemodel>0, we need the n2 tracer!"
868             call abort_physic("initracer","missing n2 tracer",1)
869           endif
870         endif
871       endif
872
873       if (scavenging) then
874       ! verify that we indeed have ccn_mass and ccn_number tracers
875         if (igcm_ccn_mass.eq.0 .and. igcm_ccnco2_mass.eq.0) then
876           write(*,*) "initracer: error !!"
877           write(*,*) "  cannot use scavenging option without ",
878     &                "a ccn_mass or ccnco2_mass tracer !"
879             call abort_physic("initracer","scavenging issue",1)
880         endif
881         if (igcm_ccn_number.eq.0 .and. igcm_ccnco2_number.eq.0 ) then
882           write(*,*) "initracer: error !!"
883           write(*,*) "  cannot use scavenging option without ",
884     &                "a ccn_number or ccnco2_number tracer !"
885             call abort_physic("initracer","scavenging issue",1)
886         endif
887       endif ! of if (scavenging)
888
889       if (photochem .or. callthermos) then
890       ! verify that we indeed have the chemistry tracers
891         if (igcm_co2.eq.0) then
892           write(*,*) "initracer: error !!"
893           write(*,*) "  cannot use chemistry option without ",
894     &                "a co2 tracer !"
895           call abort_physic("initracer","missing co2 tracer",1)
896         endif
897         if (igcm_co.eq.0) then
898           write(*,*) "initracer: error !!"
899           write(*,*) "  cannot use chemistry option without ",
900     &                "a co tracer !"
901           call abort_physic("initracer","missing co tracer",1)
902         endif
903         if (igcm_o.eq.0) then
904           write(*,*) "initracer: error !!"
905           write(*,*) "  cannot use chemistry option without ",
906     &                "a o tracer !"
907           call abort_physic("initracer","missing o tracer",1)
908         endif
909         if (igcm_o1d.eq.0) then
910           write(*,*) "initracer: error !!"
911           write(*,*) "  cannot use chemistry option without ",
912     &                "a o1d tracer !"
913           call abort_physic("initracer","missing o1d tracer",1)
914         endif
915         if (igcm_o2.eq.0) then
916           write(*,*) "initracer: error !!"
917           write(*,*) "  cannot use chemistry option without ",
918     &                "an o2 tracer !"
919           call abort_physic("initracer","missing o2 tracer",1)
920         endif
921         if (igcm_o3.eq.0) then
922           write(*,*) "initracer: error !!"
923           write(*,*) "  cannot use chemistry option without ",
924     &                "an o3 tracer !"
925           call abort_physic("initracer","missing o3 tracer",1)
926         endif
927         if (igcm_h.eq.0) then
928           write(*,*) "initracer: error !!"
929           write(*,*) "  cannot use chemistry option without ",
930     &                "a h tracer !"
931           call abort_physic("initracer","missing h tracer",1)
932         endif
933         if (igcm_h2.eq.0) then
934           write(*,*) "initracer: error !!"
935           write(*,*) "  cannot use chemistry option without ",
936     &                "a h2 tracer !"
937           call abort_physic("initracer","missing h2 tracer",1)
938         endif
939         if (igcm_oh.eq.0) then
940           write(*,*) "initracer: error !!"
941           write(*,*) "  cannot use chemistry option without ",
942     &                "an oh tracer !"
943           call abort_physic("initracer","missing oh tracer",1)
944         endif
945         if (igcm_ho2.eq.0) then
946           write(*,*) "initracer: error !!"
947           write(*,*) "  cannot use chemistry option without ",
948     &                "a ho2 tracer !"
949           call abort_physic("initracer","missing ho2 tracer",1)
950      endif
951         if (igcm_h2o2.eq.0) then
952           write(*,*) "initracer: error !!"
953           write(*,*) "  cannot use chemistry option without ",
954     &                "a h2o2 tracer !"
955           call abort_physic("initracer","missing h2o2 tracer",1)
956         endif
957         if (igcm_n2.eq.0) then
958           write(*,*) "initracer: error !!"
959           write(*,*) "  cannot use chemistry option without ",
960     &                "a n2 tracer !"
961           call abort_physic("initracer","missing n2 tracer",1)
962         endif
963         if (igcm_ar.eq.0) then
964           write(*,*) "initracer: error !!"
965           write(*,*) "  cannot use chemistry option without ",
966     &                "an ar tracer !"
967           call abort_physic("initracer","missing ar tracer",1)
968         endif
969       endif ! of if (photochem .or. callthermos)
970
971      end
Note: See TracBrowser for help on using the repository browser.