MODULE mod_chem_spack_rodas3_dyndt USE mem_chem1, ONLY: & chem1_vars ! Type USE mem_spack, ONLY: & spack_type, & ! Type spack_type_2d ! Type USE Extras, ONLY: & ext3d ! Type ! USE solve_sparse, ONLY: & ! Solve_linear & ! Subroutine, OUT: spack_2d(ijk,inob)%DLk1 USE mod_chem_spack_ros, ONLY: & get_always_zeros, & ! Subroutine get_number_nonzeros ! Subroutine USE mod_chem_spack_jacdchemdc, ONLY: & jacdchemdc ! Subroutine USE mod_chem_spack_kinetic, ONLY: & kinetic ! Subroutine USE mod_chem_spack_fexchem, ONLY: & fexchem ! Subroutine IMPLICIT NONE PRIVATE PUBLIC chem_rodas3_dyndt CONTAINS !======================================================================================== SUBROUTINE chem_rodas3_dyndt(m1,m2,m3,dtlt,pp,pi0,theta,rv,dn0,cosz,rcp, & nspecies,nr,nr_photo,weight,PhotojMethod,maxnspecies, & nspecies_chem_transported,nspecies_chem_no_transported, & transp_chem_index,no_transp_chem_index,chem1_g,nob,maxblock_size, & block_end,indexk,indexi,indexj,kij_index,last_accepted_dt, & spack,spack_2d,atol,rtol,get_non_zeros,& N_DYN_CHEM,split_method,cp,cpor,p00,na_extra3d,extra3d,& jphoto,att) use mem_grid, only: time IMPLICIT NONE INTEGER , INTENT(IN) :: m1 INTEGER , INTENT(IN) :: m2 INTEGER , INTENT(IN) :: m3 ! mem_grid REAL , INTENT(IN) :: dtlt ! mem_basic REAL , INTENT(IN) :: pp(m1,m2,m3) REAL , INTENT(IN) :: pi0(m1,m2,m3) REAL , INTENT(IN) :: theta(m1,m2,m3) REAL , INTENT(IN) :: rv(m1,m2,m3) REAL , INTENT(IN) :: dn0(m1,m2,m3) ! mem_radiate REAL , INTENT(IN) :: cosz(m2,m3) ! mem_micro REAL , INTENT(IN) :: rcp(m1,m2,m3) ! chem1_list INTEGER , INTENT(IN) :: nspecies INTEGER , INTENT(IN) :: nr_photo INTEGER , INTENT(IN) :: nr REAL , INTENT(IN) :: weight(nspecies) CHARACTER(LEN=10), INTENT(IN) :: PhotojMethod ! mem_chem1 INTEGER , INTENT(IN) :: maxnspecies INTEGER , INTENT(IN) :: nspecies_chem_transported INTEGER , INTENT(IN) :: nspecies_chem_no_transported INTEGER , INTENT(IN) :: transp_chem_index(maxnspecies) INTEGER , INTENT(IN) :: no_transp_chem_index(maxnspecies) CHARACTER(LEN=20), INTENT(IN) :: split_method INTEGER , INTENT(IN) :: N_DYN_CHEM ! mem_chem1 TYPE (chem1_vars), INTENT(INOUT) :: chem1_g(nspecies) ! spack_utils INTEGER , INTENT(IN) :: nob INTEGER , INTENT(IN) :: maxblock_size INTEGER , INTENT(IN) :: block_end(:) INTEGER , INTENT(IN) :: indexk(:,:) INTEGER , INTENT(IN) :: indexi(:,:) INTEGER , INTENT(IN) :: indexj(:,:) INTEGER , INTENT(IN) :: kij_index(:,:) DOUBLE PRECISION, INTENT(INOUT) :: last_accepted_dt(:) ! INTEGER , INTENT(IN) :: block_end(nob) ! INTEGER , INTENT(IN) :: indexk(maxblock_size,nob) ! INTEGER , INTENT(IN) :: indexi(maxblock_size,nob) ! INTEGER , INTENT(IN) :: indexj(maxblock_size,nob) ! INTEGER , INTENT(IN) :: kij_index(maxblock_size,nob) ! DOUBLE PRECISION, INTENT(INOUT) :: last_accepted_dt(nob) ! mem_spack TYPE(spack_type) , INTENT(INOUT) :: spack(1) ! TYPE(spack_type_2d), INTENT(INOUT) :: spack_2d(maxblock_size,1) TYPE(spack_type_2d), INTENT(INOUT) :: spack_2d(:,:) DOUBLE PRECISION , INTENT(IN) :: atol(nspecies) DOUBLE PRECISION , INTENT(IN) :: rtol(nspecies) LOGICAL , INTENT(INOUT) :: get_non_zeros ! rconstants REAL, INTENT(IN) :: cp REAL, INTENT(IN) :: cpor REAL, INTENT(IN) :: p00 ! extra INTEGER, INTENT(IN) :: na_extra3d TYPE(ext3d), INTENT(INOUT) :: extra3d(na_extra3d) ! Output JNO2 ! FastJX DOUBLE PRECISION, INTENT(IN) :: jphoto(nr_photo,m1,m2,m3) ! uv_atten REAL, INTENT(IN) :: att(m2,m3) ! INTEGER,PARAMETER :: block_size=2*32 ! INTEGER,PARAMETER :: block_size=4 ! INTEGER,PARAMETER :: block_size=1 DOUBLE PRECISION,PARAMETER :: pmar=28.96D0 DOUBLE PRECISION,PARAMETER :: Threshold1= 0.0D0 DOUBLE PRECISION,PARAMETER :: Threshold2=-1.D-1 DOUBLE PRECISION,PARAMETER :: Igamma = 0.5D0 DOUBLE PRECISION,PARAMETER :: c43 = 4.D0/ 3.D0 DOUBLE PRECISION,PARAMETER :: c83 = 8.D0/ 3.D0 DOUBLE PRECISION,PARAMETER :: c56 = 5.D0/ 6.D0 DOUBLE PRECISION,PARAMETER :: c16 = 1.D0/ 6.D0 DOUBLE PRECISION,PARAMETER :: c112 = 1.D0/ 12.D0 DOUBLE PRECISION dble_dtlt_i,fxc, Igamma_dtstep,dt_chem,dt_chem_i INTEGER :: i,ijk,n,j,k,ispc,Ji,Jj,k_,i_,j_,kij_,kij !integer,parameter :: NOB_MEM=1 ! if 1 : alloc just one block (uses less memory) ! ! if 0 : alloc all blocks (uses more memory) !- default number of allocatable blocks (re-use scratch arrays spack()% ) integer,parameter :: inob=1 !- parameters for dynamic timestep control integer all_accepted double precision :: dt_min,dt_max,dt_actual,dt_new double precision :: time_f,time_c double precision, parameter :: FacMin =0.2d0 ! lower bound on step decrease factor (default=0.2) double precision, parameter :: FacMax =6.0D0 ! upper bound on step increase factor (default=6) double precision, parameter :: FacRej =0.1D0 ! step decrease factor after multiple rejections double precision, parameter :: FacSafe=0.9D0 ! step by which the new step is slightly smaller ! than the predicted value (default=0.9) double precision, parameter :: uround = 1.d-15, elo = 3.0D0, Roundoff = 1.d-8 double precision Fac, Tol,err1,max_err1 double precision, allocatable :: DLmat(:,:,:) integer, allocatable :: ipos(:) integer, allocatable :: jpos(:) integer :: NumberOfNonZeros, nz integer :: offset, offsetDLmat, offsetnz integer :: blocksize, sizeOfMatrix, maxnonzeros, blocknonzeros real :: start, finish real :: elapsed_time_solver,elapsed_time,elapsed_time_alloc,elapsed_time_dealloc,elapsed_time_copy INTEGER(KIND=8) :: matrix_Id INTEGER :: error INTEGER,PARAMETER :: complex_t=0 INTEGER(KIND=8),ALLOCATABLE,DIMENSION(:) :: element INTEGER(kind=8), EXTERNAL :: sfCreate_solve INTEGER(kind=8), EXTERNAL :: sfGetElement INTEGER, EXTERNAL :: sfFactor !! if(.not. spack_alloc) then !! !- Determining number of blocks and masks, allocating spack type !! call allocIndex(block_size,N_DYN_CHEM) !! call alloc_spack(CHEMISTRY) !! end if ! IF(.NOT. get_non_zeros) then !- Get always zero elements of Jacobian matrix to save !- computational time at solver: ! non-opt !get_non_zeros = .true. !call Prepare(nspecies) ! opt !!call get_always_zeros(nr_photo,nr,nspecies,get_non_zeros,spack(1)%rk(1,1:nr), & ! spack(1)%jphoto(1,1:nr_photo),p00) call get_number_nonzeros(nr_photo,nr,nspecies,spack(1)%rk(1,1:nr), & spack(1)%jphoto(1,1:nr_photo),p00,maxnonzeros) sizeOfMatrix = nspecies ALLOCATE(ipos(maxnonzeros));ipos=1 ALLOCATE(jpos(maxnonzeros));jpos=1 ! ENDIF DO i=1,nob !- loop over all blocks !- copying structure from input to internal DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i)!- kij brams tendency index k_=indexk(ijk,i) ! index_g%indexk(ijk,i) ! k brams index i_=indexi(ijk,i) ! index_g%indexi(ijk,i) ! i brams index j_=indexj(ijk,i) ! index_g%indexj(ijk,i) ! j brams index !- Exner function divided by Cp spack(inob)%press(ijk) = ( pp (k_,i_,j_) + pi0(k_,i_,j_) )/cp !- Air temperature (K) spack(inob)%temp(ijk)= theta(k_,i_,j_) * spack(inob)%press(ijk) !- transform from Exner function to pressure spack(inob)%press(ijk)= spack(inob)%press(ijk)**cpor*p00 !- Water vapor spack(inob)%vapp(ijk) = rv(k_,i_,j_)!& ! * dn0(k_,i_,j_)*6.02D20/18.D0 !- volmol e' usado para converter de ppbm para molecule/cm^3 (fator fxc incluido) ! spack(inob)%volmol(ijk)=(spack(inob)%press(ijk))/(8.314e0*spack(inob)%temp(ijk)) ! fxc = (6.02e23*1e-15*spack(inob)%volmol(ijk))*pmar ! spack(inob)%volmol(ijk)=(6.02D23*1D-15*pmar)*(spack(inob)%press(ijk))/(8.314D0*spack(inob)%temp(ijk)) !- inverse of volmol to get back to ppbm at the end of routine spack(inob)%volmol_I(ijk)=1.0d0 / spack(inob)%volmol(ijk) !- get liquid water content spack(inob)%xlw(ijk) = rcp(k_,i_,j_) * dn0(k_,i_,j_)*1.D-3 !- convert from brams chem (ppbm) arrays to spack (molec/cm3) !- no transported species section DO ispc=1,nspecies_chem_no_transported !- map the species to NO transported ones n=no_transp_chem_index(ispc) !- initialize no-transported species (don't need to convert, because these !- species are already saved using molecule/cm^3 units) spack(inob)%sc_p(ijk,n) = chem1_g(n)%sc_p(k_,i_,j_) ! spack(inob)%sc_p_4(ijk,n) = chem1_g(n)%sc_p(k_,i_,j_) END DO END DO !- convert from brams chem (ppbm) arrays to spack (molec/cm3) !- transported species section IF(split_method == 'PARALLEL' .and. N_DYN_CHEM > 1) then DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i)!- kij brams tendency index k_=indexk(ijk,i) !index_g%indexk(ijk,i) ! k brams index i_=indexi(ijk,i) !index_g%indexi(ijk,i) ! i brams index j_=indexj(ijk,i) !index_g%indexj(ijk,i) ! j brams index DO ispc=1,nspecies_chem_transported !- map the species to transported ones n=transp_chem_index(ispc) !- conversion from ppbm to molecule/cm^3 !- get back the concentration at the begin of the chemistry timestep integration: spack(inob)%sc_p(ijk,n) = ( chem1_g(n)%sc_p(k_,i_,j_) - & ! updated mixing ratio chem1_g(n)%sc_t_dyn(kij_)*N_DYN_CHEM*dtlt )& ! accumulated tendency * spack(inob)%volmol(ijk)/weight(n) !- for testing !spack(inob)%sc_p(ijk,n) = max(0.,spack(inob)%sc_p(ijk,n)) !spack(inob)%sc_p_4(ijk,n) =(chem1_g(n)%sc_p(k_,i_,j_) - beta*chem1_g(n)%sc_t_dyn(kij_)*N_DYN_CHEM*dtlt)& !dtlt \E9 o dinamico ! * spack(inob)%volmol(ijk)/weight(n) END DO END DO ELSE DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i)!- kij brams tendency index k_=indexk(ijk,i) !index_g%indexk(ijk,i) ! k brams index i_=indexi(ijk,i) !index_g%indexi(ijk,i) ! i brams index j_=indexj(ijk,i) !index_g%indexj(ijk,i) ! j brams index DO ispc=1,nspecies_chem_transported !- map the species to transported ones n=transp_chem_index(ispc) !- conversion from ppbm to molecule/cm^3 spack(inob)%sc_p(ijk,n) = chem1_g(n)%sc_p(k_,i_,j_) * spack(inob)%volmol(ijk)/weight(n) spack(inob)%sc_p(ijk,n) = max(0.D0,spack(inob)%sc_p(ijk,n)) END DO END DO ENDIF !- Photolysis section IF(trim(PhotojMethod) == 'FAST-JX' .or. trim(PhotojMethod) == 'FAST-TUV' ) THEN DO ijk=1,block_end(i) !index_g%block_end(i) k_=indexk(ijk,i) !index_g%indexk(ijk,i) ! k brams index i_=indexi(ijk,i) !index_g%indexi(ijk,i) ! i brams index j_=indexj(ijk,i) !index_g%indexj(ijk,i) ! j brams index DO n=1,nr_photo spack(inob)%jphoto(ijk,n)= jphoto(n,k_,i_,j_) !fast_JX_g%jphoto(n,k_,i_,j_) END DO ENDDO ELSEIF(trim(PhotojMethod) == 'LUT') then DO ijk=1,block_end(i) !index_g%block_end(i) i_=indexi(ijk,i) !index_g%indexi(ijk,i) ! i brams index j_=indexj(ijk,i) !index_g%indexj(ijk,i) ! j brams index !- UV attenuation un function AOT spack(inob)%att(ijk)= att(i_,j_) !uv_atten_g%att(i_,j_) !- get zenital angle (for LUT PhotojMethod) spack(inob)%cosz(ijk)= cosz(i_,j_) ENDDO ENDIF !- compute kinetical and photochemical reactions (array: spack(inob)%rk) CALL kinetic(nr_photo,spack(inob)%Jphoto & ,spack(inob)%rk & ,spack(inob)%temp & ,spack(inob)%vapp & ,spack(inob)%Press & ,spack(inob)%cosz & ,spack(inob)%att & ,1 & ,block_end(i) & !index_g%block_end(i), & ,maxblock_size,nr) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !tmp ! ! !- srf: to use these arrays, be sure NA_EXTRAD3D RAMSIN is at least 8 ! DO ijk=1,block_end(i) !index_g%block_end(i) ! k_=indexk(ijk,i) !index_g%indexk(ijk,i) ! k brams index ! i_=indexi(ijk,i) !index_g%indexi(ijk,i) ! i brams index ! j_=indexj(ijk,i) !index_g%indexj(ijk,i) ! j brams index ! extra3d(7)%d3(k_,i_,j_) =spack(inob)%rk(ijk,1) !extra3d(8)%d3(k_,i_,j_) =spack(inob)%rk(ijk,2) !extra3d(9)%d3(k_,i_,j_) =spack(inob)%rk(ijk,3) !extra3d(10)%d3(k_,i_,j_)=spack(inob)%rk(ijk,4) !extra3d(11)%d3(k_,i_,j_)=spack(inob)%rk(ijk,5) !extra3d(12)%d3(k_,i_,j_)=spack(inob)%rk(ijk,6) !extra3d(13)%d3(k_,i_,j_)=spack(inob)%rk(ijk,7) !extra3d(14)%d3(k_,i_,j_)=spack(inob)%rk(ijk,8) !extra3d(15)%d3(k_,i_,j_)=spack(inob)%rk(ijk,9) !extra3d(16)%d3(k_,i_,j_)=spack(inob)%rk(ijk,10) !extra3d(17)%d3(k_,i_,j_)=spack(inob)%rk(ijk,11) !extra3d(18)%d3(k_,i_,j_)=spack(inob)%rk(ijk,12) !extra3d(19)%d3(k_,i_,j_)=spack(inob)%rk(ijk,13) !extra3d(20)%d3(k_,i_,j_)=spack(inob)%rk(ijk,14) !extra3d(21)%d3(k_,i_,j_)=spack(inob)%rk(ijk,16) !extra3d(22)%d3(k_,i_,j_)=spack(inob)%rk(ijk,17) !extra3d(23)%d3(k_,i_,j_)=spack(inob)%rk(ijk,18) !extra3d(24)%d3(k_,i_,j_)=spack(inob)%rk(ijk,19) !extra3d(25)%d3(k_,i_,j_)=spack(inob)%rk(ijk,20) !extra3d(26)%d3(k_,i_,j_)=spack(inob)%rk(ijk,21) ! extra3d(27)%d3(k_,i_,j_)=spack(inob)%rk(ijk,22) ! extra3d(28)%d3(k_,i_,j_)=spack(inob)%rk(ijk,23) ! ENDDO !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !tmp !- ROSENBROCK METHOD ---------------------------------------------------------------------- !- kml: Inicio do equivalente a roschem !- srf: extending to RODAS 3 (4 stages, order 3) dt_chem = last_accepted_dt(i) !index_g%last_accepted_dt(i) ! dble(dtlt) dt_min = max(1.0D0,1.D-2*dble(dtlt*N_DYN_CHEM)) dt_max = dble(dtlt*N_DYN_CHEM) dt_new = 0.0d0 time_c = 0.0d0 time_f = dble(dtlt*N_DYN_CHEM) run_until_integr_ends: DO WHILE (time_c + Roundoff < time_f) !- Compute the Jacobian (DLRDC). CALL jacdchemdc (spack(inob)%sc_p & ! CALL jacdchemdc (spack(inob)%sc_p_4 & ,spack(inob)%rk & ,spack(inob)%DLdrdc& ! Jacobian matrix ,nspecies,1 & ,block_end(i) & !index_g%block_end(i), & ,maxblock_size,nr) !- Compute chemical net production terms (actually, P-L) at initial time (array spack(inob)%DLr) CALL fexchem (spack(inob)%sc_p & ! CALL fexchem (spack(inob)%sc_p_4 & ,spack(inob)%rk & ,spack(inob)%DLr & !production term ,nspecies,1 & ,block_end(i) & !index_g%block_end(i), & ,maxblock_size,nr) DO Ji=1,nspecies DO ijk=1,block_end(i) !index_g%block_end(i) spack_2d(ijk,inob)%DLb1(Ji) = spack(inob)%DLr(ijk,Ji) ENDDO ENDDO !- compute matrix (1/(Igamma*dt) - Jacobian) where Jacobian = DLdrdc !- fill DLMAT with non-diagonal (fixed in the timestep) Jacobian DO Jj=1,nspecies DO Ji=1,nspecies DO ijk=1,block_end(i) !index_g%block_end(i) spack_2d(ijk,inob)%DLmat(Ji,Jj) = - spack(inob)%DLdrdc(ijk,Ji,Jj) ENDDO ENDDO ENDDO UntilAccepted: DO !- fill DLMAT with diagonal Jacobian (which changes in time, because of dt_chem) Igamma_dtstep = 1.0D0/(Igamma * dt_chem) DO Jj=1,nspecies DO ijk=1,block_end(i) !index_g%block_end(i) spack_2d(ijk,inob)%DLmat(Jj,Jj) = Igamma_dtstep - spack(inob)%DLdrdc(ijk,Jj,Jj) ENDDO ENDDO if ( (i .eq. 1) ) then blocknonzeros = 0 DO Ji=1,nspecies DO Jj=1,nspecies !IF(spack_2d(1,inob)%DLmat(Ji,Jj)<=0 .AND. spack_2d(1,inob)%DLmat(Ji,Jj)>=(-1)*0) CYCLE if ( (spack_2d(1,inob)%DLmat(Ji,Jj) .ne. 0.0d+00) ) then blocknonzeros = blocknonzeros + 1 ipos(blocknonzeros)=Ji jpos(blocknonzeros)=Jj endif ENDDO ENDDO NumberOfNonZeros = blocknonzeros !create matrix matrix_Id = sfCreate_solve(sizeOfMatrix,complex_t,error) IF(ALLOCATED(element)) DEALLOCATE (element) ALLOCATE(element(NumberOfNonZeros)) DO nz=1,NumberOfNonZeros Ji = ipos(nz) Jj = jpos(nz) element(nz)=sfGetElement(matrix_Id,Ji,Jj) !element(nz)=sfGetElement(matrix_Id,INT(DLmat(1,nz,1)),INT(DLmat(1,nz,2))) END DO CALL sfZero(matrix_Id) DO nz=1,NumberOfNonZeros Ji = ipos(nz) Jj = jpos(nz) CALL sfAdd1Real(element(nz),spack_2d(1,inob)%DLmat(Ji,Jj)) !CALL sfAdd1Real(element(nz),DLmat(1,nz,3)) END DO error = sfFactor(matrix_Id) ENDIF !@LNCC: begin points loop DO ijk=1,block_end(i) !index_g%block_end(i) CALL sfZero(matrix_Id) DO nz=1,NumberOfNonZeros Ji = ipos(nz) Jj = jpos(nz) CALL sfAdd1Real(element(nz),spack_2d(ijk,inob)%DLmat(Ji,Jj)) !CALL sfAdd1Real(element(nz),DLmat(ijk,nz,3)) END DO error = sfFactor(matrix_Id) !--------------------------------------------------------------------------------------------------------------------- ! 1- First step !- Compute DLk1 by Solving (1/(Igamma*dt) - DLRDC) DLk1=DLR !- solver sparse call sfSolve(matrix_Id, spack_2d(ijk,inob)%DLb1, spack_2d(ijk,inob)%DLk1) ! !--------------------------------------------------------------------------------------------------------------------- ! 2- Second step ! compute K2 by solving (1/0.5 h JAC)K2 = 4/h * K1 + F(Yn) ! compute DLK2 by solving (1/Igama*dt - DLRDC)DLK2 = 4/h * DLK1 + DLb1 dt_chem_i = 1.0d0/dt_chem DO Ji=1,nspecies spack_2d(ijk,inob)%DLb2(Ji) = (4.0D0*dt_chem_i)*spack_2d(ijk,inob)%DLk1(Ji) + & spack_2d(ijk,inob)%DLb1(Ji) ENDDO call sfSolve(matrix_Id, spack_2d(ijk,inob)%DLb2, spack_2d(ijk,inob)%DLk2) !--------------------------------------------------------------------------------------------------------------------- ! 3- Third step ! a) update concentrations !dt_chem_i = 1.0d0/dt_chem DO Ji=1,nspecies spack(inob)%sc_p_new(ijk,Ji) = spack(inob)%sc_p(ijk,Ji) + 2.0D0* spack_2d(ijk,inob)%DLk1(Ji) IF (spack(inob)%sc_p_new(ijk,Ji) .LT. Threshold1) THEN spack(inob)%sc_p_new(ijk,Ji) = Threshold1 spack_2d(ijk,inob)%DLk1(Ji) = 0.5D0*(spack(inob)%sc_p_new(ijk,Ji) - spack(inob)%sc_p(ijk,Ji)) ENDIF ENDDO ! ! b) update the net production term (DLr= P-L = F(Y3)) at this stage with the first-order ! approximation with the new concentration ! CALL fexchem (spack(inob)%sc_p_new & ,spack(inob)%rk & ,spack(inob)%DLr & ,nspecies,ijk & ,ijk & !index_g%block_end(i), & ,maxblock_size,nr) ! c) compute K3 by solving (1 /(0.5 h) - JAC)K3 = F(Y3) + 0.5 (K1-K2) DO Ji=1,nspecies spack_2d(ijk,inob)%DLb3(Ji) = spack(inob)%DLr(ijk,Ji) + dt_chem_i* & ( spack_2d(ijk,inob)%DLk1(Ji) - spack_2d(ijk,inob)%DLk2(Ji)) ENDDO call sfSolve(matrix_Id, spack_2d(ijk,inob)%DLb3, spack_2d(ijk,inob)%DLk3) !--------------------------------------------------------------------------------------------------------------------- ! 4- Fourth step ! a) update concentrations ! Y4 = Yn + 2 * k1 + K3 dt_chem_i = 1.0d0/dt_chem DO Ji=1,nspecies spack(inob)%sc_p_new(ijk,Ji) = spack(inob)%sc_p(ijk,Ji) + & 2.0D0 * spack_2d(ijk,inob)%DLk1(Ji) + & spack_2d(ijk,inob)%DLk3(Ji) IF (spack(inob)%sc_p_new(ijk,Ji) .LT. Threshold1) THEN spack(inob)%sc_p_new(ijk,Ji) = Threshold1 spack_2d(ijk,inob)%DLk3(Ji) = (spack(inob)%sc_p_new(ijk,Ji) - spack(inob)%sc_p(ijk,Ji)) & -2.0D0 * spack_2d(ijk,inob)%DLk1(Ji) ENDIF ENDDO ! b) update the net production term (DLr= P-L = F(Y4) ) at this stage with the 3rd-order ! approximation with the new concentration ! CALL fexchem (spack(inob)%sc_p_new & ! Y4 ,spack(inob)%rk & ,spack(inob)%DLr & ! F(Y4) ,nspecies,ijk & ,ijk & !index_g%block_end(i) ,maxblock_size,nr) ! c) compute K4 by solving (1/(0.5 h)- JAC)K4 = F(Y3) + K1/h -K2/h -8/3 K3/h DO Ji=1,nspecies spack_2d(ijk,inob)%DLb4(Ji) = spack(inob)%DLr (ijk,Ji) + dt_chem_i* & ! F(Y4) ( spack_2d(ijk,inob)%DLk1(Ji) & - spack_2d(ijk,inob)%DLk2(Ji) & - c83 * spack_2d(ijk,inob)%DLk3(Ji) ) ENDDO call sfSolve(matrix_Id, spack_2d(ijk,inob)%DLb4, spack_2d(ijk,inob)%DLk4) ENDDO !@LNCC: end points loop !--------------------------------------------------------------------------------------------------------------------- ! - the solution dt_chem_i = 1.0d0/dt_chem DO Ji=1,nspecies DO ijk=1,block_end(i) !index_g%block_end(i) spack(inob)%sc_p_new(ijk,Ji) = spack(inob)%sc_p(ijk,Ji) + & 2.0D0 * spack_2d(ijk,inob)%DLk1(Ji) & + spack_2d(ijk,inob)%DLk3(Ji) & + spack_2d(ijk,inob)%DLk4(Ji) spack(inob)%sc_p_new(ijk,Ji) = max (spack(inob)%sc_p_new(ijk,Ji),Threshold1) ! IF (spack(inob)%sc_p_new(ijk,Ji) .LT. Threshold1) THEN ! spack(inob)%sc_p_new(ijk,Ji) = Threshold1 ! spack_2d(ijk,inob)%DLk3(Ji) = (spack(inob)%sc_p_new(ijk,Ji) - spack(inob)%sc_p(ijk,Ji)) * dt_chem_i ! ENDIF !- keep non-negative values for concentrations ! spack(inob)%sc_p_new(ijk,Ji) = max(spack(inob)%sc_p_new(ijk,Ji),Threshold1) !- update DLK2, in case sc_p_new is changed by the statement above ! spack_2d(ijk,inob)%DLk2(Ji) = 2.0D0*( (spack(inob)%sc_p_new(ijk,Ji)-spack(inob)%sc_p(ijk,Ji)) & ! * dt_chem_i - 1.5D0 * spack_2d(ijk,inob)%DLk1(Ji) ) ENDDO ENDDO ! !- Compute the error estimation : spack(inob)%err DO ijk=1,block_end(i) !index_g%block_end(i) spack(inob)%err(ijk) = 0.0d0 DO Ji=1,nspecies Tol = ATol(Ji) + RTol(Ji)*DMAX1( DABS(spack(inob)%sc_p(ijk,Ji)),DABS(spack(inob)%sc_p_new(ijk,Ji)) ) err1= spack_2d(ijk,inob)%DLk4(Ji) spack(inob)%err(ijk) = spack(inob)%err(ijk) + (err1/Tol)**2.0D0 ENDDO spack(inob)%err(ijk) = DMAX1( uround, DSQRT( spack(inob)%err(ijk)/nspecies ) ) ENDDO all_accepted = 1 DO ijk=1,block_end(i) !index_g%block_end(i) if( spack(inob)%err(ijk) - Roundoff > 1.0D0) then all_accepted = 0 ; exit endif ENDDO !- find the maximum error occurred max_err1 = maxval( spack(inob)%err(1:block_end(i))) !index_g%block_end(i) ) ) !- use it to determine the new time step for all block elements !- new step size is bounded by FacMin <= Hnew/H <= FacMax Fac = MIN(FacMax,MAX(FacMin,FacSafe/max_err1**(1.0D0/elo))) !- possible new timestep dt_new = dt_chem * Fac dt_new = max(dt_min,min(dt_max,dt_new)) !- to reset the timestep resizing in function of the error estimation, use the statements below: ! all_accepted = 1; dt_new=dt_max if( all_accepted == 0 .and. dt_new > dt_min) then ! current solution is not accepted !- resize the timestep and try again dt_chem = dt_new else ! current solution is accepted !- go ahead, updating spack(inob)%sc_p with the solution (spack(inob)%sc_p_new) !- next time time_c = time_c + dt_chem !- next timestep (dt_new but limited by the time_f-time_c, the rest of time integration interval) dt_chem = min(dt_new, time_f-time_c) !- save the accepted timestep for the next integration interval if(time_c < time_f) last_accepted_dt(i) = dt_new !index_g%last_accepted_dt(i) = dt_new !- pointer (does not work yet) ! spack(inob)%sc_p=>spack(inob)%sc_p_new ! POINTER !- copy DO Ji=1,nspecies DO ijk=1,block_end(i) !index_g%block_end(i) spack(inob)%sc_p(ijk,Ji) = spack(inob)%sc_p_new(ijk,Ji) ! spack(inob)%sc_p_4(ijk,Ji) = spack(inob)%sc_p_new(ijk,Ji) ENDDO ENDDO EXIT UntilAccepted endif END DO UntilAccepted ENDDO run_until_integr_ends ! time-spliting !-------------------------------------------------------------------------------------------- !- Restoring species tendencies OR updated mixing ratios from internal to brams structure !- transported species section if(split_method == 'PARALLEL' .and. N_DYN_CHEM > 1) then DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i) k_ =indexk(ijk,i) !index_g%indexk(ijk,i) i_ =indexi(ijk,i) !index_g%indexi(ijk,i) j_ =indexj(ijk,i) !index_g%indexj(ijk,i) DO ispc=1,nspecies_chem_transported !- map the species to transported ones n=transp_chem_index(ispc) !- include the chemical tendency at total tendency (convert to unit: ppbm/s) chem1_g(n)%sc_p(k_,i_,j_ ) = chem1_g(n)%sc_t_dyn(kij_)*N_DYN_CHEM*dtlt +& spack(inob)%sc_p(ijk,n)*weight(n)*spack(inob)%volmol_I(ijk) chem1_g(n)%sc_p(k_,i_,j_ ) = max(0.,chem1_g(n)%sc_p(k_,i_,j_ )) END DO ENDDO ELSEIF(split_method == 'PARALLEL' .and. N_DYN_CHEM == 1) then dble_dtlt_i=1.0d0 / dble( dtlt) DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i) k_ =indexk(ijk,i) !index_g%indexk(ijk,i) i_ =indexi(ijk,i) !index_g%indexi(ijk,i) j_ =indexj(ijk,i) !index_g%indexj(ijk,i) DO ispc=1,nspecies_chem_transported !- map the species to transported ones n=transp_chem_index(ispc) !- include the chemical tendency at total tendency (convert to unit: ppbm/s) ! chem1_g(n)%sc_t(kij_) = + &! use this for update only chemistry (No dyn/emissions) chem1_g(n)%sc_t(kij_) = chem1_g(n)%sc_t(kij_) + &! previous tendency ( spack(inob)%sc_p(ijk,n)*weight(n)*spack(inob)%volmol_I(ijk) - &! new mixing ratio chem1_g(n)%sc_p(k_,i_,j_ ) ) &! old mixing ratio * dble_dtlt_i ! inverse of timestep END DO END DO ELSE DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i) k_ =indexk(ijk,i) !index_g%indexk(ijk,i) i_ =indexi(ijk,i) !index_g%indexi(ijk,i) j_ =indexj(ijk,i) !index_g%indexj(ijk,i) DO ispc=1,nspecies_chem_transported !- map the species to transported ones n=transp_chem_index(ispc) chem1_g(n)%sc_p(k_,i_,j_ ) = spack(inob)%sc_p(ijk,n)*weight(n)*spack(inob)%volmol_I(ijk) chem1_g(n)%sc_p(k_,i_,j_ ) = max(0.,chem1_g(n)%sc_p(k_,i_,j_ )) END DO END DO ENDIF !- no transported species section DO ijk=1,block_end(i) !index_g%block_end(i) kij_ =kij_index(ijk,i) !index_g%kij_index(ijk,i) k_ =indexk(ijk,i) !index_g%indexk(ijk,i) i_ =indexi(ijk,i) !index_g%indexi(ijk,i) j_ =indexj(ijk,i) !index_g%indexj(ijk,i) DO ispc=1,nspecies_chem_no_transported !- map the species to no transported ones n=no_transp_chem_index(ispc) !- save no-transported species (keep current unit : molec/cm3) chem1_g(n)%sc_p(k_,i_,j_ ) = max(0., real ( spack(inob)%sc_p(ijk,n) )) END DO END DO END DO ! enddo loop over all blocks !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX !DO n=1,nspecies !print*,'max1=',spc_name(n),maxval(chem1_g(n)%sc_p(:,:,:)),maxloc(chem1_g(n)%sc_p(:,:,:)) !print*,'max2=',n,spc_name(n),maxval(spack(:)%sc_p(:,n)),maxloc(spack(:)%sc_p(:,n)) !call flush(6) !END DO !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX END SUBROUTINE chem_rodas3_dyndt !--------------------------------------------------------------------------------------------------- END MODULE mod_chem_spack_rodas3_dyndt