PROGRAM TROPICS2_OCEAN
C
	PARAMETER(NX=43,NY=43,NY1=NY+1,NT=4000)
	REAL U1(NX,NY),U2(NX,NY),U3(NX,NY),V1(NX,NY1),V2(NX,NY1)
      REAL V3(NX,NY1),CLRAD1(NX,NY),CLRAD2(NX,NY),CLRAD3(NX,NY)
	REAL S1(NX,NY),S2(NX,NY),S3(NX,NY),SB1(NX,NY),SB2(NX,NY)
	REAL S01(NX,NY),S02(NX,NY),S03(NX,NY),SOLAR(NX,NY),HMIX
	REAL SB3(NX,NY),S0(NX,NY),S0I(NX,NY),W(NX,NY),X(NX),Y(NY),TI(NT)
	REAL SM1(NX,NY),SM2(NX,NY),SM3(NX,NY)
      REAL M1(NX,NY),M2(NX,NY),M3(NX,NY), DUM(NX)
	REAL WHOV(NX,NT),UHOV(NX,NT),VHOV(NX,NT)
	REAL WHOVA(NX,NT), UHOVA(NX,NT),VHOVA(NX,NT)
      REAL MEQ,H,HM,N2,LX,LY,LD(NX),YM(NX),MF
	REAL MMOVIE(NX,NY,NT),UMOVIE(NX,NY,NT)
	REAL VMOVIE(NX,NY,NT),THMOVIE(NX,NY,NT)
	INTEGER*4 IRAN
	CHARACTER*60 ANAME
	CHARACTER*10 KNAME
	CHARACTER*1 READOUTPUT,EPINT,IOCEAN
C
C   Read in parameters from text file
C
	OPEN(UNIT=9,FILE='input.txt', STATUS='OLD')
	READ(9,5)
	READ(9,5)
	READ(9,*)CD
	READ(9,*)CK
	READ(9,*)UMEAN
	READ(9,*)USTAR
	READ(9,*)H
	READ(9,*)N2
	READ(9,*)EP0
	READ(9,*)RAD0
	READ(9,*)RFEED
	READ(9,*)DTHE
	READ(9,*)GAMMA
	READ(9,*)RADMOIST
	READ(9,*)EPINT
	READ(9,5)
	READ(9,*)TAU
	READ(9,*)TAUCL
	READ(9,*)TAUM
	READ(9,5)
	READ(9,*)LX
	READ(9,*)LY
	READ(9,*)SSTMAX
	READ(9,*)SSTMIN
	READ(9,*)DTEMP
	READ(9,*)YM0
	READ(9,*)YM1
	READ(9,*)SSTA
	READ(9,*)SSTY
	READ(9,*)SSTDY
	READ(9,5)
	READ(9,*)IOCEAN
	READ(9,*)HMIX
	READ(9,*)TAUOCEAN
	READ(9,5)
	READ(9,*)DT
	READ(9,*)ENDTIME
	READ(9,*)PRINTTIME
	READ(9,*)TMOVIE
	READ(9,*)PRINTMOVIE
	READ(9,*)READOUTPUT
    5 FORMAT(2X)
C
C  Define some thermodynamics and other constants
C
	CP=1005.
	ALV=2.5E6
	A=6.38E6
	PI=2.*ASIN(1.0)
	IRAN=9973
C
C  Normalize some parameters and define others
C
	CD=CD*1.0E-3
	CK=CK*1.0E-3
	H=H*1000.
	AIH=1./H
	GFAC=GAMMA*AIH
	CDH=CD*AIH
	CKH=CK*AIH
	LX=2.*PI*A*LX/360.0
	LY=2.*PI*A*LY/360.0
	YM0=YM0*2.*PI*A/360.0
	YM1=YM1*2.*PI*A/360.0
	SSTY=SSTY*2.*PI*A/360.0
	SSTDY=SSTDY*2.*PI*A/360.0
	BETA=2.*PI/(12.*3600.*A)
	CHID=CP*DTEMP*LOG(350./(350.-DTHE))
	RAD0=RAD0*CP*DTEMP*2./(270.*24.*3600.)
c	HM=2.*CHID*(EP0+GAMMA)/(N2*EP0)
	TAU=TAU*3600.
	TAUI=1./MAX(TAU,1.0)
	TAUM=TAUM*24.*3600.
	TAUMI=1./TAUM
	TAUCL=TAUCL*3600.
	TAUCLI=1./TAUCL
	TAUOCEAN=TAUOCEAN*24.*3600.
	TAUOCEANI=1./TAUOCEAN
	US2=USTAR**2
	ENDTIME=ENDTIME*24.*3600.
	PRINTTIME=PRINTTIME*3600.
	PRINTMOVIE=PRINTMOVIE*24.*3600.
	TMOVIE=TMOVIE*24.*3600.
	NY2=NY/2	
	NX2=NX/2
	FNYI=1./(2.*FLOAT(NY2-1))
	DX=LX/FLOAT(NX-1)
	DY=2.*LY/FLOAT(NY-1)
	DXI=1./DX
	DYI=1./DY
C
	GDGM=3.0
	SFAC=0.001/(4.167*HMIX)
	CLFAC=GDGM*1.E6*5.E4*100./(10.*287.*1.E3*HMIX*4167.0*300.0)
C
C  This is the numerical damping used in the Asselin Filter of the leap frog time scheme
C
	DAMP=0.1
C
C  Create arrays of x and y variables
C
	DO I=1,NX
	 X(I)=DX*FLOAT(I-1)
c	 LD(I)=1000.*(1000.+800.*COS(2.*PI*X(I)/LX))
c	 YM(I)=1000.*(200.-200.*COS(2.*PI*X(I)/LX))
	 LD(I)=YM1
	 YM(I)=YM0
	END DO
	DO J=1,NY
	 Y(J)=DY*FLOAT(J-1-NY2)
	END DO
C
C	Initial conditions
C
	DO I=1,NX
	 DO K=1,NT
	  UHOV(I,K)=0.0
	  VHOV(I,K)=0.0
	  WHOV(I,K)=0.0
	  UHOVA(I,K)=0.0
	  VHOVA(I,K)=0.0
	  WHOVA(I,K)=0.0
	 END DO
	 DO J=1,NY
C
C  Here we prescribe the mean sea surface temperature distribution and add anomalies to it
C
	  SST=SSTMIN+(SSTMAX-SSTMIN)*EXP(-((Y(J)-YM(I))/LD(I))**2)
	  SST=SST+SSTA*EXP(-((Y(J)-SSTY)/SSTDY)**2)*COS(2.*PI*X(I)/LX)
C
	  SSTABS=SST+273.15
	  ES=6.112*EXP(17.67*SST/(243.5+SST))
	  QS=0.622*ES/(1013.-ES)
	  S0I(I,J)=DTEMP*(CP*LOG(SSTABS)+ALV*QS/SSTABS)	  
	  S1(I,J)=S0I(I,J)-80.*DTEMP
	  S2(I,J)=S1(I,J)
	  SOLAR(I,J)=CK*SQRT(UMEAN**2+US2)*80.*DTEMP
	  SB1(I,J)=S1(I,J)
	  SB2(I,J)=SB1(I,J)
	  SM1(I,J)=S1(I,J)-CHID
	  SM2(I,J)=SM1(I,J)
	  V1(I,J)=0.3*(RAN(IRA)-0.5)
	  V2(I,J)=V1(I,J)
	  U1(I,J)=0.5*(RAN(IRA)-0.5)
	  U2(I,J)=U1(I,J)
	  M1(I,J)=0.0
	  M2(I,J)=M1(I,J)
	  CLRAD1(I,J)=0.0
	  CLRAD2(I,J)=CLRAD1(I,J)
	  S01(I,J)=0.0
	  S02(I,J)=0.0
	 END DO
	 V1(I,NY1)=0.0
	 V2(I,NY1)=0.0
	END DO
c	U1(NX2,NY2)=-0.5
c	U2(NX2,NY2)=U1(NX2,NY2)
c	V1(NX2+2,NY2)=0.6
c	V2(NX2+2,NY2)=V1(NX2+2,NY2)
C
C  Restart from previous output, if desired
C
	IF(READOUTPUT.EQ.'y')THEN
	OPEN(UNIT=13,FILE='output/u.out',STATUS='OLD')
	OPEN(UNIT=14,FILE='output/v.out',STATUS='OLD')
	OPEN(UNIT=17,FILE='output/s.out',STATUS='OLD')
	OPEN(UNIT=18,FILE='output/sb.out',STATUS='OLD')
	OPEN(UNIT=19,FILE='output/sm.out',STATUS='OLD')
	OPEN(UNIT=20,FILE='output/m.out',STATUS='OLD')
	OPEN(UNIT=21,FILE='output/s0.out',STATUS='OLD')
C
	DO J=2,NY-1
	 READ(13,200)(U2(I,J),I=2,NX-1)
	 READ(14,200)(V2(I,J),I=2,NX-1)
	 READ(20,200)(DUM(I),I=2,NX-1)
	 DO I=2,NX-1
	  M2(I,J)=0.01*DUM(I)
	 END DO
	 READ(17,200)(DUM(I),I=2,NX-1)
	 DO I=2,NX-1
	  S2(I,J)=DTEMP*CP*LOG(DUM(I))
	 END DO
	 READ(18,200)(DUM(I),I=2,NX-1)
	 DO I=2,NX-1
	  SB2(I,J)=DTEMP*CP*LOG(DUM(I))
	 END DO
	 READ(19,200)(DUM(I),I=2,NX-1)
	 DO I=2,NX-1
	  SM2(I,J)=DTEMP*CP*LOG(DUM(I))
	 END DO
	 READ(21,200)(DUM(I),I=2,NX-1)
	 DO I=2,NX-1
	  S01(I,J)=DTEMP*CP*LOG(DUM(I))-S0I(I,J)
	  S02(I,J)=S01(I,J)
	 END DO
	 DO I=2,NX-1
	  S1(I,J)=S2(I,J)
	  SB1(I,J)=SB2(I,J)
	  U1(I,J)=U2(I,J)
	  V1(I,J)=V2(I,J)
	  SM1(I,J)=SM2(I,J)
	  M1(I,J)=M2(I,J)
	 END DO
	END DO
	CLOSE(13)
	CLOSE(14)
	CLOSE(17)
	CLOSE(18)
	CLOSE(19)
	CLOSE(20)
	CLOSE(21)
	END IF
C
C
	TIME=0.0
	TPRINT=0.0
	IPRINT=0
	TPRINTM=0.0
	IPRINTM=0
C
C	Time marching begins here
C
  100	CONTINUE
	TIME=TIME+DT
	TPRINT=TPRINT+DT
	TPRINTM=TPRINTM+DT
	IF(TIME.GE.ENDTIME)THEN
	 GOTO 300
	END IF
C
	YMM=2.*PI*A*20./360.
	DO I=1,NX
	 DO J=1,NY
C	  IF(X(I).LE.0.25*LX.OR.X(I).GE.0.75*LX)THEN
c	  S0(I,J)=S0I(I,J)-DTEMP*35.*(1.+0.4*COS(2.*PI*TIME/(40.*24.*3600.
c     1   )))*SIN(2.*PI*X(I)/LX)*EXP(-((Y(J)-YMM)/(600.*1000.))**2)
C	  ELSE
	   S0(I,J)=S0I(I,J)
C	  END IF
	 END DO
	END DO
C
C	Boundary conditions
C
	DO I=1,NX
	 V2(I,1)=0.0
	 V2(I,2)=0.0
	 V2(I,NY)=0.0
	 S2(I,1)=S2(I,2)
	 S2(I,NY)=S2(I,NY-1)
	 SM1(I,1)=SM1(I,2)
	 SM1(I,NY)=SM1(I,NY-1)
	 U2(I,1)=0.0
	 U2(I,NY)=0.0
	END DO
C
	DO J=1,NY
	 U2(1,J)=U2(NX-1,J)
	 U2(NX,J)=U2(2,J)
	 V2(1,J)=V2(NX-1,J)
	 V2(NX,J)=V2(2,J)
	 S2(1,J)=S2(NX-1,J)
	 S2(NX,J)=S2(2,J)
	 SB2(1,J)=SB2(NX-1,J)
	 SB2(NX,J)=SB2(2,J)
	 SM2(1,J)=SM2(NX-1,J)
	 SM2(NX,J)=SM2(2,J)
	 CLRAD1(1,J)=CLRAD1(NX-1,J)
	 CLRAD1(NX,J)=CLRAD1(2,J)
	END DO	 
C
C  Calculate surface wind speed at various points and forcings for time evolution of variables
C
	DO I=2,NX-1
	 DO J=2,NY-1
	  VBAR=0.25*(V2(I-1,J)+V2(I,J)+V2(I,J+1)+V2(I-1,J+1))
	  VS=SQRT((U2(I,J)+UMEAN)**2+VBAR**2+US2)
C
	  FU=DXI*(S2(I,J)-S2(I-1,J))+BETA*Y(J)*VBAR-CDH*VS*U1(I,J)
C
	  UU=MAX(UMEAN,0.)
	  UD=MIN(UMEAN,0.)
	  UADV=-DXI*(UU*(U2(I,J)-U2(I-1,J))+UD*(U2(I+1,J)-U2(I,J)))
	  FU=FU+UADV
C
	  UBAR=0.25*(U2(I+1,J)+U2(I,J)+U2(I+1,J-1)+U2(I,J-1))
	  VS=SQRT(V2(I,J)**2+(UBAR+UMEAN)**2+US2)
C
	  FV=DYI*(S2(I,J)-S2(I,J-1))-BETA*(Y(J)-0.5*DY)*UBAR-
     1    CDH*VS*V1(I,J)
C
	  VADV=-DXI*(UU*(V2(I,J)-V2(I-1,J))+UD*(V2(I+1,J)-V2(I,J)))
	  FV=FV+VADV
C
	  W(I,J)=-H*(DXI*(U2(I+1,J)-U2(I,J))+DYI*(V2(I,J+1)-V2(I,J)))
C
	  UBAR=0.5*(U2(I+1,J)+U2(I,J))+UMEAN
	  VBAR=0.5*(V2(I,J+1)+V2(I,J))
	  VS=SQRT(UBAR**2+VBAR**2+US2)
	  IF(EPINT.EQ.'y')THEN
	   EP=1.-(1.-EP0)*(SB2(I,J)-SM2(I,J))/CHID
	   EP=MIN(EP,0.9)
c	   EP=MAX(EP,EP0)
         EP=MAX(EP,0.1)
        ELSE
	   EP=EP0
	  END IF
	  DENOM=MAX(10.0,(SB1(I,J)-SM1(I,J)))
C
	  SADV=-DXI*(UU*(SB2(I,J)-SB2(I-1,J))+UD*(SB2(I+1,J)-SB2(I,J)))
C
	  MEQ=W(I,J)+CK*VS*(S0I(I,J)+S01(I,J)-SB1(I,J)+H*SADV)/DENOM
	  IF(SB2(I,J).LT.(S2(I,J)-0.1))MEQ=0.0
C
	  IF(TAU.EQ.0.0)THEN
	   MF=0.0
	   M1(I,J)=MEQ
	   M2(I,J)=MEQ
        ELSE
	   MF=TAUI*(MEQ-M1(I,J))
	  END IF
C
	  RADEQ=RFEED*N2*(-W(I,J)+EP*(M2(I,J)-RAD0/(EP*N2)))+
     1   RAD0*(1.0-RADMOIST+RADMOIST*(SB1(I,J)-SM1(I,J))/CHID)
C
	  CLRADF=TAUCLI*(RADEQ-CLRAD1(I,J))
c	  CLRADF=CLRADF-5.0*(CLRAD1(I,J)-CLRAD1(I-1,J))*DXI
C
	  FS=N2*(-W(I,J)+EP*M2(I,J))-CLRAD1(I,J)
C
	  AMFAC=MIN(0.0,(W(I,J)-M2(I,J)))
C
	  FSB=AIH*(CK*VS*(S0I(I,J)+S01(I,J)-SB1(I,J))+
     1    AMFAC*(SB1(I,J)-SM1(I,J)))+SADV
C
	  FS0=0.0
	  IF(IOCEAN.EQ.'y')THEN
	   FS0=(-CK*VS*(S0I(I,J)+S01(I,J)-SB1(I,J))+SOLAR(I,J))*SFAC+
     1    CLFAC*(CLRAD1(I,J)-RAD0)-TAUOCEANI*S01(I,J)
	  END IF
C
C	  FAC=GAMMA*(M2(I,J)-RAD0/(EP*N2))
C	  FAC=MAX(FAC,0.0)
C
c	  FSM=(SB1(I,J)-SM1(I,J))*FAC/H
	  FSM=GFAC*(M2(I,J)*(SB1(I,J)-SM1(I,J))-RAD0*CHID/(EP*N2))
	  FSM=FSM+TAUMI*(SB1(I,J)-SM1(I,J)-CHID)
	  FSM=FSM+0.2*TAUMI*(SM1(I,J+1)+SM1(I,J-1)-2.*SM1(I,J))
	  SMADV=-DXI*(UU*(SM2(I,J)-SM2(I-1,J))+UD*(SM2(I+1,J)-SM2(I,J)))
	  FSM=FSM+SMADV
c	  FSM=FSM+7.0*(SM1(I+1,J)-SM1(I,J))*DXI
C	     
C	Advance in time
C
	  U3(I,J)=U1(I,J)+2.*DT*FU
	  V3(I,J)=V1(I,J)+2.*DT*FV
	  S3(I,J)=S1(I,J)+2.*DT*FS
	  M3(I,J)=M1(I,J)+2.*DT*MF
	  M3(I,J)=MAX(M3(I,J),0.0)
	  SB3(I,J)=SB1(I,J)+2.*DT*FSB
	  SB3(I,J)=MIN(SB3(I,J),S3(I,J))
	  S03(I,J)=S01(I,J)+2.*DT*FS0
	  SM3(I,J)=SM1(I,J)+2.*DT*FSM
	  SM3(I,J)=MIN(SM3(I,J),SB3(I,J))
C	  SM3(I,J)=MAX(SM3(I,J),(SB3(I,J)-2.*CHID))
	  CLRAD3(I,J)=CLRAD1(I,J)+2.*DT*CLRADF
C
	 END DO
	END DO
C
C	Apply Asselin Filter
C
	DO I=1,NX
	 DO J=1,NY
	  U1(I,J)=U2(I,J)+DAMP*(U1(I,J)+U3(I,J)-2.*U2(I,J))
	  V1(I,J)=V2(I,J)+DAMP*(V1(I,J)+V3(I,J)-2.*V2(I,J))
	  S1(I,J)=S2(I,J)+DAMP*(S1(I,J)+S3(I,J)-2.*S2(I,J))
	  S01(I,J)=S02(I,J)+DAMP*(S01(I,J)+S03(I,J)-2.*S02(I,J))
	  M1(I,J)=M2(I,J)+DAMP*(M1(I,J)+M3(I,J)-2.*M2(I,J))
	  SB1(I,J)=SB2(I,J)+DAMP*(SB1(I,J)+SB3(I,J)-2.*SB2(I,J))
	  SM1(I,J)=SM2(I,J)+DAMP*(SM1(I,J)+SM3(I,J)-2.*SM2(I,J)) 
	  CLRAD1(I,J)=CLRAD2(I,J)+DAMP*(CLRAD1(I,J)+CLRAD3(I,J)-
     1	  2.*CLRAD2(I,J))
C
	  U2(I,J)=U3(I,J)
	  V2(I,J)=V3(I,J)
	  S2(I,J)=S3(I,J)
	  S02(I,J)=S03(I,J)
	  M2(I,J)=M3(I,J)
	  SB2(I,J)=SB3(I,J)
	  SM2(I,J)=SM3(I,J)
	  CLRAD2(I,J)=CLRAD3(I,J)
	 END DO
	END DO
C
C	Create time series
C
	IF(TPRINT.GE.PRINTTIME)THEN
	 TPRINT=0.0
	 IPRINT=IPRINT+1
	 TI(IPRINT)=TIME/(3600.*24.)
	 DO I=1,NX
	  DO K=1,NY2-1
	   K2=NY-K
	   UHOV(I,IPRINT)=UHOV(I,IPRINT)+(U2(I,K)+U2(I,K2))
	   UHOVA(I,IPRINT)=UHOVA(I,IPRINT)+(U2(I,K)-U2(I,K2))
	   VHOV(I,IPRINT)=VHOV(I,IPRINT)+(V2(I,K)+V2(I,K2))
	   VHOVA(I,IPRINT)=VHOVA(I,IPRINT)+(V2(I,K)-V2(I,K2))
	   WHOV(I,IPRINT)=WHOV(I,IPRINT)+100.*(W(I,K)+W(I,K2))
	   WHOVA(I,IPRINT)=WHOVA(I,IPRINT)+100.*(W(I,K)-W(I,K2))
	  END DO
	  UHOV(I,IPRINT)=UHOV(I,IPRINT)*FNYI
	  UHOVA(I,IPRINT)=UHOVA(I,IPRINT)*FNYI
	  VHOV(I,IPRINT)=VHOV(I,IPRINT)*FNYI
	  VHOVA(I,IPRINT)=VHOVA(I,IPRINT)*FNYI
	  WHOV(I,IPRINT)=WHOV(I,IPRINT)*FNYI
	  WHOVA(I,IPRINT)=WHOVA(I,IPRINT)*FNYI
c	  WHOV(I,IPRINT)=100.*W(I,NY2)
c	  UHOV(I,IPRINT)=U2(I,NY2)
c	  VHOV(I,IPRINT)=V2(I,NY2)
	 END DO
	END IF
C
C  Create movie
C
	IF(TIME.GT.(ENDTIME-TMOVIE-10.))THEN
	IF(TPRINTM.GE.PRINTMOVIE)THEN
	 TPRINTM=0.0
	 IPRINTM=IPRINTM+1
	 DO I=1,NX
	  DO J=1,NY
	   MMOVIE(I,J,IPRINTM)=100.*M2(I,J)
	   UMOVIE(I,J,IPRINTM)=U2(I,J)
	   VMOVIE(I,J,IPRINTM)=V2(I,J)
	   THMOVIE(I,J,IPRINTM)=EXP(S2(I,J)/(CP*DTEMP))
	  END DO
	 END DO
	END IF
	END IF
C
	GOTO 100
C
C	End of time integration
C
  300	CONTINUE
C
C Open output files
C
	OPEN(UNIT=10,FILE='output/t.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=11,FILE='output/x.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=12,FILE='output/y.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=13,FILE='output/u.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=14,FILE='output/v.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=15,FILE='output/w.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=16,FILE='output/m.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=17,FILE='output/s.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=18,FILE='output/sb.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=19,FILE='output/s0.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=20,FILE='output/sm.out',
     1   STATUS='UNKNOWN')
C
C  Print
C
	DO K=1,IPRINT
	 WRITE(10,200)TI(K)
	END DO
	WRITE(11,150)(0.001*X(I),I=2,NX-1)
	WRITE(12,150)(0.001*Y(J),J=2,NY-1)
C
	DO J=2,NY-1
	 WRITE(13,200)(U2(I,J),I=2,NX-1)
	 WRITE(14,200)(V2(I,J),I=2,NX-1)
	 WRITE(15,200)(100.*W(I,J),I=2,NX-1)
	 WRITE(16,200)(100.*M2(I,J),I=2,NX-1)
	 DO I=1,NX
	  DUM(I)=EXP(S2(I,J)/(CP*DTEMP))
	 END DO
	 WRITE(17,200)(DUM(I),I=2,NX-1)
	 DO I=1,NX
	  DUM(I)=EXP(SB2(I,J)/(CP*DTEMP))
	 END DO
	 WRITE(18,200)(DUM(I),I=2,NX-1)
	 DO I=1,NX
	  DUM(I)=EXP((S0(I,J)+S02(I,J))/(CP*DTEMP))
	 END DO
	 WRITE(19,200)(DUM(I),I=2,NX-1)
	 DO I=1,NX
	  DUM(I)=EXP(SM2(I,J)/(CP*DTEMP))
	 END DO
	 WRITE(20,200)(DUM(I),I=2,NX-1)
	END DO
C
  150	FORMAT(2X,81(F7.0,1X))
  200	FORMAT(2X,81(F7.2,1X))
  220 FORMAT(I2.2)
  221 FORMAT(I1)
C
	OPEN(UNIT=21,FILE='output/whov.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=22,FILE='output/uhov.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=23,FILE='output/vhov.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=24,FILE='output/uhova.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=25,FILE='output/vhova.out',
     1   STATUS='UNKNOWN')
	OPEN(UNIT=26,FILE='output/whova.out',
     1   STATUS='UNKNOWN')
	DO K=1,IPRINT
	 WRITE(21,200)(WHOV(I,K),I=2,NX-1)
	 WRITE(22,200)(UHOV(I,K),I=2,NX-1)
	 WRITE(23,200)(VHOV(I,K),I=2,NX-1)
	 WRITE(24,200)(UHOVA(I,K),I=2,NX-1)
	 WRITE(25,200)(VHOVA(I,K),I=2,NX-1)
	 WRITE(26,200)(WHOVA(I,K),I=2,NX-1)
	END DO
	OPEN(UNIT=30,FILE='output/jlength.out',STATUS='UNKNOWN')
	WRITE(30,*) IPRINTM
	CLOSE(21)
	CLOSE(22)
	CLOSE(23)
	CLOSE(24)
	CLOSE(25)
	CLOSE(26)
C
	DO K=1,IPRINTM
	 K2=K+40
	 WRITE(KNAME,220)K
	 ANAME='output/mmov'//KNAME
	 OPEN(UNIT=K2,FILE=ANAME,STATUS='UNKNOWN')
	 DO J=2,NY-1
	  WRITE(K2,200)(MMOVIE(I,J,K),I=2,NX-1)
	 END DO
	 CLOSE(UNIT=K2)
	 ANAME='output/umov'//KNAME
	 OPEN(UNIT=K2,FILE=ANAME,STATUS='UNKNOWN')
	 DO J=2,NY-1
	  WRITE(K2,200)(UMOVIE(I,J,K),I=2,NX-1)
	 END DO
	 CLOSE(UNIT=K2)
	 ANAME='output/vmov'//KNAME
	 OPEN(UNIT=K2,FILE=ANAME,STATUS='UNKNOWN')
	 DO J=2,NY-1
	  WRITE(K2,200)(VMOVIE(I,J,K),I=2,NX-1)
	 END DO
	 CLOSE(UNIT=K2)
	 ANAME='output/tmov'//KNAME
	 OPEN(UNIT=K2,FILE=ANAME,STATUS='UNKNOWN')
	 DO J=2,NY-1
	  WRITE(K2,200)(THMOVIE(I,J,K),I=2,NX-1)
	 END DO
	 CLOSE(UNIT=K2)
	END DO
C
	STOP
	END