source: trunk/LMDZ.GENERIC/libf/phystd/initracer.F @ 1834

Last change on this file since 1834 was 1801, checked in by bclmd, 7 years ago

Adding photochemistry to LMDZ Generic
File size: 16.4 KB
Line 
1      SUBROUTINE initracer(ngrid,nq,nametrac)
2
3      use surfdat_h
4      USE tracer_h
5      USE callkeys_mod, only: water
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 COMMON tracer in tracer.h, using tracer names provided
17c   by the argument nametrac
18c
19c   author: F.Forget
20c   ------
21c            Ehouarn Millour (oct. 2008) identify tracers by their names
22c=======================================================================
23
24      integer :: ngrid,nq
25
26!      real qsurf(ngrid,nq)       ! tracer on surface (e.g.  kg.m-2)
27!      real co2ice(ngrid)           ! co2 ice mass on surface (e.g.  kg.m-2)
28      character(len=20) :: txt ! to store some text
29      integer iq,ig,count
30      real r0_lift , reff_lift
31!      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
32
33      character*20 nametrac(nq)   ! name of the tracer from dynamics
34
35
36c-----------------------------------------------------------------------
37c  radius(nq)      ! aerosol particle radius (m)
38c  rho_q(nq)       ! tracer densities (kg.m-3)
39c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
40c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
41c  alpha_devil(nq) ! lifting coeeficient by dust devil
42c  rho_dust          ! Mars dust density
43c  rho_ice           ! Water ice density
44c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
45c  varian            ! Characteristic variance of log-normal distribution
46c-----------------------------------------------------------------------
47
48       nqtot=nq
49       !! we allocate once for all arrays in common in tracer_h.F90
50       !! (supposedly those are not used before call to initracer)
51       IF (.NOT.ALLOCATED(noms))         ALLOCATE(noms(nq))
52       IF (.NOT.ALLOCATED(mmol))         ALLOCATE(mmol(nq))
53       IF (.NOT.ALLOCATED(radius))       ALLOCATE(radius(nq))
54       IF (.NOT.ALLOCATED(rho_q))        ALLOCATE(rho_q(nq))
55       IF (.NOT.ALLOCATED(qext))         ALLOCATE(qext(nq))
56       IF (.NOT.ALLOCATED(alpha_lift))   ALLOCATE(alpha_lift(nq))
57       IF (.NOT.ALLOCATED(alpha_devil))  ALLOCATE(alpha_devil(nq))
58       IF (.NOT.ALLOCATED(qextrhor))     ALLOCATE(qextrhor(nq))
59       IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
60       !! initialization
61       alpha_lift(:)=0.
62       alpha_devil(:)=0.
63       
64       ! Added by JVO 2017 : these arrays are handled later
65       ! -> initialization is the least we can do, please !!!
66       radius(:)=0.
67       qext(:)=0.
68
69! Initialization: get tracer names from the dynamics and check if we are
70!                 using 'old' tracer convention ('q01',q02',...)
71!                 or new convention (full tracer names)
72      ! check if tracers have 'old' names
73
74! copy tracer names from dynamics
75        do iq=1,nq
76          noms(iq)=nametrac(iq)
77        enddo
78
79
80! Identify tracers by their names: (and set corresponding values of mmol)
81      ! 0. initialize tracer indexes to zero:
82      ! NB: igcm_* indexes are commons in 'tracer.h'
83      do iq=1,nq
84        igcm_dustbin(iq)=0
85      enddo
86      igcm_dust_mass=0
87      igcm_dust_number=0
88      igcm_h2o_vap=0
89      igcm_h2o_ice=0
90      igcm_co2=0
91      igcm_co=0
92      igcm_o=0
93      igcm_o1d=0
94      igcm_o2=0
95      igcm_o3=0
96      igcm_h=0
97      igcm_h2=0
98      igcm_oh=0
99      igcm_ho2=0
100      igcm_h2o2=0
101      igcm_n2=0
102      igcm_n=0
103      igcm_n2d=0
104      igcm_no=0
105      igcm_no2=0
106      igcm_ar=0
107      igcm_ar_n2=0
108      igcm_co2_ice=0
109
110      igcm_ch4=0
111      igcm_ch3=0
112      igcm_ch=0
113      igcm_3ch2=0
114      igcm_1ch2=0
115      igcm_cho=0
116      igcm_ch2o=0
117      igcm_ch3o=0
118      igcm_c=0
119      igcm_c2=0
120      igcm_c2h=0
121      igcm_c2h2=0
122      igcm_c2h3=0
123      igcm_c2h4=0
124      igcm_c2h6=0
125      igcm_ch2co=0
126      igcm_ch3co=0
127      igcm_hcaer=0
128
129
130      write(*,*) 'initracer: noms() ', noms
131
132
133      !print*,'Setting dustbin = 0 in initracer.F'
134      !dustbin=0
135
136      ! 1. find dust tracers
137      count=0
138!      if (dustbin.gt.0) then
139!        do iq=1,nq
140!          txt=" "
141!          write(txt,'(a4,i2.2)')'dust',count+1   
142!          if (noms(iq).eq.txt) then
143!            count=count+1
144!            igcm_dustbin(count)=iq
145!            mmol(iq)=100.
146!          endif
147!        enddo !do iq=1,nq
148!      endif ! of if (dustbin.gt.0)
149
150
151!      if (doubleq) then
152!        do iq=1,nq
153!          if (noms(iq).eq."dust_mass") then
154!            igcm_dust_mass=iq
155!            count=count+1
156!          endif
157!          if (noms(iq).eq."dust_number") then
158!            igcm_dust_number=iq
159!            count=count+1
160!          endif
161!        enddo
162!      endif ! of if (doubleq)
163      ! 2. find chemistry and water tracers
164      do iq=1,nq
165        if (noms(iq).eq."co2") then
166          igcm_co2=iq
167          mmol(igcm_co2)=44.
168          count=count+1
169!          write(*,*) 'co2: count=',count
170        endif
171        if (noms(iq).eq."co2_ice") then
172          igcm_co2_ice=iq
173          mmol(igcm_co2_ice)=44.
174          count=count+1
175!          write(*,*) 'co2_ice: count=',count
176        endif
177        if (noms(iq).eq."h2o_vap") then
178          igcm_h2o_vap=iq
179          mmol(igcm_h2o_vap)=18.
180          count=count+1
181!          write(*,*) 'h2o_vap: count=',count
182        endif
183        if (noms(iq).eq."h2o_ice") then
184          igcm_h2o_ice=iq
185          mmol(igcm_h2o_ice)=18.
186          count=count+1
187!          write(*,*) 'h2o_ice: count=',count
188        endif
189        if (noms(iq).eq."co") then
190          igcm_co=iq
191          mmol(igcm_co)=28.
192          count=count+1
193        endif
194        if (noms(iq).eq."o") then
195          igcm_o=iq
196          mmol(igcm_o)=16.
197          count=count+1
198        endif
199        if (noms(iq).eq."o1d") then
200          igcm_o1d=iq
201          mmol(igcm_o1d)=16.
202          count=count+1
203        endif
204        if (noms(iq).eq."o2") then
205          igcm_o2=iq
206          mmol(igcm_o2)=32.
207          count=count+1
208        endif
209        if (noms(iq).eq."o3") then
210          igcm_o3=iq
211          mmol(igcm_o3)=48.
212          count=count+1
213        endif
214        if (noms(iq).eq."h") then
215          igcm_h=iq
216          mmol(igcm_h)=1.
217          count=count+1
218        endif
219        if (noms(iq).eq."h2") then
220          igcm_h2=iq
221          mmol(igcm_h2)=2.
222          count=count+1
223        endif
224        if (noms(iq).eq."oh") then
225          igcm_oh=iq
226          mmol(igcm_oh)=17.
227          count=count+1
228        endif
229        if (noms(iq).eq."ho2") then
230          igcm_ho2=iq
231          mmol(igcm_ho2)=33.
232          count=count+1
233        endif
234        if (noms(iq).eq."h2o2") then
235          igcm_h2o2=iq
236          mmol(igcm_h2o2)=34.
237          count=count+1
238        endif
239        if (noms(iq).eq."n2") then
240          igcm_n2=iq
241          mmol(igcm_n2)=28.
242          count=count+1
243        endif
244        if (noms(iq).eq."ch4") then
245          igcm_ch4=iq
246          mmol(igcm_ch4)=16.
247          count=count+1
248        endif
249        if (noms(iq).eq."ar") then
250          igcm_ar=iq
251          mmol(igcm_ar)=40.
252          count=count+1
253        endif
254        if (noms(iq).eq."n") then
255          igcm_n=iq
256          mmol(igcm_n)=14.
257          count=count+1
258        endif
259        if (noms(iq).eq."no") then
260          igcm_no=iq
261          mmol(igcm_no)=30.
262          count=count+1
263        endif
264        if (noms(iq).eq."no2") then
265          igcm_no2=iq
266          mmol(igcm_no2)=46.
267          count=count+1
268        endif
269        if (noms(iq).eq."n2d") then
270          igcm_n2d=iq
271          mmol(igcm_n2d)=28.
272          count=count+1
273        endif
274        if (noms(iq).eq."ch3") then
275          igcm_ch3=iq
276          mmol(igcm_ch3)=15.
277          count=count+1
278        endif
279        if (noms(iq).eq."ch") then
280          igcm_ch=iq
281          mmol(igcm_ch)=13.
282          count=count+1
283        endif
284        if (noms(iq).eq."3ch2") then
285          igcm_3ch2=iq
286          mmol(igcm_3ch2)=14.
287          count=count+1
288        endif
289        if (noms(iq).eq."1ch2") then
290          igcm_1ch2=iq
291          mmol(igcm_1ch2)=14.
292          count=count+1
293        endif
294        if (noms(iq).eq."cho") then
295          igcm_cho=iq
296          mmol(igcm_cho)=29.
297          count=count+1
298        endif
299        if (noms(iq).eq."ch2o") then
300          igcm_ch2o=iq
301          mmol(igcm_ch2o)=30.
302          count=count+1
303        endif
304        if (noms(iq).eq."ch3o") then
305          igcm_ch3o=iq
306          mmol(igcm_ch3o)=31.
307          count=count+1
308        endif
309        if (noms(iq).eq."c") then
310          igcm_c=iq
311          mmol(igcm_c)=12.
312          count=count+1
313        endif
314        if (noms(iq).eq."c2") then
315          igcm_c2=iq
316          mmol(igcm_c2)=24.
317          count=count+1
318        endif
319        if (noms(iq).eq."c2h") then
320          igcm_c2h=iq
321          mmol(igcm_c2h)=25.
322          count=count+1
323        endif
324        if (noms(iq).eq."c2h2") then
325          igcm_c2h2=iq
326          mmol(igcm_c2h2)=26.
327          count=count+1
328        endif
329        if (noms(iq).eq."c2h3") then
330          igcm_c2h3=iq
331          mmol(igcm_c2h3)=27.
332          count=count+1
333        endif
334        if (noms(iq).eq."c2h4") then
335          igcm_c2h4=iq
336          mmol(igcm_c2h4)=28.
337          count=count+1
338        endif
339        if (noms(iq).eq."c2h6") then
340          igcm_c2h6=iq
341          mmol(igcm_c2h6)=30.
342          count=count+1
343        endif
344        if (noms(iq).eq."ch2co") then
345          igcm_ch2co=iq
346          mmol(igcm_ch2co)=42.
347          count=count+1
348        endif
349        if (noms(iq).eq."ch3co") then
350          igcm_ch3co=iq
351          mmol(igcm_ch3co)=43.
352          count=count+1
353        endif
354        if (noms(iq).eq."hcaer") then
355          igcm_hcaer=iq
356          mmol(igcm_hcaer)=50.
357          count=count+1
358        endif
359      enddo ! of do iq=1,nq
360     
361      ! check that we identified all tracers:
362      if (count.ne.nq) then
363        write(*,*) "initracer: found only ",count," tracers"
364        write(*,*) "               expected ",nq
365        do iq=1,count
366          write(*,*)'      ',iq,' ',trim(noms(iq))
367        enddo
368!        stop
369      else
370        write(*,*) "initracer: found all expected tracers, namely:"
371        do iq=1,nq
372          write(*,*)'      ',iq,' ',trim(noms(iq))
373        enddo
374      endif
375
376
377c------------------------------------------------------------
378c     Initialisation tracers ....
379c------------------------------------------------------------
380      call zerophys(nq,rho_q)
381
382      rho_dust=2500.  ! Mars dust density (kg.m-3)
383      rho_ice=920.    ! Water ice density (kg.m-3)
384      rho_co2=1620.   ! CO2 ice density (kg.m-3)
385
386
387
388c$$$      if (doubleq) then
389c$$$c       "doubleq" technique
390c$$$c       -------------------
391c$$$c      (transport of mass and number mixing ratio)
392c$$$c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
393c$$$
394c$$$        if( (nq.lt.2).or.(water.and.(nq.lt.3)) ) then
395c$$$          write(*,*)'initracer: nq is too low : nq=', nq
396c$$$          write(*,*)'water= ',water,' doubleq= ',doubleq   
397c$$$        end if
398c$$$
399c$$$        varian=0.637    ! Characteristic variance   
400c$$$        qext(igcm_dust_mass)=3.04   ! reference extinction at 0.67 um for ref dust
401c$$$        qext(igcm_dust_number)=3.04 ! reference extinction at 0.67 um for ref dust
402c$$$        rho_q(igcm_dust_mass)=rho_dust
403c$$$        rho_q(igcm_dust_number)=rho_dust
404c$$$
405c$$$c       Intermediate calcul for computing geometric mean radius r0
406c$$$c       as a function of mass and number mixing ratio Q and N
407c$$$c       (r0 = (r3n_q * Q/ N)
408c$$$        r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
409c$$$
410c$$$c       Intermediate calcul for computing effective radius reff
411c$$$c       from geometric mean radius r0
412c$$$c       (reff = ref_r0 * r0)
413c$$$        ref_r0 = exp(2.5*varian**2)
414c$$$       
415c$$$c       lifted dust :
416c$$$c       '''''''''''
417c$$$        reff_lift = 3.e-6      !  Effective radius of lifted dust (m)
418c$$$        alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
419c$$$        alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff
420c$$$
421c$$$        r0_lift = reff_lift/ref_r0
422c$$$        alpha_devil(igcm_dust_number)=r3n_q*
423c$$$     &                        alpha_devil(igcm_dust_mass)/r0_lift**3
424c$$$        alpha_lift(igcm_dust_number)=r3n_q*
425c$$$     &                        alpha_lift(igcm_dust_mass)/r0_lift**3
426c$$$
427c$$$c       Not used:
428c$$$        radius(igcm_dust_mass) = 0.
429c$$$        radius(igcm_dust_number) = 0.
430c$$$
431c$$$      else
432
433
434c$$$       if (dustbin.gt.1) then
435c$$$        print*,'ATTENTION:',
436c$$$     $   ' properties of dust need input in initracer !!!'
437c$$$        stop
438c$$$
439c$$$       else if (dustbin.eq.1) then
440c$$$
441c$$$c       This will be used for 1 dust particle size:
442c$$$c       ------------------------------------------
443c$$$        radius(igcm_dustbin(1))=3.e-6
444c$$$        Qext(igcm_dustbin(1))=3.04
445c$$$        alpha_lift(igcm_dustbin(1))=0.0e-6
446c$$$        alpha_devil(igcm_dustbin(1))=7.65e-9
447c$$$        qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1))
448c$$$     &                         /(rho_dust*radius(igcm_dustbin(1)))
449c$$$        rho_q(igcm_dustbin(1))=rho_dust
450c$$$
451c$$$       endif
452c$$$!      end if    ! (doubleq)
453
454c     Initialization for water vapor
455c     ------------------------------
456      if(water) then
457         radius(igcm_h2o_vap)=0.
458         Qext(igcm_h2o_vap)=0.
459         alpha_lift(igcm_h2o_vap) =0.
460         alpha_devil(igcm_h2o_vap)=0.
461         qextrhor(igcm_h2o_vap)= 0.
462
463c       "Dryness coefficient" controlling the evaporation and
464c        sublimation from the ground water ice (close to 1)
465c        HERE, the goal is to correct for the fact
466c        that the simulated permanent water ice polar caps
467c        is larger than the actual cap and the atmospheric
468c        opacity not always realistic.
469
470
471!         if(ngrid.eq.1)
472
473
474!     to be modified for BC+ version?
475
476         !! this is defined in surfdat_h.F90
477         IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
478         IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
479
480         do ig=1,ngrid
481           if (ngrid.ne.1) watercaptag(ig)=.false.
482           dryness(ig) = 1.
483         enddo
484
485
486
487
488!         IF (caps) THEN
489c Perennial H20 north cap defined by watercaptag=true (allows surface to be
490c hollowed by sublimation in vdifc).
491!         do ig=1,ngrid
492!           if (lati(ig)*180./pi.gt.84) then
493!             if (ngrid.ne.1) watercaptag(ig)=.true.
494!             dryness(ig) = 1.
495c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
496c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
497c               if (ngrid.ne.1) watercaptag(ig)=.true.
498c               dryness(ig) = 1.
499c             endif
500c             if (lati(ig)*180./pi.ge.85) then
501c               if (ngrid.ne.1) watercaptag(ig)=.true.
502c               dryness(ig) = 1.
503c             endif
504!           endif  ! (lati>80 deg)
505!         end do ! (ngrid)
506!        ENDIF ! (caps)
507
508!         if(iceparty.and.(nq.ge.2)) then
509
510           radius(igcm_h2o_ice)=3.e-6
511           rho_q(igcm_h2o_ice)=rho_ice
512           Qext(igcm_h2o_ice)=0.
513!           alpha_lift(igcm_h2o_ice) =0.
514!           alpha_devil(igcm_h2o_ice)=0.
515           qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice)
516     $       / (rho_ice*radius(igcm_h2o_ice))
517
518
519
520!         elseif(iceparty.and.(nq.lt.2)) then
521!            write(*,*) 'nq is too low : nq=', nq
522!            write(*,*) 'water= ',water,' iceparty= ',iceparty
523!         endif
524
525      end if  ! (water)
526
527
528!
529!     some extra (possibly redundant) sanity checks for tracers:
530!     ---------------------------------------------------------
531       if (water) then
532       ! verify that we indeed have h2o_vap and h2o_ice tracers
533         if (igcm_h2o_vap.eq.0) then
534           write(*,*) "initracer: error !!"
535           write(*,*) "  cannot use water option without ",
536     &                "an h2o_vap tracer !"
537           stop
538         endif
539         if (igcm_h2o_ice.eq.0) then
540           write(*,*) "initracer: error !!"
541           write(*,*) "  cannot use water option without ",
542     &                "an h2o_ice tracer !"
543           stop
544         endif
545       endif
546
547
548c     Output for records:
549c     ~~~~~~~~~~~~~~~~~~
550      write(*,*)
551      Write(*,*) '******** initracer : dust transport parameters :'
552      write(*,*) 'alpha_lift = ', alpha_lift
553      write(*,*) 'alpha_devil = ', alpha_devil
554      write(*,*) 'radius  = ', radius
555!      if(doubleq) then
556!        write(*,*) 'reff_lift (um) =  ', reff_lift
557!        write(*,*) 'size distribution variance  = ', varian
558!        write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
559!      end if
560      write(*,*) 'Qext  = ', qext
561      write(*,*)
562
563      end
Note: See TracBrowser for help on using the repository browser.