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

Last change on this file since 2546 was 2543, checked in by aslmd, 4 years ago

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ
File size: 18.3 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      USE recombin_corrk_mod, ONLY: ini_recombin
7      IMPLICIT NONE
8c=======================================================================
9c   subject:
10c   --------
11c   Initialization related to tracer
12c   (transported dust, water, chemical species, ice...)
13c
14c   Name of the tracer
15c
16c   Test of dimension :
17c   Initialize COMMON tracer in tracer.h, using tracer names provided
18c   by the argument nametrac
19c
20c   author: F.Forget
21c   ------
22c            Ehouarn Millour (oct. 2008) identify tracers by their names
23c            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
24c=======================================================================
25
26      integer,intent(in) :: ngrid,nq
27      character(len=30),intent(in) :: nametrac(nq) ! name of the tracer from dynamics
28
29      character(len=500) :: tracline   ! to read traceur.def lines
30      character(len=30)  :: txt        ! to store some text
31      integer iq,ig,count,ierr
32
33
34
35c-----------------------------------------------------------------------
36c  radius(nq)      ! aerosol particle radius (m)
37c  rho_q(nq)       ! tracer densities (kg.m-3)
38c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
39c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
40c  alpha_devil(nq) ! lifting coeeficient by dust devil
41c  rho_dust          ! Mars dust density
42c  rho_ice           ! Water ice density
43c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
44c  varian            ! Characteristic variance of log-normal distribution
45c-----------------------------------------------------------------------
46
47      moderntracdef=.false. ! For modern traceur.def (default false, old type)
48
49      open(407, form = 'formatted', status = 'old',
50     $          file = 'traceur.def', iostat=ierr)
51      if (ierr /=0) then
52        call abort_physic('initracer',
53     $  'Problem in opening traceur.def',1)
54      end if
55!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
56       READ(407,'(A)') tracline
57       IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
58          READ(tracline,*) nqtot ! Try standard traceur.def
59       ELSE
60         moderntracdef = .true.
61         DO
62           READ(407,'(A)',iostat=ierr) tracline
63           IF (ierr==0) THEN
64             IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
65               READ(tracline,*) nqtot
66               ! Temporary not implemented solution
67               if (nqtot/=nq) then
68                 call abort_physic('initracer','Different number of '//
69     & 'tracers in dynamics and physics not managed yet',1)
70               endif
71               EXIT
72             ENDIF
73           ELSE ! If pb, or if reached EOF without having found nqtot
74             call abort_physic('initracer','Unable to read numbers '//
75     & 'of tracers in traceur.def',1)
76           ENDIF
77         ENDDO
78       ENDIF ! if modern or standard traceur.def
79!! -----------------------------------------------------------------------
80       !! For the moment number of tracers in dynamics and physics are equal
81       nqtot=nq
82       !! we allocate once for all arrays in common in tracer_h.F90
83       !! (supposedly those are not used before call to initracer)
84       IF (.NOT.ALLOCATED(noms))           ALLOCATE(noms(nq))
85       IF (.NOT.ALLOCATED(mmol))           ALLOCATE(mmol(nq))
86       IF (.NOT.ALLOCATED(aki))            ALLOCATE(aki(nqtot))
87       IF (.NOT.ALLOCATED(cpi))            ALLOCATE(cpi(nqtot))
88       IF (.NOT.ALLOCATED(radius))         ALLOCATE(radius(nq))
89       IF (.NOT.ALLOCATED(rho_q))          ALLOCATE(rho_q(nq))
90       IF (.NOT.ALLOCATED(qext))           ALLOCATE(qext(nq))
91       IF (.NOT.ALLOCATED(alpha_lift))     ALLOCATE(alpha_lift(nq))
92       IF (.NOT.ALLOCATED(alpha_devil))    ALLOCATE(alpha_devil(nq))
93       IF (.NOT.ALLOCATED(qextrhor))       ALLOCATE(qextrhor(nq))
94       IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
95       IF (.NOT.ALLOCATED(is_chim))        ALLOCATE(is_chim(nqtot))
96       IF (.NOT.ALLOCATED(is_rad))         ALLOCATE(is_rad(nqtot))
97       IF (.NOT.ALLOCATED(is_recomb))      ALLOCATE(is_recomb(nqtot))
98       IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
99         ALLOCATE(is_recomb_qset(nqtot))
100       ENDIF
101       IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
102         ALLOCATE(is_recomb_qotf(nqtot))
103       ENDIF
104       !! initialization
105       alpha_lift(:)     = 0.
106       alpha_devil(:)    = 0.
107       mmol(:)           = 0.
108       aki(:)            = 0.
109       cpi(:)            = 0.
110       is_chim(:)        = 0
111       is_rad(:)         = 0
112       is_recomb(:)      = 0
113       is_recomb_qset(:) = 0
114       is_recomb_qotf(:) = 0
115       
116       ! Added by JVO 2017 : these arrays are handled later
117       ! -> initialization is the least we can do, please !!!
118       radius(:)=0.
119       qext(:)=0.
120
121
122! Initialization: copy tracer names from dynamics
123        do iq=1,nq
124          noms(iq)=nametrac(iq)
125          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
126        enddo
127
128
129! Identify tracers by their names: (and set corresponding values of mmol)
130      ! 0. initialize tracer indexes to zero:
131      ! NB: igcm_* indexes are commons in 'tracer.h'
132      do iq=1,nq
133        igcm_dustbin(iq)=0
134      enddo
135      igcm_dust_mass=0
136      igcm_dust_number=0
137      igcm_h2o_vap=0
138      igcm_h2o_ice=0
139      igcm_co2=0
140      igcm_co=0
141      igcm_o=0
142      igcm_o1d=0
143      igcm_o2=0
144      igcm_o3=0
145      igcm_h=0
146      igcm_h2=0
147      igcm_oh=0
148      igcm_ho2=0
149      igcm_h2o2=0
150      igcm_n2=0
151      igcm_n=0
152      igcm_n2d=0
153      igcm_no=0
154      igcm_no2=0
155      igcm_ar=0
156      igcm_ar_n2=0
157      igcm_co2_ice=0
158
159      igcm_ch4=0
160      igcm_ch3=0
161      igcm_ch=0
162      igcm_3ch2=0
163      igcm_1ch2=0
164      igcm_cho=0
165      igcm_ch2o=0
166      igcm_ch3o=0
167      igcm_c=0
168      igcm_c2=0
169      igcm_c2h=0
170      igcm_c2h2=0
171      igcm_c2h3=0
172      igcm_c2h4=0
173      igcm_c2h6=0
174      igcm_ch2co=0
175      igcm_ch3co=0
176      igcm_hcaer=0
177
178
179
180      ! 1. find dust tracers
181      count=0
182
183      ! 2. find chemistry and water tracers
184      do iq=1,nq
185        if (noms(iq).eq."co2") then
186          igcm_co2=iq
187          mmol(igcm_co2)=44.
188          count=count+1
189!          write(*,*) 'co2: count=',count
190        endif
191        if (noms(iq).eq."co2_ice") then
192          igcm_co2_ice=iq
193          mmol(igcm_co2_ice)=44.
194          count=count+1
195!          write(*,*) 'co2_ice: count=',count
196        endif
197        if (noms(iq).eq."h2o_vap") then
198          igcm_h2o_vap=iq
199          mmol(igcm_h2o_vap)=18.
200          count=count+1
201!          write(*,*) 'h2o_vap: count=',count
202        endif
203        if (noms(iq).eq."h2o_ice") then
204          igcm_h2o_ice=iq
205          mmol(igcm_h2o_ice)=18.
206          count=count+1
207!          write(*,*) 'h2o_ice: count=',count
208        endif
209        if (noms(iq).eq."co") then
210          igcm_co=iq
211          mmol(igcm_co)=28.
212          count=count+1
213        endif
214        if (noms(iq).eq."o") then
215          igcm_o=iq
216          mmol(igcm_o)=16.
217          count=count+1
218        endif
219        if (noms(iq).eq."o1d") then
220          igcm_o1d=iq
221          mmol(igcm_o1d)=16.
222          count=count+1
223        endif
224        if (noms(iq).eq."o2") then
225          igcm_o2=iq
226          mmol(igcm_o2)=32.
227          count=count+1
228        endif
229        if (noms(iq).eq."o3") then
230          igcm_o3=iq
231          mmol(igcm_o3)=48.
232          count=count+1
233        endif
234        if (noms(iq).eq."h") then
235          igcm_h=iq
236          mmol(igcm_h)=1.
237          count=count+1
238        endif
239        if (noms(iq).eq."h2") then
240          igcm_h2=iq
241          mmol(igcm_h2)=2.
242          count=count+1
243        endif
244        if (noms(iq).eq."oh") then
245          igcm_oh=iq
246          mmol(igcm_oh)=17.
247          count=count+1
248        endif
249        if (noms(iq).eq."ho2") then
250          igcm_ho2=iq
251          mmol(igcm_ho2)=33.
252          count=count+1
253        endif
254        if (noms(iq).eq."h2o2") then
255          igcm_h2o2=iq
256          mmol(igcm_h2o2)=34.
257          count=count+1
258        endif
259        if (noms(iq).eq."n2") then
260          igcm_n2=iq
261          mmol(igcm_n2)=28.
262          count=count+1
263        endif
264        if (noms(iq).eq."ch4") then
265          igcm_ch4=iq
266          mmol(igcm_ch4)=16.
267          count=count+1
268        endif
269        if (noms(iq).eq."ar") then
270          igcm_ar=iq
271          mmol(igcm_ar)=40.
272          count=count+1
273        endif
274        if (noms(iq).eq."n") then
275          igcm_n=iq
276          mmol(igcm_n)=14.
277          count=count+1
278        endif
279        if (noms(iq).eq."no") then
280          igcm_no=iq
281          mmol(igcm_no)=30.
282          count=count+1
283        endif
284        if (noms(iq).eq."no2") then
285          igcm_no2=iq
286          mmol(igcm_no2)=46.
287          count=count+1
288        endif
289        if (noms(iq).eq."n2d") then
290          igcm_n2d=iq
291          mmol(igcm_n2d)=28.
292          count=count+1
293        endif
294        if (noms(iq).eq."ch3") then
295          igcm_ch3=iq
296          mmol(igcm_ch3)=15.
297          count=count+1
298        endif
299        if (noms(iq).eq."ch") then
300          igcm_ch=iq
301          mmol(igcm_ch)=13.
302          count=count+1
303        endif
304        if (noms(iq).eq."3ch2") then
305          igcm_3ch2=iq
306          mmol(igcm_3ch2)=14.
307          count=count+1
308        endif
309        if (noms(iq).eq."1ch2") then
310          igcm_1ch2=iq
311          mmol(igcm_1ch2)=14.
312          count=count+1
313        endif
314        if (noms(iq).eq."cho") then
315          igcm_cho=iq
316          mmol(igcm_cho)=29.
317          count=count+1
318        endif
319        if (noms(iq).eq."ch2o") then
320          igcm_ch2o=iq
321          mmol(igcm_ch2o)=30.
322          count=count+1
323        endif
324        if (noms(iq).eq."ch3o") then
325          igcm_ch3o=iq
326          mmol(igcm_ch3o)=31.
327          count=count+1
328        endif
329        if (noms(iq).eq."c") then
330          igcm_c=iq
331          mmol(igcm_c)=12.
332          count=count+1
333        endif
334        if (noms(iq).eq."c2") then
335          igcm_c2=iq
336          mmol(igcm_c2)=24.
337          count=count+1
338        endif
339        if (noms(iq).eq."c2h") then
340          igcm_c2h=iq
341          mmol(igcm_c2h)=25.
342          count=count+1
343        endif
344        if (noms(iq).eq."c2h2") then
345          igcm_c2h2=iq
346          mmol(igcm_c2h2)=26.
347          count=count+1
348        endif
349        if (noms(iq).eq."c2h3") then
350          igcm_c2h3=iq
351          mmol(igcm_c2h3)=27.
352          count=count+1
353        endif
354        if (noms(iq).eq."c2h4") then
355          igcm_c2h4=iq
356          mmol(igcm_c2h4)=28.
357          count=count+1
358        endif
359        if (noms(iq).eq."c2h6") then
360          igcm_c2h6=iq
361          mmol(igcm_c2h6)=30.
362          count=count+1
363        endif
364        if (noms(iq).eq."ch2co") then
365          igcm_ch2co=iq
366          mmol(igcm_ch2co)=42.
367          count=count+1
368        endif
369        if (noms(iq).eq."ch3co") then
370          igcm_ch3co=iq
371          mmol(igcm_ch3co)=43.
372          count=count+1
373        endif
374        if (noms(iq).eq."hcaer") then
375          igcm_hcaer=iq
376          mmol(igcm_hcaer)=50.
377          count=count+1
378        endif
379      enddo ! of do iq=1,nq
380     
381      ! check that we identified all tracers:
382      if (count.ne.nq) then
383        write(*,*) "initracer: found only ",count," tracers"
384        write(*,*) "               expected ",nq
385        do iq=1,count
386          write(*,*)'      ',iq,' ',trim(noms(iq))
387        enddo
388!        stop
389      else
390        write(*,*) "initracer: found all expected tracers, namely:"
391        do iq=1,nq
392          write(*,*)'      ',iq,' ',trim(noms(iq))
393        enddo
394      endif
395
396      ! Get data of tracers
397      do iq=1,nqtot
398        read(407,'(A)') tracline
399        call get_tracdat(iq, tracline)
400      enddo
401
402      close(407)
403
404      ! Calculate number of species in the chemistry
405      nesp = sum(is_chim)
406      write(*,*) 'Number of species in the chemistry nesp = ',nesp
407
408!     Processing modern traceur options
409      if(moderntracdef) then
410        call ini_recombin
411      endif
412
413c------------------------------------------------------------
414c     Initialisation tracers ....
415c------------------------------------------------------------
416      rho_q(1:nq)=0
417
418      rho_dust=2500.  ! Mars dust density (kg.m-3)
419      rho_ice=920.    ! Water ice density (kg.m-3)
420      rho_co2=1620.   ! CO2 ice density (kg.m-3)
421
422c     Initialization for water vapor
423c     ------------------------------
424      if(water) then
425         radius(igcm_h2o_vap)=0.
426         Qext(igcm_h2o_vap)=0.
427         alpha_lift(igcm_h2o_vap) =0.
428         alpha_devil(igcm_h2o_vap)=0.
429         qextrhor(igcm_h2o_vap)= 0.
430
431         !! this is defined in surfdat_h.F90
432         IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
433         IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
434
435         do ig=1,ngrid
436           if (ngrid.ne.1) watercaptag(ig)=.false.
437           dryness(ig) = 1.
438         enddo
439
440
441           radius(igcm_h2o_ice)=3.e-6
442           rho_q(igcm_h2o_ice)=rho_ice
443           Qext(igcm_h2o_ice)=0.
444!           alpha_lift(igcm_h2o_ice) =0.
445!           alpha_devil(igcm_h2o_ice)=0.
446           qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice)
447     $       / (rho_ice*radius(igcm_h2o_ice))
448
449
450      end if  ! (water)
451
452
453!
454!     some extra (possibly redundant) sanity checks for tracers:
455!     ---------------------------------------------------------
456       if (water) then
457       ! verify that we indeed have h2o_vap and h2o_ice tracers
458         if (igcm_h2o_vap.eq.0) then
459           write(*,*) "initracer: error !!"
460           write(*,*) "  cannot use water option without ",
461     &                "an h2o_vap tracer !"
462           stop
463         endif
464         if (igcm_h2o_ice.eq.0) then
465           write(*,*) "initracer: error !!"
466           write(*,*) "  cannot use water option without ",
467     &                "an h2o_ice tracer !"
468           stop
469         endif
470       endif
471
472
473c     Output for records:
474c     ~~~~~~~~~~~~~~~~~~
475      write(*,*)
476      Write(*,*) '******** initracer : dust transport parameters :'
477      write(*,*) 'alpha_lift = ', alpha_lift
478      write(*,*) 'alpha_devil = ', alpha_devil
479      write(*,*) 'radius  = ', radius
480      write(*,*) 'Qext  = ', qext
481      write(*,*)
482
483      contains
484
485      subroutine get_tracdat(iq,tracline)
486        !-------------------ADDING NEW OPTIONS-------------------
487        ! Duplicate if sentence for adding new options
488        ! Do not forget to add the new variables in tracer_h.F90
489        ! Do not forget to allocate and initialize the new variables above
490        ! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
491        !-------------------------------------------------------
492          use tracer_h
493          implicit none
494          integer,           intent(in) :: iq       ! tracer index
495          character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
496 
497          read(tracline,*) noms(iq)
498          ! JVO 20 : We should add a sanity check aborting when duplicates in names !
499          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
500          ! option mmol
501          if (index(tracline,'mmol='   ) /= 0) then
502              read(tracline(index(tracline,'mmol=')+len('mmol='):),*)
503     $             mmol(iq)
504              write(*,*) ' Parameter value (traceur.def) : mmol=',
505     $             mmol(iq)
506          else
507              write(*,*) ' Parameter value (default)     : mmol=',
508     $             mmol(iq)
509          end if
510          ! option aki
511          if (index(tracline,'aki='   ) /= 0) then
512              read(tracline(index(tracline,'aki=')+len('aki='):),*)
513     $             aki(iq)
514              write(*,*) ' Parameter value (traceur.def) : aki=',
515     $             aki(iq)
516          else
517              write(*,*) ' Parameter value (default)     : aki=',
518     $             aki(iq)
519          end if
520          ! option cpi
521          if (index(tracline,'cpi='   ) /= 0) then
522              read(tracline(index(tracline,'cpi=')+len('cpi='):),*)
523     $             cpi(iq)
524              write(*,*) ' Parameter value (traceur.def) : cpi=',
525     $             cpi(iq)
526          else
527              write(*,*) ' Parameter value (default)     : cpi=',
528     $             cpi(iq)
529          end if
530          ! option is_chim
531          if (index(tracline,'is_chim='   ) /= 0) then
532          read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*)
533     $             is_chim(iq)
534              write(*,*) ' Parameter value (traceur.def) : is_chim=',
535     $             is_chim(iq)
536          else
537              write(*,*) ' Parameter value (default)     : is_chim=',
538     $             is_chim(iq)
539          end if
540          ! option is_rad
541          if (index(tracline,'is_rad=') /= 0) then
542              read(tracline(index(tracline,'is_rad=')
543     $         +len('is_rad='):),*) is_rad(iq)
544              write(*,*) ' Parameter value (traceur.def) : is_rad=',
545     $         is_rad(iq)
546          else
547              write(*,*) ' Parameter value (default)     : is_rad=',
548     $         is_rad(iq)
549          end if
550          ! option is_recomb
551          if (index(tracline,'is_recomb=') /= 0) then
552              read(tracline(index(tracline,'is_recomb=')
553     $         +len('is_recomb='):),*) is_recomb(iq)
554              write(*,*) ' Parameter value (traceur.def) : is_recomb=',
555     $         is_recomb(iq)
556          else
557              write(*,*) ' Parameter value (default)     : is_recomb=',
558     $         is_recomb(iq)
559          end if
560          ! option is_recomb_qset
561          if (index(tracline,'is_recomb_qset=') /= 0) then
562              read(tracline(index(tracline,'is_recomb_qset=')
563     $         +len('is_recomb_qset='):),*) is_recomb_qset(iq)
564              write(*,*) ' Parameter value (traceur.def) :'//
565     $         ' is_recomb_qset=',
566     $         is_recomb_qset(iq)
567          else
568              write(*,*) ' Parameter value (default)     :'//
569     $         ' is_recomb_qset=',
570     $         is_recomb_qset(iq)
571          endif
572          ! option is_recomb_qotf
573          if (index(tracline,'is_recomb_qotf=') /= 0) then
574              read(tracline(index(tracline,'is_recomb_qotf=')
575     $         +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
576              write(*,*) ' Parameter value (traceur.def) :'//
577     $         ' is_recomb_qotf=',
578     $         is_recomb_qotf(iq)
579          else
580              write(*,*) ' Parameter value (default)     :'//
581     $         ' is_recomb_qotf=',
582     $         is_recomb_qotf(iq)
583          end if
584      end subroutine get_tracdat
585
586      end
587
Note: See TracBrowser for help on using the repository browser.