source: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 @ 3936

Last change on this file since 3936 was 3936, checked in by tbertrand, 2 months ago

PLUTO PCM : correcting a bug in hazecloud (wrong lyman alpha fluxes due to mu0 being negative during nighttime) + cleaning routines
TB

File size: 19.2 KB
Line 
1      SUBROUTINE initracer(ngrid,nq)
2
3      use surfdat_h, ONLY: dryness
4      USE tracer_h
5      USE callkeys_mod, only: optichaze,nb_monomer,haze,fractal,fasthaze,rad_haze,callmufi
6      USE recombin_corrk_mod, ONLY: ini_recombin
7      USE mod_phys_lmdz_para, only: is_master, bcast
8      use aerosol_mod, only: iaero_haze,i_haze
9      IMPLICIT NONE
10!=======================================================================
11!   subject:
12!   --------
13!   Initialization related to tracer
14!   (transported dust, water, chemical species, ice...)
15!
16!   Name of the tracer
17!
18!   Test of dimension :
19!   Initialize COMMON tracer in tracer.h, using tracer names provided
20!   by the argument nametrac
21!
22!   author: F.Forget
23!   ------
24!            Ehouarn Millour (oct. 2008): identify tracers by their names
25!            Y. Jaziri & J. Vatant d'Ollone (2020) : modern traceur.def
26!            B. de Batz de Trenquelléon (2024): specific microphysical tracers   
27!=======================================================================
28
29      integer,intent(in) :: ngrid,nq
30
31      character(len=500) :: tracline   ! to read traceur.def lines
32      integer :: blank      !to store the index of 1st blank when reading tracers names
33      integer iq,ig,count,ierr
34      real r0_lift , reff_lift, rho_haze
35      integer nqhaze(nq)               ! to store haze tracers
36      integer i, ia, block, j
37      character(len=20) :: txt ! to store some text
38      character(LEN=20) :: tracername ! to temporarily store text
39      character(LEN=20) :: str
40
41!-----------------------------------------------------------------------
42!  radius(nq)      ! aerosol particle radius (m)
43!  rho_q(nq)       ! tracer densities (kg.m-3)
44!  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
45!  rho_dust          ! Mars dust density
46!  rho_ice           ! Water ice density
47!  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
48!  varian            ! Characteristic variance of log-normal distribution
49!-----------------------------------------------------------------------
50
51      if (is_master) then ! only the master proc/thread needs do this
52
53        moderntracdef=.false. ! For modern traceur.def (default false, old type)
54
55        open(407, form = 'formatted', status = 'old', &
56             file = 'traceur.def', iostat=ierr)
57        if (ierr /=0) then
58          ! call abort_physic('initracer',&
59          ! 'Problem in opening traceur.def',1)
60          print*,'Problem in opening traceur.def'
61          stop
62        end if
63!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
64        READ(407,'(A)') tracline
65        IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
66          READ(tracline,*) nqtot ! Try standard traceur.def
67        ELSE
68         moderntracdef = .true.
69         DO
70           READ(407,'(A)',iostat=ierr) tracline
71           IF (ierr==0) THEN
72             IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
73               READ(tracline,*) nqtot
74               ! Temporary not implemented solution
75               if (nqtot/=nq) then
76      !            call abort_physic('initracer','Different number of '// &
77      ! 'tracers in dynamics and physics not managed yet',1)
78                print*,'!= nbr oftracers in dynamics and physics not managed yet'
79                stop
80               endif
81               EXIT
82             ENDIF
83           ELSE ! If pb, or if reached EOF without having found nqtot
84      !        call abort_physic('initracer','Unable to read numbers '// &
85      !  'of tracers in traceur.def',1)
86            print*,"unable to read numbers of tracer in tracer.def"
87            stop
88           ENDIF
89         ENDDO
90        ENDIF ! if modern or standard traceur.def
91
92       endif ! of if (is_master)
93
94       ! share the information with other procs/threads (if any)
95       CALL bcast(nqtot)
96       CALL bcast(moderntracdef)
97
98!! -----------------------------------------------------------------------
99       !! For the moment number of tracers in dynamics and physics are equal
100       nqtot=nq
101       !! we allocate once for all arrays in common in tracer_h.F90
102       !! (supposedly those are not used before call to initracer)
103       IF (.NOT.ALLOCATED(noms))           ALLOCATE(noms(nq))
104       IF (.NOT.ALLOCATED(mmol))           ALLOCATE(mmol(nq))
105       IF (.NOT.ALLOCATED(aki))            ALLOCATE(aki(nqtot))
106       IF (.NOT.ALLOCATED(cpi))            ALLOCATE(cpi(nqtot))
107       IF (.NOT.ALLOCATED(radius))         ALLOCATE(radius(nq))
108       IF (.NOT.ALLOCATED(rho_q))          ALLOCATE(rho_q(nq))
109       IF (.NOT.ALLOCATED(qext))           ALLOCATE(qext(nq))
110       IF (.NOT.ALLOCATED(qextrhor))       ALLOCATE(qextrhor(nq))
111      !  IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
112       IF (.NOT.ALLOCATED(is_rad))         ALLOCATE(is_rad(nqtot))
113       IF (.NOT.ALLOCATED(is_recomb))      ALLOCATE(is_recomb(nqtot))
114       IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
115         ALLOCATE(is_recomb_qset(nqtot))
116       ENDIF
117       IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
118         ALLOCATE(is_recomb_qotf(nqtot))
119       ENDIF
120
121       !! initialization
122       mmol(:)           = 0.
123       aki(:)            = 0.
124       cpi(:)            = 0.
125       is_rad(:)         = 0
126       is_recomb(:)      = 0
127       is_recomb_qset(:) = 0
128       is_recomb_qotf(:) = 0
129
130       ! Added by JVO 2017 : these arrays are handled later
131       ! -> initialization is the least we can do, please !!!
132       radius(:)=0.
133       qext(:)=0.
134
135       rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat
136
137
138! Initialization: Read tracers names from traceur.def
139        do iq=1,nq
140          if (is_master) read(407,'(A)') tracline
141          call bcast(tracline)
142          blank = index(tracline,' ') ! Find position of 1st blank in tracline
143          noms(iq) = tracline(1:blank) !ensure that in Modern-trac case, noms is actually the name and not all properties
144        enddo
145
146! Identify tracers by their names: (and set corresponding values of mmol)
147      ! 0. initialize tracer indexes to zero:
148      ! NB: igcm_* indexes are commons in 'tracer.h'
149      igcm_n2=0
150      igcm_ch4_gas=0
151      igcm_ch4_ice=0
152      igcm_prec_haze=0
153      igcm_co_gas=0
154      igcm_co_ice=0
155
156      nqhaze(:)=0
157      i=0
158      DO iq=1,nq
159         txt=noms(iq)
160         IF (txt(1:4).eq."haze") THEN
161            i=i+1
162            nqhaze(i)=iq
163         ENDIF
164      ENDDO
165      if ((haze.or.fasthaze).and.i==0) then
166         print*, 'Haze active but no haze tracer in traceur.def'
167         stop
168      endif
169
170      igcm_haze=0
171      igcm_haze10=0
172      igcm_haze30=0
173      igcm_haze50=0
174      igcm_haze100=0
175
176!     Eddy diffusion tracers
177      igcm_eddy1e6=0
178      igcm_eddy1e7=0
179      igcm_eddy5e7=0
180      igcm_eddy1e8=0
181      igcm_eddy5e8=0
182      write(*,*) 'initracer: noms() ', noms
183      rho_n2=1000        ! n2 ice
184      rho_ch4_ice=520.       ! rho ch4 ice
185      rho_co_ice=520.       ! rho ch4 ice
186      rho_haze=800.     ! haze
187
188      ! 1. find dust tracers
189      count=0
190
191      ! 2. find chemistry and water tracers
192      do iq=1,nq
193        if (noms(iq).eq."n2") then
194          igcm_n2=iq
195          mmol(igcm_n2)=28.
196          count=count+1
197          write(*,*) 'Tracer ',count,' = n2'
198        endif
199        if (noms(iq).eq."ch4_gas") then
200          igcm_ch4_gas=iq
201          mmol(igcm_ch4_gas)=16.
202          count=count+1
203          write(*,*) 'Tracer ',count,' = ch4 gas'
204        endif
205        if (noms(iq).eq."ch4_ice") then
206          igcm_ch4_ice=iq
207          mmol(igcm_ch4_ice)=16.
208          radius(iq)=3.e-6
209          rho_q(iq)=rho_ch4_ice
210          count=count+1
211          write(*,*) 'Tracer ',count,' = ch4 ice'
212        endif
213        if (noms(iq).eq."co_gas") then
214          igcm_co_gas=iq
215          mmol(igcm_co_gas)=28.
216          count=count+1
217          write(*,*) 'Tracer ',count,' = co gas'
218        endif
219        if (noms(iq).eq."co_ice") then
220          igcm_co_ice=iq
221          mmol(igcm_co_ice)=28.
222          radius(iq)=3.e-6
223          rho_q(iq)=rho_co_ice
224          count=count+1
225          write(*,*) 'Tracer ',count,' = co ice'
226        endif
227        if (noms(iq).eq."prec_haze") then
228          igcm_prec_haze=iq
229          count=count+1
230          write(*,*) 'Tracer ',count,' = prec haze'
231        endif
232        if (noms(iq).eq."haze") then
233          igcm_haze=iq
234          count=count+1
235          radius(iq)=rad_haze
236          rho_q(iq)=rho_haze
237          write(*,*) 'Tracer ',count,' = haze'
238        endif
239        if (noms(iq).eq."haze10") then
240          igcm_haze10=iq
241          count=count+1
242          radius(iq)=10.e-9
243          rho_q(iq)=rho_haze
244          write(*,*) 'Tracer ',count,' = haze10'
245        endif
246        if (noms(iq).eq."haze30") then
247          igcm_haze30=iq
248          count=count+1
249          radius(iq)=30.e-9
250          rho_q(iq)=rho_haze
251          write(*,*) 'Tracer ',count,' = haze30'
252        endif
253        if (noms(iq).eq."haze50") then
254          igcm_haze50=iq
255          count=count+1
256          radius(iq)=50.e-9
257          rho_q(iq)=rho_haze
258          write(*,*) 'Tracer ',count,' = haze50'
259        endif
260        if (noms(iq).eq."haze100") then
261          igcm_haze100=iq
262          count=count+1
263          radius(iq)=100.e-9
264          rho_q(iq)=rho_haze
265          write(*,*) 'Tracer ',count,' = haze100'
266        endif
267!       Eddy diffusion tracers
268        if (noms(iq).eq."eddy1e6") then
269          igcm_eddy1e6=iq
270          count=count+1
271          write(*,*) 'Tracer ',count,' = eddy1e6'
272        endif
273        if (noms(iq).eq."eddy1e7") then
274          igcm_eddy1e7=iq
275          count=count+1
276          write(*,*) 'Tracer ',count,' = eddy1e7'
277        endif
278        if (noms(iq).eq."eddy5e7") then
279          igcm_eddy5e7=iq
280          count=count+1
281          write(*,*) 'Tracer ',count,' = eddy5e7'
282        endif
283        if (noms(iq).eq."eddy1e8") then
284          igcm_eddy1e8=iq
285          count=count+1
286          write(*,*) 'Tracer ',count,' = eddy1e8'
287        endif
288        if (noms(iq).eq."eddy5e8") then
289          igcm_eddy5e8=iq
290          count=count+1
291          write(*,*) 'Tracer ',count,' = eddy5e8'
292        endif
293      enddo ! of do iq=1,nq
294
295      ! ! 3. find condensable traceurs different from h2o and n2
296      ! do iq=1,nq
297      !   if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
298      !     count=count+1
299      !   endif
300      !   if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
301      !     count=count+1
302      !   endif
303
304      ! enddo ! of do iq=1,nq
305
306      ! check that we identified all tracers:
307      if (count.ne.nq) then
308        write(*,*) "initracer: found only ",count," tracers"
309        write(*,*) "               expected ",nq
310        do iq=1,count
311          write(*,*)'      ',iq,' ',trim(noms(iq))
312        enddo
313!        stop
314      else
315        write(*,*) "initracer: found all expected tracers, namely:"
316        do iq=1,nq
317          write(*,*)'      ',iq,' ',trim(noms(iq))
318        enddo
319      endif
320     
321      ! Compute number of microphysics tracers:
322      ! By convention they all have the prefix "mu_" (case sensitive !)
323      nmicro = 0
324      IF (callmufi) THEN
325         DO iq=1,nq
326            str = noms(iq)
327            IF (str(1:3) == "mu_") THEN
328               nmicro = nmicro+1
329               count = count+1
330            ENDIF
331         ENDDO
332
333         ! Checking the expected number of tracers:
334         ! Microphysics moment model: nmicro = 4
335         IF (nmicro < 4) THEN
336            WRITE(*,*) "initracer:error:"," Inconsistent number of microphysical tracers"
337            WRITE(*,*) "expected at least 4 tracers,", nmicro, " given"
338            CALL abort
339         ELSE IF (nmicro > 4) THEN
340            WRITE(*,*) "!!! WARNING !!! initracer: I was expecting only four tracers, you gave me more."
341            CALL abort
342         ENDIF
343
344         ! microphysics indexes share the same values than original tracname.
345         IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
346         j = 1
347         DO i=1,nq
348            str = noms(i)
349            IF (str(1:3) == "mu_") THEN
350               micro_indx(j) = i
351               j=j+1
352            ENDIF
353         ENDDO
354     
355      ELSE
356         IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
357     
358      ENDIF ! end of callmufi
359
360      ! Get data of tracers. Need to rewind traceur.def first
361      if (is_master) then
362       rewind(407)
363       do
364        read(407,'(A)') tracline
365        if (index(tracline,"#") .ne. 1) then
366          exit
367        endif
368       enddo
369      endif
370      do iq=1,nqtot
371        if (is_master) read(407,'(A)') tracline
372        call bcast(tracline)
373        call get_tracdat(iq, tracline)
374      enddo
375
376      if (is_master) close(407)
377
378      ! Calculate number of microphysical tracer
379      write(*,*) 'Number of microphysical tracer nmicro = ',nmicro
380      IF (callmufi) THEN
381         call dumptracers(micro_indx)
382      ENDIF
383
384!     Processing modern traceur options
385      if(moderntracdef) then
386        call ini_recombin
387      endif
388
389!------------------------------------------------------------
390!     Initialisation tracers ....
391!------------------------------------------------------------
392      ! rho_q(1:nq)=0
393
394      rho_n2=1000.   ! N2 ice density (kg.m-3)
395
396      lw_co=274000.
397      lw_ch4=586700.
398      lw_n2=2.5e5
399      write(*,*) "lw_n2 = ", lw_n2
400
401      if (callmufi) then
402        if (optichaze) then
403          iaero_haze = 2
404          write(*,*) 'Microphysical moment model'
405          write(*,*) '--- number of haze aerosol = ', iaero_haze
406        endif ! end optichaze
407     
408      elseif (haze) then
409        ! the sedimentation radius remains radius(igcm_haze)
410        if (fractal) then
411          nmono=nb_monomer
412        else
413          nmono=1
414        endif ! end fractal
415       
416        ia = 0
417        if (optichaze) then
418          ia = ia + 1
419          iaero_haze = ia
420          write(*,*) '--- number of haze aerosol = ', iaero_haze
421          block = 0 ! Only one type of haze is active : the first one set in traceur.def
422          do iq=1,nq
423            tracername=noms(iq)
424            write(*,*) "--> tracername ",iq,'/',nq,' = ',tracername
425            if (tracername(1:4).eq."haze".and.block.eq.0) then
426              i_haze=iq
427              block=1
428              write(*,*) "i_haze=",i_haze
429              write(*,*) "Careful: if you set many haze traceurs in &
430              traceur.def,only ",tracername," will be radiatively active &
431              (first one in traceur.def)"
432            endif
433            enddo
434        endif ! end optichaze
435      endif ! end callmufi or haze
436
437!     Output for records:
438!     ~~~~~~~~~~~~~~~~~~
439      write(*,*)
440      Write(*,*) '******** initracer : dust transport parameters :'
441      write(*,*) 'radius  = ', radius
442      write(*,*) 'Qext  = ', qext
443      write(*,*)
444
445      contains
446
447      subroutine get_tracdat(iq,tracline)
448        !-------------------ADDING NEW OPTIONS-------------------
449        ! Duplicate if sentence for adding new options
450        ! Do not forget to add the new variables in tracer_h.F90
451        ! Do not forget to allocate and initialize the new variables above
452        ! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
453        !-------------------------------------------------------
454          use tracer_h
455          implicit none
456          integer,           intent(in) :: iq       ! tracer index
457          character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
458
459          read(tracline,*) noms(iq)
460          ! JVO 20 : We should add a sanity check aborting when duplicates in names !
461          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
462          ! option mmol
463          if (index(tracline,'mmol='   ) /= 0) then
464              read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
465                  mmol(iq)
466              write(*,*) ' Parameter value (traceur.def) : mmol=', &
467                  mmol(iq)
468          else
469              write(*,*) ' Parameter value (default)     : mmol=', &
470                 mmol(iq)
471          end if
472          ! option aki
473          if (index(tracline,'aki='   ) /= 0) then
474              read(tracline(index(tracline,'aki=')+len('aki='):),*) &
475                  aki(iq)
476              write(*,*) ' Parameter value (traceur.def) : aki=', &
477                  aki(iq)
478          else
479              write(*,*) ' Parameter value (default)     : aki=', &
480                  aki(iq)
481          end if
482          ! option cpi
483          if (index(tracline,'cpi='   ) /= 0) then
484              read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
485                  cpi(iq)
486              write(*,*) ' Parameter value (traceur.def) : cpi=', &
487                  cpi(iq)
488          else
489              write(*,*) ' Parameter value (default)     : cpi=', &
490                  cpi(iq)
491          end if
492          ! option is_rad
493          if (index(tracline,'is_rad=') /= 0) then
494              read(tracline(index(tracline,'is_rad=') &
495              +len('is_rad='):),*) is_rad(iq)
496              write(*,*) ' Parameter value (traceur.def) : is_rad=', &
497              is_rad(iq)
498          else
499              write(*,*) ' Parameter value (default)     : is_rad=',  &
500              is_rad(iq)
501          end if
502          ! option is_recomb
503          if (index(tracline,'is_recomb=') /= 0) then
504              read(tracline(index(tracline,'is_recomb=') &
505              +len('is_recomb='):),*) is_recomb(iq)
506              write(*,*) ' Parameter value (traceur.def) : is_recomb=', &
507              is_recomb(iq)
508          else
509              write(*,*) ' Parameter value (default)     : is_recomb=', &
510              is_recomb(iq)
511          end if
512          ! option is_recomb_qset
513          if (index(tracline,'is_recomb_qset=') /= 0) then
514              read(tracline(index(tracline,'is_recomb_qset=') &
515              +len('is_recomb_qset='):),*) is_recomb_qset(iq)
516              write(*,*) ' Parameter value (traceur.def) :'// &
517              ' is_recomb_qset=', &
518              is_recomb_qset(iq)
519          else
520              write(*,*) ' Parameter value (default)     :'// &
521              ' is_recomb_qset=', &
522              is_recomb_qset(iq)
523          endif
524          ! option is_recomb_qotf
525          if (index(tracline,'is_recomb_qotf=') /= 0) then
526              read(tracline(index(tracline,'is_recomb_qotf=') &
527              +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
528              write(*,*) ' Parameter value (traceur.def) :'// &
529              ' is_recomb_qotf=', &
530              is_recomb_qotf(iq)
531          else
532              write(*,*) ' Parameter value (default)     :'// &
533              ' is_recomb_qotf=',  &
534              is_recomb_qotf(iq)
535          end if
536          !option radius
537          if (index(tracline,'radius=') .ne. 0) then
538            read(tracline(index(tracline,'radius=') &
539              +len('radius='):),*) radius(iq)
540            write(*,*)'Parameter value (traceur.def) :'// &
541              "radius=",radius(iq), "m"
542          else
543            write(*,*) ' Parameter value (default)     :'// &
544            ' radius=', radius(iq)," m"
545          endif
546          !option rho
547          if (index(tracline,'rho=') .ne. 0) then
548            read(tracline(index(tracline,'rho=') &
549              +len('rho='):),*) rho_q(iq)
550            write(*,*)'Parameter value (traceur.def) :'// &
551              "rho=",rho_q(iq)
552          else
553            write(*,*) ' Parameter value (default)     :'// &
554              ' rho=', rho_q(iq)
555          endif
556      end subroutine get_tracdat
557
558      end
559
Note: See TracBrowser for help on using the repository browser.