SUBROUTINE initracer(ngrid,nq) use surfdat_h, ONLY: dryness USE tracer_h USE callkeys_mod, only: aerohaze,nb_monomer,haze,fractal,fasthaze,rad_haze USE recombin_corrk_mod, ONLY: ini_recombin USE mod_phys_lmdz_para, only: is_master, bcast use generic_cloud_common_h use aerosol_mod, only: iaero_haze,i_haze IMPLICIT NONE !======================================================================= ! subject: ! -------- ! Initialization related to tracer ! (transported dust, water, chemical species, ice...) ! ! Name of the tracer ! ! Test of dimension : ! Initialize COMMON tracer in tracer.h, using tracer names provided ! by the argument nametrac ! ! author: F.Forget ! ------ ! Ehouarn Millour (oct. 2008) identify tracers by their names ! Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def ! L Teinturier (2022): Tracer names are now read here instead of ! inside interfaces !======================================================================= integer,intent(in) :: ngrid,nq character(len=500) :: tracline ! to read traceur.def lines integer :: blank !to store the index of 1st blank when reading tracers names integer iq,ig,count,ierr real r0_lift , reff_lift, rho_haze integer nqhaze(nq) ! to store haze tracers integer i, ia, block character(len=20) :: txt ! to store some text CHARACTER(LEN=20) :: tracername ! to temporarily store text !----------------------------------------------------------------------- ! radius(nq) ! aerosol particle radius (m) ! rho_q(nq) ! tracer densities (kg.m-3) ! qext(nq) ! Single Scat. Extinction coeff at 0.67 um ! alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1) ! alpha_devil(nq) ! lifting coeeficient by dust devil ! rho_dust ! Mars dust density ! rho_ice ! Water ice density ! doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio ! varian ! Characteristic variance of log-normal distribution !----------------------------------------------------------------------- if (is_master) then ! only the master proc/thread needs do this moderntracdef=.false. ! For modern traceur.def (default false, old type) open(407, form = 'formatted', status = 'old', & file = 'traceur.def', iostat=ierr) if (ierr /=0) then ! call abort_physic('initracer',& ! 'Problem in opening traceur.def',1) print*,'Problem in opening traceur.def' stop end if !! - Modif. by JVO and YJ for modern planetary traceur.def --------------- READ(407,'(A)') tracline IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def READ(tracline,*) nqtot ! Try standard traceur.def ELSE moderntracdef = .true. DO READ(407,'(A)',iostat=ierr) tracline IF (ierr==0) THEN IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header READ(tracline,*) nqtot ! Temporary not implemented solution if (nqtot/=nq) then ! call abort_physic('initracer','Different number of '// & ! 'tracers in dynamics and physics not managed yet',1) print*,'!= nbr oftracers in dynamics and physics not managed yet' stop endif EXIT ENDIF ELSE ! If pb, or if reached EOF without having found nqtot ! call abort_physic('initracer','Unable to read numbers '// & ! 'of tracers in traceur.def',1) print*,"unable to read numbers of tracer in tracer.def" stop ENDIF ENDDO ENDIF ! if modern or standard traceur.def endif ! of if (is_master) ! share the information with other procs/threads (if any) CALL bcast(nqtot) CALL bcast(moderntracdef) !! ----------------------------------------------------------------------- !! For the moment number of tracers in dynamics and physics are equal nqtot=nq !! we allocate once for all arrays in common in tracer_h.F90 !! (supposedly those are not used before call to initracer) IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) IF (.NOT.ALLOCATED(mmol)) ALLOCATE(mmol(nq)) IF (.NOT.ALLOCATED(aki)) ALLOCATE(aki(nqtot)) IF (.NOT.ALLOCATED(cpi)) ALLOCATE(cpi(nqtot)) IF (.NOT.ALLOCATED(radius)) ALLOCATE(radius(nq)) IF (.NOT.ALLOCATED(rho_q)) ALLOCATE(rho_q(nq)) IF (.NOT.ALLOCATED(qext)) ALLOCATE(qext(nq)) IF (.NOT.ALLOCATED(alpha_lift)) ALLOCATE(alpha_lift(nq)) IF (.NOT.ALLOCATED(alpha_devil)) ALLOCATE(alpha_devil(nq)) IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq)) ! IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq)) IF (.NOT.ALLOCATED(is_chim)) ALLOCATE(is_chim(nqtot)) IF (.NOT.ALLOCATED(is_rad)) ALLOCATE(is_rad(nqtot)) IF (.NOT.ALLOCATED(is_recomb)) ALLOCATE(is_recomb(nqtot)) IF (.NOT.ALLOCATED(is_recomb_qset)) THEN ALLOCATE(is_recomb_qset(nqtot)) ENDIF IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN ALLOCATE(is_recomb_qotf(nqtot)) ENDIF IF (.NOT. allocated(is_condensable)) allocate(is_condensable(nq)) !LT IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq)) IF (.NOT. allocated(constants_delta_gasH)) allocate(constants_delta_gasH(nq)) IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq)) IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq)) IF (.NOT. allocated(constants_epsi_generic)) allocate(constants_epsi_generic(nq)) IF (.NOT. allocated(constants_RLVTT_generic)) allocate(constants_RLVTT_generic(nq)) IF (.NOT. allocated(constants_metallicity_coeff)) allocate(constants_metallicity_coeff(nq)) IF (.NOT. allocated(constants_RCPV_generic)) allocate(constants_RCPV_generic(nq)) !! initialization alpha_lift(:) = 0. alpha_devil(:) = 0. mmol(:) = 0. aki(:) = 0. cpi(:) = 0. is_chim(:) = 0 is_rad(:) = 0 is_recomb(:) = 0 is_recomb_qset(:) = 0 is_recomb_qotf(:) = 0 ! Added by JVO 2017 : these arrays are handled later ! -> initialization is the least we can do, please !!! radius(:)=0. qext(:)=0. ! For condensable tracers, by Lucas Teinturier and Noé Clément (2022) is_condensable(:)= 0 is_rgcs(:) = 0 constants_mass(:)=0 constants_delta_gasH(:)=0 constants_Tref(:)=0 constants_Pref(:)=0 constants_epsi_generic(:)=0 constants_RLVTT_generic(:)=0 constants_metallicity_coeff(:)=0 constants_RCPV_generic(:)=0 rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat ! Initialization: Read tracers names from traceur.def do iq=1,nq if (is_master) read(407,'(A)') tracline call bcast(tracline) blank = index(tracline,' ') ! Find position of 1st blank in tracline noms(iq) = tracline(1:blank) !ensure that in Modern-trac case, noms is actually the name and not all properties enddo ! Identify tracers by their names: (and set corresponding values of mmol) ! 0. initialize tracer indexes to zero: ! NB: igcm_* indexes are commons in 'tracer.h' igcm_n2=0 igcm_ch4_gas=0 igcm_ch4_ice=0 igcm_prec_haze=0 igcm_co_gas=0 igcm_co_ice=0 nqhaze(:)=0 i=0 DO iq=1,nq txt=noms(iq) IF (txt(1:4).eq."haze") THEN i=i+1 nqhaze(i)=iq ENDIF ENDDO if ((haze.or.fasthaze).and.i==0) then print*, 'Haze active but no haze tracer in traceur.def' stop endif igcm_haze=0 igcm_haze10=0 igcm_haze30=0 igcm_haze50=0 igcm_haze100=0 ! Eddy diffusion tracers igcm_eddy1e6=0 igcm_eddy1e7=0 igcm_eddy5e7=0 igcm_eddy1e8=0 igcm_eddy5e8=0 write(*,*) 'initracer: noms() ', noms rho_n2=1000 ! n2 ice rho_ch4_ice=520. ! rho ch4 ice rho_co_ice=520. ! rho ch4 ice rho_haze=800. ! haze ! 1. find dust tracers count=0 ! 2. find chemistry and water tracers do iq=1,nq if (noms(iq).eq."n2") then igcm_n2=iq mmol(igcm_n2)=28. count=count+1 write(*,*) 'Tracer ',count,' = n2' endif if (noms(iq).eq."ch4_gas") then igcm_ch4_gas=iq mmol(igcm_ch4_gas)=16. count=count+1 write(*,*) 'Tracer ',count,' = ch4 gas' endif if (noms(iq).eq."ch4_ice") then igcm_ch4_ice=iq mmol(igcm_ch4_ice)=16. radius(iq)=3.e-6 rho_q(iq)=rho_ch4_ice count=count+1 write(*,*) 'Tracer ',count,' = ch4 ice' endif if (noms(iq).eq."co_gas") then igcm_co_gas=iq mmol(igcm_co_gas)=28. count=count+1 write(*,*) 'Tracer ',count,' = co gas' endif if (noms(iq).eq."co_ice") then igcm_co_ice=iq mmol(igcm_co_ice)=28. radius(iq)=3.e-6 rho_q(iq)=rho_co_ice count=count+1 write(*,*) 'Tracer ',count,' = co ice' endif if (noms(iq).eq."prec_haze") then igcm_prec_haze=iq count=count+1 write(*,*) 'Tracer ',count,' = prec haze' endif if (noms(iq).eq."haze") then igcm_haze=iq count=count+1 radius(iq)=rad_haze rho_q(iq)=rho_haze write(*,*) 'Tracer ',count,' = haze' endif if (noms(iq).eq."haze10") then igcm_haze10=iq count=count+1 radius(iq)=10.e-9 rho_q(iq)=rho_haze write(*,*) 'Tracer ',count,' = haze10' endif if (noms(iq).eq."haze30") then igcm_haze30=iq count=count+1 radius(iq)=30.e-9 rho_q(iq)=rho_haze write(*,*) 'Tracer ',count,' = haze30' endif if (noms(iq).eq."haze50") then igcm_haze50=iq count=count+1 radius(iq)=50.e-9 rho_q(iq)=rho_haze write(*,*) 'Tracer ',count,' = haze50' endif if (noms(iq).eq."haze100") then igcm_haze100=iq count=count+1 radius(iq)=100.e-9 rho_q(iq)=rho_haze write(*,*) 'Tracer ',count,' = haze100' endif ! Eddy diffusion tracers if (noms(iq).eq."eddy1e6") then igcm_eddy1e6=iq count=count+1 write(*,*) 'Tracer ',count,' = eddy1e6' endif if (noms(iq).eq."eddy1e7") then igcm_eddy1e7=iq count=count+1 write(*,*) 'Tracer ',count,' = eddy1e7' endif if (noms(iq).eq."eddy5e7") then igcm_eddy5e7=iq count=count+1 write(*,*) 'Tracer ',count,' = eddy5e7' endif if (noms(iq).eq."eddy1e8") then igcm_eddy1e8=iq count=count+1 write(*,*) 'Tracer ',count,' = eddy1e8' endif if (noms(iq).eq."eddy5e8") then igcm_eddy5e8=iq count=count+1 write(*,*) 'Tracer ',count,' = eddy5e8' endif enddo ! of do iq=1,nq ! ! 3. find condensable traceurs different from h2o and n2 ! do iq=1,nq ! if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then ! count=count+1 ! endif ! if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then ! count=count+1 ! endif ! enddo ! of do iq=1,nq ! check that we identified all tracers: if (count.ne.nq) then write(*,*) "initracer: found only ",count," tracers" write(*,*) " expected ",nq do iq=1,count write(*,*)' ',iq,' ',trim(noms(iq)) enddo ! stop else write(*,*) "initracer: found all expected tracers, namely:" do iq=1,nq write(*,*)' ',iq,' ',trim(noms(iq)) enddo endif ! Get data of tracers. Need to rewind traceur.def first if (is_master) then rewind(407) do read(407,'(A)') tracline if (index(tracline,"#") .ne. 1) then exit endif enddo endif do iq=1,nqtot if (is_master) read(407,'(A)') tracline call bcast(tracline) call get_tracdat(iq, tracline) enddo if (is_master) close(407) ! Get specific data of condensable tracers do iq=1,nq if((is_condensable(iq)==1)) then write(*,*) "There is a specie which is condensable, for generic condensation : ", noms(iq) write(*,*) 'looking specie parameters for : ',noms(iq)(1:len(trim(noms(iq)))-4) ! call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4)) ! constants_mass(iq)=m ! constants_delta_gasH(iq)=delta_gasH ! constants_Tref(iq)=Tref ! constants_Pref(iq)=Pref ! constants_epsi_generic(iq)=epsi_generic ! constants_RLVTT_generic(iq)=RLVTT_generic ! constants_metallicity_coeff(iq)=metallicity_coeff ! constants_RCPV_generic(iq)=RCPV_generic else write(*,*) "This tracer is not condensable, for generic condensation : : ", noms(iq) write(*,*) "We keep condensable constants at zero" endif !(is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0)) enddo ! iq=1,nq ! Calculate number of species in the chemistry nesp = sum(is_chim) write(*,*) 'Number of species in the chemistry nesp = ',nesp ! Calculate number of generic tracers ngt = sum(is_condensable) write(*,*) 'Number of generic tracer is ngt = ',ngt ! Calculate number of radiative generic condensable species n_rgcs = sum(is_rgcs) write(*,*)'Number of Radiative Generic Condensable Species is n_rgcs = ',n_rgcs if (n_rgcs> ngt/2) then write(*,*) 'You have more Radiative Generic Condensable Species than Generic Condensable Species' write(*,*)'This is not possible: check your Modern traceur.def' call abort_physic("initracer, issue with # of RGCS and GCS") endif ! Processing modern traceur options if(moderntracdef) then call ini_recombin endif !------------------------------------------------------------ ! Initialisation tracers .... !------------------------------------------------------------ ! rho_q(1:nq)=0 rho_n2=1000. ! N2 ice density (kg.m-3) lw_co=274000. lw_ch4=586700. lw_n2=2.5e5 write(*,*) "lw_n2 = ", lw_n2 if (haze) then ! the sedimentation radius remains radius(igcm_haze) if (fractal) then nmono=nb_monomer else nmono=1 endif ia=0 if (aerohaze) then ia=ia+1 iaero_haze=ia write(*,*) '--- number of haze aerosol = ', iaero_haze block=0 ! Only one type of haze is active : the first one set in traceur.def do iq=1,nq tracername=noms(iq) write(*,*) "--> tracername ",iq,'/',nq,' = ',tracername if (tracername(1:4).eq."haze".and.block.eq.0) then i_haze=iq block=1 write(*,*) "i_haze=",i_haze write(*,*) "Careful: if you set many haze traceurs in& traceur.def,only ",tracername," will be radiatively active& (first one in traceur.def)" endif enddo endif endif ! Initialization for water vapor !AF24: removed ! Output for records: ! ~~~~~~~~~~~~~~~~~~ write(*,*) Write(*,*) '******** initracer : dust transport parameters :' write(*,*) 'alpha_lift = ', alpha_lift write(*,*) 'alpha_devil = ', alpha_devil write(*,*) 'radius = ', radius write(*,*) 'Qext = ', qext write(*,*) contains subroutine get_tracdat(iq,tracline) !-------------------ADDING NEW OPTIONS------------------- ! Duplicate if sentence for adding new options ! Do not forget to add the new variables in tracer_h.F90 ! Do not forget to allocate and initialize the new variables above ! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern" !------------------------------------------------------- use tracer_h implicit none integer, intent(in) :: iq ! tracer index character(len=500),intent(in) :: tracline ! traceur.def lines with parameters read(tracline,*) noms(iq) ! JVO 20 : We should add a sanity check aborting when duplicates in names ! write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq)) ! option mmol if (index(tracline,'mmol=' ) /= 0) then read(tracline(index(tracline,'mmol=')+len('mmol='):),*)& mmol(iq) write(*,*) ' Parameter value (traceur.def) : mmol=', & mmol(iq) else write(*,*) ' Parameter value (default) : mmol=', & mmol(iq) end if ! option aki if (index(tracline,'aki=' ) /= 0) then read(tracline(index(tracline,'aki=')+len('aki='):),*) & aki(iq) write(*,*) ' Parameter value (traceur.def) : aki=', & aki(iq) else write(*,*) ' Parameter value (default) : aki=', & aki(iq) end if ! option cpi if (index(tracline,'cpi=' ) /= 0) then read(tracline(index(tracline,'cpi=')+len('cpi='):),*) & cpi(iq) write(*,*) ' Parameter value (traceur.def) : cpi=', & cpi(iq) else write(*,*) ' Parameter value (default) : cpi=', & cpi(iq) end if ! option is_chim if (index(tracline,'is_chim=' ) /= 0) then read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) & is_chim(iq) write(*,*) ' Parameter value (traceur.def) : is_chim=', & is_chim(iq) else write(*,*) ' Parameter value (default) : is_chim=', & is_chim(iq) end if ! option is_rad if (index(tracline,'is_rad=') /= 0) then read(tracline(index(tracline,'is_rad=') & +len('is_rad='):),*) is_rad(iq) write(*,*) ' Parameter value (traceur.def) : is_rad=', & is_rad(iq) else write(*,*) ' Parameter value (default) : is_rad=', & is_rad(iq) end if ! option is_recomb if (index(tracline,'is_recomb=') /= 0) then read(tracline(index(tracline,'is_recomb=') & +len('is_recomb='):),*) is_recomb(iq) write(*,*) ' Parameter value (traceur.def) : is_recomb=', & is_recomb(iq) else write(*,*) ' Parameter value (default) : is_recomb=', & is_recomb(iq) end if ! option is_recomb_qset if (index(tracline,'is_recomb_qset=') /= 0) then read(tracline(index(tracline,'is_recomb_qset=') & +len('is_recomb_qset='):),*) is_recomb_qset(iq) write(*,*) ' Parameter value (traceur.def) :'// & ' is_recomb_qset=', & is_recomb_qset(iq) else write(*,*) ' Parameter value (default) :'// & ' is_recomb_qset=', & is_recomb_qset(iq) endif ! option is_recomb_qotf if (index(tracline,'is_recomb_qotf=') /= 0) then read(tracline(index(tracline,'is_recomb_qotf=') & +len('is_recomb_qotf='):),*) is_recomb_qotf(iq) write(*,*) ' Parameter value (traceur.def) :'// & ' is_recomb_qotf=', & is_recomb_qotf(iq) else write(*,*) ' Parameter value (default) :'// & ' is_recomb_qotf=', & is_recomb_qotf(iq) end if !option is_condensable (LT) if (index(tracline,'is_condensable=') /=0) then read(tracline(index(tracline,'is_condensable=') & +len('is_condensable='):),*) is_condensable(iq) write(*,*) ' Parameter value (traceur.def) :'// & ' is_condensable=', is_condensable(iq) else write(*,*) ' Parameter value (default) :'// & ' is_condensable=', is_condensable(iq) endif !option radius if (index(tracline,'radius=') .ne. 0) then read(tracline(index(tracline,'radius=') & +len('radius='):),*) radius(iq) write(*,*)'Parameter value (traceur.def) :'// & "radius=",radius(iq), "m" else write(*,*) ' Parameter value (default) :'// & ' radius=', radius(iq)," m" endif !option rho if (index(tracline,'rho=') .ne. 0) then read(tracline(index(tracline,'rho=') & +len('rho='):),*) rho_q(iq) write(*,*)'Parameter value (traceur.def) :'// & "rho=",rho_q(iq) else write(*,*) ' Parameter value (default) :'// & ' rho=', rho_q(iq) endif !option is_rgcs if (index(tracline,'is_rgcs') .ne. 0) then read(tracline(index(tracline,'is_rgcs=') & +len('is_rgcs='):),*) is_rgcs(iq) write(*,*)'Parameter value (traceur.def) :'// & 'is_rgcs=',is_rgcs(iq) else write(*,*)'Parameter value (default) : '// & 'is_rgcs = ',is_rgcs(iq) endif end subroutine get_tracdat end