Changeset 1496 for trunk

Timestamp:

Nov 10, 2015, 2:36:29 PM (10 years ago)

Author:

flefevre

Message:

Optimisation du solveur chimique ASIS

File:

• trunk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90 (modified) (14 diffs)

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

@@ -99 +99 @@
 real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3) 
  
-!real (kind = 8), dimension(ngridmx,nlayer,nesp), save :: cold3d ! cold saved from previous timestep
-!real (kind = 8), dimension(ngridmx,nlayer), save      :: dt3d   ! saved timestep
-
 ! reaction rates
  
@@ -111 +108 @@
 ! matrix
 
-real (kind = 8), dimension(nesp,nesp) :: mat
+real (kind = 8), dimension(nesp,nesp) :: mat, mat1
 integer, dimension(nesp)              :: indx
 integer                               :: code
+
+! production and loss terms (for first-guess solution only)
+
+real (kind = 8), dimension(nesp) :: prod, loss
 
 ! curvatures
@@ -128 +129 @@
                i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
                i_n, i_n2d, i_no, i_no2, i_n2)
-!  dt3d(:,:) = 600.
    firstcall = .false.
 end if
@@ -199 +199 @@
    dt_guess = ctimestep
    cold(:) = c(ilev,:)
-!  dt_guess = dt3d(ig,ilev)
-!  cold(:) = cold3d(ig,ilev,:)
 
 !  internal loop for the chemistry
@@ -206 +204 @@
    do while (time < ptimestep)
 
-   iter(ilev) = iter(ilev) + 1
+   iter(ilev) = iter(ilev) + 1    ! iteration counter
   
 !  first-guess: fill matrix
 
-   call fill_matrix(ilev, nesp, c, nlayer, v_phot, v_3, v_4, dt_guess, mat)
+   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
+
+! first guess : form the matrix identity + mat*dt_guess
+
+   mat(:,:) = mat1(:,:)*dt_guess
+   do iesp = 1,nesp
+      mat(iesp,iesp) = 1. + mat(iesp,iesp)
+   end do
 
 !  first-guess: solve linear system
@@ -266 +271 @@
 
    if (time + dt_corrected > ptimestep) then
-!     dt3d(ig,ilev) = dt_corrected
       dt_corrected = ptimestep - time
    end if
@@ -272 +276 @@
    if (dt_corrected /= dt_guess) then  ! the timestep has been modified
 
-!  fill matrix
-
-      call fill_matrix(ilev, nesp, c, nlayer, v_phot, v_3, v_4, dt_corrected, mat)
+!  form the matrix identity + mat*dt_guess
+
+      mat(:,:) = mat1(:,:)*dt_corrected
+      do iesp = 1,nesp
+         mat(iesp,iesp) = 1. + mat(iesp,iesp)
+      end do
 
 !  solve linear system
@@ -956 +963 @@
 !======================================================================
 
- subroutine fill_matrix(ilev, nesp, c, nlayer, v_phot, v_3, v_4,  &
-                        dtx, mat)
+ subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
 
 !======================================================================
@@ -974 +980 @@
 integer             :: nesp    ! number of species in the chemistry
 integer, intent(in) :: nlayer  ! number of atmospheric layers
-
-real                :: dtx     ! timestep (s)
-real(kind = jprb)   :: eps, eps_4  ! implicit/explicit coefficient
 
 real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
@@ -985 +988 @@
 ! output
 
-real (kind = 8), dimension(nesp,nesp)     :: mat  ! matrix
+real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
+real (kind = 8), dimension(nesp), intent(out)      :: prod, loss
 
 ! local
@@ -995 +999 @@
 integer :: iphot,i3,i4
 
-! initialisation 
+real(kind = jprb) :: eps, eps_4  ! implicit/explicit coefficient
+
+! initialisations 
 
 mat(:,:) = 0.
-do iesp = 1,nesp
-   mat(iesp,iesp) = 1.
-end do  
+prod(:)  = 0.
+loss(:)  = 0.
 
 ! photodissociations
@@ -1012 +1017 @@
   ind_phot_6 = indice_phot(iphot)%z6
 
-  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)*dtx
-  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)*dtx
-  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)*dtx
+  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
+  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
+  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
+
+  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
+  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
 
 end do
@@ -1026 +1035 @@
   ind_3_6 = indice_3(i3)%z6
 
-  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)*dtx
-  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*dtx
-  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*dtx
+  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
+  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)
+  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)
+
+  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
+  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
 
 end do
@@ -1046 +1059 @@
   eps_4 = min(eps_4,1.0_jprb)
 
-  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*dtx*(1. - eps_4)*c(ilev,ind_4_4) 
-  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*dtx*eps_4*c(ilev,ind_4_2)
-  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*dtx*(1. - eps_4)*c(ilev,ind_4_4)
-  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*dtx*eps_4*c(ilev,ind_4_2)   
-  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*dtx*(1. - eps_4)*c(ilev,ind_4_4)
-  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*dtx*eps_4*c(ilev,ind_4_2)
-  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*dtx*(1. - eps_4)*c(ilev,ind_4_4)
-  mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*dtx*eps_4*c(ilev,ind_4_2)
+  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4) 
+  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
+  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)   
+  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
+  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
+
+  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
+  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
+  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
 
 end do