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

Last change on this file since 2613 was 2550, checked in by emillour, 3 years ago

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