Index: trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/calchim_asis_early_earth.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/calchim_asis_early_earth.F90 (revision 1807) +++ trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/calchim_asis_early_earth.F90 (revision 1810) @@ -70,5 +70,5 @@ !======================================================================= -#include "chimiedata.h" +#include "chimiedata_early_earth.h" ! input: @@ -159,5 +159,5 @@ 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,16 +543,16 @@ ! 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 Index: trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/chimiedata_early_earth.h =================================================================== --- trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/chimiedata_early_earth.h (revision 1807) +++ trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/chimiedata_early_earth.h (revision 1810) @@ -19,5 +19,5 @@ 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) Index: trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/photochemistry_asis_early_earth.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/photochemistry_asis_early_earth.F90 (revision 1807) +++ trunk/LMDZ.GENERIC/libf/aeronostd/chemistry_early_earth/photochemistry_asis_early_earth.F90 (revision 1810) @@ -19,5 +19,5 @@ implicit none -#include "chimiedata.h" +#include "chimiedata_early_earth.h" !=================================================================== @@ -178,5 +178,5 @@ !=================================================================== -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,21 +249,35 @@ 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,75 +286,5 @@ 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,4 +300,7 @@ end do ! ilev + + + !=================================================================== @@ -370,4 +317,125 @@ 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,5 +462,5 @@ implicit none -#include "chimiedata.h" +#include "chimiedata_early_earth.h" !---------------------------------------------------------------------- @@ -1796,5 +1864,5 @@ implicit none -#include "chimiedata.h" +#include "chimiedata_early_earth.h" ! input @@ -1925,5 +1993,5 @@ implicit none -#include "chimiedata.h" +#include "chimiedata_early_earth.h" ! input @@ -3249,4 +3317,7 @@ (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,4 +3707,4 @@ end subroutine chimtogcm -end subroutine photochemistry_asis - +end subroutine photochemistry_asis_early_earth +