      PROGRAM WINDOW
C ===============================
C David E. Dietrich, September, 1995
C ===========================================================
C MAX-lat = 89., MIN-lat = -90., MAX-lon = 359., MIN-lon = 0.
C ===========================================================
C
C Northern hemisphere,
C Southernmost cell centered at 0.5 deg latitude
C Northernmost cell centered at 80.5 deg latitude
C
      INTEGER*2 IN
      DIMENSION T(360,180,33),S(360,180,33)
C TK and IN include ghost zones
      DIMENSION TK(362,182),SK(362,182),IN(362,182)
C ===============
C Input data file
      OPEN(7,file='newlev_sum.dat',form='unformatted',status='old')
      REWIND 7
      READ(7) T,S
      CLOSE(7)
C Output data files
      OPEN(10,file='inview',status='unknown')
      OPEN(11,file='landfills',status='unknown')
      REWIND 10
      REWIND 11
      OPEN(70,file='modelev_sum.dat',form='unformatted')
      REWIND 70
C ===================
C Dietrich, May, 1995
C Use Poisson equation to get smooth fields over land.
C This avoids difficulties with slight mismatches between Levitus
C  and etopo5 dataset land boundaries.
C This can also be used to relax out bad Levitus data by setting
C  IN(I,J)=1 at bad data locations (but should be done only at the
C  bad K-levels)
C Dietrich, May, 1995
C ===================
C LPER=1 for periodic b.c.'s in longitudinal (I) direction
C We solve T globally, even though we write out only the model region,
C so we ALWAYS use the periodic option
      LPER=1
C Show western hemisphere only
C Show model working region only (106 E to 121 E)
      DO 30 K=1,33
      SMIN=1.E10
      SMAX=-1.E10
      sbar0=0
      tbar0=0
      DO 1000 J=91,180
      sbar0=sbar
      tbar0=tbar
      SBAR=0.
      TBAR=0.
      SUMIN=0.
      DO 600 I=1,360
      IN(I+1,J+1)=1
      IF (T(I,J,K).LT.-10..OR.S(I,J,K).LT.-10.) IN(I+1,J+1)=0
      IF (IN(I+1,J+1).EQ.1) SMIN=MIN(SMIN,S(I,J,K))
      IF (IN(I+1,J+1).EQ.1) SMAX=MAX(SMAX,S(I,J,K))
      SBAR=SBAR+IN(I+1,J+1)*S(I,J,K)
      TBAR=TBAR+IN(I+1,J+1)*T(I,J,K)
 600  SUMIN=SUMIN+IN(I+1,J+1)
      IF (SUMIN.NE.0.) SBAR=SBAR/SUMIN
      IF (SUMIN.NE.0.) TBAR=TBAR/SUMIN
      if (sbar.eq.0.0) sbar=sbar0
      if (tbar.eq.0.0) tbar=tbar0
C     WRITE(10,7) K,J,SBAR,TBAR,SUMIN
      WRITE(*,7) K,J,SBAR,TBAR,SUMIN
      DO 900 I=1,360
      SK(I+1,J+1)=IN(I+1,J+1)*S(I,J,K)+(1-IN(I+1,J+1))*SBAR
 900  TK(I+1,J+1)=IN(I+1,J+1)*T(I,J,K)+(1-IN(I+1,J+1))*TBAR
C This allows user to see region shade map
 1000 WRITE(10,11) (IN(I+1,J+1),I=106,121)
      sbar0=0
      tbar0=0
      DO 10 J=90,1,-1
      sbar0=sbar
      tbar0=tbar
      SBAR=0.
      TBAR=0.
      SUMIN=0.
      DO 6 I=1,360
      IN(I+1,J+1)=1
      IF (T(I,J,K).LT.-10..OR.S(I,J,K).LT.-10.) IN(I+1,J+1)=0
      IF (IN(I+1,J+1).EQ.1) SMIN=MIN(SMIN,S(I,J,K))
      IF (IN(I+1,J+1).EQ.1) SMAX=MAX(SMAX,S(I,J,K))
      SBAR=SBAR+IN(I+1,J+1)*S(I,J,K)
      TBAR=TBAR+IN(I+1,J+1)*T(I,J,K)
 6    SUMIN=SUMIN+IN(I+1,J+1)
      IF (SUMIN.NE.0.) SBAR=SBAR/SUMIN
      IF (SUMIN.NE.0.) TBAR=TBAR/SUMIN
      if (sbar.eq.0.0) sbar=sbar0
      if (tbar.eq.0.0) tbar=tbar0
C     WRITE(10,7) K,J,SBAR,TBAR,SUMIN
      WRITE(*,7) K,J,SBAR,TBAR,SUMIN
 7    FORMAT('At K=',I3,',J=',I4,
     1 ': SBAR=',F6.2,', TBAR=',F6.2,', SUMIN=',F6.2)
      DO 9 I=1,360
      SK(I+1,J+1)=IN(I+1,J+1)*S(I,J,K)+(1-IN(I+1,J+1))*SBAR
 9    TK(I+1,J+1)=IN(I+1,J+1)*T(I,J,K)+(1-IN(I+1,J+1))*TBAR
C This allows user to see region shade map
 10   WRITE(10,11) (IN(I+1,J+1),I=106,121)
      WRITE(*,2222) K,SMIN,SMAX
 2222 FORMAT('K=',I3,', SMIN,SMAX=',2F10.5)
      WRITE(10,12) K
 11   FORMAT(180I1)
 12   FORMAT('level ',I2,': land points are zeros')
      WRITE(*,14) K
 14   FORMAT(/'Level ',I2,' Poisson land fills')
      WRITE(11,13)
 13   FORMAT('Fields before POISOL land fills')
      DO 15 J=2,181
      IF (K.LT.3) WRITE(11,16) J,(SK(I,J),I=2,361)
 15   IF (K.LT.3) WRITE(11,17) J,(TK(I,J),I=2,361)
 16   FORMAT(/'salinity at J=',I3/(20F6.2))
 17   FORMAT(/'temperature at J=',I3/(20F6.2))
      WRITE(*,18)
 18   FORMAT(/'salinity')
      CALL POISOL(SK,0.01,2,361,2,181,IN,LPER)
      WRITE(*,19)
 19   FORMAT(/'temperature')
      CALL POISOL(TK,0.01,2,361,2,181,IN,LPER)
      IF (K.LT.3) WRITE(11,20)
 20   FORMAT(/'Fields after POISOL land fills')
      DO 25 J=2,181
      IF (K.LT.3) WRITE(11,16) J,(SK(I,J),I=2,361)
 25   IF (K.LT.3) WRITE(11,17) J,(TK(I,J),I=2,361)
C Include South China Sea window only in output T, but with 1/6 deg res
C Western most point at 106.0 east, eastern most point at 121.0 east
C TK(106,J) is at 106 deg east
C T(1,J,K) is at 106.0 deg east, T(2,J,K) is at 106.16667 east
 30   WRITE(70) ((TK(I,J),I=2,361),J=2,181),((SK(I,J),I=2,361),J=2,181)
      END
C ----------------------------------------------------------------------
      SUBROUTINE POISOL(A,EPS,IL,IH,JL,JH,IN,LPER)
C ----------------------------------------------------------------------
      INTEGER*2 IN
      DIMENSION A(362,182),IN(362,182)
      DIMENSION B(362)
C Output every NCHECK iterations
      NCHECK=20
C Relaxation factor (overrelaxation if RLX>1)
      RLX=1.
      RNORM=(IH-IL+1)*(JH-JL+1)
      DO 50 J=JL,JH
      DO 50 I=IL,IH
 50   RNORM=RNORM-IN(I,J)
      RNORM=1./RNORM
      N=0
 100  N=N+1
      IF (LPER.EQ.0) GO TO 151
C Periodic option
      DO 150 J=JL,JH
      A(IL-1,J)=A(IH,J)
 150  A(IH+1,J)=A(IL,J)
      GO TO 159
C Neuman option
 151  DO 152 J=JL,JH
      A(IL-1,J)=A(IL,J)
 152  A(IH+1,J)=A(IH,J)
 159  DO 160 I=IL,IH
      A(I,JL-1)=A(I,JL)
 160  A(I,JH+1)=A(I,JH)
      RMSD=0.
      DO 250 J=JL,JH
      DO 230 I=IL,IH
 230  B(I)=(1-IN(I,J))*.25*(A(I-1,J)+A(I+1,J)+A(I,J+1)+A(I,J-1))
     1 +IN(I,J)*A(I,J)
      DO 240 I=IL,IH
      TMP=RLX*(B(I)-A(I,J))
      A(I,J)=A(I,J)+TMP
 240  RMSD=RMSD+TMP**2
 250  CONTINUE
      RMSD=SQRT(RMSD*RNORM)
      IF (MOD(N,NCHECK).NE.1.AND.RMSD.GT.EPS) GO TO 600
      TMIN=1.E10
      TMAX=-1.E10
      DO 280 J=JL,JH
      DO 280 I=IL,IH
      TMIN=MIN(TMIN,A(I,J))
      TMAX=MAX(TMAX,A(I,J))
 280  CONTINUE
      WRITE(*,284) N,TMIN,TMAX,RMSD
      IF (RMSD.GT.EPS) GO TO 100
      WRITE(10,284) N,TMIN,TMAX,RMSD
 284  FORMAT('iteration #',I5,
     1 1P,', min=',E12.5,', max=',1E12.5,', RMSchange= ',E8.1)
 600  IF (RMSD.GT.EPS) GO TO 100
      END
