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

Last change on this file since 2176 was 1980, checked in by emillour, 6 years ago

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