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
RevLine 
[787]1      SUBROUTINE initracer(ngrid,nq,nametrac)
[135]2
[1980]3      use surfdat_h, ONLY: dryness, watercaptag
[787]4      USE tracer_h
[1397]5      USE callkeys_mod, only: water
[2543]6      USE recombin_corrk_mod, ONLY: ini_recombin
[135]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
[787]18c   by the argument nametrac
[135]19c
20c   author: F.Forget
21c   ------
22c            Ehouarn Millour (oct. 2008) identify tracers by their names
[2436]23c            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
[135]24c=======================================================================
25
[1980]26      integer,intent(in) :: ngrid,nq
27      character(len=30),intent(in) :: nametrac(nq) ! name of the tracer from dynamics
[787]28
[2436]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
[135]32
33
[787]34
[135]35c-----------------------------------------------------------------------
[787]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
[135]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
[2436]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
[2542]68                 call abort_physic('initracer','Different number of '//
69     & 'tracers in dynamics and physics not managed yet',1)
[2436]70               endif
71               EXIT
72             ENDIF
73           ELSE ! If pb, or if reached EOF without having found nqtot
[2542]74             call abort_physic('initracer','Unable to read numbers '//
75     & 'of tracers in traceur.def',1)
[2436]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
[1621]81       nqtot=nq
[787]82       !! we allocate once for all arrays in common in tracer_h.F90
83       !! (supposedly those are not used before call to initracer)
[2543]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
[787]104       !! initialization
[2543]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
[1764]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.
[787]120
[135]121
[1980]122! Initialization: copy tracer names from dynamics
[787]123        do iq=1,nq
124          noms(iq)=nametrac(iq)
[1980]125          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
[135]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'
[787]132      do iq=1,nq
[135]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
[1801]151      igcm_n=0
152      igcm_n2d=0
153      igcm_no=0
154      igcm_no2=0
[135]155      igcm_ar=0
156      igcm_ar_n2=0
157      igcm_co2_ice=0
158
[1801]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
[135]179
180      ! 1. find dust tracers
181      count=0
182
183      ! 2. find chemistry and water tracers
[787]184      do iq=1,nq
[135]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
[1801]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
[787]379      enddo ! of do iq=1,nq
[135]380     
381      ! check that we identified all tracers:
[787]382      if (count.ne.nq) then
[135]383        write(*,*) "initracer: found only ",count," tracers"
[787]384        write(*,*) "               expected ",nq
[135]385        do iq=1,count
386          write(*,*)'      ',iq,' ',trim(noms(iq))
387        enddo
[1297]388!        stop
[135]389      else
390        write(*,*) "initracer: found all expected tracers, namely:"
[787]391        do iq=1,nq
[135]392          write(*,*)'      ',iq,' ',trim(noms(iq))
393        enddo
394      endif
395
[2436]396      ! Get data of tracers
397      do iq=1,nqtot
398        read(407,'(A)') tracline
399        call get_tracdat(iq, tracline)
400      enddo
[135]401
[2436]402      close(407)
403
[2542]404      ! Calculate number of species in the chemistry
405      nesp = sum(is_chim)
406      write(*,*) 'Number of species in the chemistry nesp = ',nesp
407
[2543]408!     Processing modern traceur options
409      if(moderntracdef) then
410        call ini_recombin
411      endif
412
[135]413c------------------------------------------------------------
414c     Initialisation tracers ....
415c------------------------------------------------------------
[2278]416      rho_q(1:nq)=0
[135]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
[787]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.
[135]437           dryness(ig) = 1.
438         enddo
439
[253]440
[135]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
[1801]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
[135]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
[2436]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))
[2542]500          ! option mmol
[2436]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
[2542]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
[2543]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
[2436]584      end subroutine get_tracdat
585
[135]586      end
[2436]587
Note: See TracBrowser for help on using the repository browser.