      program CIP2Dnew
      implicit none
*******************************************************************
******二次元CIPのプログラム，これは自由空間中を波が伝搬していく．**
*****************************TE波**********************************
      
      integer NX,NY
      integer NT,I,J,K
      parameter (NX=160,NY=160,NT=700)

      integer xxx,BOUND(NX,NY)

      real*8 dx,dy,dt,pai,c0,eps1,eps2
      real*8 theta,eps0,mu0,z1,z2,R0,R
      
      parameter (eps1=1.0d0,eps2=36.0d0)
      parameter (dx=0.5d-3, dy=dx, dt=0.1d0*dx/3.0d8)
      parameter (R0=NX/8.0d0)
      
****** 移流方程式の変数（p→中身が正，m→中身が負）

****** 例えば，fp=Ez+Hy*Z
****** sが語尾についているものはStep1が終了した時点のモノ
****** nを現在の時刻として，未来の時刻をn+1としたら，
****** fn→fs→fn+1となる．
      real*8 FP(NX,NY),GPX(NX,NY),GPY(NX,NY)
      real*8 FM(NX,NY),GMX(NX,NY),GMY(NX,NY)
      real*8 FPS(NX,NY),GPXS(NX,NY),GPYS(NX,NY)
      real*8 FMS(NX,NY),GMXS(NX,NY),GMYS(NX,NY)

****** 電界，磁界，光速，インピーダンス
      real*8 EZ(NX,NY),EZDX(NX,NY),EZDY(NX,NY)
      real*8 HX(NX,NY),HXDX(NX,NY),HXDY(NX,NY)
      real*8 HY(NX,NY),HYDX(NX,NY),HYDY(NX,NY),dypp
      real*8 C(NX,NY),Z(NX,NY),dxpp
****** integer cd,year,month,day,hour,min,sec,ssec
****** character fnum*4, cyear*4, cmonth*2, cday*2
****** character chour*2, cmin*2, csec*2, date*14
      
****** π，μ，ε，インピーダンス
      pai=3.14159265358979324d0
      mu0=4.0d-7*pai
      eps0=8.85418782d-12
      c0=1.0d0/DSQRT(eps0*mu0)
      z1=DSQRT(mu0/eps0/eps1)
      z2=DSQRT(mu0/eps0/eps2)
      
      
**********結果格納用ディレクトリ生成*********************************
******         call getdat(year, month, day)
******         call gettim(hour, min, sec, ssec)
******         write(cyear,'(I4.4)') year
******         write(cmonth,'(I2.2)') month
******         write(cday,'(I2.2)') day
******         write(chour,'(I2.2)') hour
******         write(cmin,'(I2.2)') min
******         write(csec,'(I2.2)') sec
******         write(fnum,'(I4.4)') NT
******         date = cyear // cmonth // cday // chour // cmin // csec
        
******         mkdir = SYSTEM ('mkdir '//date//'_NT_'//fnum)
******         cd = chdir(date//'_NT_'//fnum)
*********************************************************************
      
****** 境界生成
      open(10,file='bound.dat')
       DO 7020 I=1,NX
          DO 7030 J=1,NY
**             IF ((I.GE.NX/2).AND.(J.EQ.NY/2)) THEN
**                 BOUND(I,J)=1
**             ELSE IF((J.GE.NY/2).AND.(I.EQ.NX/2)) THEN
**                 BOUND(I,J)=1
**             IF(I.EQ.NX/2) then
**              BOUND(I,J)=1
**             ELSE
**                 BOUND(I,J)=0
**             ENDIF

**     半分が誘電体
      IF(I.EQ.NX/2)then
            BOUND(I,J)=1
      else
            BOUND(I,J)=0
      end if


**           円柱の境界生成
**            R=INT(DSQRT(DFLOAT((I-NX/5*3)**2+(J-NY/2)**2)))
**            IF(R.EQ.R0) THEN
**                  BOUND(I,J)=1
**            ELSE
**                  BOUND(I,J)=0
**            END IF
             write (10,*) I, J, BOUND(I,J) 
7030       CONTINUE
7020     CONTINUE
      close(10)
      
**     空間の作成
      open(11,file='z.dat')        
      DO 7000 J=1,NY
           DO 7010 I=1,NX
           
**              IF(I.GT.NX/2) THEN
**               IF ((I.GT.NX/2).AND.(J.GT.NY/2)) THEN
**                C(I,J)=c0/DSQRT(eps2)
**                Z(I,J)=z2
**              ELSE IF(I.EQ.NX/2) THEN
**               else if (BOUND(I,J).EQ.1) THEN
**                C(I,J)=0.0
**                Z(I,J)=0.0
**              ELSE
**                C(I,J)=c0/DSQRT(eps1)
**                Z(I,J)=z1
**              END IF
            
**     半分から右が誘電体
      if(I.GT.NX/2) then
            C(I,J)=c0/DSQRT(eps2)
            Z(I,J)=z2
      else if(I.EQ.NX/2) then
            C(I,J)=((c0/DSQRT(eps1))+(c0/DSQRT(eps1)))/2.0d0
            Z(I,J)=(z1+z2)/2.0d0
      else
            C(I,J)=c0/DSQRT(eps1)
            Z(I,J)=z1
      end if
      
**            円柱の中身の作成
**            R=INT(DSQRT(DFLOAT((I-NX/5*3)**2+(J-NY/2)**2)))
**            IF(R.LT.R0) THEN
**             C(I,J)=c0/DSQRT(eps2)
**             Z(I,J)=z2
**            ELSE IF(BOUND(I,J).EQ.1) THEN
**             C(I,J)=((c0/DSQRT(eps1))+(c0/DSQRT(eps1)))/2.0d0
**             Z(I,J)=(z1+z2)/2.0d0
**             C(I,J)=0.0d0
**             Z(I,J)=0.0d0
**            ELSE
**             C(I,J)=c0/DSQRT(eps1)
**             Z(I,J)=z1
**            END IF
             write (11,*) I, J, Z(I,J) 
7010       CONTINUE
7000    CONTINUE

      close(11)

**     初期化
      do 1030 I=1,NX
            do 1040 J=1,NY
              EZ(I,J)=0.0d0
              EZDX(I,J)=0.0d0
              EZDY(I,J)=0.0d0
              
              HX(I,J)=0.0d0
              HXDX(I,J)=0.0d0
              HXDY(I,J)=0.0d0
              
              HY(I,J)=0.0d0
              HYDX(I,J)=0.0d0
              HYDY(I,J)=0.0d0
1040        continue
1030  continue

**     初期電界
      do 1060 I=1,NX
           do 1070 J=1,NY
              theta=DFLOAT(I-NX/10*3)
              EZ(I,J)=EXP(-0.01d0*theta**2)
              EZDX(I,J)=-2.0d0*0.01d0*theta*EXP(-0.01d0*theta**2)
              EZDY(I,J)=0.0d0

            if (BOUND(I,J).EQ.1) then
              HY(I,J)=0.0d0
              HYDX(I,J)=0.0d0
              HYDY(I,J)=0.0d0
            else
              HY(I,J)=-EZ(I,J)/Z(I,J)
              HYDX(I,J)=2.0d0*0.01d0*theta*EXP(-0.01d0*theta**2)/Z(I,J)
              HYDY(I,J)=0.0d0
            end if
            
1070      CONTINUE
1060    CONTINUE

      
      OPEN(8050,FILE='init_ez.dat')
      do 7080 J=1,NX
           do 7090 I=1,NY
            WRITE(8050,*) I, J, EZ(I,J)
7090       CONTINUE
7080	CONTINUE
      CLOSE(8050)
      
**     ここからタイムループ
      do 1095 K=1,NT
      
      write(*,*) 'Time Step K=',K
      
**     FとGの計算
      do 1100 I=1,NX
            do 1110 J=1,NY
*************************************************************
****** 和，差の移流方程式とその微分値
      if ((BOUND(I,J).EQ.1).and.((I.NE.NX).or.(I.NE.1))) then
       FP(I,J)=EZ(I,J)+HY(I,J)*Z(I-1,J)
       FM(I,J)=EZ(I,J)-HY(I,J)*Z(I+1,J)
       GPX(I,J)=EZDX(I,J)+HYDX(I,J)*Z(I-1,J)
       GMX(I,J)=EZDX(I,J)-HYDX(I,J)*Z(I+1,J)
       GPY(I,J)=EZDY(I,J)+HYDY(I,J)*Z(I-1,J)
       GMY(I,J)=EZDY(I,J)-HYDY(I,J)*Z(I+1,J)
      else
       FP(I,J)=EZ(I,J)+HY(I,J)*Z(I,J)
       FM(I,J)=EZ(I,J)-HY(I,J)*Z(I,J)
       GPX(I,J)=EZDX(I,J)+HYDX(I,J)*Z(I,J)
       GMX(I,J)=EZDX(I,J)-HYDX(I,J)*Z(I,J)
       GPY(I,J)=EZDY(I,J)+HYDY(I,J)*Z(I,J)
       GMY(I,J)=EZDY(I,J)-HYDY(I,J)*Z(I,J)
      end if
      
*************************************************************
1110        continue
1100  continue

**     Step1の計算～その1

**      OPEN(8110,FILE='xxx1.dat')
      do 1200 I=1,NX
            do 1250 J=1,NY      
            dxpp=0    
            

****** 左へ移流する処理，領域の外側から２セルは通常の処理とする．
****** 配列の問題回避のため
      if((I.LE.2).OR.(I.GE.NX-2))then
******      通常の処理
            call STEP1(I,J,NX,NY,-1.0,C,FP,GPX,GPY,FPS,GPXS,GPYS,dx,dt)
            
      else
******      境界の直前 I-1
            if ((BOUND(I-1,J).EQ.1).and.(BOUND(I-2,J).NE.1)) then
                  dxpp=C(I,J)*dt
                  xxx=1
            call STEP1(I,J,NX,NY,-1.0,C,FP,GPX,GPY,FPS,GPXS,GPYS,
     1       dxpp,dt)
            
******      境界上
            else if(BOUND(I,J).EQ.1) then
            
            
            xxx=2
******      後で処理

******      境界の直後 I+1
            else if((BOUND(I+1,J).EQ.1).and.(BOUND(I+2,J).NE.1)) then
                  dxpp=C(I,J)*dt
                  xxx=3
            call STEP1(I,J,NX,NY,-1.0,C,FP,GPX,GPY,FPS,GPXS,GPYS,
     1       dxpp,dt)
                  
******      境界の後 I+2
            else if((BOUND(I+2,J).EQ.1).and.((BOUND(I+3,J).NE.1)
     1       .and.(BOUND(I+1,J).NE.1))) then
                  dxpp=dx-C(I,J)*dt
                  xxx=4
            call STEP1(I,J,NX,NY,-1.0,C,FP,GPX,GPY,FPS,GPXS,GPYS,
     1       dxpp,dt)
     
******      通常の処理
            else
            xxx=5
            call STEP1(I,J,NX,NY,-1.0,C,FP,GPX,GPY,FPS,GPXS,GPYS,
     1       dx,dt)
     
            end if
      
      end if
******      write (8110,*) I, xxx
1250  continue
1200  continue
******      close(8110)


******      Step1の計算～その2


******      OPEN(8120,FILE='xxx2.dat')
      do 1300 I=1,NX
            do 1350 J=1,NY
            
******      右へ移流する処理
      if((I.LE.2).OR.(I.GE.NX-2))then
      
******      通常の処理
            call STEP1(I,J,NX,NY,1.0,C,FM,GMX,GMY,FMS,GMXS,GMYS,dx,dt)
      else
      
******      境界の直前 I+1
            if ((BOUND(I+1,J).EQ.1).and.(BOUND(I+2,J).NE.1)) then
                  dxpp=C(I,J)*dt
                  xxx=1
            call STEP1(I,J,NX,NY,1.0,C,FM,GMX,GMY,FMS,GMXS,GMYS,dxpp,dt)
            
******      境界上
            else if(BOUND(I,J).EQ.1) then

******      後で処理
            xxx=2
            
            
******      境界の直後 I-1
            else if((BOUND(I-1,J).EQ.1).and.(BOUND(I-2,J).NE.1)) then
                  dxpp=C(I,J)*dt
                  xxx=3
                  
            call STEP1(I,J,NX,NY,1.0,C,FM,GMX,GMY,FMS,GMXS,GMYS,dxpp,dt)

******      境界の後 I-2
            else if((BOUND(I-2,J).EQ.1).and.((BOUND(I-3,J).NE.1)
     1       .and.(BOUND(I+1,J).NE.1))) then
                  dxpp=dx-C(I,J)*dt
                  xxx=4
                  
            call STEP1(I,J,NX,NY,1.0,C,FM,GMX,GMY,FMS,GMXS,GMYS,dxpp,dt)
            
******      通常の処理
            else
            call STEP1(I,J,NX,NY,1.0,C,FM,GMX,GMY,FMS,GMXS,GMYS,dx,dt)
            xxx=5
            end if
******      write(8120,*) I, xxx         
      end if
            
1350        continue
1300  continue
******      close(8120)

****** 境界上の処理（境界の前後が境界でないときのみ実行）
      DO 7180 I=2,NX-1
         DO 7190 J=1,NY
            if((BOUND(I,J).EQ.1).and.((BOUND(I-1,J).NE.1)
     1       .and.(BOUND(I+1,J).NE.1))) then   
                  call STEP_BOUND_X(I,J,NX,NY,Z,FP,GPX,FM,GMX,
     1                  FPS,GPXS,FMS,GMXS,EZ,HY,EZDX,HYDX)
            end if
7190  continue
7180  continue   


******      中間値の計算　ここが中間地点
      do 1400 I=1,NX
            do 1450 J=1,NY
            if (BOUND(I,J).NE.1) then
     
              EZ(I,J)=(FPS(I,J)+FMS(I,J))/2.0d0
              HY(I,J)=(FPS(I,J)-FMS(I,J))/Z(I,J)/2.0d0
              EZDX(I,J)=(GPXS(I,J)+GMXS(I,J))/2.0d0
              HYDX(I,J)=(GPXS(I,J)-GMXS(I,J))/Z(I,J)/2.0d0
              EZDY(I,J)=(GPYS(I,J)+GMYS(I,J))/2.0d0              
              HYDY(I,J)=(GPYS(I,J)-GMYS(I,J))/Z(I,J)/2.0d0
            end if

1450        continue
1400  continue
      do 1500 I=1,NX
            do 1550 J=1,NY
*******************************************************************+
****** 和，差の移流方程式とその微分値
            if ((BOUND(I,J).EQ.1).and.((J.NE.NY).or.(J.NE.1)))then
             FP(I,J)=EZ(I,J)+HX(I,J)*Z(I,J+1)
             FM(I,J)=EZ(I,J)-HX(I,J)*Z(I,J-1)
             GPXS(I,J)=EZDX(I,J)+HXDX(I,J)*Z(I,J+1)
             GMXS(I,J)=EZDX(I,J)-HXDX(I,J)*Z(I,J-1)
             GPYS(I,J)=EZDY(I,J)+HXDY(I,J)*Z(I,J+1)
             GMYS(I,J)=EZDY(I,J)-HXDY(I,J)*Z(I,J-1)
             
            else
             FP(I,J)=EZ(I,J)+HX(I,J)*Z(I,J)
             FM(I,J)=EZ(I,J)-HX(I,J)*Z(I,J)
             GPXS(I,J)=EZDX(I,J)+HXDX(I,J)*Z(I,J)
             GMXS(I,J)=EZDX(I,J)-HXDX(I,J)*Z(I,J)
             GPYS(I,J)=EZDY(I,J)+HXDY(I,J)*Z(I,J)
             GMYS(I,J)=EZDY(I,J)-HXDY(I,J)*Z(I,J)
            end if
*******************************************************************
1550        continue
1500  continue
      
****** Step2の計算～その1
      do 1301 I=1,NX
            do 1351 J=1,NY
            
****** 奥へ移流する処理
      if((J.LE.2).OR.(J.GE.NX-2))then
****** 通常の処理
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,dy,dt)
            
      else
            if ((J.EQ.1).OR.(J.EQ.NY)) then
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,
     1       dy,dt)

******      境界の直前 J+1
            else if ((BOUND(I,J+1).EQ.1).and.(BOUND(I,J+2).NE.1)) then
                  dypp=C(I,J)*dt
                  xxx=1
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,
     1       dypp,dt)
            
******      境界上・何もしない
            else if(BOUND(I,J).EQ.1) then

******      境界の直後 J-1
            else if((BOUND(I,J-1).EQ.1).and.(BOUND(I,J-2).NE.1)) then
                  dypp=C(I,J)*dt
                  xxx=3
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,
     1       dypp,dt)

******      境界の後 J-2
            else if((BOUND(I,J-2).EQ.1).and.((BOUND(I,J-3).NE.1)
     1       .and.(BOUND(I,J+1).NE.1))) then
                  dypp=dx-C(I,J)*dt
                  xxx=4
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,
     1       dypp,dt)

******      通常の処理
            else
            call STEP2(I,J,NX,NY,1.0,C,FP,GPXS,GPYS,FPS,GPX,GPY,
     1       dy,dt)
            xxx=5
            end if      
      end if
1351        continue
1301  continue


       
****** Step2の計算～その2
      do 1201 I=1,NX
            do 1251 J=1,NY      
            dypp=0.d0
                  
****** 手前へ移流する処理
      if((J.LE.2).OR.(J.GE.NX-2))then
****** 通常の処理
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,dy,dt)
      else
            if ((J.EQ.1).OR.(J.EQ.NY)) then
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,
     1       dy,dt)
******      境界の直前 J-1
            else if ((BOUND(I,J-1).EQ.1).and.(BOUND(I,J-2).NE.1)) then
                  dypp=C(I,J)*dt
                  xxx=1
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,
     1       dypp,dt)

            xxx=2

******      境界上・何もしない（境界の前後も境界）
            else if(BOUND(I,J).EQ.1) then


******      境界の直後 J+1
            else if((BOUND(I,J+1).EQ.1).and.(BOUND(I,J+2).NE.1)) then
                  dypp=C(I,J)*dt
                  xxx=3
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,
     1       dypp,dt)
                  
******      境界の後 J+2
            else if((BOUND(I,J+2).EQ.1).and.((BOUND(I,J+3).NE.1)
     1       .and.(BOUND(I,J+1).NE.1))) then
                  dypp=dx-C(I,J)*dt
                  xxx=4
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,
     1       dypp,dt)
******      通常の処理
            else
            xxx=5
            call STEP2(I,J,NX,NY,-1.0,C,FM,GMXS,GMYS,FMS,GMX,GMY,
     1       dy,dt)
            end if
      end if
1251  continue
1201  continue

****** 境界上の処理（境界の前後が境界でないときのみ実行）
      DO 7181 I=1,NX
         DO 7191 J=2,NY-1
            if ((BOUND(I,J).EQ.1).and.((BOUND(I,J-1).NE.1)
     1       .and.(BOUND(I,J+1).NE.1))) then   
                  call STEP_BOUND_Y(I,J,NX,NY,Z,FP,GPYS,FM,GMYS,
     1                  FPS,GPY,FMS,GMY,EZ,HX,EZDY,HXDY)
            end if
            
7191  continue
7181  continue   


      
****** 電磁界の計算
      do 7200 I=1,NX
            do 7210 J=1,NY
            
            if (BOUND(I,J).NE.1) then
              EZ(I,J)=(FPS(I,J)+FMS(I,J))/2.0d0
              HX(I,J)=(FPS(I,J)-FMS(I,J))/Z(I,J)/2.0d0
              EZDY(I,J)=(GPY(I,J)+GMY(I,J))/2.0d0
              HXDY(I,J)=(GPY(I,J)-GMY(I,J))/Z(I,J)/2.0d0
              EZDX(I,J)=(GPX(I,J)+GMX(I,J))/2.0d0
              HXDX(I,J)=(GPX(I,J)-GMX(I,J))/Z(I,J)/2.0d0
            end if
7210        continue 
7200  continue

****** タイムループの終端
1095  continue
      
      WRITE(*,*) 'Saving Data.'
	OPEN(8145,FILE='fin_ez.dat')

	DO 7340 I=1,NX
	   DO 7350 J=1,NY
	      WRITE(8145,77) I,J,EZ(I,J)
7350	   CONTINUE
7340	CONTINUE

	CLOSE(8145)
77    format(i7,1x,i7,1x,f12.7)
      stop
      end
      
*************************************************************************
****** 「CIP法」P.72
****** FS,GXS,GYSは戻り値

        SUBROUTINE STEP1(I0,J0,NX,NY,LAMBDA,C,F,GX,GY,FS,GXS,GYS,DX,DT)
        INTEGER I0,J0,ISG,IUP,NX,NY
        real LAMBDA
        real*8 XX,FDIF,XAM1,XBM1,DT,DX
        real*8 F(NX,NY),GX(NX,NY),GY(NX,NY),C(NX,NY),YU(NX,NY)
        real*8 FS(NX,NY),GXS(NX,NY),GYS(NX,NY)
      
	IF (LAMBDA.GT.0) THEN
	   ISG=-1
	ELSE
	   ISG=1
	ENDIF
        IUP=I0+ISG
        

        YU(I0,J0)=LAMBDA*C(I0,J0)
      
        XX=-1.0d0*YU(I0,J0)*DT
        
        IF (IUP.LE.0) THEN
           IUP=1
        ELSE IF (IUP.GT.NX) THEN
           IUP=NX
	ENDIF

        FDIF=(F(IUP,J0)-F(I0,J0))/DX*DFLOAT(ISG)        
        XAM1=(GX(I0,J0)+GX(IUP,J0)-2.0d0*FDIF)/(DX*DX)
        XBM1=(3.0d0*FDIF-2.0d0*GX(I0,J0)-GX(IUP,J0))/DX*DFLOAT(ISG)
        FS(I0,J0)=((XAM1*XX+XBM1)*XX+GX(I0,J0))*XX+F(I0,J0)
        GXS(I0,J0)=(3.0d0*XAM1*XX+2.0d0*XBM1)*XX+GX(I0,J0)
        GYS(I0,J0)=GY(I0,J0)
     1                   -XX/DX*DFLOAT(ISG)*(GY(I0,J0)-GY(IUP,J0))
	RETURN
	END

***********************************************************************
      SUBROUTINE STEP2(I0,J0,NX,NY,LAMBDA,C,FS,GXS,GYS,FN,GXN,GYN,DY,DT)
      INTEGER I0,J0,JSG,JUP,NX,NY
      real LAMBDA
      real*8 YY,FDIF,XAM1,XBM1,DT,DY
	real*8 FS(NX,NY),GYS(NX,NY),C(NX,NY),YU(NX,NY)
	real*8 FN(NX,NY),GYN(NX,NY),GXN(NX,NY),GXS(NX,NY)

	IF (LAMBDA.GT.0) THEN
	   JSG=-1
	ELSE
	   JSG=1
	ENDIF

        JUP=J0+JSG

        YU(I0,J0)=LAMBDA*C(I0,J0)

        YY=-1.0d0*YU(I0,J0)*DT

        IF (JUP.LE.0) THEN
           JUP=1
        ELSE IF (JUP.GT.NY) THEN
           JUP=NY
        END IF

        FDIF=(FS(I0,JUP)-FS(I0,J0))/DY*DFLOAT(JSG)
        XAM1=(GYS(I0,J0)+GYS(I0,JUP)-2.0d0*FDIF)/(DY*DY)
        XBM1=(3.0d0*FDIF-2.0d0*GYS(I0,J0)-GYS(I0,JUP))/DY*DFLOAT(JSG)
        FN(I0,J0)=((XAM1*YY+XBM1)*YY+GYS(I0,J0))*YY+FS(I0,J0)
        GYN(I0,J0)=(3.0d0*XAM1*YY+2.0d0*XBM1)*YY+GYS(I0,J0)
        GXN(I0,J0)=GXS(I0,J0)-YY/DY*DFLOAT(JSG)*(GXS(I0,J0)-GXS(I0,JUP))
        
	RETURN
	END
***********************************************************
      SUBROUTINE STEP_BOUND_X(I0,J0,NX,NY,Z,FP,GP,FM,GM,
     1                         FPN,GPN,FMN,GMN,EZ,HY,EZDX,HYDX)    
      integer I0,J0,NX,NY
      real*8 Z(NX,NY)
      real*8 FP(NX,NY),GP(NX,NY),FM(NX,NY),GM(NX,NY)
      real*8 FPN(NX,NY),GPN(NX,NY),FMN(NX,NY),GMN(NX,NY)
      real*8 EZ(NX,NY),HY(NX,NY),EZDX(NX,NY),HYDX(NX,NY)
       
******電磁界の計算
      EZ(I0,J0)=Z(I0-1,J0)*FP(I0+1,J0)+Z(I0+1,J0)*FM(I0-1,J0)
      EZ(I0,J0)=EZ(I0,J0)/(Z(I0-1,J0)+Z(I0+1,J0))
      
      HY(I0,J0)=(FP(I0+1,J0)-FM(I0-1,J0))/(Z(I0-1,J0)+Z(I0+1,J0))
      
      EZDX(I0,J0)=Z(I0-1,J0)*GP(I0+1,J0)+Z(I0+1,J0)*GM(I0-1,J0)
      EZDX(I0,J0)=EZDX(I0,J0)/(Z(I0-1,J0)+Z(I0+1,J0))
      
      HYDX(I0,J0)=(GP(I0+1,J0)-GM(I0-1,J0))/(Z(I0-1,J0)+Z(I0+1,J0))
      
******fnとｇnの計算（現在の電界と磁界から求める）
      FPN(I0,J0)=EZ(I0,J0)+HY(I0,J0)*Z(I0-1,J0)
      GPN(I0,J0)=EZDX(I0,J0)+HYDX(I0,J0)*Z(I0-1,J0)
      FMN(I0,J0)=EZ(I0,J0)-HY(I0,J0)*Z(I0+1,J0)
      GMN(I0,J0)=EZDX(I0,J0)-HYDX(I0,J0)*Z(I0+1,J0)

      RETURN
      END
***************************************************************************
      SUBROUTINE STEP_BOUND_Y(I0,J0,NX,NY,Z,FP,GP,FM,GM,
     1                         FPN,GPN,FMN,GMN,EZ,HY,EZDX,HYDX)

      integer I0,J0,NX,NY
      real*8 Z(NX,NY)
      real*8 FP(NX,NY),GP(NX,NY),FM(NX,NY),GM(NX,NY)
      real*8 FPN(NX,NY),GPN(NX,NY),FMN(NX,NY),GMN(NX,NY)
      real*8 EZ(NX,NY),HY(NX,NY),EZDX(NX,NY),HYDX(NX,NY)

****** 電磁界の計算
      EZ(I0,J0)=Z(I0,J0+1)*FP(I0,J0-1)+Z(I0,J0-1)*FM(I0,J0+1)
      EZ(I0,J0)=EZ(I0,J0)/(Z(I0,J0-1)+Z(I0,J0+1))
      
      HY(I0,J0)=(FP(I0,J0-1)-FM(I0,J0+1))/(Z(I0,J0-1)+Z(I0,J0+1))
      
      EZDX(I0,J0)=Z(I0,J0+1)*GP(I0,J0-1)+Z(I0,J0-1)*GM(I0,J0+1)
      EZDX(I0,J0)=EZDX(I0,J0)/(Z(I0,J0-1)+Z(I0,J0+1))
      
      HYDX(I0,J0)=(GP(I0,J0-1)-GM(I0,J0+1))/(Z(I0,J0-1)+Z(I0,J0+1)) 
      
****** fnとｇnの計算（現在の電界と磁界から求める）
      FPN(I0,J0)=EZ(I0,J0)+HY(I0,J0)*Z(I0,J0+1)
      GPN(I0,J0)=EZDX(I0,J0)+HYDX(I0,J0)*Z(I0,J0+1)
      FMN(I0,J0)=EZ(I0,J0)-HY(I0,J0)*Z(I0,J0-1)
      GMN(I0,J0)=EZDX(I0,J0)-HYDX(I0,J0)*Z(I0,J0-1)

      RETURN
      END
