Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3288)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3289)
@@ -4583,2 +4583,10 @@
 when the first (at every time step) variable sent to wstats() was not part of 
 those listed in stats.def
+
+== 21/03/2024 == Jliu
+Fix several bugs in photochemistry related to the submit-crash-resubmit
+problem: Some of the integers in the related routines are not set in
+privatethreads. Consequently, the additives such as n=n+1 are accumulated with
+threads, causing systermatic problems in the routine. This update fixed the
+submit-crash-resubmit problem. The simulated results are same with previous
+simulations as tested.
Index: trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 3288)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 3289)
@@ -32,5 +32,5 @@
       use comcstfi_h, only: pi
       use chemistrydata_mod, only: read_phototable
-      use photolysis_mod, only: init_photolysis, nphot
+      use photolysis_mod, only: init_photolysis !, nphot
       use iono_h, only: temp_elect
       use wstats_mod, only: wstats
@@ -168,12 +168,19 @@
       integer :: ig_vl1
 
-      integer :: nb_reaction_3_max     ! number of quadratic reactions
-      integer :: nb_reaction_4_max     ! number of bimolecular reactions
-      integer :: nquench               ! number of quenching + heterogeneous reactions
-      integer :: nphotion              ! number of photoionizations
-      integer :: nb_reaction_4_ion     ! quadratic reactions for ionosphere
-      integer :: nb_reaction_4_deut    ! quadratic reactions for deuterium chem
-      integer :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
-
+      integer, save :: nb_reaction_3_max     ! number of quadratic reactions
+      integer, save :: nb_reaction_4_max     ! number of bimolecular reactions
+      integer, parameter :: nb_reaction_3_max0 = 6    ! number of quadratic reactions
+      integer, parameter :: nb_reaction_4_max0 = 31   ! number of bimolecular reactions
+      integer, parameter :: nquench =  9            ! number of quenching + heterogeneous reactions
+      integer, save :: nphotion            ! number of photoionizations
+      integer, parameter :: nphotion0 = 0            ! number of photoionizations
+      integer, save :: nphot                       ! number of photolysis 
+      integer, parameter :: nphot0 = 13           ! number of photolysis 
+                                                  ! is incremented by +2 in calchim if deuterium chemisty)
+      integer, parameter :: nb_reaction_4_ion = 64     ! quadratic reactions for ionosphere
+      integer, parameter :: nb_reaction_4_deut = 35    ! quadratic reactions for deuterium chem
+      integer, save :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
+      integer, parameter :: nb_phot_max0 = 0 ! total number of photolysis+photoionizations+quenching reactions
+!$OMP THREADPRIVATE(nb_reaction_3_max,nb_reaction_4_max,nphotion,nphot,nb_phot_max) 
 
       real    :: latvl1, lonvl1
@@ -757,10 +764,11 @@
          write(*,*) 'calchim: tracer indices=',niq(:)
 
-
+         nphot = nphot0         ! number of photolysis 
+                                ! is incremented by +2 in calchim if deuterium chemisty)
          if (photochem) then
             if (jonline) then
                print*,'calchim: Read UV absorption cross-sections'
                !Add two photodissociations in deuterium chemistry included
-               if(deutchem) nphot = nphot + 2
+               if(deutchem) nphot = nphot0 + 2
                call init_photolysis     ! on-line photolysis
             else
@@ -843,17 +851,20 @@
          if (photochem) then
             ! set number of reactions, depending on ion chemistry or not
-            nb_reaction_4_ion  = 64
-            nb_reaction_4_deut = 35
+           ! nb_reaction_4_ion  = 64
+           ! nb_reaction_4_deut = 35
 
             !Default numbers if no ion and no deuterium chemistry included
 
-            nb_reaction_4_max = 31     ! set number of bimolecular reactions
-            nb_reaction_3_max = 6      ! set number of quadratic reactions
-            nquench           = 9      ! set number of quenching + heterogeneous
-            nphotion          = 0      ! set number of photoionizations
-
+           ! nb_reaction_4_max = 31     ! set number of bimolecular reactions
+           ! nb_reaction_3_max = 6      ! set number of quadratic reactions
+           ! nquench           = 9      ! set number of quenching + heterogeneous
+           ! nphotion          = 0      ! set number of photoionizations
+            nb_reaction_4_max = nb_reaction_4_max0     ! set number of bimolecular reactions
+            nb_reaction_3_max = nb_reaction_3_max0      ! set number of quadratic reactions
+            nphotion          = nphotion0      ! set number of photoionizations
+            nb_phot_max = nb_phot_max0
             if (ionchem) then
-               nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion  
-               nphotion          = 18   ! set number of photoionizations
+               nb_reaction_4_max = nb_reaction_4_max0 + nb_reaction_4_ion  
+               nphotion          = nphotion0 + 18   ! set number of photoionizations
             endif
             if(deutchem) then
@@ -864,6 +875,8 @@
 !        numerically as a photolysis:
 
-            nb_phot_max = nphot + nphotion + nquench
-
+            nb_phot_max = nb_phot_max0 + nphot + nphotion + nquench
+!          print*, 'nb_phot_max = ', nb_phot_max
+!	  print*, 'nb_reaction_4_max=', nb_reaction_4_max 
+!          print*, 'nb_reaction_3_max=', nb_reaction_3_max 
 !        call main photochemistry routine
 
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 3288)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 3289)
@@ -161,6 +161,6 @@
 !integer,parameter :: i_hdo2    = 38
 
-integer :: i_last
-
+integer, save :: i_last
+!$OMP THREADPRIVATE(i_last)
 !Tracer indexes for photionization coeffs
 
@@ -180,7 +180,7 @@
 
 integer :: ilay
-integer :: ind_norec
-integer :: ind_orec
-
+integer, save :: ind_norec
+integer, save :: ind_orec
+!$OMP THREADPRIVATE(ind_norec, ind_orec)
 real :: ctimestep           ! standard timestep for the chemistry (s) 
 real :: dt_guess            ! first-guess timestep (s) 
@@ -503,5 +503,6 @@
 em_no(:)=c(:,i_o)*c(:,i_n)*v_4(:,ind_norec)   !2.8e-17*(300./temp(:)))**0.5
 em_o2(:)=0.75*c(:,i_o)*c(:,i_o)*c(:,i_co2)*v_3(:,ind_orec)   !2.5*9.46e-34*exp(485./temp(:))*dens(:)
-
+!print*, "ind_norec= ", ind_norec
+!print*, "ind_orec= ", ind_orec
 !===================================================================
 !     save chemical species for the gcm       
@@ -659,5 +660,5 @@
 
 use comcstfi_h
-use photolysis_mod, only : nphot
+!use photolysis_mod, only : nphot
 
 implicit none
@@ -675,4 +676,5 @@
 logical, intent(in)     :: deutchem
 integer                 :: lswitch           ! interface level between lower
+integer,parameter :: nphot0 =13
                                              ! atmosphere and thermosphere chemistries
 real, dimension(nlayer) :: dens              ! total number density (molecule.cm-3)
@@ -695,6 +697,6 @@
 real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
 real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
-integer :: ind_norec
-integer :: ind_orec
+integer, intent(out) :: ind_norec
+integer, intent(out) :: ind_orec
 
 !----------------------------------------------------------------------
@@ -703,5 +705,6 @@
 
 integer          :: ilev
-integer          :: nb_phot, nb_reaction_3, nb_reaction_4
+integer, save    :: nb_phot, nb_reaction_3, nb_reaction_4
+!$OMP THREADPRIVATE(nb_phot, nb_reaction_3, nb_reaction_4)
 real :: ak0, ak1, xpo, rate, rate1, rate2
 real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
@@ -737,6 +740,6 @@
 !     initialisation
 !----------------------------------------------------------------------
-
-      nb_phot       = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
+      nb_phot       = 0
+      nb_phot       = nphot0 + nphotion ! initialised to the number of photolysis + number of photoionization rates
       nb_reaction_3 = 0
       nb_reaction_4 = 0
@@ -776,5 +779,5 @@
       v_3(:,nb_reaction_3) = a002(:)
       ind_orec=nb_reaction_3
-
+!      print*,"ind_orec",ind_orec
 !---  a003: o + o3 -> o2 + o2
 
@@ -1200,5 +1203,5 @@
       v_4(:,nb_reaction_4) = d007(:)
       ind_norec = nb_reaction_4
-
+!      print*,"ind_norec=",ind_norec
 !---  d008: n + ho2 -> no + oh
 
