Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/calchim.F (revision 616) +++ (revision ) @@ -1,502 +1,0 @@ - subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, - $ zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, - $ dqscloud,tauref,co2ice, - $ pu,pdu,pv,pdv,surfdust,surfice) -c - implicit none -c -c======================================================================= -c -c subject: -c -------- -c -c Prepare the call for the photochemical module, and send back the -c tendencies from photochemistry in the chemical species mass mixing ratios -c -c Author: Sebastien Lebonnois (08/11/2002) -c ------- -c update 12/06/2003 for water ice clouds and compatibility with dust -c update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll) -c update 03/05/2005 cosmetic changes (Franck Lefevre) -c update sept. 2008 identify tracers by their names (Ehouarn Millour) -c update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre) -c update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre) -c -c Arguments: -c ---------- -c -c Input: -c -c ptimestep timestep (s) -c pplay(ngridmx,nlayermx) Pressure at the middle of the layers (Pa) -c pplev(ngridmx,nlayermx+1) Intermediate pressure levels (Pa) -c pt(ngridmx,nlayermx) Temperature (K) -c pdt(ngridmx,nlayermx) Temperature tendency (K) -c pu(ngridmx,nlayermx) u component of the wind (ms-1) -c pdu(ngridmx,nlayermx) u component tendency (K) -c pv(ngridmx,nlayermx) v component of the wind (ms-1) -c pdv(ngridmx,nlayermx) v component tendency (K) -c dist_sol distance of the sun (AU) -c mu0(ngridmx) cos of solar zenith angle (=1 when sun at zenith) -c pq(ngridmx,nlayermx,nqmx) Advected fields, ie chemical species here -c pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq -c tauref(ngridmx) Optical depth at 7 hPa -c co2ice(ngridmx) co2 ice surface layer (kg.m-2) -c surfdust(ngridmx,nlayermx) dust surface area (m2/m3) -c surfice(ngridmx,nlayermx) ice surface area (m2/m3) -c -c Output: -c -c dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry -c dqschim(ngridmx,nqmx) ! tendencies on qsurf -c -c======================================================================= - -c Declarations : -c -------------- - -#include "dimensions.h" -#include "dimphys.h" -#include "chimiedata.h" -#include "tracer.h" -#include "comcstfi.h" -#include "callkeys.h" -#include "conc.h" - -c Arguments : -c ----------- - -c inputs: -c ------- - - real ptimestep - real pplay(ngridmx,nlayermx) ! pressure at the middle of the layers - real zzlay(ngridmx,nlayermx) ! pressure at the middle of the layers - real pplev(ngridmx,nlayermx+1) ! intermediate pressure levels - real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries - real pt(ngridmx,nlayermx) ! temperature - real pdt(ngridmx,nlayermx) ! temperature tendency - real pu(ngridmx,nlayermx) ! u component of the wind (m.s-1) - real pdu(ngridmx,nlayermx) ! u component tendency - real pv(ngridmx,nlayermx) ! v component of the wind (m.s-1) - real pdv(ngridmx,nlayermx) ! v component tendency - real dist_sol ! distance of the sun (AU) - real mu0(ngridmx) ! cos of solar zenith angle (=1 when sun at zenith) - real pq(ngridmx,nlayermx,nqmx) ! tracers mass mixing ratio - real pdq(ngridmx,nlayermx,nqmx) ! previous tendencies - real zday ! date (time since Ls=0, in martian days) - real tauref(ngridmx) ! optical depth at 7 hPa - real co2ice(ngridmx) ! co2 ice surface layer (kg.m-2) - real surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3) - real surfice(ngridmx,nlayermx) ! ice surface area (m2/m3) - -c outputs: -c -------- - - real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry - real dqschim(ngridmx,nqmx) ! tendencies on qsurf - real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation - real dqscloud(ngridmx,nqmx) ! tendencies on qsurf - -c Local variables : -c ----------------- - - character*20 str20 - integer ig,l,i,iq - integer foundswitch, lswitch - integer ig_vl1 - real latvl1, lonvl1 - real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry - ! new mole fraction after - real zt(ngridmx,nlayermx) ! temperature - real zu(ngridmx,nlayermx) ! u component of the wind - real zv(ngridmx,nlayermx) ! v component of the wind - real taucol ! optical depth at 7 hPa - - logical depos ! switch for dry deposition - parameter (depos=.false.) -c -c for each column of atmosphere: -c - real zpress(nlayermx) ! Pressure (mbar) - real zdens(nlayermx) ! Density (cm-3) - real ztemp(nlayermx) ! Temperature (K) - real zlocal(nlayermx) ! Altitude (km) - real zycol(nlayermx,nqmx) ! Composition (mole fractions) - real szacol ! Solar zenith angle - real surfice1d(nlayermx) ! Ice surface area (cm2/cm3) - real surfdust1d(nlayermx) ! Dust surface area (cm2/cm3) - real jo3(nlayermx) ! Photodissociation rate O3->O1D (s-1) -c -c for output: -c - real jo3_3d(ngridmx,nlayermx) ! Photodissociation rate O3->O1D (s-1) - - logical output ! to issue calls to writediagfi and stats - parameter (output=.true.) ! see at end of routine - - logical,save :: firstcall=.true. - integer,save :: nbq ! number of tracers used in the chemistry - integer,save :: niq(nqmx) ! array storing the indexes of the tracers - -! index of tracers: - - integer,save :: i_co2=0 - integer,save :: i_co=0 - integer,save :: i_o=0 - integer,save :: i_o1d=0 - integer,save :: i_o2=0 - integer,save :: i_o3=0 - integer,save :: i_h=0 - integer,save :: i_h2=0 - integer,save :: i_oh=0 - integer,save :: i_ho2=0 - integer,save :: i_h2o2=0 - integer,save :: i_ch4=0 - integer,save :: i_n2=0 - integer,save :: i_ar=0 - integer,save :: i_ice=0 ! water ice - integer,save :: i_h2o=0 ! water vapour -c -c======================================================================= -c initialization of the chemistry (first call only) -c======================================================================= -c - if (firstcall) then -c - if (photochem) then - print*,'calchim: Read photolysis lookup table' - call read_phototable - end if - - ! find index of chemical tracers to use - nbq=0 ! to count number of tracers - i_co2=igcm_co2 - if (i_co2.eq.0) then - write(*,*) "calchim: Error; no CO2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_co2 - endif - i_co=igcm_co - if (i_co.eq.0) then - write(*,*) "calchim: Error; no CO tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_co - endif - i_o=igcm_o - if (i_o.eq.0) then - write(*,*) "calchim: Error; no O tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_o - endif - i_o1d=igcm_o1d - if (i_o1d.eq.0) then - write(*,*) "calchim: Error; no O1D tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_o1d - endif - i_o2=igcm_o2 - if (i_o2.eq.0) then - write(*,*) "calchim: Error; no O2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_o2 - endif - i_o3=igcm_o3 - if (i_o3.eq.0) then - write(*,*) "calchim: Error; no O3 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_o3 - endif - i_h=igcm_h - if (i_h.eq.0) then - write(*,*) "calchim: Error; no H tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_h - endif - i_h2=igcm_h2 - if (i_h2.eq.0) then - write(*,*) "calchim: Error; no H2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_h2 - endif - i_oh=igcm_oh - if (i_oh.eq.0) then - write(*,*) "calchim: Error; no OH tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_oh - endif - i_ho2=igcm_ho2 - if (i_ho2.eq.0) then - write(*,*) "calchim: Error; no HO2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_ho2 - endif - i_h2o2=igcm_h2o2 - if (i_h2o2.eq.0) then - write(*,*) "calchim: Error; no H2O2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_h2o2 - endif - i_ch4=igcm_ch4 - if (i_ch4.eq.0) then - write(*,*) "calchim: Error; no CH4 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_ch4 - endif - i_n2=igcm_n2 - if (i_n2.eq.0) then - write(*,*) "calchim: Error; no N2 tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_n2 - endif - i_ar=igcm_ar - if (i_ar.eq.0) then - write(*,*) "calchim: Error; no AR tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_ar - endif - i_ice=igcm_h2o_ice - if (i_ice.eq.0) then - write(*,*) "calchim: Error; no water ice tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_ice - endif - i_h2o=igcm_h2o_vap - if (i_h2o.eq.0) then - write(*,*) "calchim: Error; no water vapor tracer !!!" - stop - else - nbq=nbq+1 - niq(nbq)=i_h2o - endif - - write(*,*) 'calchim: found nbq = ',nbq,' tracers' - write(*,*) ' i_co2 = ',i_co2 - write(*,*) ' i_co = ',i_co - write(*,*) ' i_o = ',i_o - write(*,*) ' i_o1d = ',i_o1d - write(*,*) ' i_o2 = ',i_o2 - write(*,*) ' i_o3 = ',i_o3 - write(*,*) ' i_h = ',i_h - write(*,*) ' i_h2 = ',i_h2 - write(*,*) ' i_oh = ',i_oh - write(*,*) ' i_ho2 = ',i_ho2 - write(*,*) ' i_h2o2 = ',i_h2o2 - write(*,*) ' i_ch4 = ',i_ch4 - write(*,*) ' i_n2 = ',i_n2 - write(*,*) ' i_ar = ',i_ar - write(*,*) ' i_ice = ',i_ice - write(*,*) ' i_h2o = ',i_h2o - - firstcall = .false. - end if ! if (firstcall) - -! Initialize output tendencies to zero (to handle case of tracers which -! are not used in the chemistry (e.g. dust)) - - dqchim(:,:,:) = 0 - dqschim(:,:) = 0 -c -! latvl1= 22.27 -! lonvl1= -47.94 -! ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + -! $ int(1.5+(lonvl1+180)*iim/360.) -c -c======================================================================= -c loop over grid -c======================================================================= -c - do ig = 1,ngridmx -c -c local updates -c - foundswitch = 0 - do l = 1,nlayermx - do i = 1,nbq - iq = niq(i) ! get tracer index - zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep - zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq) - 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 -c - 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. -c -c surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3 -c - surfdust1d(l) = surfdust(ig,l)*1.e-2 - surfice1d(l) = surfice(ig,l)*1.e-2 -c -c search for switch index between regions -c - if (photochem .and. thermochem) then - if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then - lswitch = l - foundswitch=1 - end if - end if - if ( .not. photochem) then - lswitch = 22 - end if - if (.not. thermochem) then - lswitch = min(50,nlayermx+1) ! FL (original value: 33) - end if -c - end do ! of do l=1,nlayermx -c - szacol = acos(mu0(ig))*180./pi - taucol = tauref(ig) -c -c======================================================================= -c call chemical subroutine -c======================================================================= -c -c chemistry in lower atmosphere -c - if (photochem) then - call photochemistry(lswitch,zycol,szacol,ptimestep, - $ zpress,ztemp,zdens,dist_sol, - $ surfdust1d,surfice1d,jo3,taucol) - end if -c -c chemistry in upper atmosphere -c - if (thermochem) then - call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress, - $ zlocal,szacol,ptimestep,zday) - end if -c -c dry deposition -c - if (depos) then - call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev, - $ zu, zv, zt, zycol, ptimestep, co2ice) - end if -c -c======================================================================= -c tendencies -c======================================================================= -c -c must be 0. for water ice: -c - if (water) then - do l = 1,nlayermx - dqchim(ig,l,i_ice) = 0. - end do - end if -c -c tendency for CO2 = - sum of others for lower atmosphere -c tendency for O = - sum of others for upper atmosphere -c - do l = 1,nlayermx - if (l .lt. lswitch) then - do i=1,nbq - iq=niq(i) ! get tracer index - if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then - dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) - $ - zq(ig,l,iq))/ptimestep - else if (iq.eq.i_co2) then - dqchim(ig,l,iq) = 0. - end if - dqschim(ig,iq) = 0. - end do ! of do i=1,nbq - - do i=1,nbq - iq=niq(i) ! get tracer index - if (iq.ne.i_co2) then - dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2) - $ - dqchim(ig,l,iq) - end if - end do - else if (l .ge. lswitch) then - do i=1,nbq - iq=niq(i) ! get tracer index - if ((iq.ne.i_o).and.(iq.ne.i_ice)) then - dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq) - $ /mmean(ig,l) - $ - zq(ig,l,iq))/ptimestep - else if (iq.eq.i_o) then - dqchim(ig,l,iq) = 0. - end if - enddo - - do i=1,nbq - iq=niq(i) ! get tracer index - if (iq.ne.i_o) then - dqchim(ig,l,i_o) = dqchim(ig,l,i_o) - $ - dqchim(ig,l,iq) - end if - end do - end if ! of if (l.lt.lswitch) else if (l.ge.lswitch) - end do ! of do l = 1,nlayermx -c -c condensation of h2o2 -c - call perosat(ig,ptimestep,pplev,pplay, - $ ztemp,zycol,dqcloud,dqscloud) -c -c density and j(o3->o1d), for outputs -c - do l = 1,nlayermx - jo3_3d(ig,l) = jo3(l) - end do -c -c======================================================================= -c end of loop over grid -c======================================================================= -c - end do ! of do ig=1,ngridmx -c -c======================================================================= -c write outputs -c======================================================================= -c -! value of parameter 'output' to trigger writting of outputs -! is set above at the declaration of the variable. - - if (output) then - if (ngridmx .gt. 1) then - call writediagfi(ngridmx,'jo3','j o3->o1d', - $ 's-1',3,jo3_3d(1,1)) - if (callstats) then - call wstats(ngridmx,'jo3','j o3->o1d', - $ 's-1',3,jo3_3d(1,1)) - endif - end if ! of if (ngridmx.gt.1) - endif ! of if (output) -c - end Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90 (revision 618) +++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90 (revision 618) @@ -0,0 +1,420 @@ + subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, & + zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, & + dqscloud,tauref,co2ice, & + pu,pdu,pv,pdv,surfdust,surfice) + + 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) +! +! Arguments: +! ---------- +! +! Input: +! +! ptimestep timestep (s) +! pplay(ngridmx,nlayermx) Pressure at the middle of the layers (Pa) +! pplev(ngridmx,nlayermx+1) Intermediate pressure levels (Pa) +! pt(ngridmx,nlayermx) Temperature (K) +! pdt(ngridmx,nlayermx) Temperature tendency (K) +! pu(ngridmx,nlayermx) u component of the wind (ms-1) +! pdu(ngridmx,nlayermx) u component tendency (K) +! pv(ngridmx,nlayermx) v component of the wind (ms-1) +! pdv(ngridmx,nlayermx) v component tendency (K) +! dist_sol distance of the sun (AU) +! mu0(ngridmx) cos of solar zenith angle (=1 when sun at zenith) +! pq(ngridmx,nlayermx,nqmx) Advected fields, ie chemical species here +! pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq +! tauref(ngridmx) Optical depth at 7 hPa +! co2ice(ngridmx) co2 ice surface layer (kg.m-2) +! surfdust(ngridmx,nlayermx) dust surface area (m2/m3) +! surfice(ngridmx,nlayermx) ice surface area (m2/m3) +! +! Output: +! +! dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry +! dqschim(ngridmx,nqmx) ! tendencies on qsurf +! +!======================================================================= + +#include "dimensions.h" +#include "dimphys.h" +#include "chimiedata.h" +#include "tracer.h" +#include "comcstfi.h" +#include "callkeys.h" +#include "conc.h" + +! input: + + real :: ptimestep + real :: pplay(ngridmx,nlayermx) ! pressure at the middle of the layers + real :: zzlay(ngridmx,nlayermx) ! pressure at the middle of the layers + real :: pplev(ngridmx,nlayermx+1) ! intermediate pressure levels + real :: zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries + real :: pt(ngridmx,nlayermx) ! temperature + real :: pdt(ngridmx,nlayermx) ! temperature tendency + real :: pu(ngridmx,nlayermx) ! u component of the wind (m.s-1) + real :: pdu(ngridmx,nlayermx) ! u component tendency + real :: pv(ngridmx,nlayermx) ! v component of the wind (m.s-1) + real :: pdv(ngridmx,nlayermx) ! v component tendency + real :: dist_sol ! distance of the sun (AU) + real :: mu0(ngridmx) ! cos of solar zenith angle (=1 when sun at zenith) + real :: pq(ngridmx,nlayermx,nqmx) ! tracers mass mixing ratio + real :: pdq(ngridmx,nlayermx,nqmx) ! previous tendencies + real :: zday ! date (time since Ls=0, in martian days) + real :: tauref(ngridmx) ! optical depth at 7 hPa + real :: co2ice(ngridmx) ! co2 ice surface layer (kg.m-2) + real :: surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3) + real :: surfice(ngridmx,nlayermx) ! ice surface area (m2/m3) + +! output: + + real :: dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry + real :: dqschim(ngridmx,nqmx) ! tendencies on qsurf + real :: dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation + real :: dqscloud(ngridmx,nqmx) ! tendencies on qsurf + +! local variables: + + integer,save :: nbq ! number of tracers used in the chemistry + integer,save :: niq(nqmx) ! array storing the indexes of the tracers + integer :: major(nlayermx) ! index of major species + integer :: ig,l,i,iq,iqmax + integer :: foundswitch, lswitch + + integer,save :: i_co2 = 0 + integer,save :: i_co = 0 + integer,save :: i_o = 0 + integer,save :: i_o1d = 0 + integer,save :: i_o2 = 0 + integer,save :: i_o3 = 0 + integer,save :: i_h = 0 + integer,save :: i_h2 = 0 + integer,save :: i_oh = 0 + integer,save :: i_ho2 = 0 + integer,save :: i_h2o2 = 0 + integer,save :: i_ch4 = 0 + integer,save :: i_n2 = 0 + integer,save :: i_h2o = 0 + + integer :: ig_vl1 + + real :: latvl1, lonvl1 + real :: zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry + ! new mole fraction after + real :: zt(ngridmx,nlayermx) ! temperature + real :: zu(ngridmx,nlayermx) ! u component of the wind + real :: zv(ngridmx,nlayermx) ! v component of the wind + real :: taucol ! optical depth at 7 hPa + + logical,save :: firstcall = .true. + logical,save :: depos = .false. ! switch for dry deposition + +! for each column of atmosphere: + + real :: zpress(nlayermx) ! Pressure (mbar) + real :: zdens(nlayermx) ! Density (cm-3) + real :: ztemp(nlayermx) ! Temperature (K) + real :: zlocal(nlayermx) ! Altitude (km) + real :: zycol(nlayermx,nqmx) ! Composition (mole fractions) + real :: szacol ! Solar zenith angle + real :: surfice1d(nlayermx) ! Ice surface area (cm2/cm3) + real :: surfdust1d(nlayermx) ! Dust surface area (cm2/cm3) + real :: jo3(nlayermx) ! Photodissociation rate O3->O1D (s-1) + +! for output: + + logical :: output ! to issue calls to writediagfi and stats + parameter (output = .true.) + real :: jo3_3d(ngridmx,nlayermx) ! Photodissociation rate O3->O1D (s-1) + +!======================================================================= +! initialization of the chemistry (first call only) +!======================================================================= + + if (firstcall) then + + if (photochem) then + print*,'calchim: Read photolysis lookup table' + call read_phototable + end if + + ! find index of chemical tracers to use + nbq = 0 ! to count number of tracers + i_co2 = igcm_co2 + if (i_co2 == 0) then + write(*,*) "calchim: Error; no CO2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_co2 + end if + i_co = igcm_co + if (i_co == 0) then + write(*,*) "calchim: Error; no CO tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_co + end if + i_o = igcm_o + if (i_o == 0) then + write(*,*) "calchim: Error; no O tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_o + end if + i_o1d = igcm_o1d + if (i_o1d == 0) then + write(*,*) "calchim: Error; no O1D tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_o1d + end if + i_o2 = igcm_o2 + if (i_o2 == 0) then + write(*,*) "calchim: Error; no O2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_o2 + end if + i_o3 = igcm_o3 + if (i_o3 == 0) then + write(*,*) "calchim: Error; no O3 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_o3 + end if + i_h = igcm_h + if (i_h == 0) then + write(*,*) "calchim: Error; no H tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_h + end if + i_h2 = igcm_h2 + if (i_h2 == 0) then + write(*,*) "calchim: Error; no H2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_h2 + end if + i_oh = igcm_oh + if (i_oh == 0) then + write(*,*) "calchim: Error; no OH tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_oh + end if + i_ho2 = igcm_ho2 + if (i_ho2 == 0) then + write(*,*) "calchim: Error; no HO2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_ho2 + end if + i_h2o2 = igcm_h2o2 + if (i_h2o2 == 0) then + write(*,*) "calchim: Error; no H2O2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_h2o2 + end if + i_ch4 = igcm_ch4 + if (i_ch4 == 0) then + write(*,*) "calchim: no CH4 tracer !!!" + write(*,*) "CH4 will be ignored in the chemistry" + else + nbq = nbq + 1 + niq(nbq) = i_ch4 + end if + i_n2 = igcm_n2 + if (i_n2 == 0) then + write(*,*) "calchim: Error; no N2 tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_n2 + end if + i_h2o = igcm_h2o_vap + if (i_h2o == 0) then + write(*,*) "calchim: Error; no water vapor tracer !!!" + stop + else + nbq = nbq + 1 + niq(nbq) = i_h2o + end if + + write(*,*) 'calchim: found nbq = ',nbq,' tracers' + + firstcall = .false. + end if ! if (firstcall) + +! Initializations + + zycol(:,:) = 0. + dqchim(:,:,:) = 0 + dqschim(:,:) = 0 + +! latvl1= 22.27 +! lonvl1= -47.94 +! ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + & +! int(1.5+(lonvl1+180)*iim/360.) + +!======================================================================= +! loop over grid +!======================================================================= + + do ig = 1,ngridmx + + foundswitch = 0 + do l = 1,nlayermx + do i = 1,nbq + iq = niq(i) ! get tracer index + zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep + zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq) + 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. + +! 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-3) then + lswitch = l + foundswitch = 1 + end if + end if + if (.not. photochem) then + lswitch = 22 + end if + if (.not. thermochem) then + lswitch = min(50,nlayermx+1) + end if + + end do ! of do l=1,nlayermx + + szacol = acos(mu0(ig))*180./pi + taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars + +!======================================================================= +! call chemical subroutines +!======================================================================= + +! chemistry in lower atmosphere + + if (photochem) then + call photochemistry(lswitch,zycol,szacol,ptimestep, & + zpress,ztemp,zdens,dist_sol, & + surfdust1d,surfice1d,jo3,taucol) + +! ozone photolysis, for output + + do l = 1,nlayermx + jo3_3d(ig,l) = jo3(l) + end do + +! condensation of h2o2 + + call perosat(ig,ptimestep,pplev,pplay, & + ztemp,zycol,dqcloud,dqscloud) + end if + +! chemistry in upper atmosphere + + if (thermochem) then + call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress, & + zlocal,szacol,ptimestep,zday) + end if + +! dry deposition + + if (depos) then + call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,& + zu, zv, zt, zycol, ptimestep, co2ice) + 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,nlayermx + iqmax = major(l) + do i = 1,nbq + iq = niq(i) ! get tracer index + if (iq /= iqmax) then + dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) & + - 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,nlayermx + +!======================================================================= +! end of loop over grid +!======================================================================= + + end do ! of do ig=1,ngridmx + +!======================================================================= +! write outputs +!======================================================================= + +! value of parameter 'output' to trigger writting of outputs +! is set above at the declaration of the variable. + + if (photochem .and. output) then + if (ngridmx > 1) then + call writediagfi(ngridmx,'jo3','j o3->o1d', & + 's-1',3,jo3_3d(1,1)) + if (callstats) then + call wstats(ngridmx,'jo3','j o3->o1d', & + 's-1',3,jo3_3d(1,1)) + endif + end if ! of if (ngridmx.gt.1) + end if ! of if (output) + + return + end Index: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/concentrations.F (revision 616) +++ trunk/LMDZ.MARS/libf/aeronomars/concentrations.F (revision 618) @@ -3,17 +3,16 @@ implicit none -c======================================================================= -c CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R -c -c mmean(ngridmx,nlayermx) amu -c cpnew(ngridmx,nlayermx) J/kg/K -c rnew(ngridmx,nlayermx) J/kg/K -c akknew(ngridmx,nlayermx) coefficient of thermal concduction -c -c version: March 2011 - Franck Lefevre -c======================================================================= +!======================================================================= +! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R +! +! mmean(ngridmx,nlayermx) amu +! cpnew(ngridmx,nlayermx) J/kg/K +! rnew(ngridmx,nlayermx) J/kg/K +! akknew(ngridmx,nlayermx) coefficient of thermal concduction +! +! version: April 2012 - Franck Lefevre +!======================================================================= -c Declarations -c ------------ +! declarations #include "dimensions.h" @@ -26,6 +25,5 @@ #include "conc.h" -c Input/Output -c ------------ +! input/output real pplay(ngridmx,nlayermx) @@ -36,109 +34,109 @@ real ptimestep -c Local variables -c --------------- +! local variables - integer :: l, ig, n, k - integer, save :: gind(ncomp) + integer :: i, l, ig, iq + integer, save :: nbq, niq(nqmx) real :: ni(nqmx), ntot - real :: zq(ngridmx,nlayermx,ncomp) - real :: zt(ngridmx,nlayermx) - real, save :: aki(ncomp) - real, save :: cpi(ncomp) + real :: zq(ngridmx, nlayermx, nqmx) + real :: zt(ngridmx, nlayermx) + real, save :: aki(nqmx) + real, save :: cpi(nqmx) logical, save :: firstcall = .true. -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c tracer numbering for the thermal conduction and -c specific heat coefficients -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - integer,parameter :: i_co2 = 1 - integer,parameter :: i_co = 2 - integer,parameter :: i_o = 3 - integer,parameter :: i_o1d = 4 - integer,parameter :: i_o2 = 5 - integer,parameter :: i_o3 = 6 - integer,parameter :: i_h = 7 - integer,parameter :: i_h2 = 8 - integer,parameter :: i_oh = 9 - integer,parameter :: i_ho2 = 10 - integer,parameter :: i_h2o2 = 11 - integer,parameter :: i_ch4 = 12 - integer,parameter :: i_n2 = 13 - integer,parameter :: i_ar = 14 - integer,parameter :: i_h2o = 15 - if (firstcall) then -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c initializations at first call: -c fill local array of tracer indexes -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc +! find index of chemical tracers to use +! initialize thermal conductivity and specific heat coefficients +! !? values are estimated - gind(i_co2) = igcm_co2 ! co2 - gind(i_co) = igcm_co ! co - gind(i_o) = igcm_o ! o - gind(i_o1d) = igcm_o1d ! o1d - gind(i_o2) = igcm_o2 ! o2 - gind(i_o3) = igcm_o3 ! o3 - gind(i_h) = igcm_h ! h - gind(i_h2) = igcm_h2 ! h2 - gind(i_oh) = igcm_oh ! oh - gind(i_ho2) = igcm_ho2 ! ho2 - gind(i_h2o2) = igcm_h2o2 ! h2o2 - gind(i_ch4) = igcm_ch4 ! ch4 - gind(i_n2) = igcm_n2 ! n2 - gind(i_ar) = igcm_ar ! ar - gind(i_h2o) = igcm_h2o_vap ! h2o + nbq = 0 ! to count number of tracers used in this subroutine -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c Thermal conductivity and specific heat coefficients -cccccccccccccccccccccccccccccccccccccccccccccccccccccccc + if (igcm_co2 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_co2 + aki(nbq) = 3.072e-4 + cpi(nbq) = 0.834e3 + end if + if (igcm_co /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_co + aki(nbq) = 4.87e-4 + cpi(nbq) = 1.034e3 + end if + if (igcm_o /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_o + aki(nbq) = 7.59e-4 + cpi(nbq) = 1.3e3 + end if + if (igcm_o1d /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_o1d + aki(nbq) = 7.59e-4 !? + cpi(nbq) = 1.3e3 !? + end if + if (igcm_o2 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_o2 + aki(nbq) = 5.68e-4 + cpi(nbq) = 0.9194e3 + end if + if (igcm_o3 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_o3 + aki(nbq) = 3.00e-4 !? + cpi(nbq) = 0.800e3 !? + end if + if (igcm_h /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_h + aki(nbq) = 0.0 + cpi(nbq) = 20.780e3 + end if + if (igcm_h2 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_h2 + aki(nbq) = 36.314e-4 + cpi(nbq) = 14.266e3 + end if + if (igcm_oh /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_oh + aki(nbq) = 7.00e-4 !? + cpi(nbq) = 1.045e3 + end if + if (igcm_ho2 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_ho2 + aki(nbq) = 0.0 + cpi(nbq) = 1.065e3 !? + end if + if (igcm_n2 /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_n2 + aki(nbq) = 5.6e-4 + cpi(nbq) = 1.034e3 + end if + if (igcm_ar /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_ar + aki(nbq) = 0.0 !? + cpi(nbq) = 1.000e3 !? + end if + if (igcm_h2o_vap /= 0) then + nbq = nbq + 1 + niq(nbq) = igcm_h2o_vap + aki(nbq) = 0.0 + cpi(nbq) = 1.870e3 + end if - aki(i_co2) = 3.072e-4 - aki(i_co) = 4.87e-4 - aki(i_o) = 7.59e-4 - aki(i_o1d) = 7.59e-4 !? - aki(i_o2) = 5.68e-4 - aki(i_o3) = 3.00e-4 !? - aki(i_h) = 0.0 - aki(i_h2) = 36.314e-4 - aki(i_oh) = 7.00e-4 !? - aki(i_ho2) = 0.0 - aki(i_h2o2) = 0.0 - aki(i_ch4) = 0.0 !? - aki(i_n2) = 5.6e-4 - aki(i_ar) = 0.0 !? - aki(i_h2o) = 0.0 + firstcall = .false. - cpi(i_co2) = 0.834e3 - cpi(i_co) = 1.034e3 - cpi(i_o) = 1.3e3 - cpi(i_o1d) = 1.3e3 !? - cpi(i_o2) = 0.9194e3 - cpi(i_o3) = 0.800e3 !? - cpi(i_h) = 20.780e3 - cpi(i_h2) = 14.266e3 - cpi(i_oh) = 1.045e3 - cpi(i_ho2) = 1.065e3 !? - cpi(i_h2o2) = 1.000e3 !? - cpi(i_ch4) = 1.000e3 !? - cpi(i_n2) = 1.034e3 - cpi(i_ar) = 1.000e3 !? - cpi(i_h2o) = 1.870e3 + end if ! if (firstcall) - firstcall=.false. +! update temperature - end if ! of if (firstcall) -c -c initializations -c - mmean(:,:) = 0. - cpnew(:,:) = 0. - akknew(:,:) = 0. -c -c update temperature -c do l = 1,nlayermx do ig = 1,ngridmx @@ -146,23 +144,27 @@ end do end do -c -c update tracers -c + +! update tracers + do l = 1,nlayermx do ig = 1,ngridmx - do n = 1,ncomp - zq(ig,l,n) = max(1.e-30, pq(ig,l,gind(n)) - $ + pdq(ig,l,gind(n))*ptimestep) + do i = 1,nbq + iq = niq(i) + zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq) + $ + pdq(ig,l,iq)*ptimestep) end do end do end do -c -c mmean : mean molecular mass -c rnew : specific gas constant -c + +! mmean : mean molecular mass +! rnew : specific gas constant + + mmean(:,:) = 0. + do l = 1,nlayermx do ig = 1,ngridmx - do n = 1, ncomp - mmean(ig,l) = mmean(ig,l) + zq(ig,l,n)/mmol(gind(n)) + do i = 1,nbq + iq = niq(i) + mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq) end do mmean(ig,l) = 1./mmean(ig,l) @@ -170,20 +172,24 @@ end do end do -c -c cpnew : specicic heat -c akknew : thermal conductivity cofficient -c + +! cpnew : specicic heat +! akknew : thermal conductivity cofficient + + cpnew(:,:) = 0. + akknew(:,:) = 0. + do l = 1,nlayermx do ig = 1,ngridmx ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6 ! in #/cm3 - do n = 1,ncomp - ni(n) = ntot*zq(ig,l,n)*mmean(ig,l)/mmol(gind(n)) - cpnew(ig,l) = cpnew(ig,l) + ni(n)*cpi(n) - akknew(ig,l) = akknew(ig,l) + ni(n)*aki(n) + do i = 1,nbq + iq = niq(i) + ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq) + cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(iq) + akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(iq) end do cpnew(ig,l) = cpnew(ig,l)/ntot akknew(ig,l) = akknew(ig,l)/ntot end do -c print*, l, mmean(1,l), cpnew(1,l), rnew(1,l) +! print*, l, mmean(1,l), cpnew(1,l), rnew(1,l) end do Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F (revision 616) +++ (revision ) @@ -1,644 +1,0 @@ - subroutine inichim_newstart(pq,ps,sig,nq,latfi,lonfi,airefi, - $ thermo,qsurf) - - implicit none - -c======================================================================= -c -c subject: -c -------- -c -c Initialization of the composition for use in the building of a new start file -c used by program newstart.F -c also used by program testphys1d.F -c -c Author: Sebastien Lebonnois (08/11/2002) -c ------- -c Revised 07/2003 by Monica Angelats-i-Coll to use input files -c Modified 10/2008 Identify tracers by their names Ehouarn Millour -c Modified 11/2011 Addition of methane (Franck Lefevre) -c -c Arguments: -c ---------- -c -c pq(iip1,jjp1,llm,nqmx) Advected fields, ie chemical species here -c sig = sigma vertical coordinate (interface of layers) -c nq: number of tracers to initialize (used to evaluate if only -c chemistry species are to be initialized or if water vapour -c and water ice should also be initialized. -c NB: number of chemistry tracers is ncomp (set in chimiedata.h) -c======================================================================= - -c Declarations : -c -------------- - -#include "dimensions.h" -#include "dimphys.h" -#include "paramet.h" -#include "chimiedata.h" -#include "tracer.h" -#include "comcstfi.h" -#include "comdiurn.h" -#include "callkeys.h" -#include "temps.h" -#include "datafile.h" - -c Arguments : -c ----------- - - real ps(iip1,jjp1) - real pq(iip1,jjp1,llm,nqmx) ! Advected fields, ie chemical species - real sig(llm+1) ! vertical coordinate (interface of layers) - integer nq ! =nqmx - ! or =nqmx-1 if H2O kept - ! or =nqmx-2 if H2O and ice kept - REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) - integer thermo ! flag for thermosphere initialisation only - real :: qsurf(ngridmx,nqmx) ! surface tracers - -c Local variables : -c ----------------- - - integer iq,i,j,l,n, nspecini,inil,iqmax,iiq - INTEGER aa(nqmx) - REAL qtot(iip1,jjp1,llm) - real densconc,zgcm,zfile(252),vmrint,nt, mmean - real nxf(252),zfilemin(252),sigfile(252) - real vmrf(252,14),nf(252,14) - real tfile(252),zzfile(252) - real ni(14),niold(14) - real mmolf(14) ! molecular mass in amu (species in files) - data mmolf/44.,40.,28.,32.,28.,16.,2., ! majors - & 1.,17.,33.,18.,34.,16.,48./ ! minors - character*3 tmp ! temporary variable - integer ierr - - logical :: oldnames ! =.true. if old tracer naming convention (q01,...) - character(len=20) :: txt ! to store some text - integer :: nqchem_start,count - integer :: nqchem(nqmx) ! local chemistry tracer index array - integer :: nbqchem ! total number of chemical species (water included) - -c----------------------------------------------------------------------- -c Input/Output -c ------------ -! read in 'callphys.def' - call inichim_readcallphys() - - ! check if tracers have 'old' names - oldnames=.false. - count=0 - do iq=1,nqmx - txt=" " - write(txt,'(a1,i2.2)') 'q',iq - if (txt.eq.noms(iq)) then - count=count+1 - endif - enddo ! of do iq=1,nqmx - - if (count.eq.nqmx) then - write(*,*) "inichim_newstart: tracers seem to follow old ", - & "naming convention (q01,q02,...)" - write(*,*) " => will work for now ... " - write(*,*) " but you should rename them" - oldnames=.true. - endif - - ! copy/set tracer names - if (oldnames) then ! old names (derived from iq & option values) - ! 1. dust: - if (dustbin.ne.0) then ! transported dust - do iq=1,dustbin - txt=" " - write(txt,'(a4,i2.2)') 'dust',iq - noms(iq)=txt - mmol(iq)=100. - enddo - ! special case if doubleq - if (doubleq) then - if (dustbin.ne.2) then - write(*,*) 'initracer: error expected dustbin=2' - else -! noms(1)='dust_mass' ! dust mass mixing ratio -! noms(2)='dust_number' ! dust number mixing ratio -! same as above, but avoid explict possible out-of-bounds on array - noms(1)='dust_mass' ! dust mass mixing ratio - do iq=2,2 - noms(iq)='dust_number' ! dust number mixing ratio - enddo - endif - endif - endif - ! 2. water & ice - if (water) then - noms(nqmx)='h2o_vap' - mmol(nqmx)=18. -! noms(nqmx-1)='h2o_ice' -! mmol(nqmx-1)=18. - !"loop" to avoid potential out-of-bounds in array - do iq=nqmx-1,nqmx-1 - noms(iq)='h2o_ice' - mmol(iq)=18. - enddo - endif - ! 3. Chemistry - if (photochem .or. callthermos) then - if (doubleq) then - nqchem_start=3 - else - nqchem_start=dustbin+1 - end if - endif ! of if (photochem .or. callthermos) - do iq=nqchem_start,nqchem_start+ncomp-2 - noms(iq)=nomchem(iq-nqchem_start+1) - mmol(iq)=mmolchem(iq-nqchem_start+1) - enddo - ! 4. Other tracers ???? - if ((dustbin.eq.0).and.(.not.water)) then - noms(1)='co2' - mmol(1)=44 - if (nqmx.eq.2) then - noms(nqmx)='Ar_N2' - mmol(nqmx)=30 - endif - endif - else - ! noms(:) already contain tracer names - ! mmol(:) array is initialized later (see below) - endif ! of if (oldnames) - - ! special modification when using 'old' tracers: - ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour - ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice - if (oldnames.and.water) then - write(*,*)'inichim_newstart: moving surface water ice to index ' - & ,nqmx-1 - ! "loop" to avoid potential out-of-bounds on arrays - do iq=nqmx-1,nqmx-1 - qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1) - enddo - qsurf(1:ngridmx,nqmx)=0 - endif - -c Dimension test: -c --------------- - - if (water) then -! if ((nqchem_min+ncomp+1).ne.nqmx) then - if ((dustbin+ncomp+2).ne.nqmx) then -! print*,'********* Dimension problem! ********' -! print*,"nqchem_min+ncomp+1).ne.nqmx" - print*,'inichim_newstart: tracer dimension problem:' - print*,"(dustbin+ncomp+2).ne.nqmx" - print*,"ncomp: ",ncomp -! print*,"nqchem_min: ",nqchem_min - print*,"nqmx: ",nqmx - print*,'Change ncomp in chimiedata.h' - STOP - endif - else -! if ((nqchem_min+ncomp).ne.nqmx) then - if ((dustbin+ncomp+1).ne.nqmx) then -! print*,'********* Dimension problem! ********' -! print*,"nqchem_min+ncomp).ne.nqmx" - print*,'initracer: tracer dimension problem:' - print*,"(dustbin+ncomp+1).ne.nqmx" - print*,"ncomp: ",ncomp -! print*,"nqchem_min: ",nqchem_min - print*,"nqmx: ",nqmx - print*,'Change ncomp in chimiedata.h' - STOP - endif - endif - -c noms and mmol vectors: -c ---------------------- -! -! if (iceparty) then -! do iq=nqchem_min,nqmx-2 -! noms(iq) = nomchem(iq-nqchem_min+1) -! mmol(iq) = mmolchem(iq-nqchem_min+1) -! enddo -! noms(nqmx-1) = 'ice' -! mmol(nqmx-1) = 18. -! noms(nqmx) = 'h2o' -! mmol(nqmx) = 18. -! else -! do iq=nqchem_min,nqmx-1 -! noms(iq) = nomchem(iq-nqchem_min+1) -! mmol(iq) = mmolchem(iq-nqchem_min+1) -! enddo -! noms(nqmx) = 'h2o' -! mmol(nqmx) = 18. -! endif -! if (nqchem_min.gt.1) then -! do iq=1,nqchem_min-1 -! noms(iq) = 'dust' -! mmol(iq) = 100. -! enddo -! endif - - -! Identify tracers by their names: (and set corresponding values of mmol) - ! 0. initialize tracer indexes to zero: - do iq=1,nqmx - igcm_dustbin(iq)=0 - enddo - igcm_dust_mass=0 - igcm_dust_number=0 - igcm_h2o_vap=0 - igcm_h2o_ice=0 - igcm_co2=0 - igcm_co=0 - igcm_o=0 - igcm_o1d=0 - igcm_o2=0 - igcm_o3=0 - igcm_h=0 - igcm_h2=0 - igcm_oh=0 - igcm_ho2=0 - igcm_h2o2=0 - igcm_ch4=0 - igcm_n2=0 - igcm_ar=0 - igcm_ar_n2=0 - - ! 1. find dust tracers - count=0 - if (dustbin.gt.0) then - do iq=1,nqmx - txt=" " - write(txt,'(a4,i2.2)')'dust',count+1 - if (noms(iq).eq.txt) then - count=count+1 - igcm_dustbin(count)=iq - mmol(iq)=100. - endif - enddo !do iq=1,nqmx - endif ! of if (dustbin.gt.0) - if (doubleq) then - do iq=1,nqmx - if (noms(iq).eq."dust_mass") then - igcm_dust_mass=iq - count=count+1 - endif - if (noms(iq).eq."dust_number") then - igcm_dust_number=iq - count=count+1 - endif - enddo - endif ! of if (doubleq) - ! 2. find chemistry and water tracers - nbqchem=0 - do iq=1,nqmx - if (noms(iq).eq."co2") then - igcm_co2=iq - mmol(igcm_co2)=44. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."co") then - igcm_co=iq - mmol(igcm_co)=28. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."o") then - igcm_o=iq - mmol(igcm_o)=16. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."o1d") then - igcm_o1d=iq - mmol(igcm_o1d)=16. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."o2") then - igcm_o2=iq - mmol(igcm_o2)=32. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."o3") then - igcm_o3=iq - mmol(igcm_o3)=48. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."h") then - igcm_h=iq - mmol(igcm_h)=1. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."h2") then - igcm_h2=iq - mmol(igcm_h2)=2. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."oh") then - igcm_oh=iq - mmol(igcm_oh)=17. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."ho2") then - igcm_ho2=iq - mmol(igcm_ho2)=33. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."h2o2") then - igcm_h2o2=iq - mmol(igcm_h2o2)=34. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."ch4") then - igcm_ch4=iq - mmol(igcm_ch4)=16. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."n2") then - igcm_n2=iq - mmol(igcm_n2)=28. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."ar") then - igcm_ar=iq - mmol(igcm_ar)=40. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."h2o_vap") then - igcm_h2o_vap=iq - mmol(igcm_h2o_vap)=18. - count=count+1 - nbqchem=nbqchem+1 - endif - if (noms(iq).eq."h2o_ice") then - igcm_h2o_ice=iq - mmol(igcm_h2o_ice)=18. - count=count+1 - nbqchem=nbqchem+1 - endif - ! Other stuff: e.g. for simulations using co2 + neutral gaz - if (noms(iq).eq."Ar_N2") then - igcm_ar_n2=iq - mmol(igcm_ar_n2)=30. - count=count+1 - endif - enddo ! of do iq=1,nqmx - - ! check that we identified all tracers: - if (count.ne.nqmx) then - write(*,*) "inichim_newstart: found only ",count," tracers" - write(*,*) " expected ",nqmx - do iq=1,count - write(*,*)' ',iq,' ',trim(noms(iq)) - enddo - stop - else - write(*,*)"inichim_newstart: found all expected tracers, namely:" - do iq=1,nqmx - write(*,*)' ',iq,' ',trim(noms(iq)) - enddo - endif - - ! build local chemistry tracer index array nqchem(:) - ! as in the old days, water vapour is last and water ice, - ! if present, is just before water vapour - do iq=1,16 ! loop so as to avoid out-of-bounds on array - if (iq==1) nqchem(i)=igcm_co2 - if (iq==2) nqchem(i)=igcm_co - if (iq==3) nqchem(i)=igcm_o - if (iq==4) nqchem(i)=igcm_o1d - if (iq==5) nqchem(i)=igcm_o2 - if (iq==6) nqchem(i)=igcm_o3 - if (iq==7) nqchem(i)=igcm_h - if (iq==8) nqchem(i)=igcm_h2 - if (iq==9) nqchem(i)=igcm_oh - if (iq==10) nqchem(i)=igcm_ho2 - if (iq==11) nqchem(i)=igcm_h2o2 - if (iq==12) nqchem(i)=igcm_ch4 - if (iq==13) nqchem(i)=igcm_n2 - if (iq==14) nqchem(i)=igcm_ar - if (iq==15) nqchem(i)=igcm_h2o_ice - if (iq==16) nqchem(i)=igcm_h2o_vap - enddo - -! Load in chemistry data for initialization: -c---------------------------------------------------------------------- -c Order of Major species in file: -c -c 1=CO2 -c 2=AR -c 3=N2 -c 4=O2 -c 5=CO -c 6=O -c 7=H2 -c -c Order of Minor species in files: -c -c 1=H -c 2=OH -c 3=HO2 -c 4=H2O -c 5=H2O2 -c 6=O1D -c 7=O3 -c -c---------------------------------------------------------------------- -c Major species: - aa(igcm_co2) = 1 - aa(igcm_ar) = 2 - aa(igcm_n2) = 3 - aa(igcm_o2) = 4 - aa(igcm_co) = 5 - aa(igcm_o) = 6 - aa(igcm_h2) = 7 -c Minor species: - aa(igcm_h) = 8 - aa(igcm_oh) = 9 - aa(igcm_ho2) = 10 - aa(igcm_h2o_vap) = 11 - aa(igcm_h2o2) = 12 - aa(igcm_o1d) = 13 - aa(igcm_o3) = 14 -c---------------------------------------------------------------------- - open(210, iostat=ierr,file= - & trim(datafile)//'/atmosfera_LMD_may.dat') - if (ierr.ne.0) then - write(*,*)'Error : cannot open file atmosfera_LMD_may.dat ' - write(*,*)'(in aeronomars/inichim.F)' - write(*,*)'It should be in :', trim(datafile),'/' - write(*,*)'1) You can change this path in callphys.def with' - write(*,*)' datadir=/path/to/datafiles/' - write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)' - write(*,*)' can be obtained online on:' - write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile' - STOP - endif - open(220, iostat=ierr,file= - & trim(datafile)//'/atmosfera_LMD_min.dat') - if (ierr.ne.0) then - write(*,*)'Error : cannot open file atmosfera_LMD_min.dat ' - write(*,*)'(in aeronomars/inichim.F)' - write(*,*)'It should be in :', trim(datafile),'/' - write(*,*)'1) You can change this path in callphys.def with' - write(*,*)' datadir=/path/to/datafiles/' - write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)' - write(*,*)' can be obtained online on:' - write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile' - STOP - endif - - read(210,*) tmp - read(220,*) tmp - do l=1,252 - read(210,112) tfile(l),zfile(l),nxf(l),(vmrf(l,n),n=1,7) - zfile(l)=zfile(l)*100. !pressure (Pa) - sigfile(l)=zfile(l)/zfile(1) - enddo - do l=1,252 - read(220,113)zfilemin(l),(vmrf(l,n),n=8,14) - zfilemin(l)=zfilemin(l)*1000. !height (m) - do n=1,14 - nf(l,n)=vmrf(l,n)*nxf(l) - enddo - enddo - close(210) - close(220) - -c flag thermo set to 1 or 0 in newstart -c inil=33 for initialisation of thermosphere only -c inil=1 for initialisation of all layers - if (thermo.eq.1) then - inil=33 - else - inil=1 - endif - -! Initialization of tracers - - do i=1,iip1 - do j=1,jjp1 - do l=inil,llm - -c zgcm=sig(l) - zgcm=sig(l)*ps(i,j) - densconc=0. - nt=0. - - do n=1,14 - call intrplf(zgcm,vmrint,zfile,nf(1,n),252) -c call intrplf(zgcm,vmrint,sigfile,nf(1,n),252) - ni(n)=vmrint - - densconc=densconc+ni(n)*mmolf(n) - nt=nt+ni(n) - enddo - - mmean=densconc/nt ! in amu - - if (nq.ne.nbqchem) then ! don't initialize water vapour - do n=1,nq-1 - pq(i,j,l,nqchem(n))= - & ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean - if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n)) - enddo - if (water .and. (nq.gt.nbqchem-2)) then - pq(i,j,l,nqchem(nq)) = 0. - if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) - else - pq(i,j,l,nqchem(nq))= - & ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean - if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) - endif - endif ! of if (nq.ne.nbqchem) - - if (nq.eq.nbqchem) then ! also initialize water vapour - do n=1,nq-2 - pq(i,j,l,nqchem(n))= - & ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean - if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n)) - enddo - if (water) then - pq(i,j,l,nqchem(nq-1)) = 0. - if(i.eq.iip1) then - pq(i,j,l,nqchem(nq-1))=pq(1,j,l,nqchem(nq-1)) - endif - pq(i,j,l,nqchem(nq))= - & ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean - if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) - else - do n=nq-1,nq - pq(i,j,l,nqchem(nq))= - & ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean - if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) - enddo - endif ! of if (water) - endif ! of if (nq.eq.nqmx) - enddo !nlayer - enddo !ngrid_j - enddo !ngrid_i - -cccccccccccccccccccccccccccccc -c make sure that sum of q = 1 -c dominent species is = 1 - sum(all other species) -cccccccccccccccccccccccccccccc -! iqmax=nqchem_min - iqmax=1 - -! if ((nqmx-nqchem_min).gt.10) then - if (nbqchem.gt.10) then - do l=1,llm - do j=1,jjp1 - do i=1,iip1 -! do iq=nqchem_min,nqmx - do iiq=1,nbqchem - iq=nqchem(iiq) - if ( (pq(i,j,l,iq).gt.pq(i,j,l,iqmax)).and. - . (noms(iq).ne."ice") ) then - iqmax=iq - endif - enddo - pq(i,j,l,iqmax)=1. - qtot(i,j,l)=0 -! do iq=nqchem_min,nqmx - do iiq=1,nbqchem - iq=nqchem(iiq) - if ( (iq.ne.iqmax).and. - . (noms(iq).ne."ice") ) then - pq(i,j,l,iqmax)=pq(i,j,l,iqmax)-pq(i,j,l,iq) - endif - enddo !iq -! do iq=nqchem_min,nqmx - do iiq=1,nbqchem - iq=nqchem(iiq) - if (noms(iq).ne."ice") then - qtot(i,j,l)=qtot(i,j,l)+pq(i,j,l,iq) - endif -c if (i.eq.1.and.j.eq.1.and.l.eq.1) write(*,*) 'qtot(i,j,l)', -c $ qtot(i,j,l) - enddo !iq - if (i.eq.1.and.j.eq.1) write(*,*) 'inichim_newstart: ', - $ 'qtot(i,j,l)=',qtot(i,j,l) - enddo !i - enddo !j - enddo !l - endif ! of if (nbqchem.gt.10) -ccccccccccccccccccccccccccccccc - -66 format(i2,6(1x,e11.5)) -112 format(7x,f7.3,3x,e12.6,3x,e12.6,7(3x,e12.6)) -113 format(1x,f6.2,7(3x,e12.6)) - - end Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90 (revision 618) +++ trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90 (revision 618) @@ -0,0 +1,530 @@ + subroutine inichim_newstart(pq, qsurf, ps, flagh2o, flagthermo) + + implicit none + +!======================================================================= +! +! Purpose: +! -------- +! +! Initialization of the chemistry for use in the building of a new start file +! used by program newstart.F +! also used by program testphys1d.F +! +! Authors: +! ------- +! Initial version 11/2002 by Sebastien Lebonnois +! Revised 07/2003 by Monica Angelats-i-Coll to use input files +! Modified 10/2008 Identify tracers by their names Ehouarn Millour +! Modified 11/2011 Addition of methane Franck Lefevre +! Rewritten 04/2012 Franck Lefevre +! +! Arguments: +! ---------- +! +! pq(iip1,jjp1,llm,nqmx) Advected fields, ie chemical species here +! qsurf(ngridmx,nqmx) Amount of tracer on the surface (kg/m2) +! ps(iip1,jjp1) Surface pressure (Pa) +! flagh2o flag for initialisation of h2o (1: yes / 0: no) +! flagthermo flag for initialisation of thermosphere only (1: yes / 0: no) +! +!======================================================================= + +#include "dimensions.h" +#include "dimphys.h" +#include "paramet.h" +#include "tracer.h" +#include "comvert.h" +#include "callkeys.h" +#include "datafile.h" + +! inputs : + + real,intent(in) :: ps(iip1,jjp1) ! surface pressure in the gcm (Pa) + integer,intent(in) :: flagh2o ! flag for h2o initialisation + integer,intent(in) :: flagthermo ! flag for thermosphere initialisation only + +! outputs : + + real,intent(out) :: pq(iip1,jjp1,llm,nqmx) ! advected fields, ie chemical species + real,intent(out) :: qsurf(ngridmx,nqmx) ! surface values (kg/m2) of tracers + +! local : + + integer :: iq, i, j, l, n, nbqchem + integer :: count, ierr, dummy + real :: mmean(iip1,jjp1,llm) ! mean molecular mass (g) + real :: pgcm ! pressure at each layer in the gcm (Pa) + + integer, parameter :: nalt = 252 ! number of altitudes in the initialization files + integer, parameter :: nspe = 14 ! number of species in the initialization files + integer, dimension(nspe) :: niq ! local index of species in initialization files + real, dimension(nalt) :: tinit, zzfile ! temperature in initialization files + real, dimension(nalt) :: pinit ! pressure in initialization files + real, dimension(nalt) :: densinit ! total number density in initialization files + real, dimension(nalt,nspe) :: vmrinit ! mixing ratios in initialization files + real, dimension(nspe) :: vmrint ! mixing ratio interpolated onto the gcm vertical grid + real :: vmr + + character(len=20) :: txt ! to store some text + +! 1. identify tracers by their names: (and set corresponding values of mmol) + +! 1.1 initialize tracer indexes to zero: + + do iq = 1,nqmx + igcm_dustbin(iq) = 0 + end do + + igcm_dust_mass = 0 + igcm_dust_number = 0 + igcm_ccn_mass = 0 + igcm_ccn_number = 0 + igcm_h2o_vap = 0 + igcm_h2o_ice = 0 + igcm_co2 = 0 + igcm_co = 0 + igcm_o = 0 + igcm_o1d = 0 + igcm_o2 = 0 + igcm_o3 = 0 + igcm_h = 0 + igcm_h2 = 0 + igcm_oh = 0 + igcm_ho2 = 0 + igcm_h2o2 = 0 + igcm_ch4 = 0 + igcm_n2 = 0 + igcm_ar = 0 + igcm_ar_n2 = 0 + igcm_co2plus = 0 + igcm_oplus = 0 + igcm_o2plus = 0 + igcm_coplus = 0 + igcm_cplus = 0 + igcm_nplus = 0 + igcm_noplus = 0 + igcm_n2plus = 0 + igcm_hplus = 0 + igcm_elec = 0 + +! 1.2 find dust tracers + + count = 0 + + if (dustbin > 0) then + do iq = 1,nqmx + txt = " " + write(txt,'(a4,i2.2)') 'dust', count + 1 + if (noms(iq) == txt) then + count = count + 1 + igcm_dustbin(count) = iq + mmol(iq) = 100. + end if + end do !do iq=1,nqmx + end if ! of if (dustbin.gt.0) + + if (doubleq) then + do iq = 1,nqmx + if (noms(iq) == "dust_mass") then + igcm_dust_mass = iq + count = count + 1 + end if + if (noms(iq) == "dust_number") then + igcm_dust_number = iq + count = count + 1 + end if + end do + end if ! of if (doubleq) + + if (scavenging) then + do iq = 1,nqmx + if (noms(iq) == "ccn_mass") then + igcm_ccn_mass = iq + count = count + 1 + end if + if (noms(iq) == "ccn_number") then + igcm_ccn_number = iq + count = count + 1 + end if + end do + end if ! of if (scavenging) + + if (submicron) then + do iq=1,nqmx + if (noms(iq) == "dust_submicron") then + igcm_dust_submicron = iq + mmol(iq) = 100. + count = count + 1 + end if + end do + end if ! of if (submicron) + +! 1.3 find chemistry and water tracers + + nbqchem = 0 + do iq = 1,nqmx + if (noms(iq) == "co2") then + igcm_co2 = iq + mmol(igcm_co2) = 44. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "co") then + igcm_co = iq + mmol(igcm_co) = 28. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "o") then + igcm_o = iq + mmol(igcm_o) = 16. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "o1d") then + igcm_o1d = iq + mmol(igcm_o1d) = 16. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "o2") then + igcm_o2 = iq + mmol(igcm_o2) = 32. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "o3") then + igcm_o3 = iq + mmol(igcm_o3) = 48. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "h") then + igcm_h = iq + mmol(igcm_h) = 1. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "h2") then + igcm_h2 = iq + mmol(igcm_h2) = 2. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "oh") then + igcm_oh = iq + mmol(igcm_oh) = 17. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "ho2") then + igcm_ho2 = iq + mmol(igcm_ho2) = 33. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "h2o2") then + igcm_h2o2 = iq + mmol(igcm_h2o2) = 34. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "ch4") then + igcm_ch4 = iq + mmol(igcm_ch4) = 16. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "n2") then + igcm_n2 = iq + mmol(igcm_n2) = 28. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "n") then + igcm_n = iq + mmol(igcm_n) = 14. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "n2d") then + igcm_n2d = iq + mmol(igcm_n2d) = 14. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "no") then + igcm_no = iq + mmol(igcm_no) = 30. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "no2") then + igcm_no2 = iq + mmol(igcm_no2) = 46. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "ar") then + igcm_ar = iq + mmol(igcm_ar) = 40. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "h2o_vap") then + igcm_h2o_vap = iq + mmol(igcm_h2o_vap) = 18. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "h2o_ice") then + igcm_h2o_ice = iq + mmol(igcm_h2o_ice) = 18. + count = count + 1 + nbqchem = nbqchem + 1 + end if + +! 1.4 find ions + + if (noms(iq) == "co2plus") then + igcm_co2plus = iq + mmol(igcm_co2plus) = 44. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "oplus") then + igcm_oplus = iq + mmol(igcm_oplus) = 16. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "o2plus") then + igcm_o2plus = iq + mmol(igcm_o2plus) = 32. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "coplus") then + igcm_coplus = iq + mmol(igcm_coplus) = 28. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "cplus") then + igcm_cplus = iq + mmol(igcm_cplus) = 12. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "nplus") then + igcm_nplus = iq + mmol(igcm_nplus) = 14. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "noplus") then + igcm_noplus = iq + mmol(igcm_noplus) = 30. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "n2plus") then + igcm_n2plus = iq + mmol(igcm_n2plus) = 28. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "hplus") then + igcm_hplus = iq + mmol(igcm_hplus) = 1. + count = count + 1 + nbqchem = nbqchem + 1 + end if + if (noms(iq) == "elec") then + igcm_elec = iq + mmol(igcm_elec) = 1./1822.89 + count = count + 1 + end if + +! 1.5 find idealized non-condensible tracer + + if (noms(iq) == "Ar_N2") then + igcm_ar_n2 = iq + mmol(igcm_ar_n2) = 30. + count = count + 1 + end if + + end do ! of do iq=1,nqmx + +! 1.6 check that we identified all tracers: + + if (count /= nqmx) then + write(*,*) "inichim_newstart: found only ",count," tracers" + write(*,*) " expected ",nqmx + do iq = 1,count + write(*,*) ' ', iq, ' ', trim(noms(iq)) + end do + stop + else + write(*,*) "inichim_newstart: found all expected tracers" + do iq = 1,nqmx + write(*,*) ' ', iq, ' ', trim(noms(iq)) + end do + end if + +! 2. load in chemistry data for initialization: + +! order of major species in initialization file: +! +! 1: co2 +! 2: ar +! 3: n2 +! 4: o2 +! 5: co +! 6: o +! 7: h2 +! +! order of minor species in initialization file: +! +! 1: h +! 2: oh +! 3: ho2 +! 4: h2o +! 5: h2o2 +! 6: o1d +! 7: o3 + +! major species: + + niq(1) = igcm_co2 + niq(2) = igcm_ar + niq(3) = igcm_n2 + niq(4) = igcm_o2 + niq(5) = igcm_co + niq(6) = igcm_o + niq(7) = igcm_h2 + +! minor species: + + niq(8) = igcm_h + niq(9) = igcm_oh + niq(10) = igcm_ho2 + niq(11) = igcm_h2o_vap + niq(12) = igcm_h2o2 + niq(13) = igcm_o1d + niq(14) = igcm_o3 + +! 2.1 open initialization files + + open(210, iostat=ierr,file=trim(datafile)// & + '/atmosfera_LMD_may.dat') + if (ierr /= 0) then + write(*,*)'Error : cannot open file atmosfera_LMD_may.dat ' + write(*,*)'(in aeronomars/inichim.F)' + write(*,*)'It should be in :', trim(datafile),'/' + write(*,*)'1) You can change this path in callphys.def with' + write(*,*)' datadir=/path/to/datafiles/' + write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)' + write(*,*)' can be obtained online on:' + write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile' + stop + end if + open(220, iostat=ierr,file=trim(datafile)// & + '/atmosfera_LMD_min.dat') + if (ierr /= 0) then + write(*,*)'Error : cannot open file atmosfera_LMD_min.dat ' + write(*,*)'(in aeronomars/inichim.F)' + write(*,*)'It should be in :', trim(datafile),'/' + write(*,*)'1) You can change this path in callphys.def with' + write(*,*)' datadir=/path/to/datafiles/' + write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)' + write(*,*)' can be obtained online on:' + write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile' + stop + end if + +! 2.2 read initialization files + +! major species + + read(210,*) + do l = 1,nalt + read(210,*) dummy, tinit(l), pinit(l), densinit(l), & + (vmrinit(l,n), n = 1,7) + pinit(l) = pinit(l)*100. ! conversion in Pa + pinit(l) = log(pinit(l)) ! for the vertical interpolation + end do + close(210) + +! minor species + + read(220,*) + do l = 1,nalt + read(220,*) dummy, (vmrinit(l,n), n = 8,14) + end do + close(220) + +! 3. initialization of tracers + + do i = 1,iip1 + do j = 1,jjp1 + do l = 1,llm + + pgcm = aps(l) + bps(l)*ps(i,j) ! gcm pressure + pgcm = log(pgcm) ! for the vertical interpolation + mmean(i,j,l) = 0. + +! 3.1 vertical interpolation + + do n = 1,nspe + call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt) + vmrint(n) = vmr + iq = niq(n) + mmean(i,j,l) = mmean(i,j,l) + vmrint(n)*mmol(iq) + end do + +! 3.2 attribute mixing ratio: - all layers or only thermosphere +! - with our without h2o + + if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 1.e-3)) then + do n = 1,nspe + iq = niq(n) + if (iq /= igcm_h2o_vap .or. flagh2o == 1) then + pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l) + end if + end do + end if + + end do + end do + end do + + ! set surface values of chemistry tracers to zero + if (flagthermo == 0) then + ! NB: no problem for "surface water vapour" tracer which is always 0 + do n=1,nspe + iq=niq(n) + qsurf(1:ngridmx,iq)=0 + enddo + endif + + +! 3.3 initialization of tracers not contained in the initialization files + +! methane : 10 ppbv + + if (igcm_ch4 /= 0) then + vmr = 10.e-9 + do i = 1,iip1 + do j = 1,jjp1 + do l = 1,llm + pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l) + end do + end do + end do + ! set surface value to zero + qsurf(1:ngridmx,igcm_ch4)=0 + end if + + end Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_readcallphys.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/inichim_readcallphys.F (revision 616) +++ (revision ) @@ -1,520 +1,0 @@ - SUBROUTINE inichim_readcallphys -! -!======================================================================= -! -! subject: -! -------- -! -! Initialisation for the physical parametrisations of the LMD -! martian atmospheric general circulation model. -! This routine is called by inichim_newstart (which is used by programs -! newstart and testphys1d but not the GCM) -! -! author: Frederic Hourdin 15 / 10 /93 -! ------- -! modifications: Sebastien Lebonnois 11/06/2003 (new callphys.def) -! 10/2008 adapted to using tracer sby name. Ehouarn Millour -! 07/2009 use 'getin' to read callphys.def. EM -! -! arguments: -! ---------- -! -! input: -! ------ -! -!======================================================================= -! -!----------------------------------------------------------------------- -! declarations: -! ------------- -! to use 'getin' - USE ioipsl_getincom - IMPLICIT NONE -#include "dimensions.h" -#include "dimphys.h" -#include "planete.h" -#include "comcstfi.h" -#include "comsaison.h" -#include "comdiurn.h" -#include "comgeomfi.h" -#include "callkeys.h" -#include "surfdat.h" - - character*12 ch1 - integer ierr - logical chem, h2o - - -! -------------------------------------------------------------- -! Reading the "callphys.def" file controlling some key options -! -------------------------------------------------------------- - - ! check that 'callphys.def' file is around - OPEN(99,file='callphys.def',status='old',form='formatted' - & ,iostat=ierr) - CLOSE(99) - - IF(ierr.EQ.0) THEN - PRINT* - PRINT* - PRINT*,'--------------------------------------------' - PRINT*,' inichim_readcallphys: Parametres pour la physique', - & ' (callphys.def)' - PRINT*,'--------------------------------------------' - - - write(*,*) "Run with or without tracer transport ?" - tracer=.false. ! default value - call getin("tracer",tracer) - write(*,*) " tracer = ",tracer - - write(*,*) "Diurnal cycle ?" - write(*,*) "(if diurnal=False, diurnal averaged solar heating)" - diurnal=.true. ! default value - call getin("diurnal",diurnal) - write(*,*) " diurnal = ",diurnal - - write(*,*) "Seasonal cycle ?" - write(*,*) "(if season=False, Ls stays constant, to value ", - & "set in 'start'" - season=.true. ! default value - call getin("season",season) - write(*,*) " season = ",season - - write(*,*) "Write some extra output to the screen ?" - lwrite=.false. ! default value - call getin("lwrite",lwrite) - write(*,*) " lwrite = ",lwrite - - write(*,*) "Save statistics in file stats.nc ?" - callstats=.true. ! default value - call getin("callstats",callstats) - write(*,*) " callstats = ",callstats - - write(*,*) "Save EOF profiles in file 'profiles' for ", - & "Climate Database?" - calleofdump=.false. ! default value - call getin("calleofdump",calleofdump) - write(*,*) " calleofdump = ",calleofdump - - write(*,*) "Dust scenario:" - iaervar=3 ! default value - call getin("iaervar",iaervar) - write(*,*) " iaervar = ",iaervar - - write(*,*) "Dust vertical distribution:" - write(*,*) "(=1 Dust opt.deph read in startfi;", - & " =2 Viking scenario; =3 MGS scenario,", - & " =4 Mars Year 24 from TES assimilation)" - iddist=3 ! default value - call getin("iddist",iddist) - write(*,*) " iddist = ",iddist - - write(*,*) "Dust top altitude (km). (Matters only if iddist=1)" - topdustref= 90.0 ! default value - call getin("topdustref",topdustref) - write(*,*) " topdustref = ",topdustref - - - write(*,*) "call radiative transfer ?" - callrad=.true. ! default value - call getin("callrad",callrad) - write(*,*) " callrad = ",callrad - - write(*,*) "call NLTE radiative schemes ?", - & "(matters only if callrad=T)" - callnlte=.false. ! default value - call getin("callnlte",callnlte) - write(*,*) " callnlte = ",callnlte - - write(*,*) "call CO2 NIR absorption ?", - & "(matters only if callrad=T)" - callnirco2=.false. ! default value - call getin("callnirco2",callnirco2) - write(*,*) " callnirco2 = ",callnirco2 - - write(*,*) "call turbulent vertical diffusion ?" - calldifv=.true. ! default value - call getin("calldifv",calldifv) - write(*,*) " calldifv = ",calldifv - - write(*,*) "call convective adjustment ?" - calladj=.true. ! default value - call getin("calladj",calladj) - write(*,*) " calladj = ",calladj - - - write(*,*) "call CO2 condensation ?" - callcond=.true. ! default value - call getin("callcond",callcond) - write(*,*) " callcond = ",callcond - - write(*,*)"call thermal conduction in the soil ?" - callsoil=.true. ! default value - call getin("callsoil",callsoil) - write(*,*) " callsoil = ",callsoil - - - write(*,*)"call Lott's gravity wave/subgrid topography ", - & "scheme ?" - calllott=.true. ! default value - call getin("calllott",calllott) - write(*,*)" calllott = ",calllott - - - write(*,*)"rad.transfer is computed every iradia", - & " physical timestep" - iradia=1 ! default value - call getin("iradia",iradia) - write(*,*)" iradia = ",iradia - - - write(*,*)"Output of the exchange coefficient mattrix ?", - & "(for diagnostics only)" - callg2d=.false. ! default value - call getin("callg2d",callg2d) - write(*,*)" callg2d = ",callg2d - - write(*,*)"Rayleigh scattering : (should be .false. for now)" - rayleigh=.false. - call getin("rayleigh",rayleigh) - write(*,*)" rayleigh = ",rayleigh - - -! TRACERS: - - write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)" - dustbin=0 ! default value - call getin("dustbin",dustbin) - write(*,*)" dustbin = ",dustbin - - write(*,*)"Radiatively active dust ? (matters if dustbin>0)" - active=.false. ! default value - call getin("active",active) - write(*,*)" active = ",active - -! Test of incompatibility: -! if active is used, then dustbin should be > 0 - - if (active.and.(dustbin.lt.1)) then - print*,'if active is used, then dustbin should > 0' - stop - endif - - write(*,*)"use mass and number mixing ratios to predict", - & " dust size ?" - doubleq=.false. ! default value - call getin("doubleq",doubleq) - write(*,*)" doubleq = ",doubleq - -! Test of incompatibility: -! if doubleq is used, then dustbin should be at least 2 - - if (doubleq.and.(dustbin.lt.2)) then - print*,'if doubleq is used, then dustbin should be > 1' - stop - endif - - write(*,*)"dust lifted by GCM surface winds ?" - lifting=.false. ! default value - call getin("lifting",lifting) - write(*,*)" lifting = ",lifting - -! Test of incompatibility: -! if lifting is used, then dustbin should be > 0 - - if (lifting.and.(dustbin.lt.1)) then - print*,'if lifting is used, then dustbin should > 0' - stop - endif - - write(*,*)" dust lifted by dust devils ?" - callddevil=.false. !default value - call getin("callddevil",callddevil) - write(*,*)" callddevil = ",callddevil - - -! Test of incompatibility: -! if dustdevil is used, then dustbin should be > 0 - - if (callddevil.and.(dustbin.lt.1)) then - print*,'if dustdevil is used, then dustbin should > 0' - stop - endif - - write(*,*)"Dust scavenging by CO2 snowfall ?" - scavenging=.false. ! default value - call getin("scavenging",scavenging) - write(*,*)" scavenging = ",scavenging - - -! Test of incompatibility: -! if scavenging is used, then dustbin should be > 0 - - if (scavenging.and.(dustbin.lt.1)) then - print*,'if scavenging is used, then dustbin should > 0' - stop - endif - - write(*,*) "Gravitationnal sedimentation ?" - sedimentation=.true. ! default value - call getin("sedimentation",sedimentation) - write(*,*) " sedimentation = ",sedimentation - - write(*,*) "includes water ice", - & "(if true, 'water' must also be .true.)" - - write(*,*) "Radiatively active transported atmospheric ", - & "water ice ?" - activice=.false. ! default value - call getin("activice",activice) - write(*,*) " activice = ",activice - - write(*,*) "Compute water cycle ?" - water=.false. ! default value - call getin("water",water) - write(*,*) " water = ",water - write(*,*) "Permanent water caps at poles ?", - & " .true. is RECOMMENDED" - write(*,*) "(with .true., North cap is a source of water ", - & "and South pole is a cold trap)" - caps=.true. ! default value - call getin("caps",caps) - write(*,*) " caps = ",caps - -! Test of incompatibility: -! if activice is used, then water should be used too - - if (activice.and..not.water) then - print*,'if activice is used, water should be used too' - stop - endif - - write(*,*) "photochemistry: include chemical species" - photochem=.false. ! default value - call getin("photochem",photochem) - write(*,*) " photochem = ",photochem - - -! THERMOSPHERE - - write(*,*) "call thermosphere ?" - callthermos=.false. ! default value - call getin("callthermos",callthermos) - write(*,*) " callthermos = ",callthermos - - - write(*,*) " water included without cycle ", - & "(only if water=.false.)" - thermoswater=.false. ! default value - call getin("thermoswater",thermoswater) - write(*,*) " thermoswater = ",thermoswater - - write(*,*) "call thermal conduction ?", - & " (only if callthermos=.true.)" - callconduct=.false. ! default value - call getin("callconduct",callconduct) - write(*,*) " callconduct = ",callconduct - - write(*,*) "call EUV heating ?", - & " (only if callthermos=.true.)" - calleuv=.false. ! default value - call getin("calleuv",calleuv) - write(*,*) " calleuv = ",calleuv - - write(*,*) "call molecular viscosity ?", - & " (only if callthermos=.true.)" - callmolvis=.false. ! default value - call getin("callmolvis",callmolvis) - write(*,*) " callmolvis = ",callmolvis - - write(*,*) "call molecular diffusion ?", - & " (only if callthermos=.true.)" - callmoldiff=.false. ! default value - call getin("callmoldiff",callmoldiff) - write(*,*) " callmoldiff = ",callmoldiff - - - write(*,*) "call thermospheric photochemistry ?", - & " (only if callthermos=.true.)" - thermochem=.false. ! default value - call getin("thermochem",thermochem) - write(*,*) " thermochem = ",thermochem - - write(*,*) "date for solar flux calculation:", - & " (1985 < date < 2002)" - write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)" - solarcondate=1993.4 ! default value - call getin("solarcondate",solarcondate) - write(*,*) " solarcondate = ",solarcondate - - - if (.not.callthermos) then - if (thermoswater) then - print*,'if thermoswater is set, callthermos must be true' - stop - endif - if (callconduct) then - print*,'if callconduct is set, callthermos must be true' - stop - endif - if (calleuv) then - print*,'if calleuv is set, callthermos must be true' - stop - endif - if (callmolvis) then - print*,'if callmolvis is set, callthermos must be true' - stop - endif - if (callmoldiff) then - print*,'if callmoldiff is set, callthermos must be true' - stop - endif - if (thermochem) then - print*,'if thermochem is set, callthermos must be true' - stop - endif - endif - -! Test of incompatibility: -! if photochem is used, then water should be used too - - if (photochem.and..not.water) then - print*,'if photochem is used, water should be used too' - stop - endif - -! if callthermos is used, then thermoswater should be used too -! (if water not used already) - - if (callthermos .and. .not.water) then - if (callthermos .and. .not.thermoswater) then - print*,'if callthermos is used, water or thermoswater - & should be used too' - stop - endif - endif - - PRINT*,'--------------------------------------------' - PRINT* - PRINT* - ELSE - write(*,*) - write(*,*) 'Cannot read file callphys.def. Is it here ?' - stop - ENDIF - CLOSE(99) - - pi=2.*asin(1.) - -! managing the tracers, and tests: -! ------------------------------- - - if(tracer) then - -! when photochem is used, nqchem_min is the rank -! of the first chemical species - -! Ehouarn: nqchem_min is now meaningless and no longer used -! nqchem_min = 1 - if (photochem .or. callthermos) then - chem = .true. -! if (doubleq) then -! nqchem_min = 3 -! else -! nqchem_min = dustbin+1 -! end if - end if - - if (water .or. thermoswater) h2o = .true. - -c TESTS - - print*,'inichim_readcallphys: TRACERS:' - - if ((doubleq).and.(h2o).and. - $ (chem)) then -! print*,' 1: dust ; 2: dust (doubleq)' -! print*,' 3 to ',nqmx-2,': chemistry' -! print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - print*,' 2 dust tracers (doubleq)' - print*,' 1 water vapour tracer' - print*,' 1 water ice tracer' - print*,nqmx-4,' chemistry tracers' - endif - - if ((doubleq).and.(h2o).and. - $ .not.(chem)) then -! print*,' 1: dust ; 2: dust (doubleq)' -! print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - print*,' 2 dust tracers (doubleq)' - print*,' 1 water vapour tracer' - print*,' 1 water ice tracer' - if (nqmx.lt.4) then - print*,'nqmx should be at least 4 with these options.' - print*,'(or check callphys.def)' - stop - endif - endif - -! if ((doubleq).and..not.(h2o)) then -! print*,' 1: dust ; 2: dust (doubleq)' -! print*,' 2 dust tracers (doubleq)' -! if (nqmx.ne.2) then -! print*,'nqmx should be 2 with these options...' -! print*,'(or check callphys.def)' -! stop -! endif -! endif - -! if (.not.(doubleq).and.(h2o).and. -! $ (chem)) then -! if (dustbin.gt.0) then -! print*,' 1 to ',dustbin,': dust bins' -! print*,dustbin,' dust bins' -! endif -! print*,nqchem_min,' to ',nqmx-2,': chemistry' -! print*,nqmx-1,': water ice ; ',nqmx,': water vapor' -! print*,nqmx-2-dustbin,' chemistry tracers' -! print*,' 1 water vapour tracer' -! print*,' 1 water ice tracer' -! endif - -! if (.not.(doubleq).and.(h2o).and. -! $ .not.(chem)) then -! if (dustbin.gt.0) then -! print*,' 1 to ',dustbin,': dust bins' -! print*,dustbin,' dust bins' -! endif -! print*,nqmx-1,': water ice ; ',nqmx,': water vapor' -! print*,' 1 water vapour tracer' -! print*,' 1 water ice tracer' -! if (nqmx.ne.(dustbin+2)) then -! print*,'nqmx should be ',(dustbin+2), -! $ ' with these options...' -! print*,'(or check callphys.def)' -! stop -! endif -! endif - -! if (.not.(doubleq).and..not.(h2o)) then -! if (dustbin.gt.0) then -! print*,' 1 to ',dustbin,': dust bins' -! print*,dustbin,' dust bins' -! if (nqmx.ne.dustbin) then -! print*,'nqmx should be ',dustbin, -! $ ' with these options...' -! print*,'(or check callphys.def)' -! stop -! endif -! else -! print*,'dustbin=',dustbin, -! $ ': tracer should be F with these options...' -! $ ,'UNLESS you just want to move tracers around ' -! endif -! endif - - endif ! of if (tracer) - - RETURN - END Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F (revision 616) +++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F (revision 618) @@ -60,5 +60,6 @@ parameter (nesp = 16) ! number of species in the chemistry c - real ctimestep ! chemistry timestep (s) + real stimestep ! standard timestep for the chemistry (s) + real ctimestep ! real timestep for the chemistry (s) real rm(nlayermx,nesp) ! species volume mixing ratio real j(nlayermx,nd) ! interpolated photolysis rates (s-1) @@ -85,9 +86,12 @@ c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c ctimestep : chemistry timestep (s) c -c 1/3 of physical timestep c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c - phychemrat = 3 +c stimestep : standard timestep for the chemistry (s) c +c ctimestep : real timestep for the chemistry (s) c +c phychemrat : ratio physical/chemical timestep c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + stimestep = 600. ! standard value : 10 mn +c + phychemrat = nint(ptimestep/stimestep) c ctimestep = ptimestep/real(phychemrat) @@ -1175,5 +1179,4 @@ c integer l,iq - logical,save :: firstcall = .true. c tracer indexes in the chemistry: @@ -1196,73 +1199,4 @@ integer,parameter :: i_ox = 16 c -c first call initializations -c - if (firstcall) then - -c identify the indexes of the tracers we need - - if (igcm_co2.eq.0) then - write(*,*) "concentrations: Error; no CO2 tracer !!!" - stop - endif - if (igcm_co.eq.0) then - write(*,*) "concentrations: Error; no CO tracer !!!" - stop - endif - if (igcm_o.eq.0) then - write(*,*) "concentrations: Error; no O tracer !!!" - stop - endif - if (igcm_o1d.eq.0) then - write(*,*) "concentrations: Error; no O1D tracer !!!" - stop - endif - if (igcm_o2.eq.0) then - write(*,*) "concentrations: Error; no O2 tracer !!!" - stop - endif - if (igcm_o3.eq.0) then - write(*,*) "concentrations: Error; no O3 tracer !!!" - stop - endif - if (igcm_h.eq.0) then - write(*,*) "concentrations: Error; no H tracer !!!" - stop - endif - if (igcm_h2.eq.0) then - write(*,*) "concentrations: Error; no H2 tracer !!!" - stop - endif - if (igcm_oh.eq.0) then - write(*,*) "concentrations: Error; no OH tracer !!!" - stop - endif - if (igcm_ho2.eq.0) then - write(*,*) "concentrations: Error; no HO2 tracer !!!" - stop - endif - if (igcm_h2o2.eq.0) then - write(*,*) "concentrations: Error; no H2O2 tracer !!!" - stop - endif - if (igcm_ch4.eq.0) then - write(*,*) "concentrations: Error; no CH4 tracer !!!" - stop - endif - if (igcm_n2.eq.0) then - write(*,*) "concentrations: Error; no N2 tracer !!!" - stop - endif - if (igcm_ar.eq.0) then - write(*,*) "concentrations: Error; no AR tracer !!!" - stop - endif - if (igcm_h2o_vap.eq.0) then - write(*,*) "concentrations: Error; no water vapor tracer !!!" - stop - endif - firstcall = .false. - end if ! of if (firstcall) -c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c initialise chemical species @@ -1281,8 +1215,17 @@ rm(l,i_ho2) = max(zycol(l, igcm_ho2), 1.e-30) rm(l,i_h2o2) = max(zycol(l, igcm_h2o2), 1.e-30) - rm(l,i_ch4) = max(zycol(l, igcm_ch4), 1.e-30) rm(l,i_n2) = max(zycol(l, igcm_n2), 1.e-30) rm(l,i_h2o) = max(zycol(l, igcm_h2o_vap), 1.e-30) end do + + if (igcm_ch4 .eq. 0) then + do l = 1,lswitch-1 + rm(l,i_ch4) = 0. + end do + else + do l = 1,lswitch-1 + rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30) + end do + end if c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -1370,8 +1313,13 @@ zycol(l, igcm_ho2) = rm(l,i_ho2) zycol(l, igcm_h2o2) = rm(l,i_h2o2) - zycol(l, igcm_ch4) = rm(l,i_ch4) zycol(l, igcm_n2) = rm(l,i_n2) zycol(l, igcm_h2o_vap) = rm(l,i_h2o) end do + + if (igcm_ch4 .ne. 0) then + do l = 1,lswitch-1 + zycol(l,igcm_ch4) = rm(l,i_ch4) + end do + end if c return @@ -2000,3 +1948,2 @@ return end - Index: trunk/LMDZ.MARS/libf/dyn3d/newstart.F =================================================================== --- trunk/LMDZ.MARS/libf/dyn3d/newstart.F (revision 616) +++ trunk/LMDZ.MARS/libf/dyn3d/newstart.F (revision 618) @@ -139,5 +139,6 @@ real choix_1 character*80 fichnom - integer Lmodif,iq,thermo + integer Lmodif,iq + integer flagthermo, flagh2o character modif*20 real z_reel(iip1,jjp1) @@ -450,17 +451,12 @@ enddo ! of do iq=1,nqmx + ! initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...) + call initracer(qsurf,co2ice) + if (count.eq.nqmx) then write(*,*) 'Newstart: updating tracer names' - ! do things the easy but dirty way: - ! 1. call inichim_readcallphys (so that callphys.def is read) - call inichim_readcallphys() - ! 2. call initracer to set all new tracer names (in noms(:)) - call initracer(qsurf,co2ice) - ! 3. copy noms(:) to tnom(:) + ! copy noms(:) to tnom(:) to have matching tracer names in physics + ! and dynamics tnom(1:nqmx)=noms(1:nqmx) - write(*,*) 'Newstart: updated tracer names' - else - ! initialize tracer names and indexes (igcm_co2, igcm_h2o_vap, ...) - call initracer(qsurf,co2ice) endif @@ -474,36 +470,37 @@ write(*,*) write(*,*) 'List of possible changes :' - write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) - write(*,*) 'flat : no topography ("aquaplanet")' - write(*,*) 'bilball : uniform albedo and thermal inertia' - write(*,*) 'z0 : set a uniform surface roughness length' - write(*,*) 'coldspole : cold subsurface and high albedo at S.pole' - write(*,*) 'qname : change tracer name' - write(*,*) 'q=0 : ALL tracer =zero' - write(*,*) 'q=x : give a specific uniform value to one tracer' - write(*,*) 'q=profile : specify a profile for a tracer' - write(*,*) 'ini_q : tracers initialisation for chemistry, water an - $d ice ' - write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and - $ice ' - write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on - $ly ' - write(*,*) 'ini_h2osurf : reinitialize surface water ice ' - write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' - write(*,*) 'watercapn : H20 ice on permanent N polar cap ' - write(*,*) 'watercaps : H20 ice on permanent S polar cap ' - write(*,*) 'wetstart : start with a wet atmosphere' - write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero' - write(*,*) 'co2ice=0 : remove CO2 polar cap' - write(*,*) 'ptot : change total pressure' - write(*,*) 'therm_ini_s: Set soil thermal inertia to reference sur - &face values' - write(*,*) 'subsoilice_n: Put deep underground ice layer in northe - &rn hemisphere' - write(*,*) 'subsoilice_s: Put deep underground ice layer in southe - &rn hemisphere' - write(*,*) 'mons_ice: Put underground ice layer according to MONS- - &derived data' + write(*,*) 'flat : no topography ("aquaplanet")' + write(*,*) 'bilball : uniform albedo and thermal inertia' + write(*,*) 'z0 : set a uniform surface roughness length' + write(*,*) 'coldspole : cold subsurface and high albedo at + $ S.Pole' + write(*,*) 'qname : change tracer name' + write(*,*) 'q=0 : ALL tracer =zero' + write(*,*) 'q=x : give a specific uniform value to one + $ tracer' + write(*,*) 'q=profile : specify a profile for a tracer' + write(*,*) 'ini_q : tracers initialization for chemistry + $ and water vapour' + write(*,*) 'ini_q-h2o : tracers initialization for chemistry + $ only' + write(*,*) 'ini_h2osurf : reinitialize surface water ice ' + write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' + write(*,*) 'watercapn : H20 ice on permanent N polar cap ' + write(*,*) 'watercaps : H20 ice on permanent S polar cap ' + write(*,*) 'wetstart : start with a wet atmosphere' + write(*,*) 'isotherm : Isothermal Temperatures, wind set to + $ zero' + write(*,*) 'co2ice=0 : remove CO2 polar cap' + write(*,*) 'ptot : change total pressure' + write(*,*) 'therm_ini_s : set soil thermal inertia to reference + $ surface values' + write(*,*) 'subsoilice_n : put deep underground ice layer in + $ northern hemisphere' + write(*,*) 'subsoilice_s : put deep underground ice layer in + $ southern hemisphere' + write(*,*) 'mons_ice : put underground ice layer according + $ to MONS derived data' write(*,*) @@ -834,7 +831,8 @@ c ----------------------------------------------- else if (trim(modif) .eq. 'ini_q') then + flagh2o = 1 + flagthermo = 0 + yes=' ' c For more than 32 layers, possible to initiate thermosphere only - thermo=0 - yes=' ' if (llm.gt.32) then do while ((yes.ne.'y').and.(yes.ne.'n')) @@ -843,30 +841,23 @@ read(*,fmt='(a)') yes if (yes.eq.'y') then - thermo=1 + flagthermo=1 else - thermo=0 + flagthermo=0 endif enddo endif - call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi, - $ thermo,qsurf) - write(*,*) 'Chemical species initialized' - - if (thermo.eq.0) then -c mise a 0 des qsurf (traceurs a la surface) - DO iq =1, nqmx - DO ig=1,ngridmx - qsurf(ig,iq)=0. - ENDDO - ENDDO - endif - -c ini_q-H2O : as above exept for the water vapour tracer + call inichim_newstart(q, qsurf, ps, flagh2o, flagthermo) + write(*,*) 'inichim_newstart: chemical species and + $ water vapour initialised' + + +c ini_q-h2o : as above exept for the water vapour tracer c ------------------------------------------------------ - else if (trim(modif) .eq. 'ini_q-H2O') then + else if (trim(modif) .eq. 'ini_q-h2o') then + flagh2o = 0 + flagthermo = 0 + yes=' ' ! for more than 32 layers, possible to initiate thermosphere only - thermo=0 - yes=' ' if(llm.gt.32) then do while ((yes.ne.'y').and.(yes.ne.'n')) @@ -875,58 +866,14 @@ read(*,fmt='(a)') yes if (yes.eq.'y') then - thermo=1 + flagthermo=1 else - thermo=0 + flagthermo=0 endif enddo endif - call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi, - $ thermo,qsurf) - write(*,*) 'Initialized chem. species exept last (H2O)' - - if (thermo.eq.0) then -c set surface tracers to zero, except water ice - DO iq =1, nqmx - if (iq.ne.igcm_h2o_ice) then - DO ig=1,ngridmx - qsurf(ig,iq)=0. - ENDDO - endif - ENDDO - endif - -c ini_q-iceH2O : as above exept for ice et H2O -c ----------------------------------------------- - else if (trim(modif) .eq. 'ini_q-iceH2O') then -c For more than 32 layers, possible to initiate thermosphere only - thermo=0 - yes=' ' - if(llm.gt.32) then - do while ((yes.ne.'y').and.(yes.ne.'n')) - write(*,*)'', - & 'initialisation for thermosphere only? (y/n)' - read(*,fmt='(a)') yes - if (yes.eq.'y') then - thermo=1 - else - thermo=0 - endif - enddo - endif - - call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi, - $ thermo,qsurf) - write(*,*) 'Initialized chem. species exept ice and H2O' - - if (thermo.eq.0) then -c set surface tracers to zero, except water ice - DO iq =1, nqmx - if (iq.ne.igcm_h2o_ice) then - DO ig=1,ngridmx - qsurf(ig,iq)=0. - ENDDO - endif - ENDDO - endif + call inichim_newstart(q, qsurf, ps, flagh2o, flagthermo) + write(*,*) 'inichim_newstart: chemical species initialised + $ (except water vapour)' + c wetstart : wet atmosphere with a north to south gradient