subroutine calchim_asis(ngrid,nlayer,nq, & ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,fract, & zzlev,zzlay,zday,pq,pdq,dqchim,dqschim, & tauref,co2ice, & pu,pdu,pv,pdv,surfdust,surfice,icount,zdtchim) use tracer_h, only: mmol, noms, nesp, is_chim use conc_mod, only: mmean ! mean molecular mass of the atmosphere USE comcstfi_mod use callkeys_mod use time_phylmdz_mod, only: ecritphy, iphysiq ! output rate set by ecritphy use types_asis, only: v_phot_3d, jlabel, nb_phot_hv_max use chimiedata_h use radcommon_h, only: glat use wstats_mod, only: wstats implicit none !======================================================================= ! ! subject: ! -------- ! ! Prepare the call for the photochemical module, and send back the ! tendencies from photochemistry in the chemical species mass mixing ratios ! ! Author: Sebastien Lebonnois (08/11/2002) ! ------- ! update 12/06/2003 for water ice clouds and compatibility with dust ! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll) ! update 03/05/2005 cosmetic changes (Franck Lefevre) ! update sept. 2008 identify tracers by their names (Ehouarn Millour) ! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre) ! update 16/03/2012 optimization (Franck Lefevre) ! update 11/12/2013 optimization (Franck Lefevre) ! update 20/10/2017 adapted to LMDZ GENERIC+cosmetic changes (Benjamin Charnay) ! update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri) ! ! Arguments: ! ---------- ! ! Input: ! ! ptimestep timestep (s) ! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa) ! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa) ! pt(ngrid,nlayer) Temperature (K) ! pdt(ngrid,nlayer) Temperature tendency (K) ! pu(ngrid,nlayer) u component of the wind (ms-1) ! pdu(ngrid,nlayer) u component tendency (K) ! pv(ngrid,nlayer) v component of the wind (ms-1) ! pdv(ngrid,nlayer) v component tendency (K) ! dist_sol distance of the sun (AU) ! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith) ! fract(ngrid) day fraction (for diurnal=false) ! pq(ngrid,nlayer,nq) Advected fields, ie chemical species here ! pdq(ngrid,nlayer,nq) Previous tendencies on pq ! tauref(ngrid) Optical depth at 7 hPa ! co2ice(ngrid) co2 ice surface layer (kg.m-2) ! surfdust(ngrid,nlayer) dust surface area (m2/m3) ! surfice(ngrid,nlayer) ice surface area (m2/m3) ! ! Output: ! ! dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry ! dqschim(ngrid,nesp) ! tendencies on qsurf ! !======================================================================= ! input: integer,intent(in) :: ngrid ! number of atmospheric columns integer,intent(in) :: nlayer ! number of atmospheric layers integer,intent(in) :: nq ! number of tracers real :: ptimestep real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers real :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries real :: pt(ngrid,nlayer) ! temperature real :: pdt(ngrid,nlayer) ! temperature tendency real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1) real :: pdu(ngrid,nlayer) ! u component tendency real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1) real :: pdv(ngrid,nlayer) ! v component tendency real :: dist_sol ! distance of the sun (AU) real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith) real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio real :: pdq(ngrid,nlayer,nq) ! previous tendencies real :: zday ! date (time since Ls=0, in martian days) real :: tauref(ngrid) ! optical depth at 7 hPa real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2) real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3) real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3) ! output: real :: dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry real :: dqschim(ngrid,nesp) ! tendencies on qsurf ! local variables: integer :: iloc(1) ! index of major species integer :: ig,l,i,iq,iqmax,iesp integer :: foundswitch, lswitch ! to switch between photochem and thermochem ? (YJ) integer,save :: chemthermod real :: zq(ngrid,nlayer,nesp) ! pq+pdq*ptimestep before chemistry ! new mole fraction after real :: zt(ngrid,nlayer) ! temperature real :: zu(ngrid,nlayer) ! u component of the wind real :: zv(ngrid,nlayer) ! v component of the wind real :: taucol ! optical depth at 7 hPa real :: xmmol(nlayer,nesp) ! mmol/mmean but only for chemical species logical,save :: firstcall = .true. ! for each column of atmosphere: real :: zpress(nlayer) ! Pressure (mbar) real :: zdens(nlayer) ! Density (cm-3) real :: ztemp(nlayer) ! Temperature (K) real :: zlocal(nlayer) ! Altitude (km) real :: zycol(nlayer,nesp) ! Composition (mole fractions) real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1) real :: szacol ! Solar zenith angle real :: fract(ngrid) ! day fraction real :: fractcol ! day fraction real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3) real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3) integer :: iter(nlayer) ! Number of chemical iterations ! within one physical timestep integer :: icount ! for output: logical :: output ! to issue calls to writediagfi and stats parameter (output = .true.) real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations ! within one physical timestep integer :: ierr real :: zdtchim(ngrid,nlayer) ! temperature modification by chemistry real :: dEzchim(ngrid,nlayer) ! energy modification by chemistry real :: zdtchim_output(ngrid) ! flux modification by chemistry in W.m-2 !======================================================================= ! initialization of the chemistry (first call only) !======================================================================= if (firstcall) then !! Moved to routine indice in photochemistry_asis !! because nb_phot_hv_max value needed in order !! to choose if we call read_phototable or not. !! A cleaner solution need to be find. ! if (photochem .and. .not. jonline) then ! print*,'calchim: Read photolysis lookup table' ! call read_phototable ! end if if (.not.allocated(SF_mode)) allocate(SF_mode(nesp)) if (.not.allocated(SF_value)) allocate(SF_value(nesp)) if (.not.allocated(prod_rate)) allocate(prod_rate(nesp)) if (.not.allocated(surface_flux)) allocate(surface_flux(ngrid,nesp)) if (.not.allocated(surface_flux2)) allocate(surface_flux2(ngrid,nesp)) if (.not.allocated(escape)) allocate(escape(ngrid,nesp)) if (.not.allocated(chemnoms)) allocate(chemnoms(nesp)) surface_flux(:,:) = 0. surface_flux2(:,:) = 0. escape(:,:) = 0. SF_mode(:) = 2 SF_value(:) = 0. prod_rate(:) = 0. iter_3d(:,:) = 0. iter(:) = 0. call ini_tracchim ! Sanity check mmol /= 0. in chemistry do iq = 1,nq if (is_chim(iq).eq.1 .and. mmol(iq).eq.0.) then write(*,*) 'Error in calchim:' write(*,*) 'Mmol cannot be equal to 0 for chemical species' stop end if end do firstcall = .false. end if ! if (firstcall) ! Initializations zycol(:,:) = 0. dqchim(:,:,:) = 0. dqschim(:,:) = 0. !======================================================================= ! loop over grid !======================================================================= do ig = 1,ngrid foundswitch = 0 do l = 1,nlayer iesp = 0 do iq = 1,nq if (is_chim(iq).eq.1) then iesp = iesp + 1 zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep xmmol(l,iesp) = mmol(iq)/mmean(ig,l) zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp) end if end do zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep zpress(l) = pplay(ig,l)/100. ztemp(l) = zt(ig,l) zdens(l) = zpress(l)/(kb*1.e4*ztemp(l)) zlocal(l) = zzlay(ig,l)/1000. zmmean(l) = mmean(ig,l) ! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3 surfdust1d(l) = surfdust(ig,l)*1.e-2 surfice1d(l) = surfice(ig,l)*1.e-2 ! search for switch index between regions ! if (photochem .and. thermochem) then ! if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then ! lswitch = l ! foundswitch = 1 ! end if ! end if if (.not. photochem) then lswitch = 22 end if ! if (.not. thermochem) then lswitch = min(50,nlayer+1) ! end if end do ! of do l=1,nlayer szacol = acos(mu0(ig))*180./pi taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars fractcol = fract(ig) !======================================================================= ! call chemical subroutines !======================================================================= ! chemistry in lower atmosphere call photochemistry_asis(nlayer,ngrid, & ig,lswitch,zycol,szacol,fractcol,ptimestep,& zpress,zlocal,ztemp,zdens,zmmean,dist_sol, & surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:)) ! diagnostic photochemical heating zdtchim_output(ig) = 0. do l = 1,nlayer dEzchim(ig,l) = zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig) zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig) end do do l = 1,nlayer iter_3d(ig,l) = iter(l) end do ! chemistry in upper atmosphere ! if (thermochem) then ! call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,& ! zdens,zpress,zlocal,szacol,ptimestep,zday) ! end if ! dry deposition if (depos) then call deposition_source(ngrid, nlayer, nq, & ig, zzlay, zzlev, zdens, zycol, ptimestep) surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:) if (ngrid==1) then if(mod(icount,ecritphy).eq.0) then surface_flux2(ig,:) = surface_flux2(ig,:)/float(ecritphy) endif else if(mod(icount*iphysiq,ecritphy).eq.0) then surface_flux2(ig,:) = surface_flux2(ig,:)*iphysiq/float(ecritphy) endif endif end if !======================================================================= ! tendencies !======================================================================= ! index of the most abundant species at each level ! major(:) = maxloc(zycol, dim = 2) ! tendency for the most abundant species = - sum of others do l = 1,nlayer iloc=maxloc(zycol(l,:)) iqmax=iloc(1) do iq = 1,nesp if (iq /= iqmax) then dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep) dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq) end if end do end do ! of do l = 1,nlayer !======================================================================= ! end of loop over grid !======================================================================= end do ! of do ig=1,ngridbidon(ig,:) = 1 !======================================================================= ! write outputs !======================================================================= ! value of parameter 'output' to trigger writting of outputs ! is set above at the declaration of the variable. if (output) then ! photoloysis do i=1,nb_phot_hv_max call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1), & 's-1',3,v_phot_3d(1,1,i)) end do ! iter call wstats(ngrid,'iter','iterations', & ' ',3,iter_3d(1,1)) ! phothcemical heating call wstats(ngrid,'zdtchim','dT chim', & 'K.s-1',3,zdtchim) call wstats(ngrid,'dEzchim','dE chim', & 'W m-2',3,dEzchim) call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output) ! Mean molecular mass call wstats(ngrid,'mmean','mean molecular mass', & 'g.mole-1',3,mmean(1,1)) ! Chemical tendencies do iesp=1,nesp call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', & 'kg/kg s^-1',3,dqchim(1,1,iesp) ) end do if (depos) then ! Surface fluxes and escape do iesp=1,nesp call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', & 'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp))))) call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', & 'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp))))) call wstats(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape', & 'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp))))) end do endif ! end depos ! photoloysis do i=1,nb_phot_hv_max call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1), & 's-1',3,v_phot_3d(1,1,i)) end do ! iter call writediagfi(ngrid,'iter','iterations', & ' ',3,iter_3d(1,1)) ! phothcemical heating call writediagfi(ngrid,'zdtchim','dT chim', & 'K.s-1',3,zdtchim) call writediagfi(ngrid,'dEzchim','dE chim', & 'W m-2',3,dEzchim) call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output) ! Mean molecular mass call writediagfi(ngrid,'mmean','mean molecular mass', & 'g.mole-1',3,mmean(1,1)) ! Chemical tendencies do iesp=1,nesp call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', & 'kg/kg s^-1',3,dqchim(1,1,iesp) ) end do if (depos) then ! Surface fluxes and escape do iesp=1,nesp call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', & 'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp))))) call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', & 'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp))))) call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape', & 'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp))))) end do ! Restart flux average if (ngrid == 1) then if(mod(icount,ecritphy).eq.0) then surface_flux2(:,:) = 0.0 endif else if(mod(icount*iphysiq,ecritphy).eq.0) then surface_flux2(:,:) = 0.0 endif endif endif ! end depos if (jonline .and. output_writediagspecUV) then !flux at the surface call writediagspecUV(ngrid,"flux_surf","Surface flux(lon,lat,band)","photon.s-1.nm-1.cm-2",3,fluxUV(:,:,1)) endif ! end jonline end if ! of if (output) return contains !====================================================================== subroutine ini_tracchim !====================================================================== ! YJ Modern tracer.def ! Get in tracer.def chemical parameters !====================================================================== use chimiedata_h use tracer_h, only: is_chim implicit none ! local character(len=500) :: tracline ! to store a line of text integer :: nq ! number of tracers integer :: iesp, iq ! Get nq open(90,file='traceur.def',status='old',form='formatted',iostat=ierr) if (ierr.eq.0) then READ(90,'(A)') tracline IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def ELSE DO READ(90,'(A)',iostat=ierr) tracline IF (ierr.eq.0) THEN IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header READ(tracline,*,iostat=ierr) nq EXIT ENDIF ELSE ! If pb, or if reached EOF without having found number of tracer write(*,*) "calchim: error reading line of tracers" write(*,*) " (first lines of traceur.def) " stop ENDIF ENDDO ENDIF ! if modern or standard traceur.def else write(*,*) "calchim: error opening traceur.def" stop endif ! Get data of tracers iesp = 0 do iq=1,nq read(90,'(A)') tracline if (is_chim(iq).eq.1) then iesp = iesp + 1 read(tracline,*) chemnoms(iesp) write(*,*)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp)) if (index(tracline,'SF_mode=' ) /= 0) then read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp) write(*,*) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp) else write(*,*) ' Parameter value (default) : SF_mode=', SF_mode(iesp) end if if (index(tracline,'SF_value=' ) /= 0) then read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp) write(*,*) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp) else write(*,*) ' Parameter value (default) : SF_value=', SF_value(iesp) end if if (index(tracline,'prod_rate=' ) /= 0) then read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp) write(*,*) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp) else write(*,*) ' Parameter value (default) : prod_rate=', prod_rate(iesp) end if if (index(tracline,'escape=' ) /= 0) then read(tracline(index(tracline,'escape=')+len('escape='):),*) escape(:,iesp) write(*,*) ' Parameter value (traceur.def) : escape=', escape(:,iesp) else write(*,*) ' Parameter value (default) : escape=', escape(:,iesp) end if end if end do close(90) end subroutine ini_tracchim end subroutine calchim_asis