((nxPML_1 - i) / (nxPML_1 - 1.0))**m
      be_x_1(i) = EXP(-(sige_x_PML_1(i) / kappae_x_PML_1(i) +   &
                                 alphae_x_PML_1(i))*dt/epsO)
      if ((sige_x_PML_1(i) == 0.0) .and.        &
         (alphae_x_PML_1(i) == 0.0) .and. (i == nxPML_1)) then
         ce_x_1(i) = 0.0
      else
         ce_x_1(i) = sige_x_PML_1(i)*(be_x_1(i)-1.0)/       &
            (sige_x_PML_1(i)+kappae_x_PML_1(i)*alphae_x_PML_1(i)) &
            / kappae_x_PML_1(i)
      endif
   ENDDO
   DO i = 1,nxPML_1-1
      sigh_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i - 0.5)/(nxPML_1-1.0))**m
      alphah_x_PML_1(i) = alpha_x_max*((i-0.5)/(nxPML_1-1.0))**ma
      kappah_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   &
                            ((nxPML_1 - i - 0.5) / (nxPML_1 - 1.0))**m
      bh_x_1(i) = EXP(-(sigh_x_PML_1(i) / kappah_x_PML_1(i) +   &
                                 alphah_x_PML_1(i))*dt/epsO)
      ch_x_1(i) = sigh_x_PML_1(i)*(bh_x_1(i)-1.0)/      &
                  (sigh_x_PML_1(i)+kappah_x_PML_1(i)*alphah_x_PML_1(i)) &
                  / kappah_x_PML_1(i)
   ENDDO

   DO i = 1,nxPML_2
      sige_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i) / (nxPML_2 - 1.0) )**m
      alphae_x_PML_2(i) = alpha_x_max*((i-1.0)/(nxPML_2-1.0))**ma
      kappae_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   &
                                 ((nxPML_2 - i) / (nxPML_2 - 1.0))**m
      be_x_2(i) = EXP(-(sige_x_PML_2(i) / kappae_x_PML_2(i) +   &
                                 alphae_x_PML_2(i))*dt/epsO)
      if ((sige_x_PML_2(i) == 0.0) .and.        &
         (alphae_x_PML_2(i) == 0.0) .and. (i == nxPML_2)) then
         ce_x_2(i) = 0.0
      else
         ce_x_2(i) = sige_x_PML_2(i)*(be_x_2(i)-1.0)/       &
            (sige_x_PML_2(i)+kappae_x_PML_2(i)*alphae_x_PML_2(i)) &
            / kappae_x_PML_2(i)
      endif
   ENDDO
   DO i = 1,nxPML_2-1
      sigh_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i - 0.5)/(nxPML_2-1.0))**m
      alphah_x_PML_2(i) = alpha_x_max*((i-0.5)/(nxPML_2-1.0))**ma
      kappah_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   &
                            ((nxPML_2 - i - 0.5) / (nxPML_2 - 1.0))**m
      bh_x_2(i) = EXP(-(sigh_x_PML_2(i) / kappah_x_PML_2(i) +   &
                                 alphah_x_PML_2(i))*dt/epsO)
      ch_x_2(i) = sigh_x_PML_2(i)*(bh_x_2(i)-1.0)/      &
                  (sigh_x_PML_2(i)+kappah_x_PML_2(i)*alphah_x_PML_2(i)) &
                  / kappah_x_PML_2(i)
   ENDDO

   DO j = 1,nyPML_1
      sige_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j ) / (nyPML_1 - 1.0) )**m
      alphae_y_PML_1(j) = alpha_y_max*((j-1)/(nyPML_1-1.0))**ma
      kappae_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   &
                                 ((nyPML_1 - j) / (nyPML_1 - 1.0))**m
      be_y_1(j) = EXP(-(sige_y_PML_1(j) / kappae_y_PML_1(j) +   &
                                 alphae_y_PML_1(j))*dt/epsO)
      if ((sige_y_PML_1(j) == 0.0) .and.        &
         (alphae_y_PML_1(j) == 0.0) .and. (j == nyPML_1)) then
         ce_y_1(j) = 0.0
      else
         ce_y_1(j) = sige_y_PML_1(j)*(be_y_1(j)-1.0)/       &
            (sige_y_PML_1(j)+kappae_y_PML_1(j)*alphae_y_PML_1(j)) &
            / kappae_y_PML_1(j)
      endif
   ENDDO
   DO j = 1,nyPML_1-1
      sigh_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j - 0.5)/(nyPML_1-1.0))**m
      alphah_y_PML_1(j) = alpha_y_max*((j-0.5)/(nyPML_1-1.0))**ma
      kappah_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   &
                            ((nyPML_1 - j - 0.5) / (nyPML_1 - 1.0))**m
      bh_y_1(j) = EXP(-(sigh_y_PML_1(j) / kappah_y_PML_1(j) +   &
                                 alphah_y_PML_1(j))*dt/epsO)
      ch_y_1(j) = sigh_y_PML_1(j)*(bh_y_1(j)-1.0)/      &
                  (sigh_y_PML_1(j)+kappah_y_PML_1(j)*alphah_y_PML_1(j)) &
                  / kappah_y_PML_1(j)
   ENDDO
   DO j = 1,nyPML_2
      sige_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j ) / (nyPML_2 - 1.0) )**m
      alphae_y_PML_2(j) = alpha_y_max*((j-1)/(nyPML_2-1.0))**ma
      kappae_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   &
                                 ((nyPML_2 - j) / (nyPML_2 - 1.0))**m
      be_y_2(j) = EXP(-(sige_y_PML_2(j) / kappae_y_PML_2(j) +   &
                                 alphae_y_PML_2(j))*dt/epsO)
      if ((sige_y_PML_2(j) == 0.0) .and.        &
         (alphae_y_PML_2(j) == 0.0) .and. (j == nyPML_2)) then
         ce_y_2(j) = 0.0
      else
         ce_y_2(j) = sige_y_PML_2(j)*(be_y_2(j)-1.0)/       &
            (sige_y_PML_2(j)+kappae_y_PML_2(j)*alphae_y_PML_2(j)) &
            / kappae_y_PML_2(j)
      endif
   ENDDO
   DO j = 1,nyPML_2-1
      sigh_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j - 0.5)/(nyPML_2-1.0))**m
      alphah_y_PML_2(j) = alpha_y_max*((j-0.5)/(nyPML_2-1.0))**ma
      kappah_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   &
                            ((nyPML_2 - j - 0.5) / (nyPML_2 - 1.0))**m
      bh_y_2(j) = EXP(-(sigh_y_PML_2(j) / kappah_y_PML_2(j) +   &
                                 alphah_y_PML_2(j))*dt/epsO)
      ch_y_2(j) = sigh_y_PML_2(j)*(bh_y_2(j)-1.0)/      &
                  (sigh_y_PML_2(j)+kappah_y_PML_2(j)*alphah_y_PML_2(j)) &
                  / kappah_y_PML_2(j)
   ENDDO

   DO k = 1,nzPML_1
      sige_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k ) / (nzPML_1 - 1.0) )**m
      alphae_z_PML_1(k) = alpha_z_max*((k-1)/(nzPML_1-1.0))**ma
      kappae_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   &
                                 ((nzPML_1 - k) / (nzPML_1 - 1.0))**m
      be_z_1(k) = EXP(-(sige_z_PML_1(k) / kappae_z_PML_1(k) +   &
                                 alphae_z_PML_1(k))*dt/epsO)
      if ((sige_z_PML_1(k) == 0.0) .and.        &
         (alphae_z_PML_1(k) == 0.0) .and. (k == nzPML_1)) then
         ce_z_1(k) = 0.0
      else
         ce_z_1(k) = sige_z_PML_1(k)*(be_z_1(k)-1.0)/       &
            (sige_z_PML_1(k)+kappae_z_PML_1(k)*alphae_z_PML_1(k)) &
            / kappae_z_PML_1(k)
      endif
   ENDDO
   DO k = 1,nzPML_1-1
      sigh_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k - 0.5)/(nzPML_1-1.0))**m
      alphah_z_PML_1(k) = alpha_z_max*((k-0.5)/(nzPML_1-1.0))**ma
      kappah_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   &
                            ((nzPML_1 - k - 0.5) / (nzPML_1 - 1.0))**m
      bh_z_1(k) = EXP(-(sigh_z_PML_1(k) / kappah_z_PML_1(k) +   &
                                 alphah_z_PML_1(k))*dt/epsO)
      ch_z_1(k) = sigh_z_PML_1(k)*(bh_z_1(k)-1.0)/      &
                  (sigh_z_PML_1(k)+kappah_z_PML_1(k)*alphah_z_PML_1(k)) &
                  / kappah_z_PML_1(k)
   ENDDO

   DO k = 1,nzPML_2
      sige_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k ) / (nzPML_2 - 1.0) )**m
      alphae_z_PML_2(k) = alpha_z_max*((k-1)/(nzPML_2-1.0))**ma
      kappae_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   &
                                 ((nzPML_2 - k) / (nzPML_2 - 1.0))**m
      be_z_2(k) = EXP(-(sige_z_PML_2(k) / kappae_z_PML_2(k) +   &
                                 alphae_z_PML_2(k))*dt/epsO)
      if ((sige_z_PML_2(k) == 0.0) .and.        &
         (alphae_z_PML_2(k) == 0.0) .and. (k == nzPML_2)) then
         ce_z_2(k) = 0.0
      else
         ce_z_2(k) = sige_z_PML_2(k)*(be_z_2(k)-1.0)/       &
            (sige_z_PML_2(k)+kappae_z_PML_2(k)*alphae_z_PML_2(k)) &
            / kappae_z_PML_2(k)
      endif
   ENDDO
   DO k = 1,nzPML_2-1
      sigh_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k - 0.5)/(nzPML_2-1.0))**m
      alphah_z_PML_2(k) = alpha_z_max*((k-0.5)/(nzPML_2-1.0))**ma
      kappah_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   &
                            ((nzPML_2 - k - 0.5) / (nzPML_2 - 1.0))**m
      bh_z_2(k) = EXP(-(sigh_z_PML_2(k) / kappah_z_PML_2(k) +   &
                                 alphah_z_PML_2(k))*dt/epsO)
      ch_z_2(k) = sigh_z_PML_2(k)*(bh_z_2(k)-1.0)/      &
                  (sigh_z_PML_2(k)+kappah_z_PML_2(k)*alphah_z_PML_2(k)) &
                  / kappah_z_PML_2(k)
   ENDDO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  FILL IN UPDATING COEFFICIENTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DA = 1.0
   DB = (dt/(muO)) 

   DO i = 1,Imax
      DO j = 1,Jmax
         DO k = 1,Kmax
            CA(i,j,k) = (1.0 - sig(i,j,k)*dt / (2.0*eps(i,j,k))) /        &
               (1.0 + sig(i,j,k) * dt / (2.0*eps(i,j,k)))
            CB(i,j,k) = (dt/(eps(i,j,k))) /                            &
               (1.0 + sig(i,j,k)*dt / (2.0*eps(i,j,k)))
            ENDDO
	ENDDO
   ENDDO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  FILL IN DENOMINATORS FOR FIELD UPDATES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   ii = nxPML_2-1
   DO i = 1,Imax-1
      if (i <= nxPML_1-1) then
         den_hx(i) = 1.0/(kappah_x_PML_1(i)*dx)
      elseif (i >= Imax+1-nxPML_2) then
         den_hx(i) = 1.0/(kappah_x_PML_2(ii)*dx)
         ii = ii-1
      else
         den_hx(i) = 1.0/dx
      endif
   ENDDO
   jj = nyPML_2-1
   DO j = 1,Jmax-1
      if (j <= nyPML_1-1) then
         den_hy(j) = 1.0/(kappah_y_PML_1(j)*dy)
      elseif (j >= Jmax+1-nyPML_2) then
         den_hy(j) = 1.0/(kappah_y_PML_2(jj)*dy)
         jj = jj-1
      else
         den_hy(j) = 1.0/dy
      endif
   ENDDO
   kk = nzPML_2-1
   DO k = 2,Kmax-1
      if (k <= nzPML_1) then
         den_hz(k) = 1.0/(kappah_z_PML_1(k-1)*dz)
      elseif (k >= Kmax+1-nzPML_2) then
         den_hz(k) = 1.0/(kappah_z_PML_2(kk)*dz)
         kk = kk - 1
      else
         den_hz(k) = 1.0/dz
      endif
   ENDDO
   ii = nxPML_2
   DO i = 1,Imax-1
      if (i <= nxPML_1) then
         den_ex(i) = 1.0/(kappae_x_PML_1(i)*dx)
      elseif (i >= Imax+1-nxPML_2) then
         den_ex(i) = 1.0/(kappae_x_PML_2(ii)*dx)
         ii = ii-1
      else
         den_ex(i) = 1.0/dx
      endif
   ENDDO
   jj = nyPML_2
   DO j = 1,Jmax-1
      if (j <= nyPML_1) then
         den_ey(j) = 1.0/(kappae_y_PML_1(j)*dy)
      elseif (j >= Jmax+1-nyPML_2) then
         den_ey(j) = 1.0/(kappae_y_PML_2(jj)*dy)
         jj = jj-1
      else
         den_ey(j) = 1.0/dy
      endif
   ENDDO
   kk = nzPML_2
   DO k = 1,Kmax-1
      if (k <= nzPML_1) then
         den_ez(k) = 1.0/(kappae_z_PML_1(k)*dz)
      elseif (k >= Kmax-nzPML_2) then
         den_ez(k) = 1.0/(kappae_z_PML_2(kk)*dz)
         kk = kk - 1
      else
         den_ez(k) = 1.0/dz
      endif
   ENDDO


!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  BEGIN TIME STEP
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  write(*,*)"begin time-stepping"
  DO n = 1,nmax
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Hx
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 2,Kmax-1
      DO i = 1,Imax-1
	   DO j = 1,Jmax-1
	      Hx(i,j,k) = DA * Hx(i,j,k) + DB *       &
			( (Ez(i,j,k) - Ez(i,j+1,k))*den_hy(j)  +    &
			  (Ey(i,j,k) - Ey(i,j,k-1))*den_hz(k) )
	   ENDDO
	ENDDO 
      DO i = 1,Imax-1
!.....................................................................
!  PML for bottom Hx, j-direction
!.....................................................................
         DO j = 1,nyPML_1-1
  	      psi_Hxy_1(i,j,k) = bh_y_1(j)*psi_Hxy_1(i,j,k)                 &
	 			+ ch_y_1(j) *(Ez(i,j,k) - Ez(i,j+1,k))/dy
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Hx, j-direction
!.....................................................................
         jj = nyPML_2-1
         DO j = Jmax+1-nyPML_2,Jmax-1
  	      psi_Hxy_2(i,jj,k) = bh_y_2(jj)*psi_Hxy_2(i,jj,k)       &
	 			+ ch_y_2(jj) *(Ez(i,j,k) -       &
                                Ez(i,(j+1),k))/dy
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_2(i,jj,k)
            jj = jj-1
         ENDDO
	ENDDO
   ENDDO
   DO i = 1,Imax-1
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Hx, k-direction
!.....................................................................
         DO k = 2,nzPML_1
  	      psi_Hxz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hxz_1(i,j,k-1)            &
	 			+ ch_z_1(k-1) *(Ey(i,j,k) - Ey(i,j,k-1))/dz
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_1(i,j,k-1)
         ENDDO
!.....................................................................
!  PML for top Hx, k-direction
!.....................................................................
         kk = nzPML_2-1
         DO k = Kmax+1-nzPML_2,Kmax-1
  	      psi_Hxz_2(i,j,kk) = bh_z_2(kk)*psi_Hxz_2(i,j,kk)             &
	 			+ ch_z_2(kk) *(Ey(i,j,k) -       &
                                Ey(i,j,k-1))/dz
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_2(i,j,kk)
           kk = kk-1
         ENDDO
	ENDDO
   ENDDO

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Hy
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 2,Kmax-1
      DO i = 1,Imax-1
	   DO j = 1,Jmax-1
            Hy(i,j,k) = DA * Hy(i,j,k) + DB *             &
			( (Ez(i+1,j,k) - Ez(i,j,k))*den_hx(i) +      &
			  (Ex(i,j,k-1) - Ex(i,j,k))*den_hz(k) )
	   ENDDO 
      ENDDO
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Hy, i-direction