Context Navigation

← Previous Changeset
Next Changeset →

Changeset 618

Timestamp:

Apr 13, 2012, 8:40:01 AM (13 years ago)

Author:

emillour

Message:

Mars GCM:

Update of the chemistry package:

calchim.F90 :
- upgraded from calchim.F
- can run with or without CH4
- fixed the mass conservation scheme
photochemistry.F : -chemistry timestep is now independant from physics timestep -can run with or without a CH4 tracer
- removed initial tests on species, since these are already done in calchim
removed inichim_readcallphys.F (not used any more)
concentrations.F: adaptated to better handle indexes and molar masses of

tracers. "ncomp" (in chimiedata.h) no longer needed.

inichim_newstart.F90 : rewritten and cleaned (won't be compatible with

old start files where tracer names are q01,...).
Now handles hybrid levels and automaticaly adapts
depending on which tracers are available.

newstart.F : adapted to followup changes in inchim_newstart.F90 and some

cleanup around the initialization of tracer names and surface
values.

FL+EM

Location:

trunk/LMDZ.MARS

Files:

: 1 deleted
: 4 edited
: 2 moved

README (modified) (1 diff)
libf/aeronomars/calchim.F90 (moved) (moved from trunk/LMDZ.MARS/libf/aeronomars/calchim.F) (3 diffs)
libf/aeronomars/concentrations.F (modified) (5 diffs)
libf/aeronomars/inichim_newstart.F90 (moved) (moved from trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F) (1 diff)
libf/aeronomars/inichim_readcallphys.F (deleted)
libf/aeronomars/photochemistry.F (modified) (7 diffs)
libf/dyn3d/newstart.F (modified) (6 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/LMDZ.MARS/README

-                      r608
+                      r618
 == 31/03/12 == AS & FF
 >> Deleted obsolete lines in solarlong
+== 13/04/12 == FL
+>> Update of the chemistry package:
+- calchim.F90 :
+  - upgraded from calchim.F
+  - can run with or without CH4
+  - fixed the mass conservation scheme
+- photochemistry.F :
+  -chemistry timestep is now independant from physics timestep
+  -can run with or without a CH4 tracer
+  - removed initial tests on species, since these are already done in calchim
+- removed inichim_readcallphys.F (not used any more)
+- concentrations.F: adaptated to better handle indexes and molar masses of
+                    tracers. "ncomp" (in chimiedata.h) no longer needed.
+- inichim_newstart.F90 : rewritten and cleaned (won't be compatible with
+                         old start files where tracer names are q01,...).
+                         Now handles hybrid levels and automaticaly adapts
+                         depending on which tracers are available.
+- newstart.F : adapted to followup changes in inchim_newstart.F90 and some
+               cleanup around the initialization of tracer names and surface
+               values.

trunk/LMDZ.MARS/libf/aeronomars/calchim.F90

-                      r616
+                      r618
       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
+      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
+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    --------------
+!=======================================================================
+!
+!   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 "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
+!     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
+c
          if (photochem) then
             print*,'calchim: Read photolysis lookup table'
 …
          ! 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
+         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'
-         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))
+! Initializations
+      zycol(:,:)    = 0.
       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
+!     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
+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)
+               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
+            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.
+c
 c           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
+c
+!           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
+c
 c           search for switch index between regions
+c
+!           search for switch index between regions
             if (photochem .and. thermochem) then
                if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
+               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
                   lswitch = l
                   foundswitch=1
+                  foundswitch = 1
                end if
             end if
             if ( .not. photochem) then
+            if (.not. photochem) then
                lswitch = 22
             end if
             if (.not.  thermochem) then
                lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
+            if (.not. thermochem) then
+               lswitch = min(50,nlayermx+1)
             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
+         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)
+         end if
+c
+c        chemistry in upper atmosphere
+c
+            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
+c
 c        dry deposition
+c
+            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
+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.
+            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 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 l = 1,nlayermx
+!=======================================================================
+!     end of loop over grid
+!=======================================================================
       end do ! of do ig=1,ngridmx
+c
 c=======================================================================
 c     write outputs
 c=======================================================================
+c
+!=======================================================================
+!     write outputs
+!=======================================================================
 ! 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 (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))
+              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 if ! of if (output)
+      return
       end

trunk/LMDZ.MARS/libf/aeronomars/concentrations.F

-                      r370
+                      r618
       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"
 …
 #include "conc.h"
+c    Input/Output
+c    ------------
+!     input/output
       real pplay(ngridmx,nlayermx)
 …
       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
 …
          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)
 …
          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

trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90

-                      r616
+                      r618
+      subroutine inichim_newstart(pq,ps,sig,nq,latfi,lonfi,airefi,
+     $ thermo,qsurf)
+      subroutine inichim_newstart(pq, qsurf, ps, flagh2o, flagthermo)
       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    --------------
+!=======================================================================
+!
+!  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 "chimiedata.h"
 #include "tracer.h"
+#include "comcstfi.h"
+#include "comdiurn.h"
+#include "comvert.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
+! 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
+        ! 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))
+         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
+      ! 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
+      format(i2,6(1x,e11.5))
+     format(7x,f7.3,3x,e12.6,3x,e12.6,7(3x,e12.6))
+     format(1x,f6.2,7(3x,e12.6))
+! 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

trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F

-                      r576
+                      r618
       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)
 …
+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)
 …
+c
       integer l,iq
-      logical,save :: firstcall = .true.
 c     tracer indexes in the chemistry:
 …
       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
 …
          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
 …
          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
 …
       return
       end

trunk/LMDZ.MARS/libf/dyn3d/newstart.F

-                      r563
+                      r618
       real choix_1
       character*80      fichnom
+      integer Lmodif,iq,thermo
+      integer Lmodif,iq
+      integer flagthermo, flagh2o
       character modif*20
       real z_reel(iip1,jjp1)
 …
       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
 …
       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(*,*)
 …
 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'))
 …
             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'))
 …
             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

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 618

Legend:

trunk/LMDZ.MARS/README

trunk/LMDZ.MARS/libf/aeronomars/calchim.F90

trunk/LMDZ.MARS/libf/aeronomars/concentrations.F

trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90

trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F

trunk/LMDZ.MARS/libf/dyn3d/newstart.F

Download in other formats: