C
      SUBROUTINE SURFF(SSTM,WTSURF,WSSURF,WUSURF,WVSURF,IMM,JMM,ISBC)
C
C********************************************************************
C
C    An example of a subroutine to calculate surface forcing for POM
C    (windstress-WUSURF/WVSURF, heatflux-WTSURF, saltflux-WSSURF)
C
C     COADS MONTHLY SURF. FORCING FIELDS WERE ALREADY INTERPOLATED
C     ONTO MODEL GRID AND STORED IN ARRAYS THAT INCLUDE 2-MONTHS
C     AT A TIME TO ALLOW TIME INTERPOLATION EVERY TIME STEP
C     (e.g., 1/2 2/3 3/4...12/13 ; mon=13 is annual mean)
C     0<FR<1 is used for time-interpolation between months    T.E.
C********************************************************************
C 
C    SURFACE FORCING OPTIONS:
C           ZERO  CONS. (1 MON)  CONS. (ANN)   VAR. (MON)   TSURF
C    ISBC=   0  ;    1-12      ;      13     ;       14  ;  15
C
      INCLUDE 'grid'
      REAL JDAY
      COMMON /DATES/ JDAY,MON,MDAY,JYEAR
      COMMON/SINP/ UW(IM,JM,2),VW(IM,JM,2),UVW(IM,JM,2),
     2          QQ(IM,JM,2),DQDT(IM,JM,2),SST(IM,JM,2),SWR(IM,JM,2),
     3          SSS(IM,JM,2),EMP(IM,JM,2),FR,MM,JINT
      DIMENSION WUSURF(IMM,JMM),WVSURF(IMM,JMM),WTSURF(IMM,JMM),
     1           WSSURF(IMM,JMM),SWRAD(IMM,JMM),SSTM(IMM,JMM)
      DATA MM,JINT/0,0/
      JINT=JINT+1
      FROLD=FR
C
C     WRITE(6,'(''in surff '',I8,F8.3,3I8)')JINT,JDAY,MON,MDAY,JYEAR
C
         IF(MDAY.LE.15)THEN
      FR=(15.-MDAY)/28.
         ELSE
      FR=(44.-MDAY)/28.
         ENDIF
C
C     IF CONSTANT FORCING (BUT VAR. FEEDBACK TERM)
C
      IF(ISBC.EQ.0)RETURN
      IF(ISBC.LE.13)THEN
      IF(JINT.NE.1)GO TO 382
C ---- read monthly  windstress, hf, dq/dT, sst, E-P, sss
      DO 371 K=1,ISBC
      READ(60) UW,VW
      READ(65) QQ,DQDT,SST
      READ(66) EMP,SSS
      WRITE(6,'(/,A30,I5,/)')'%%%%% READ SURF. BC  MON= ',K
  371 CONTINUE
      GO TO 382
      ENDIF
C
C     IF VAR. MONTHLY FORCING
C
      MREAD=0
      IF(JINT.EQ.1)THEN
      MREAD=MON
      IF(MDAY.LE.15)MREAD=MON-1
      IF(MREAD.EQ.0)MREAD=12
          ENDIF
      IF(FR-FROLD.GT.0.9)MREAD=1
C
      DO 381 K=1,MREAD
      MM=MM+1
C                             (READ 13 is ANNUAL)
          IF(MM.EQ.13)THEN
      MM=1
      READ(60) UW,VW
      READ(65) QQ,DQDT,SST
      READ(66) EMP,SSS
      REWIND 60
      REWIND 65
      REWIND 66
           ENDIF
      M1=MM+1
      WRITE(6,'(/,A30,2I5,/)')'%%%%% READ SURF. BC  MON= ',MM,M1
      WRITE(6,'(I8,2F8.3,A10,/)')JINT,JDAY,FR
      READ(60) UW,VW
      READ(65) QQ,DQDT,SST
      READ(66) EMP,SSS
  381 CONTINUE
C
  382 CONTINUE
C
C--------------------------------- CALC WIND STRESS ---------------
      RHOA=1.226E-3
      RHOW=1.025
      RR=RHOA/RHOW
C
      DO 385 J=1,JM
      DO 385 I=1,IM
      UW1=UW(I,J,2)+FR*(UW(I,J,1)-UW(I,J,2))
      VW1=VW(I,J,2)+FR*(VW(I,J,1)-VW(I,J,2))
      SS1=SST(I,J,2)+FR*(SST(I,J,1)-SST(I,J,2))
C
      WUSURF(I,J)=-1.*UW1*RR 
      WVSURF(I,J)=-1.*VW1*RR
C
C--------------------------------------- CALC HEAT FLUX   ---------------
C (WTSURF= NET HF MINUS SHORT WAVE RAD; SWR IS USED DIRECTLY IN PROFT)
C
      QQ1=QQ(I,J,2)+FR*(QQ(I,J,1)-QQ(I,J,2))
C     QQ2=SWR(I,J,2)+FR*(SWR(I,J,1)-SWR(I,J,2))
      DQ1=DQDT(I,J,2)+FR*(DQDT(I,J,1)-DQDT(I,J,2))
C
      WTSURF(I,J)=-(QQ1+DQ1*(SSTM(I,J)-SS1))/(4.1876E+6)
      SWRAD(I,J)=0.
C     SWRAD(I,J)=-(QQ2)/(4.1876E+6)
C
C ------- TSURF as BC instead of fluxes -------------------------------
C
C     IF(ISBC.EQ.15) THEN
C         SSTM(I,J)=SS1
C       WTSURF(I,J)=0.
C        SWRAD(I,J)=0.
C     ENDIF
C
C--------------------------------------- CALC SALT FLUX   ---------------
C (E-P is in mm/3h, change to psu*m/s)
C
      EE1=EMP(I,J,2)+FR*(EMP(I,J,1)-EMP(I,J,2))
      SL1=SSS(I,J,2)+FR*(SSS(I,J,1)-SSS(I,J,2))
C
      WSSURF(I,J)=-1.*(EE1*SL1)/(1000.*3600.*3.)
C     WSSURF(I,J)=0.
C
  385 CONTINUE
C ---- if short wave radiation is used in PROFT
      WTSURF=WTSURF-SWRAD
C
      RETURN
      END
C