!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com
    
! --------------------------------Subroutine part---------------------------------------------!
subroutine Iteration
  use constantparameters
  USE CONSTANTPARAMETERS
  USE ELECTROMAGNETIC_VARIABLES
  USE RES_MODEL_PARAMETER
  USE TIME_PARAMETER
  implicit none
  real::t1,t2,t                  !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
  REAL*8 CA,CB,DELX1,DELY1,DELZ1  !ca, cb, delx1, dely1, delz1 are all middle variables used in the computation of EM field
  REAL*8 TEMP_SIG,temp_cacb       !Temp_sig and temp_cacb are middle variables used in the computation of EM field
  REAL*8 DELY2,DELZ2,delx2        !They are all middle variables as above ones.
  integer num,i,j,k,ii     !num is the number of computation fraction
  real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq
  do num=1,num_fra_com,1                                   !The outer loop which begins from the first fraction ends at the last fraction
     call cpu_time(t1)                                     !Record the cpu time at the beginning of each computing fraction
     allocate(mdelt(0:mstop(num)),meps_r(mstop(num)),mcq(mstop(num)),msource(mstop(num)))
     ! The memory of mdelt, meps_r, mcq and msource are allocated at the begining of fraction
     do ii=mstart(num),mstart(num)+mstop(num)-1,1
        mdelt(ii-mstart(num)+1)=delt(ii)
        meps_r(ii-mstart(num)+1)=eps_r(ii)
        mcq(ii-mstart(num)+1)=cq(ii)
        msource(ii-mstart(num)+1)=source(ii) !Link the local value of mdelt, meps_r, mcq and msorce to the global value of delt, eps_r, cq and source array.
     end do
     print*,'Now computing fraction:',num
     mdelt(0)=mdelt(1)
     !$acc data copy(Ex(1:nx,1:nyb,1:nzb),Ey(1:nxb,1:ny,1:nzb),Ez(1:nxb,1:nyb,1:nz))&
     !$acc  copy(Hx(1:nxb,1:ny,0:nz),Hy(1:nx,1:nyb,0:nz),Hz(1:nx,1:ny,1:nzb)),copyin(cdelx(1:nx))&
     !$acc  copyin(ccsig(1:nx,1:ny,1:nz),mdelt(0:mstop(num)),cdely(1:ny),cdelz(1:nz),mcq(1:mstop(num)),meps_r(1:mstop(num)))&
     !$acc  copyin(is_ex_in_source(1:nx,2:nyb-1),is_ey_in_source(2:nx,1:ny),msource(1:mstop(num)))
     ! OpenACC directive, copy in and out of Ex,Ey,Ez,Hx,Hy,Hz, copy in ccsig, mdelt, cdelz, mcq, meps_r, is_ex_in_source, is_ey_in_source
     do loop=1,mstop(num),1
        ! --------------------------------update the value of Ex and Ey in source area---------------------------------------!
        !$acc parallel async(1)
        !$acc loop gang
        DO J=2,NYB-1
           !$acc loop vector
           DO I=1,NX
              K=NZ/2+1
              DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
              DELZ1=CDELZ(NZ/2+1)
              TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
                   &+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
                   &+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
                   &+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
              TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
              CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*TEMP_SIG)/(2.0*Meps_r(loop)+Mdelt(LOOP-1)*TEMP_SIG)
              CB=(2.0D0*MDELT(LOOP-1))/(2.0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
              EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)-cb*Msource(loop)*is_ex_in_source(i,j)
           ENDDO
        ENDDO                  
        !$acc end parallel
        ! end of updating Ex while k=Nzs+1
        ! update the value of Ey while k=Nzs+1
        !$acc parallel async(2)
        !$acc loop gang
        DO J=1,NY
           !$acc loop vector
           DO I=2,NX
              K=NZ/2+1
              DELX1=(CDELX(I-1)+CDELX(I))/2.0
              DELZ1=CDELZ(NZ/2+1)
              TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
                   &+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
                   &+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
                   &+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
              TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
              CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
              CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
              EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)-cb*Msource(loop)*is_ey_in_source(i,j)
           ENDDO
        ENDDO                 
        !$acc end parallel
        ! end of uptating Ey while k=Nzs+1
        ! ---------------------------------------------------Ex Part-------------------------------------------------------------!
        !$acc parallel async(3)
        !$acc loop gang
        DO K=NZ/2+2,NZ
           !$acc loop worker
           DO J=2,NY
              !$acc loop vector
              DO I=1,NX
                 DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
                 DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
                 TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
                      &+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
                      &+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
                      &+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
                 TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
                 CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !$acc parallel async(4)
        !$acc loop gang
        DO K=2,NZ/2
           !$acc loop worker
           DO J=2,NY
              !$acc loop vector
              DO I=1,NX
                 DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
                 DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
                 TEMP_SIG=CCSIG(I,J-1,K-1)*CDELY(J-1)*CDELZ(K-1)&
                      &+CCSIG(I,J-1,K)*CDELY(J-1)*CDELZ(K)&
                      &+CCSIG(I,J,K-1)*CDELY(J)*CDELZ(K-1)&
                      &+CCSIG(I,J,K)*CDELY(J)*CDELZ(K)
                 TEMP_SIG=TEMP_SIG/(4.0D0*DELY1*DELZ1)
                 CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))/DELY1-(HY(I,J,K)-HY(I,J,K-1))/DELZ1)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel  
        ! ================end of updating Ex==================!
        ! -----------------------------------------update the value of Ey--------------------------------!
        !$acc parallel async(5)
        !$acc loop gang
        DO K=NZ/2+2,NZ
           !$acc loop worker
           DO J=1,NY
              !$acc loop vector
              DO I=2,NX
                 DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
                 DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
                 TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
                      &+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
                      &+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
                      &+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
                 TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
                 CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !$acc parallel async(6)
        !$acc loop gang
        DO K=2,NZ/2
           !$acc loop worker
           DO J=1,NY
              !$acc loop vector
              DO I=2,NXB-1
                 DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
                 DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
                 TEMP_SIG=CCSIG(I-1,J,K-1)*CDELX(I-1)*CDELZ(K-1)&
                      &+CCSIG(I-1,J,K)*CDELX(I-1)*CDELZ(K)&
                      &+CCSIG(I,J,K-1)*CDELX(I)*CDELZ(K-1)&
                      &+CCSIG(I,J,K)*CDELX(I)*CDELZ(K)
                 TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELZ1)
                 CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG)
                 EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))/DELZ1-(HZ(I,J,K)-HZ(I-1,J,K))/DELX1)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !===============end of updating Ey===================!
        ! -------------------------------------update the value of Ez--------------------------------------!
        !$acc parallel async(7)
        !$acc loop gang
        DO K=1,NZ
           !$acc loop worker
           DO J=2,NYB-1
              !$acc loop vector
              DO I=2,NXB-1
                 DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
                 DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
                 TEMP_SIG=CCSIG(I-1,J-1,K)*CDELX(I-1)*CDELY(J-1)&
                      &+CCSIG(I-1,J,K)*CDELX(I-1)*CDELY(J)&
                      &+CCSIG(I,J-1,K)*CDELX(I)*CDELY(J-1)&
                      &+CCSIG(I,J,K)*CDELX(I)*CDELY(J)
                 TEMP_SIG=TEMP_SIG/(4.0D0*DELX1*DELY1)
                 TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*TEMP_SIG
                 CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*TEMP_SIG)/TEMP_CACB
                 CB=(2.0D0*MDELT(LOOP-1))/TEMP_CACB
                 EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))/DELX1-(HX(I,J,K)-HX(I,J-1,K))/DELY1)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !$acc wait
        !===============end of updating Ez=========================!
        ! ------------------------------------update the value of Hx-----------------------------------------------!
        !$acc parallel async(8)
        !$acc loop gang
        DO K=1,NZ
           !$acc loop worker
           DO J=1,NY
              !$acc loop vector
              DO I=1,NXB
                 DELY2=CDELY(J)
                 DELZ2=CDELZ(K)
                 HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))/DELY2-(EY(I,J,K+1)-EY(I,J,K))/DELZ2)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !================end of updating Hx=======================!
        ! -------------------------------------update the value of Hy---------------------------------------------!
        !$acc parallel async(9)
        !$acc loop gang
        DO K=1,NZ
           !$acc loop worker
           DO J=1,NYB
              !$acc loop vector
              DO I=1,NX
                 DELZ2=CDELZ(K)
                 DELX2=CDELX(I)
                 HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))/DELZ2-(EZ(I+1,J,K)-EZ(I,J,K))/DELX2)
              ENDDO
           ENDDO
        ENDDO
        !$acc end parallel
        !$acc wait
        !===============end of updating Hy========================!
        !-------------------------------------update the value of Hz----------------------------------------------!
        !$acc kernels async(10)
        DO J=1,NY
           DO I=1,NX
              DO K=NZ,NZ/2+1,-1       !NZ,2,-1   !
                 DELX2=CDELX(I)
                 DELY2=CDELY(J)
                 DELZ2=CDELZ(K)
                 HZ(I,J,K)=HZ(I,J,K+1)+DELZ2*((HX(I+1,J,K)-HX(I,J,K))/DELX2+(HY(I,J+1,K)-HY(I,J,K))/DELY2)
              ENDDO
           ENDDO
        ENDDO
        !$acc end kernels
        !$acc kernels async(11)
        DO K=1,NZ/2-1
           DO J=1,NY
              DO I=1,NX
                 DELX2=CDELX(I)
                 DELY2=CDELY(J)
                 DELZ2=CDELZ(K)
                 HZ(I,J,K+1)=HZ(I,J,K)-DELZ2*((HX(I+1,J,K)-HX(I,J,K))/DELX2+(HY(I,J+1,K)-HY(I,J,K))/DELY2)
              ENDDO
           ENDDO
        ENDDO
        !$acc end kernels
        !$acc wait
        !===================end of updating Hz==========================!
     enddo
     !$acc end data
     call cpu_time(t2)
     t=t2-t1
     print*,'The computing time for this fraction is:', t
     deallocate(meps_r,mcq,msource,mdelt)
     call WriteRecFiles(num)
     write(*,'(1x,e20.10e3,3x,e20.10e3)')Hz(nxs,nys+2,Nzs_air(1)),Hz(Nxs,Nys+2,Nz/2+1)
     write(*,*)'Now loop is:',mstart(num)+mstop(num)-1
     print*,mstop(num),'steps have just finished'
  ENDDO
end subroutine Iteration