Changeset 1810

Timestamp:

Oct 24, 2017, 1:20:18 PM (7 years ago)

Author:

bclmd

Message:

Updating photochemical core for the early Earth

Location:

trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth

Files:

• calchim_asis_early_earth.F90 (modified) (3 diffs)
• chimiedata_early_earth.h (modified) (1 diff)
• photochemistry_asis_early_earth.F90 (modified) (11 diffs)

trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/calchim_asis_early_earth.F90

@@ -70 +70 @@
 !=======================================================================
 
-#include "chimiedata.h"
+#include "chimiedata_early_earth.h"
 
 !     input:
@@ -159 +159 @@
  
       logical,save :: firstcall = .true.
-      logical,save :: depos = .true.  ! switch for dry deposition
+      logical,save :: depos = .false.  ! switch for dry deposition
 
 !     for each column of atmosphere:
@@ -543 +543 @@
 !           search for switch index between regions
 
-            if (photochem .and. thermochem) then
-               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
-                  lswitch = l
-                  foundswitch = 1
-               end if
-            end if
+!            if (photochem .and. thermochem) then
+!               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
+!                  lswitch = l
+!                  foundswitch = 1
+!               end if
+!            end if
             if (.not. photochem) then
                lswitch = 22
             end if
-            if (.not. thermochem) then
+!            if (.not. thermochem) then
                lswitch = min(50,nlayer+1)
-            end if
+!            end if
 
          end do ! of do l=1,nlayer

trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/chimiedata_early_earth.h

@@ -19 +19 @@
       real, parameter :: kb = 1.3806e-23
 
-      common/chimiedata/jphot,colairtab,table_ozo
+      common/chimiedata_early_earth/jphot,colairtab,table_ozo
 
       real jphot(ntemp,nsza,nz,nozo,ntau,nch4,nd)

trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/photochemistry_asis_early_earth.F90

@@ -19 +19 @@
 implicit none
 
-#include "chimiedata.h" 
+#include "chimiedata_early_earth.h" 
 
 !===================================================================
@@ -178 +178 @@
 !===================================================================
 
-call photolysis_asis(nlayer, ngrid, lswitch, press, temp, sza, fractcol,tau, zmmean, dist_sol, & 
+call photolysis_asis_early_earth(nlayer, ngrid, lswitch, press, temp, sza, fractcol,tau, zmmean, dist_sol, & 
                      rm(:,i_co2), rm(:,i_o3), rm(:,i_ch4), v_phot)
 
@@ -249 +249 @@
    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, dens(ilev))
+
+   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,:)
 
-!#ifdef LAPACK
+#ifdef LAPACK
    call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
-!#else
-!   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-!   stop
-!#endif
+#else
+#   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
+#   stop
+#endif
+
+!  end if
+
 !  eliminate small values
 
@@ -272 +286 @@
    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,:)
-
-!#ifdef LAPACK
-      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
-!#else
-!   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-!   stop
-!#endif
-   end if
-
-!  eliminate small values
-
-   where (cnew(:)/dens(ilev) < 1.e-30)
-      cnew(:) = 0.
-   end where
+!  update concentrations
 
    cold(:)   = c(ilev,:)
@@ -356 +300 @@
 
 end do ! ilev
+
+
+
 
 !===================================================================
@@ -370 +317 @@
 
 contains
+
+!================================================================
+
+ subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
+                      prod, loss, dens)
+
+!================================================================
+! 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
+real                                  :: dens
+
+! output
+
+real :: dtnew
+
+! local
+
+real (kind = 8), dimension(nesp)      :: cnew
+real (kind = 8), dimension(nesp,nesp) :: mat
+real (kind = 8) :: atol, 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 :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
+real (kind = 8), parameter :: rtol    = 1./0.05   ! 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
+
+atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
+
+do iter = 1,niter
+
+if (fast_guess) then
+
+! first guess : fast semi-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
+      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
+#else
+   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
+   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
+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
+
 
 !======================================================================
@@ -394 +462 @@
 implicit none
 
-#include "chimiedata.h"
+#include "chimiedata_early_earth.h"
 
 !----------------------------------------------------------------------
@@ -1796 +1864 @@
 implicit none
 
-#include "chimiedata.h"
+#include "chimiedata_early_earth.h"
 
 ! input
@@ -1925 +1993 @@
 implicit none
 
-#include "chimiedata.h"
+#include "chimiedata_early_earth.h"
 
 ! input
@@ -3249 +3317 @@
     (nb_reaction_4 /= nb_reaction_4_max)) then
    print*, 'wrong dimensions in indice' 
+   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
    stop
 end if  
@@ -3636 +3707 @@
       end subroutine chimtogcm
 
-end subroutine photochemistry_asis
-
+end subroutine photochemistry_asis_early_earth
+