Changeset 1672

Timestamp:

Mar 8, 2017, 10:37:31 AM (8 years ago)

Author:

jvatant

Message:

Re-plug chemistry of old Titan
+ minor correction on reading Titan's startfiles and haze routine
JVO

Location:

trunk

Files:

• LMDZ.COMMON/libf/dyn3d/gcm.F90 (modified) (1 diff)
• LMDZ.COMMON/libf/dyn3dpar/gcm.F (modified) (1 diff)
• LMDZ.TITAN/libf/chimtitan/disso.c (modified) (2 diffs)
• LMDZ.TITAN/libf/phytitan/calchim.F90 (added)
• LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (modified) (2 diffs)
• LMDZ.TITAN/libf/phytitan/comcstfi_mod.F90 (modified) (1 diff)
• LMDZ.TITAN/libf/phytitan/disr_haze.F90 (modified) (5 diffs)
• LMDZ.TITAN/libf/phytitan/inicondens.F90 (added)
• LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (modified) (5 diffs)
• LMDZ.TITAN/libf/phytitan/initracer.F (modified) (3 diffs)
• LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (modified) (12 diffs)
• LMDZ.TITAN/libf/phytitan/tracer_h.F90 (modified) (1 diff)

trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90

@@ -372 +372 @@
     write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
     write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-  else if (planet_type.eq.'titan') then
-    jD_ref=1 ! Because we use old Titan starts (J.V.O 2016)
+!  else if (planet_type.eq.'titan') then
+!    jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
   else
 ! A voir pour Titan et Venus

trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F

@@ -396 +396 @@
       write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
       write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-      else if (planet_type.eq.'titan') then
-        jD_ref=1 ! Because we use old Titan starts (J.V.O 2016)
+!      else if (planet_type.eq.'titan') then
+!        jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
       else
 ! A voir pour Titan et Venus

trunk/LMDZ.TITAN/libf/chimtitan/disso.c

@@ -179 +179 @@
    double f,**flux;
    double flact;
-   char   name[60],dir[24];
+   char   name[60],dir[27];
    FILE   *fp,*out;
 
@@ -186 +186 @@
 /* lecture des flux actiniques:
    - suppose que l'executable est dans $LMDGCM/RUN/xxx/
-   - et les moyennes dans $LMDGCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
+   - et les moyennes dans $LMDGCM/datagcm/PHOT(NLAT)/Moy_(lat:1 a NLAT)
 */
-   strcpy( dir, "../../INPUT/PHOT" );
+   strcpy( dir, "../../datagcm/PHOT" );
    if( (*NLAT) < 10 )
    {

trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

@@ -4 +4 @@
       logical,save :: callrad,corrk,calldifv,UseTurbDiff
 !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
-      logical,save :: calladj,callsoil
-!$OMP THREADPRIVATE(calladj,callsoil)
+      logical,save :: calladj,callsoil,callchim
+!$OMP THREADPRIVATE(calladj,callsoil,callchim)
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
@@ -38 +38 @@
       logical,save :: oblate
 !$OMP THREADPRIVATE(nosurf,oblate)
-
+      
+      integer,save :: ichim
       integer,save :: iddist
       integer,save :: iaervar
       integer,save :: iradia
       integer,save :: startype
-!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype)
+!$OMP THREADPRIVATE(ichim,iddist,iaervar,iradia,startype)
 
       real,save :: Fat1AU

trunk/LMDZ.TITAN/libf/phytitan/comcstfi_mod.F90

@@ -11 +11 @@
       REAL,SAVE :: omeg ! planet rotation rate (rad/s)
       REAL,SAVE :: avocado ! something like 6.022e23
-!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado)
+      REAL,SAVE :: kbol ! something like 1.380649e-23
+!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado,kbol)
 
 END MODULE comcstfi_mod

trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90

@@ -59 +59 @@
    press_PL(:)=press_PL(:)*1E-2
    
-   ! Correction of the table according to Lavvas et al. 2010 (ext coeff = constant < 80 km / 20 mbar ) because of condensation
-!   do i=1,49
-!     ext_PL(i,:) = ext_PL(50,:)
-!   enddo
-
    firstcall=.false.
 endif
@@ -101 +96 @@
 ! Lin-Log interpolation : linear on wln, logarithmic on press
 
+if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
 taeros = ( ext_PL(ip,iw)*(1.-factw)   + ext_PL(ip,iw+1)  *factw ) ** (1.-factp) &
         *( ext_PL(ip+1,iw)*(1.-factw) + ext_PL(ip+1,iw+1)*factw ) ** factp
@@ -109 +105 @@
 cbar   = ( asf_PL(ip,iw)*(1.-factw)   + asf_PL(ip,iw+1)  *factw ) ** (1.-factp) &
         *( asf_PL(ip+1,iw)*(1.-factw) + asf_PL(ip+1,iw+1)*factw ) ** factp
+else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then
+taeros =  ext_PL(ip,iw) ** (1.-factp) * ext_PL(ip+1,iw) ** factp
+
+ssa    =  ssa_PL(ip,iw) ** (1.-factp) * ssa_PL(ip+1,iw) ** factp
+
+cbar   =  asf_PL(ip,iw) ** (1.-factp) * asf_PL(ip+1,iw) ** factp
 
 ! In case of vertical extension over the max of the table
@@ -114 +116 @@
 ! Arbitray threshold pressure value, just to deal with the last level press=0
 
-if(ip.eq.nblev_PL) then 
+else if(ip.eq.nblev_PL) then 
 
   tmp_p = press
@@ -122 +124 @@
   endif
 
-  fact_t = log10( ( ext_PL(ip,iw)*(1.-factw) + ext_PL(ip,iw+1)  *factw )   & 
-           / ( ext_PL(ip-5,iw)*(1.-factw) + ext_PL(ip-5,iw+1)  *factw ) )
-  fact_s = log10( ( ssa_PL(ip,iw)*(1.-factw) + ssa_PL(ip,iw+1)  *factw )   & 
-           / ( ssa_PL(ip-5,iw)*(1.-factw) + ssa_PL(ip-5,iw+1)  *factw ) )
-  fact_c = log10( ( asf_PL(ip,iw)*(1.-factw) + asf_PL(ip,iw+1)  *factw )   & 
-           / ( asf_PL(ip-5,iw)*(1.-factw) + asf_PL(ip-5,iw+1)  *factw ) )
+  if (iw.ne.nbwl_PL) then
+    fact_t = log10( ( ext_PL(ip,iw)*(1.-factw) + ext_PL(ip,iw+1)  *factw )   & 
+             / ( ext_PL(ip-5,iw)*(1.-factw) + ext_PL(ip-5,iw+1)  *factw ) )
+    fact_s = log10( ( ssa_PL(ip,iw)*(1.-factw) + ssa_PL(ip,iw+1)  *factw )   & 
+             / ( ssa_PL(ip-5,iw)*(1.-factw) + ssa_PL(ip-5,iw+1)  *factw ) )
+    fact_c = log10( ( asf_PL(ip,iw)*(1.-factw) + asf_PL(ip,iw+1)  *factw )   & 
+             / ( asf_PL(ip-5,iw)*(1.-factw) + asf_PL(ip-5,iw+1)  *factw ) )
+  else
+    fact_t = log10( ext_PL(ip,iw) / ext_PL(ip-5,iw) )
+    fact_s = log10( ssa_PL(ip,iw) / ssa_PL(ip-5,iw) )
+    fact_c = log10( asf_PL(ip,iw) / asf_PL(ip-5,iw) )
+  endif            
 
   fact_t = fact_t / log10( press_PL(ip) / press_PL(ip-5) )

trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

@@ -18 +18 @@
                               init_time, daysec, dtphys
   use comcstfi_mod, only: rad, cpp, g, r, rcp, &
-                          mugaz, pi, avocado
+                          mugaz, pi, avocado, kbol
   use planete_mod, only: nres
   use planetwide_mod, only: planetwide_sumval
@@ -84 +84 @@
   pi=2.*asin(1.)
   avocado = 6.02214179e23   ! added by RW
+  kbol = 1.38064852e-23  ! added by JVO for Titan chem
 
   ! Initialize some "temporal and calendar" related variables
@@ -162 +163 @@
      flatten = 0.0
      call getin_p("flatten", flatten)
-     write(*,*) "flatten = ", flatten
-         
+     write(*,*) "flatten = ", flatten         
 
      write(*,*) "Needed if oblate=.true.: J2"
@@ -300 +300 @@
      call getin_p("graybody",graybody)
      write(*,*)" graybody = ",graybody
+
+! Chemistry
+
+     write(*,*) "Run with or without chemistry?"
+     callchim=.false. ! default value
+     call getin_p("callchim",callchim)
+     write(*,*) " callchim = ",callchim
+
+     ! sanity check
+     if (callchim.and.(.not.tracer)) then
+       print*,"You are running chemistry without tracer"
+       print*,"Please start again with tracer =.true."
+       stop
+     endif
+
+     write(*,*)"Chemistry is computed every ichim", &
+                   " physical timestep"
+     ichim=1 ! default value
+     call getin_p("ichim",ichim)
+     write(*,*)" ichim = ",ichim
 
 ! Soil model
@@ -497 +517 @@
   PRINT*,'        or each ',iradia*dtphys,' seconds'
   PRINT*
-
+  PRINT*,'inifis: Chemistry is computed:'
+  PRINT*,'           each ',ichim,' physical time-step'
+  PRINT*,'        or each ',ichim*dtphys,' seconds'
+  PRINT*
 
 !-----------------------------------------------------------------------

trunk/LMDZ.TITAN/libf/phytitan/initracer.F

@@ -19 +19 @@
 c   ------
 c            Ehouarn Millour (oct. 2008) identify tracers by their names
+c            Jan Vatant d'Ollone (fev. 2017) chemical species for Titan
 c=======================================================================
 
@@ -26 +27 @@
       character(len=20) :: txt ! to store some text
       integer iq,ig,count
-      real r0_lift , reff_lift
 !      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
 
@@ -65 +65 @@
         enddo
 
-
 ! Identify tracers by their names: (and set corresponding values of mmol)
       ! 0. initialize tracer indexes to zero:
       ! NB: igcm_* indexes are commons in 'tracer.h'
-      igcm_n2=0
-
+      
+      igcm_h       =0
+      igcm_h2      =0
+      igcm_ch      =0
+      igcm_ch2s    =0
+      igcm_ch2     =0
+      igcm_ch3     =0
+      igcm_ch4     =0
+      igcm_c2      =0
+      igcm_c2h     =0
+      igcm_c2h2    =0
+      igcm_c2h3    =0
+      igcm_c2h4    =0
+      igcm_c2h5    =0
+      igcm_c2h6    =0
+      igcm_c3h3    =0
+      igcm_c3h5    =0
+      igcm_c3h6    =0
+      igcm_c3h7    =0
+      igcm_c4h     =0
+      igcm_c4h3    =0
+      igcm_c4h4    =0
+      igcm_c4h2s   =0
+      igcm_ch2cch2 =0
+      igcm_ch3cch  =0
+      igcm_c3h8    =0
+      igcm_c4h2    =0
+      igcm_c4h6    =0
+      igcm_c4h10   =0
+      igcm_ac6h6   =0
+      igcm_c3h2    =0
+      igcm_c4h5    =0
+      igcm_ac6h5   =0
+      igcm_n2      =0
+      igcm_n4s     =0
+      igcm_cn      =0
+      igcm_hcn     =0
+      igcm_h2cn    =0
+      igcm_chcn    =0
+      igcm_ch2cn   =0
+      igcm_ch3cn   =0
+      igcm_c3n     =0
+      igcm_hc3n    =0
+      igcm_nccn    =0
+      igcm_c4n2    =0      
+      
       write(*,*) 'initracer: noms() ', noms
 
 
       ! 1. find chemistry tracers
-      count = 0.    
-
+      count = 0.
+      
+      do iq=1,nq
+      
+        if (noms(iq).eq."h") then
+          igcm_h=iq
+          mmol(igcm_h)=1.01
+          count=count+1
+        endif
+        if (noms(iq).eq."h2") then
+          igcm_h2=iq
+          mmol(igcm_h2)=2.0158
+          count=count+1
+        endif
+        if (noms(iq).eq."ch") then
+          igcm_ch=iq
+          mmol(igcm_ch)=13.02
+          count=count+1
+        endif
+        if (noms(iq).eq."ch2s") then
+          igcm_ch2s=iq
+          mmol(igcm_ch2s)=14.03
+          count=count+1
+        endif
+        if (noms(iq).eq."ch2") then
+          igcm_ch2=iq
+          mmol(igcm_ch2)=14.03
+          count=count+1
+        endif
+        if (noms(iq).eq."ch3") then
+          igcm_ch3=iq
+          mmol(igcm_ch3)=15.03
+          count=count+1
+        endif
+        if (noms(iq).eq."ch4") then
+          igcm_ch4=iq
+          mmol(igcm_ch4)=16.04
+          count=count+1
+        endif
+        if (noms(iq).eq."c2") then
+          igcm_c2=iq
+          mmol(igcm_c2)=24.02
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h") then
+          igcm_c2h=iq
+          mmol(igcm_c2h)=25.03
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h2") then
+          igcm_c2h2=iq
+          mmol(igcm_c2h2)=26.04
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h3") then
+          igcm_c2h3=iq
+          mmol(igcm_c2h3)=27.05
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h4") then
+          igcm_c2h4=iq
+          mmol(igcm_c2h4)=28.05
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h5") then
+          igcm_c2h5=iq
+          mmol(igcm_c2h5)=29.06
+          count=count+1
+        endif
+        if (noms(iq).eq."c2h6") then
+          igcm_c2h6=iq
+          mmol(igcm_c2h6)=30.07
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h3") then
+          igcm_c3h3=iq
+          mmol(igcm_c3h3)=39.06
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h5") then
+          igcm_c3h5=iq
+          mmol(igcm_c3h5)=41.07
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h6") then
+          igcm_c3h6=iq
+          mmol(igcm_c3h6)=42.08
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h7") then
+          igcm_c3h7=iq
+          mmol(igcm_c3h7)=43.09
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h") then
+          igcm_c4h=iq
+          mmol(igcm_c4h)=49.05
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h3") then
+          igcm_c4h3=iq
+          mmol(igcm_c4h3)=51.07
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h4") then
+          igcm_c4h4=iq
+          mmol(igcm_c4h4)=52.08
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h2s") then
+          igcm_c4h2s=iq
+          mmol(igcm_c4h2s)=50.06
+          count=count+1
+        endif
+        if (noms(iq).eq."ch2cch2") then
+          igcm_ch2cch2=iq
+          mmol(igcm_ch2cch2)=40.07
+          count=count+1
+        endif
+        if (noms(iq).eq."ch3cch") then
+          igcm_ch3cch=iq
+          mmol(igcm_ch3cch)=40.07
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h8") then
+          igcm_c3h8=iq
+          mmol(igcm_c3h8)=44.11
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h2") then
+          igcm_c4h2=iq
+          mmol(igcm_c4h2)=50.06
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h6") then
+          igcm_c4h6=iq
+          mmol(igcm_c4h6)=54.09
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h10") then
+          igcm_c4h10=iq
+          mmol(igcm_c4h10)=58.13
+          count=count+1
+        endif
+        if (noms(iq).eq."ac6h6") then
+          igcm_ac6h6=iq
+          mmol(igcm_ac6h6)=78.1136
+          count=count+1
+        endif
+        if (noms(iq).eq."c3h2") then
+          igcm_c3h2=iq
+          mmol(igcm_c3h2)=38.05
+          count=count+1
+        endif
+        if (noms(iq).eq."c4h5") then
+          igcm_c4h5=iq
+          mmol(igcm_c4h5)=53.07
+          count=count+1
+        endif
+        if (noms(iq).eq."ac6h5") then
+          igcm_ac6h5=iq
+          mmol(igcm_ac6h5)=77.1136
+          count=count+1
+        endif
+        if (noms(iq).eq."n2") then
+          igcm_n2=iq
+          mmol(igcm_n2)=28.0134
+          count=count+1
+        endif
+        if (noms(iq).eq."n4s") then
+          igcm_n4s=iq
+          mmol(igcm_n4s)=14.01
+          count=count+1
+        endif
+        if (noms(iq).eq."cn") then
+          igcm_cn=iq
+          mmol(igcm_cn)=26.02
+          count=count+1
+        endif
+        if (noms(iq).eq."hcn") then
+          igcm_hcn=iq
+          mmol(igcm_hcn)=27.04
+          count=count+1
+        endif
+        if (noms(iq).eq."h2cn") then
+          igcm_h2cn=iq
+          mmol(igcm_h2cn)=28.05
+          count=count+1
+        endif
+        if (noms(iq).eq."chcn") then
+          igcm_chcn=iq
+          mmol(igcm_chcn)=39.05
+          count=count+1
+        endif
+        if (noms(iq).eq."ch2cn") then
+          igcm_ch2cn=iq
+          mmol(igcm_cn)=40.04
+          count=count+1
+        endif
+        if (noms(iq).eq."ch3cn") then
+          igcm_ch3cn=iq
+          mmol(igcm_ch3cn)=41.05
+          count=count+1
+        endif
+        if (noms(iq).eq."c3n") then
+          igcm_c3n=iq
+          mmol(igcm_c3n)=50.04
+          count=count+1
+        endif
+        if (noms(iq).eq."hc3n") then
+          igcm_hc3n=iq
+          mmol(igcm_hc3n)=51.05
+          count=count+1
+        endif
+        if (noms(iq).eq."nccn") then
+          igcm_nccn=iq
+          mmol(igcm_nccn)=52.04
+          count=count+1
+        endif
+        if (noms(iq).eq."c4n2") then
+          igcm_c4n2=iq
+          mmol(igcm_c4n2)=76.1
+          count=count+1
+        endif
+        
+
+      enddo ! of do iq=1,nq
  
       ! check that we identified all tracers:

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -33 +33 @@
       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
       use time_phylmdz_mod, only: daysec
+      use logic_mod, only: moyzon_ch
+      use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
+                            zplaybar, zzlevbar, zzlaybar, ztfibar
       use callkeys_mod
       use vertical_layers_mod, only: presnivs, pseudoalt
@@ -70 +73 @@
 !
 !      V. Tracers
-!         V.1. Aerosols and particles.
-!         V.2. Updates (pressure variations, surface budget).
-!         V.3. Surface Tracer Update.
+!         V.1. Chemistry
+!         V.2. Aerosols and particles.
+!         V.3. Updates (pressure variations, surface budget).
+!         V.4. Surface Tracer Update.
 !
 !      VI. Surface and sub-surface soil temperature.
@@ -140 +144 @@
 !           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
 !           Purge of the code : M. Turbet (2015)
-!           Fork for Titan and clean of all too-generic (ocean, water, co2 ...) routines : J. Vatant d'Ollone (2017)
+!           Fork for Titan : J. Vatant d'Ollone (2017)
+!                + clean of all too-generic (ocean, water, co2 ...) routines
+!                + Titan's chemistry
 !============================================================================================
 
@@ -277 +283 @@
       real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
       real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
+      
+      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
+      real zdqmph(ngrid,nlayer,nq)    ! Microphysical tendency ( condensation routine only for now, no microphysical routines ).
                         
       ! For Winds : (m/s/s)
@@ -357 +366 @@
       real reffcol(ngrid,naerkind)
 
-
+! Local variables for Titan chemistry and microphysics (JVO 2017)
+! ----------------------------------------------------------------------------
+     
+      integer,parameter  :: nmicro=0 ! Temporary ! To be put in start/def
+
+      real ctimestep ! Chemistry timestep (s)
+ 
+      ! Grandeurs en moyennes zonales ------------------------
+      real temp_eq(nlayer), press_eq(nlayer)
+      real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
+!      real zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
+      real ztemp(ngrid,nlayer)
+      
+      real ychim(ngrid,nlayer,nq-nmicro)   
+
+      ! 2D vmr tendencies ( chemistry and condensation )
+      real,dimension(:,:,:),allocatable,save :: dycchi
+      ! Must be saved since chemistry is not called every step
+      
+      real dycmph(ngrid,nlayer,nq-nmicro)       
+      ! ----------------------------------------------------------
+         
+      real,dimension(:,:),allocatable,save       :: qysat
+      
+      character*10,dimension(:),allocatable,save :: nomqy
+!$OMP THREADPRIVATE(dycchi,qysat,nomqy)
+
+!-----------------------------------------------------------------------------
+      
 !==================================================================================================
 
@@ -401 +438 @@
         ALLOCATE(zdtsw(ngrid,nlayer))
         ALLOCATE(tau_col(ngrid))
-
+        ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
+        ALLOCATE(qysat(nlayer,nq-nmicro))
+        ALLOCATE(nomqy(nq-nmicro+1))
+        
         ! This is defined in comsaison_h
         ALLOCATE(mu0(ngrid))
@@ -433 +473 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0(startphy_file,                                 &
-                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
+         call phyetat0(startphy_file,ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
+                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)                       
          if (.not.startphy_file) then
            ! additionnal "academic" initialization of physics
@@ -514 +553 @@
                           ptimestep,pday+nday,time_phys,cell_area,          &
                           albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
+         endif
+
+         ! Sanity check for chemistry
+         if ( ((moyzon_ch).and.(.not.callchim)) .or. ((.not.moyzon_ch).and.(callchim)) ) then
+            print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
+            print *, "This is not compatible..."
+            stop
+         endif
+
+         ! Initialize names, timestep and saturation profiles for chemistry
+
+         if ( callchim .and. (nq.gt.nmicro) ) then
+
+            ctimestep = ptimestep*REAL(ichim)
+
+            do iq=nmicro+1,nq
+               nomqy(iq-nmicro) = nametrac(iq)
+            enddo
+
+            nomqy(nq-nmicro+1) = "HV"
+            
+            ! qysat is taken at the equator ( small variations of t,p)
+            temp_eq(:)  = tmoy(:)
+            press_eq(:) = playmoy(:)/100. ! en mbar
+            
+            call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
+          
          endif
          
@@ -601 +667 @@
       enddo     
 
+      ! -------------------------------Taken from old Titan --------------------------
+      ! zonal averages needed
+      if (moyzon_ch) then
+         
+         zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
+         ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
+         !       zzlaybar(1,:)=RG*RA*RA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
+         zzlevbar(1,1)=zphisbar(1)/g
+         DO l=2,nlayer
+            z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplevbar(1,l-1)-zplevbar(1,l))
+            z2=(zplevbar(1,l)  +zplaybar(1,l))/(zplevbar(1,l)  -zplaybar(1,l))
+            zzlevbar(1,l)=(z1*zzlaybar(1,l-1)+z2*zzlaybar(1,l))/(z1+z2)
+         ENDDO
+         zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
+
+         DO ig=2,ngrid
+            if (latitude(ig).ne.latitude(ig-1)) then
+               DO l=1,nlayer
+                  zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
+                  ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
+                  !zzlaybar(ig,l)=RG*RA*RA/(RG*RA-(zphibar(ig,l)+zphisbar(ig)))-RA
+               ENDDO
+               zzlevbar(ig,1)=zphisbar(ig)/g
+               DO l=2,nlayer
+                  z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplevbar(ig,l-1)-zplevbar(ig,l))
+                  z2=(zplevbar(ig,l)  +zplaybar(ig,l))/(zplevbar(ig,l)  -zplaybar(ig,l))
+                  zzlevbar(ig,l)=(z1*zzlaybar(ig,l-1)+z2*zzlaybar(ig,l))/(z1+z2)
+               ENDDO
+               zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
+            else
+               zzlaybar(ig,:)=zzlaybar(ig-1,:)
+               zzlevbar(ig,:)=zzlevbar(ig-1,:)
+            endif
+         ENDDO
+
+      endif  ! moyzon
+      ! -------------------------------------------------------------------------------------
+      
       ! Compute potential temperature
       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
@@ -896 +1000 @@
       if (tracer) then
       
-      
   ! -------------------------
-  !   V.1. Aerosol particles
+  !   V.1. Chemistry
+  ! -------------------------
+         if (callchim) then
+
+            zplev(:,:) = zplevbar(:,:)
+            zplay(:,:) = zplaybar(:,:)
+            zzlev(:,:) = zzlevbar(:,:)
+            zzlay(:,:) = zzlaybar(:,:)
+            ztemp(:,:) = ztfibar(:,:)
+
+            if (nq.gt.nmicro) then
+               do iq = nmicro+1,nq
+                  ychim(:,:,iq-nmicro) = pq(:,:,iq) 
+               enddo
+            endif
+
+            ! Condensation tendency after the transport
+            do iq=1,nq-nmicro
+               do l=1,nlayer
+                  do ig=1,ngrid
+                     if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
+                        dycmph(ig,l,nmicro+iq)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
+                     endif
+                  enddo
+               enddo
+            enddo
+
+            if( mod(icount-1,ichim).eq.0.) then
+               
+               print *, "On passe dans la chimie..."
+               
+               call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
+                           ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
+
+            ! JVO 2017 : NLEV = nlayer+70, en accord avec le C. Quid si nlay=/ 55 ?
+               
+            endif
+            
+            if (nq.gt.nmicro) then
+               ! We convert tendencies back to 3D and mass mixing ratio
+               do iq=nmicro+1,nq
+                  zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro)*pq(:,:,iq) / ychim(:,:,iq-nmicro)
+                  zdqmph(:,:,iq) = dycmph(:,:,iq-nmicro)*pq(:,:,iq) / ychim(:,:,iq-nmicro)
+               enddo
+
+               pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
+                    zdqchi(1:ngrid,1:nlayer,1:nq) + zdqmph(1:ngrid,1:nlayer,1:nq)
+               
+            endif
+            
+         endif
+      
+  ! -------------------------
+  !   V.2. Aerosol particles
   ! -------------------------
 
@@ -919 +1075 @@
 
   ! ---------------
-  !   V.2. Updates
+  !   V.3 Updates
   ! ---------------
 
@@ -951 +1107 @@
 
   ! -----------------------------
-  !   V.3. Surface Tracer Update
+  !   V.4. Surface Tracer Update
   ! -----------------------------
 

trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90

@@ -21 +21 @@
 ! tracer indexes: these are initialized in initracer and should be 0 if the
 !                 corresponding tracer does not exist
+
       ! chemistry:
+      
+      integer :: igcm_h
+      integer :: igcm_h2
+      integer :: igcm_ch
+      integer :: igcm_ch2s
+      integer :: igcm_ch2
+      integer :: igcm_ch3
+      integer :: igcm_ch4
+      integer :: igcm_c2
+      integer :: igcm_c2h
+      integer :: igcm_c2h2
+      integer :: igcm_c2h3
+      integer :: igcm_c2h4
+      integer :: igcm_c2h5
+      integer :: igcm_c2h6
+      integer :: igcm_c3h3
+      integer :: igcm_c3h5
+      integer :: igcm_c3h6
+      integer :: igcm_c3h7
+      integer :: igcm_c4h
+      integer :: igcm_c4h3
+      integer :: igcm_c4h4
+      integer :: igcm_c4h2s
+      integer :: igcm_ch2cch2
+      integer :: igcm_ch3cch
+      integer :: igcm_c3h8
+      integer :: igcm_c4h2
+      integer :: igcm_c4h6
+      integer :: igcm_c4h10
+      integer :: igcm_ac6h6
+      integer :: igcm_c3h2
+      integer :: igcm_c4h5
+      integer :: igcm_ac6h5
       integer :: igcm_n2
-!$OMP THREADPRIVATE(igcm_n2)
+      integer :: igcm_n4s
+      integer :: igcm_cn
+      integer :: igcm_hcn
+      integer :: igcm_h2cn
+      integer :: igcm_chcn
+      integer :: igcm_ch2cn
+      integer :: igcm_ch3cn
+      integer :: igcm_c3n
+      integer :: igcm_hc3n
+      integer :: igcm_nccn
+      integer :: igcm_c4n2 
+      
+      
+!$OMP THREADPRIVATE(igcm_h,igcm_h2,igcm_ch,igcm_ch2s,igcm_ch2,igcm_ch3,igcm_ch4,   &
+      !$OMP igcm_c2,igcm_c2h,igcm_c2h2,igcm_c2h3,igcm_c2h4,igcm_c2h5,igcm_c2h6,    &
+      !$OMP igcm_c3h3,igcm_c3h5,igcm_c3h6,igcm_c3h7,igcm_c4h,igcm_c4h3,igcm_c4h4,  &
+      !$OMP igcm_c4h2s,igcm_ch2cch2,igcm_ch3cch,igcm_c3h8,igcm_c4h2,igcm_c4h6,     &
+      !$OMP igcm_c4h10,igcm_ac6h6,igcm_c3h2,igcm_c4h5,igcm_ac6h5,igcm_n2,igcm_n4s, &
+      !$OMP igcm_cn,igcm_hcn,igcm_h2cn,igcm_chcn,igcm_ch2cn,igcm_ch3cn,igcm_c3n,   &
+      !$OMP igcm_hc3n,igcm_nccn,igcm_c4n2)
 
        end module tracer_h