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

Last change on this file since 2618 was 2550, checked in by emillour, 4 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
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
[2550]7      USE mod_phys_lmdz_para, only: is_master, bcast
[135]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
[787]19c   by the argument nametrac
[135]20c
21c   author: F.Forget
22c   ------
23c            Ehouarn Millour (oct. 2008) identify tracers by their names
[2436]24c            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
[135]25c=======================================================================
26
[1980]27      integer,intent(in) :: ngrid,nq
28      character(len=30),intent(in) :: nametrac(nq) ! name of the tracer from dynamics
[787]29
[2436]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
[135]33
34
[787]35
[135]36c-----------------------------------------------------------------------
[787]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
[135]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
[2550]48      if (is_master) then ! only the master proc/thread needs do this
[2436]49
[2550]50        moderntracdef=.false. ! For modern traceur.def (default false, old type)
51
52        open(407, form = 'formatted', status = 'old',
[2436]53     $          file = 'traceur.def', iostat=ierr)
[2550]54        if (ierr /=0) then
55          call abort_physic('initracer',
56     $    'Problem in opening traceur.def',1)
57        end if
[2436]58!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
[2550]59        READ(407,'(A)') tracline
60        IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
[2436]61          READ(tracline,*) nqtot ! Try standard traceur.def
[2550]62        ELSE
[2436]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
[2542]71                 call abort_physic('initracer','Different number of '//
72     & 'tracers in dynamics and physics not managed yet',1)
[2436]73               endif
74               EXIT
75             ENDIF
76           ELSE ! If pb, or if reached EOF without having found nqtot
[2542]77             call abort_physic('initracer','Unable to read numbers '//
78     & 'of tracers in traceur.def',1)
[2436]79           ENDIF
80         ENDDO
[2550]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       
[2436]89!! -----------------------------------------------------------------------
90       !! For the moment number of tracers in dynamics and physics are equal
[1621]91       nqtot=nq
[787]92       !! we allocate once for all arrays in common in tracer_h.F90
93       !! (supposedly those are not used before call to initracer)
[2543]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
[787]114       !! initialization
[2543]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
[1764]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.
[787]130
[135]131
[1980]132! Initialization: copy tracer names from dynamics
[787]133        do iq=1,nq
134          noms(iq)=nametrac(iq)
[1980]135          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
[135]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'
[787]142      do iq=1,nq
[135]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
[1801]161      igcm_n=0
162      igcm_n2d=0
163      igcm_no=0
164      igcm_no2=0
[135]165      igcm_ar=0
166      igcm_ar_n2=0
167      igcm_co2_ice=0
168
[1801]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
[135]189
190      ! 1. find dust tracers
191      count=0
192
193      ! 2. find chemistry and water tracers
[787]194      do iq=1,nq
[135]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
[1801]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
[787]389      enddo ! of do iq=1,nq
[135]390     
391      ! check that we identified all tracers:
[787]392      if (count.ne.nq) then
[135]393        write(*,*) "initracer: found only ",count," tracers"
[787]394        write(*,*) "               expected ",nq
[135]395        do iq=1,count
396          write(*,*)'      ',iq,' ',trim(noms(iq))
397        enddo
[1297]398!        stop
[135]399      else
400        write(*,*) "initracer: found all expected tracers, namely:"
[787]401        do iq=1,nq
[135]402          write(*,*)'      ',iq,' ',trim(noms(iq))
403        enddo
404      endif
405
[2436]406      ! Get data of tracers
407      do iq=1,nqtot
[2550]408        if (is_master) read(407,'(A)') tracline
409        call bcast(tracline)
[2436]410        call get_tracdat(iq, tracline)
411      enddo
[135]412
[2550]413      if (is_master) close(407)
[2436]414
[2542]415      ! Calculate number of species in the chemistry
416      nesp = sum(is_chim)
417      write(*,*) 'Number of species in the chemistry nesp = ',nesp
418
[2543]419!     Processing modern traceur options
420      if(moderntracdef) then
421        call ini_recombin
422      endif
423
[135]424c------------------------------------------------------------
425c     Initialisation tracers ....
426c------------------------------------------------------------
[2278]427      rho_q(1:nq)=0
[135]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
[787]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.
[135]448           dryness(ig) = 1.
449         enddo
450
[253]451
[135]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
[1801]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
[135]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
[2436]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))
[2542]511          ! option mmol
[2436]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
[2542]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
[2543]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
[2436]595      end subroutine get_tracdat
596
[135]597      end
[2436]598
Note: See TracBrowser for help on using the repository browser.