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

Last change on this file since 2499 was 2436, checked in by yjaziri, 5 years ago

Generic GCM:
Implementation of an option for a new reading process of "traceur.def"
Use "#ModernTrac-v1" flag as first line of "traceur.def" to use this option
Further details in "LMDZ.GENERIC/deftank/traceur.def.modern"
YJ + JVO
File size: 14.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            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
23c=======================================================================
24
25      integer,intent(in) :: ngrid,nq
26      character(len=30),intent(in) :: nametrac(nq) ! name of the tracer from dynamics
27
28      character(len=500) :: tracline   ! to read traceur.def lines
29      character(len=30)  :: txt        ! to store some text
30      integer iq,ig,count,ierr
31
32
33
34c-----------------------------------------------------------------------
35c  radius(nq)      ! aerosol particle radius (m)
36c  rho_q(nq)       ! tracer densities (kg.m-3)
37c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
38c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
39c  alpha_devil(nq) ! lifting coeeficient by dust devil
40c  rho_dust          ! Mars dust density
41c  rho_ice           ! Water ice density
42c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
43c  varian            ! Characteristic variance of log-normal distribution
44c-----------------------------------------------------------------------
45
46      moderntracdef=.false. ! For modern traceur.def (default false, old type)
47
48      open(407, form = 'formatted', status = 'old',
49     $          file = 'traceur.def', iostat=ierr)
50      if (ierr /=0) then
51        call abort_physic('initracer',
52     $  'Problem in opening traceur.def',1)
53      end if
54!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
55       READ(407,'(A)') tracline
56       IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
57          READ(tracline,*) nqtot ! Try standard traceur.def
58       ELSE
59         moderntracdef = .true.
60         DO
61           READ(407,'(A)',iostat=ierr) tracline
62           IF (ierr==0) THEN
63             IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
64               READ(tracline,*) nqtot
65               ! Temporary not implemented solution
66               if (nqtot/=nq) then
67                 call abort_physic('initracer','Different number of&
68     & tracers in dynamics and physics not managed yet',1)
69               endif
70               EXIT
71             ENDIF
72           ELSE ! If pb, or if reached EOF without having found nqtot
73             call abort_physic('initracer','Unable to read numbers&
74     & of tracers in traceur.def',1)
75           ENDIF
76         ENDDO
77       ENDIF ! if modern or standard traceur.def
78!! -----------------------------------------------------------------------
79       !! For the moment number of tracers in dynamics and physics are equal
80       nqtot=nq
81       !! we allocate once for all arrays in common in tracer_h.F90
82       !! (supposedly those are not used before call to initracer)
83       IF (.NOT.ALLOCATED(noms))         ALLOCATE(noms(nq))
84       IF (.NOT.ALLOCATED(mmol))         ALLOCATE(mmol(nq))
85       IF (.NOT.ALLOCATED(radius))       ALLOCATE(radius(nq))
86       IF (.NOT.ALLOCATED(rho_q))        ALLOCATE(rho_q(nq))
87       IF (.NOT.ALLOCATED(qext))         ALLOCATE(qext(nq))
88       IF (.NOT.ALLOCATED(alpha_lift))   ALLOCATE(alpha_lift(nq))
89       IF (.NOT.ALLOCATED(alpha_devil))  ALLOCATE(alpha_devil(nq))
90       IF (.NOT.ALLOCATED(qextrhor))     ALLOCATE(qextrhor(nq))
91       IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
92       !! initialization
93       alpha_lift(:)=0.
94       alpha_devil(:)=0.
95       
96       ! Added by JVO 2017 : these arrays are handled later
97       ! -> initialization is the least we can do, please !!!
98       radius(:)=0.
99       qext(:)=0.
100
101
102! Initialization: copy tracer names from dynamics
103        do iq=1,nq
104          noms(iq)=nametrac(iq)
105          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
106        enddo
107
108
109! Identify tracers by their names: (and set corresponding values of mmol)
110      ! 0. initialize tracer indexes to zero:
111      ! NB: igcm_* indexes are commons in 'tracer.h'
112      do iq=1,nq
113        igcm_dustbin(iq)=0
114      enddo
115      igcm_dust_mass=0
116      igcm_dust_number=0
117      igcm_h2o_vap=0
118      igcm_h2o_ice=0
119      igcm_co2=0
120      igcm_co=0
121      igcm_o=0
122      igcm_o1d=0
123      igcm_o2=0
124      igcm_o3=0
125      igcm_h=0
126      igcm_h2=0
127      igcm_oh=0
128      igcm_ho2=0
129      igcm_h2o2=0
130      igcm_n2=0
131      igcm_n=0
132      igcm_n2d=0
133      igcm_no=0
134      igcm_no2=0
135      igcm_ar=0
136      igcm_ar_n2=0
137      igcm_co2_ice=0
138
139      igcm_ch4=0
140      igcm_ch3=0
141      igcm_ch=0
142      igcm_3ch2=0
143      igcm_1ch2=0
144      igcm_cho=0
145      igcm_ch2o=0
146      igcm_ch3o=0
147      igcm_c=0
148      igcm_c2=0
149      igcm_c2h=0
150      igcm_c2h2=0
151      igcm_c2h3=0
152      igcm_c2h4=0
153      igcm_c2h6=0
154      igcm_ch2co=0
155      igcm_ch3co=0
156      igcm_hcaer=0
157
158
159
160      ! 1. find dust tracers
161      count=0
162
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      ! Get data of tracers
377      do iq=1,nqtot
378        read(407,'(A)') tracline
379        call get_tracdat(iq, tracline)
380      enddo
381
382      close(407)
383
384c------------------------------------------------------------
385c     Initialisation tracers ....
386c------------------------------------------------------------
387      rho_q(1:nq)=0
388
389      rho_dust=2500.  ! Mars dust density (kg.m-3)
390      rho_ice=920.    ! Water ice density (kg.m-3)
391      rho_co2=1620.   ! CO2 ice density (kg.m-3)
392
393c     Initialization for water vapor
394c     ------------------------------
395      if(water) then
396         radius(igcm_h2o_vap)=0.
397         Qext(igcm_h2o_vap)=0.
398         alpha_lift(igcm_h2o_vap) =0.
399         alpha_devil(igcm_h2o_vap)=0.
400         qextrhor(igcm_h2o_vap)= 0.
401
402         !! this is defined in surfdat_h.F90
403         IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
404         IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
405
406         do ig=1,ngrid
407           if (ngrid.ne.1) watercaptag(ig)=.false.
408           dryness(ig) = 1.
409         enddo
410
411
412           radius(igcm_h2o_ice)=3.e-6
413           rho_q(igcm_h2o_ice)=rho_ice
414           Qext(igcm_h2o_ice)=0.
415!           alpha_lift(igcm_h2o_ice) =0.
416!           alpha_devil(igcm_h2o_ice)=0.
417           qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice)
418     $       / (rho_ice*radius(igcm_h2o_ice))
419
420
421      end if  ! (water)
422
423
424!
425!     some extra (possibly redundant) sanity checks for tracers:
426!     ---------------------------------------------------------
427       if (water) then
428       ! verify that we indeed have h2o_vap and h2o_ice tracers
429         if (igcm_h2o_vap.eq.0) then
430           write(*,*) "initracer: error !!"
431           write(*,*) "  cannot use water option without ",
432     &                "an h2o_vap tracer !"
433           stop
434         endif
435         if (igcm_h2o_ice.eq.0) then
436           write(*,*) "initracer: error !!"
437           write(*,*) "  cannot use water option without ",
438     &                "an h2o_ice tracer !"
439           stop
440         endif
441       endif
442
443
444c     Output for records:
445c     ~~~~~~~~~~~~~~~~~~
446      write(*,*)
447      Write(*,*) '******** initracer : dust transport parameters :'
448      write(*,*) 'alpha_lift = ', alpha_lift
449      write(*,*) 'alpha_devil = ', alpha_devil
450      write(*,*) 'radius  = ', radius
451      write(*,*) 'Qext  = ', qext
452      write(*,*)
453
454      contains
455
456      subroutine get_tracdat(iq,tracline)
457        !-------------------ADDING NEW OPTIONS-------------------
458        ! Duplicate if sentence for adding new options
459        ! Do not forget to add the new variables in tracer_h.F90
460        ! Do not forget to allocate and initialize the new variables above
461        ! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
462        !-------------------------------------------------------
463          use tracer_h
464          implicit none
465          integer,           intent(in) :: iq       ! tracer index
466          character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
467 
468          read(tracline,*) noms(iq)
469          ! JVO 20 : We should add a sanity check aborting when duplicates in names !
470          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
471          if (index(tracline,'mmol='   ) /= 0) then
472              read(tracline(index(tracline,'mmol=')+len('mmol='):),*)
473     $             mmol(iq)
474              write(*,*) ' Parameter value (traceur.def) : mmol=',
475     $             mmol(iq)
476          else
477              write(*,*) ' Parameter value (default)     : mmol=',
478     $             mmol(iq)
479          end if
480      end subroutine get_tracdat
481
482      end
483
Note: See TracBrowser for help on using the repository browser.