Changeset 1569 for trunk/LMDZ.MARS/libf/aeronomars

Timestamp:

Jun 29, 2016, 6:08:29 PM (9 years ago)

Author:

flefevre

Message:

Nouvelle version du solveur chimique ASIS:

• subroutine separee de calcul du pas de temps variable
• calcul du first-guess de la solution en semi-implicite
• balayage de toutes les especes pour le calcul de la courbure de la solution
• introduction des tolerances absolue et relative sur l'erreur de la solution (reglages provisoires)

File:

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

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

@@ -85 +85 @@
 integer,parameter :: i_n2   = 17
 
-real :: stimestep           ! standard timestep for the chemistry (s) 
-real :: ctimestep           ! real timestep for the chemistry (s) 
+real :: ctimestep           ! standard timestep for the chemistry (s) 
 real :: dt_guess            ! first-guess timestep (s) 
 real :: dt_corrected        ! corrected timestep (s) 
-real :: dt_min = 1.         ! minimum allowed timestep (s) 
-real :: dtg                 ! correction factor for the timestep (s) 
 real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
 real :: time                ! internal time (between 0 and ptimestep, in s)
@@ -115 +112 @@
 
 real (kind = 8), dimension(nesp) :: prod, loss
-
-! curvatures
-
-real :: ratio, curv, e, e1, e2, e3
 
 !===================================================================
@@ -169 +162 @@
 
 !===================================================================
-!     stimestep : standard chemical timestep (s)                           
-!     ctimestep : real chemical timestep (s),
-!                 taking into account the physical timestep                           
-!===================================================================
-
-stimestep = 600. ! standard value : 10 mn
-
-phychemrat = nint(ptimestep/stimestep)
+!     ctimestep : standard chemical timestep (s), defined as 
+!                 the fraction phychemrat of the physical timestep                           
+!===================================================================
+
 phychemrat = 1
 
 ctimestep = ptimestep/real(phychemrat)
 
-!print*, "stimestep  = ", stimestep
 !print*, "ptimestep  = ", ptimestep
 !print*, "phychemrat = ", phychemrat
@@ -210 +198 @@
    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
+!  adaptative evaluation of the sub time step
+
+   call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:),  &
+                  mat1, prod, loss)
+
+   if (time + dt_corrected > ptimestep) then
+      dt_corrected = ptimestep - time
+   end if
+
+!  if (dt_corrected /= dt_guess) then  ! the timestep has been modified
+
+!  form the matrix identity + mat*dt_corrected
+
+   mat(:,:) = mat1(:,:)*dt_corrected
    do iesp = 1,nesp
       mat(iesp,iesp) = 1. + mat(iesp,iesp)
    end do
 
-!  first-guess: solve linear system
+!  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
      
    cnew(:) = c(ilev,:)
@@ -227 +226 @@
    stop
 #endif
+
+!  end if
+
 !  eliminate small values
 
@@ -233 +235 @@
    end where
 
-!  curvature of the solution and weighted error norm for O
-!
-!  cold : concentrations at time t - dt 
-!  c    : concentrations at time t
-!  cnew : concentrations guessed at time t + dt
-
-!  ratio = dtold/dtnew
-!  here ratio = 1 since the timestep for the guess is equal to the previous timestep
-
-   ratio = 1.
-
-   curv = 2.*(ratio*cnew(i_o) - (1. + ratio)*c(ilev,i_o) + cold(i_o)) &
-          /(ratio*(1. + ratio))
-   e1 = (curv/(cnew(i_o) + cnew(i_o3)))*100.
-   e1 = abs(e1)
-
-!  curvature of the solution and weighted error norm for HO2
-
-   curv = 2.*(ratio*cnew(i_ho2) - (1. + ratio)*c(ilev,i_ho2) + cold(i_ho2)) &
-          /(ratio*(1. + ratio))
-   e2 = (curv/(cnew(i_h) + cnew(i_oh) + cnew(i_ho2) + 2.*cnew(i_h2o2)))*100.
-   e2 = abs(e2)
-
-!  curvature of the solution and weighted error norm for NO2
-
-!  curv = 2.*(ratio*cnew(i_no2) - (1. + ratio)*c(ilev,i_no2) + cold(i_no2)) &
-!         /(ratio*(1. + ratio))
-!  e3 = (curv/(cnew(i_no) + cnew(i_no2)))*100.
-!  e3 = abs(e3)
-   e3 = 0.
-
-!  timestep correction
-
-   e = max(e1, e2)
-   e = max(e,  e3)
-   e = max(e, 0.1)
-   dtg = max(0.001, min(2.5,0.8/sqrt(e)))
-!  dtg = 1.                                     ! uncomment this line to turn off the variable timestep scheme
-   dt_corrected = max(dt_min,dtg*dt_guess)      ! minimal timestep: dt_min
-   dt_corrected = min(dt_corrected,ctimestep)   ! maximal timestep: ctimestep
-
-   if (time + dt_corrected > ptimestep) then
-      dt_corrected = ptimestep - time
-   end if
-
-   if (dt_corrected /= dt_guess) then  ! the timestep has been modified
-
-!  form the matrix identity + mat*dt_corrected
-
-      mat(:,:) = mat1(:,:)*dt_corrected
-      do iesp = 1,nesp
-         mat(iesp,iesp) = 1. + mat(iesp,iesp)
-      end do
-
-!  solve linear system
-     
-      cnew(:) = c(ilev,:)
+!  update concentrations
+
+   cold(:)   = c(ilev,:)
+   c(ilev,:) = cnew(:)
+   cnew(:)   = 0.
+
+!  increment internal time
+
+   time = time + dt_corrected
+   dt_guess = dt_corrected     ! first-guess timestep for next iteration
+
+   end do ! while (time < ptimestep)
+
+end do ! ilev
+
+!===================================================================
+!     save chemical species for the gcm       
+!===================================================================
+
+call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
+               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
+               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
+               i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
+
+contains
+
+!================================================================
+
+ subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, prod, loss)
+
+!================================================================
+! iterative evaluation of the appropriate time step dtnew
+! according to curvature criterion based on
+! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
+! with r = (tn - tn-1)/(tn+1 - tn)
+!================================================================
+
+implicit none
+
+! input
+
+integer :: nesp  ! number of species in the chemistry
+
+real :: dtold, ctimestep
+real (kind = 8), dimension(nesp)      :: cold, ccur
+real (kind = 8), dimension(nesp,nesp) :: mat1
+real (kind = 8), dimension(nesp)      :: prod, loss
+
+! output
+
+real :: dtnew
+
+! local
+
+real (kind = 8), dimension(nesp)      :: cnew
+real (kind = 8), dimension(nesp,nesp) :: mat
+real (kind = 8) :: ratio, e, es, coef
+
+integer                  :: code, iesp, iter
+integer, dimension(nesp) :: indx
+
+real :: dttest
+
+! parameters
+
+real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
+real (kind = 8), parameter :: atol    = 1.e4     ! recommended value : 1.e4
+real (kind = 8), parameter :: rtol    = 1./0.1   ! 1/rtol recommended value : 0.1-0.02
+integer,         parameter :: niter   = 3        ! number of iterations
+real (kind = 8), parameter :: coefmax = 2.
+real (kind = 8), parameter :: coefmin = 0.1 
+logical                    :: fast_guess = .true.
+
+dttest = dtold   ! dttest = dtold = dt_guess
+
+do iter = 1,niter
+
+if (fast_guess) then
+
+! first guess : fast implicit method
+
+   do iesp = 1, nesp
+      cnew(iesp) = (ccur(iesp) + prod(iesp)*dttest)/(1. + loss(iesp)*dttest)
+   end do
+
+else
+
+! first guess : form the matrix identity + mat*dt_guess
+
+   mat(:,:) = mat1(:,:)*dttest
+   do iesp = 1,nesp
+      mat(iesp,iesp) = 1. + mat(iesp,iesp)
+   end do
+
+! form right-hand side (RHS) of the system
+
+   cnew(:) = ccur(:)
+
+! solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
 
 #ifdef LAPACK
@@ -297 +341 @@
    stop
 #endif
+
+end if
+
+! ratio old/new subtimestep
+
+ratio = dtold/dttest
+
+! e : local error indicatocitr
+
+e = 0.
+
+do iesp = 1,nesp
+   es = 2.*abs((ratio*cnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp))   &
+         /(1. + ratio)/max(ccur(iesp),atol))
+
+   if (es > e) then
+      e = es
    end if
-
-!  eliminate small values
-
-   where (cnew(:)/dens(ilev) < 1.e-30)
-      cnew(:) = 0.
-   end where
-
-   cold(:)   = c(ilev,:)
-   c(ilev,:) = cnew(:)
-   cnew(:)   = 0.
-
-!  increment internal time
-
-   time = time + dt_corrected
-   dt_guess = dt_corrected     ! first-guess timestep for next iteration
-
-   end do ! while (time < ptimestep)
-
-end do ! ilev
-
-!===================================================================
-!     save chemical species for the gcm       
-!===================================================================
-
-call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
-               i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
-
-contains
+end do
+e = rtol*e
+
+! timestep correction
+
+coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
+
+dttest = max(dtmin,dttest*coef)
+dttest = min(ctimestep,dttest)
+
+end do ! iter
+
+! new timestep
+
+dtnew = dttest
+
+end subroutine define_dt
 
 !======================================================================
@@ -973 +1021 @@
 integer :: iphot,i3,i4
 
-real(kind = jprb) :: eps, eps_4  ! implicit/explicit coefficient
+real(kind = 8) :: eps, eps_4  ! implicit/explicit coefficient
 
 ! initialisations 
@@ -1031 +1079 @@
 
   eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
-  eps_4 = min(eps_4,1.0_jprb)
+  eps_4 = min(eps_4,1.0)
 
   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)