!        Generated by TAPENADE     (INRIA, Tropics team)
!  Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!  
! $Id: ice_flux_in.F,v 1.14 2004/03/01 16:56:23 eclare Exp $
!=======================================================================
!
!BOP
!
! !MODULE: ice_flux_in - reads and interpolates input forcing data
!
! !DESCRIPTION:
!
! Reads and interpolates forcing data for atmospheric and oceanic quantities.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !INTERFACE:
!
MODULE ICE_FLUX_IN
  USE ice_calendar
  USE ice_constants
  USE ice_diagnostics
  USE ice_domain
  USE ice_fileunits
  USE ice_flux
  USE ice_grid
  USE ice_history
  USE ice_kinds_mod
  USE ice_model_size
  USE ice_ocean
  USE ice_read_write
  USE ice_work
  CHARACTER(len=char_len_long) :: atm_data_dir, ocn_data_dir
  REAL(KIND=DBL_KIND) :: cldf(ilo:ihi, jlo:jhi)
  REAL(KIND=DBL_KIND) :: c1intp, c2intp, ftime
  INTEGER(KIND=INT_KIND) :: fyear
  INTEGER(KIND=INT_KIND) :: oldrecnum=0
  CHARACTER(len=char_len_long) :: flw_file, fsw_file, height_file, &
&  humid_file, pott_file, rain_file, rhoa_file, sss_file, sst_file, &
&  tair_file, uwind_file, vwind_file
  REAL(KIND=DBL_KIND) :: cldf_data(ilo:ihi, jlo:jhi, 2), flw_data(ilo:&
&  ihi, jlo:jhi, 2), fsnow_data(ilo:ihi, jlo:jhi, 2), fsw_data(ilo:ihi, &
&  jlo:jhi, 2), pott_data(ilo:ihi, jlo:jhi, 2), qa_data(ilo:ihi, jlo:jhi&
&  , 2), rhoa_data(ilo:ihi, jlo:jhi, 2), sss_data(ilo:ihi, jlo:jhi, 2), &
&  sst_data(ilo:ihi, jlo:jhi, 2), tair_data(ilo:ihi, jlo:jhi, 2), &
&  uatm_data(ilo:ihi, jlo:jhi, 2), vatm_data(ilo:ihi, jlo:jhi, 2), &
&  zlvl_data(ilo:ihi, jlo:jhi, 2)
  INTEGER(KIND=INT_KIND) :: fyear_final, fyear_init, ycycle

CONTAINS
  SUBROUTINE INIT_GETFLUX()
    IMPLICIT NONE
    INTEGER :: arg1
    INTRINSIC MOD
!
! !DESCRIPTION:
!
! Determines the current and final years of the forcing cycle based on
! namelist input, initializes the forcing data filenames, and initializes
! surface temperature and salinity from data. \\
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! last year in forcing cycle
    fyear_final = fyear_init + ycycle - 1
! current year
    arg1 = nyr - year_init
    fyear = fyear_init + MOD(arg1, ycycle)
!
    IF (my_task .EQ. master_task) THEN
      WRITE(nu_diag, *) ' Initial forcing data year = ', fyear_init
      WRITE(nu_diag, *) ' Current forcing data year = ', fyear
      WRITE(nu_diag, *) ' Final   forcing data year = ', fyear_final
    END IF
!
!ech-later      call Tair_limit
!
    CALL NCAR_FILES(fyear)
  END SUBROUTINE INIT_GETFLUX
  SUBROUTINE GETFLUX()
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    INTEGER :: arg1
    INTRINSIC MOD
!
! !DESCRIPTION:
!
! Get forcing data and interpolate as necessary! 
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! current year
    arg1 = nyr - year_init
    fyear = fyear_init + MOD(arg1, ycycle)
!
! forcing time
    ftime = time
! for restarting
!-------------------------------------------------------------------
! Read and interpolate annual climatologies of SSS and SST.
! Restore model SST to data.
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Read and interpolate atmospheric data
!-------------------------------------------------------------------
    time_forc = ftime
!
!
!      call sss_sst_restore
!
    CALL NCAR_BULK_DAT()
    CALL PREPARE_FORCING()
  END SUBROUTINE GETFLUX
  SUBROUTINE NCAR_FILES(yr)
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND),INTENT(IN) :: yr
!
! !DESCRIPTION:
!
! This subroutine is based on the LANL file naming conventions.
! Edit for other directory structures or filenames.
! Note: The year number in these filenames does not matter, because
!       subroutine file\_year will insert the correct year.
! 
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
! current forcing year
!
!EOP
!
    fsw_file = '/home/jkim/atm/swdn.1997.dat'
    CALL FILE_YEAR(fsw_file, yr)
!
    flw_file = '/home/jkim/atm/cldf.1997.dat'
    CALL FILE_YEAR(flw_file, yr)
!
    rain_file = '/home/jkim/atm/prec.1997.dat'
    CALL FILE_YEAR(rain_file, yr)
!
    uwind_file = '/home/jkim/atm/u_10.1997.dat'
    CALL FILE_YEAR(uwind_file, yr)
!
    vwind_file = '/home/jkim/atm/v_10.1997.dat'
    CALL FILE_YEAR(vwind_file, yr)
!
    tair_file = '/home/jkim/atm/t_10.1997.dat'
    CALL FILE_YEAR(tair_file, yr)
!
    humid_file = '/home/jkim/atm/q_10.1997.dat'
    CALL FILE_YEAR(humid_file, yr)
!
    rhoa_file = '/home/jkim/atm/dn10.1997.dat'
    CALL FILE_YEAR(rhoa_file, yr)
! master_task
!
    IF (my_task .EQ. master_task) THEN
      WRITE(nu_diag, *) ''
      WRITE(nu_diag, *) 'Forcing data year = ', fyear
      WRITE(nu_diag, *) 'Atmospheric data files:'
      WRITE(nu_diag, *) fsw_file
      WRITE(nu_diag, *) flw_file
      WRITE(nu_diag, *) rain_file
      WRITE(nu_diag, *) uwind_file
      WRITE(nu_diag, *) vwind_file
      WRITE(nu_diag, *) tair_file
      WRITE(nu_diag, *) humid_file
      WRITE(nu_diag, *) rhoa_file
    END IF
  END SUBROUTINE NCAR_FILES
  SUBROUTINE NCAR_BULK_DAT()
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    INTEGER :: arg1
    DOUBLE PRECISION :: arg10
    INTEGER(KIND=INT_KIND) :: i, imx, ipx, ixx, j, maxrec, midmonth, &
&    recnum, recslot
    LOGICAL(KIND=LOG_KIND) :: read6, readm
    REAL(KIND=DBL_KIND) :: sec6hr
    INTRINSIC MOD, INT, DBLE
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! record numbers for neighboring months
! record number
! maximum record number
! spline slot for current record
! middle day of month
!
! number of seconds in 6 hours
!
!-------------------------------------------------------------------
! monthly data 
!
! Assume that monthly data values are located in the middle of the 
! month.
!-------------------------------------------------------------------
!
! data is given on 15th of every month
! Compute record numbers for surrounding months
    midmonth = 15
!      midmonth = fix(p5 * dble(daymo(month)))  ! exact middle
!
    maxrec = 12
    arg1 = month + maxrec - 2
    imx = MOD(arg1, maxrec) + 1
    ipx = MOD(month, maxrec) + 1
! other two points will be used
    IF (mday .GE. midmonth) imx = 99
! Determine whether interpolation will use values 1:2 or 2:3
! recslot = 2 means we use values 1:2, with the current value (2)
!  in the second slot
! recslot = 1 means we use values 2:3, with the current value (2)
!  in the first slot
    IF (mday .LT. midmonth) ipx = 99
!
! latter half of month
    recslot = 1
! first half of month
! Find interpolation coefficients
    IF (mday .LT. midmonth) recslot = 2
!
    CALL INTERP_COEFF_MONTHLY(recslot)
! Read 2 monthly values 
!
    readm = .false.
    IF (istep .EQ. 1 .OR. mday .EQ. midmonth .AND. sec .EQ. 0) readm = &
&        .true.
!
    CALL READ_DATA(readm, 0, fyear, imx, month, ipx, maxrec, fsw_file, &
&             fsw_data)
    CALL READ_DATA(readm, 0, fyear, imx, month, ipx, maxrec, flw_file, &
&             cldf_data)
    CALL READ_DATA(readm, 0, fyear, imx, month, ipx, maxrec, rain_file, &
&             fsnow_data)
!
    CALL INTERPOLATE_DATA(fsw_data, fsw)
    CALL INTERPOLATE_DATA(cldf_data, cldf)
    CALL INTERPOLATE_DATA(fsnow_data, fsnow)
!-------------------------------------------------------------------
! 6-hourly data
! 
! Assume that the 6-hourly value is located at the end of the
!  6-hour period.  This is the convention for NCEP reanalysis data.
!  E.g. record 1 gives conditions at 6 am GMT on 1 January.
!-------------------------------------------------------------------
!
!
! seconds in 6 hours
    sec6hr = secday/c4
! 365*4
! current record number
    maxrec = 1460
!
! Compute record numbers for surrounding data (2 on each side)
    arg10 = DBLE(sec)/sec6hr
    recnum = 4*INT(yday) - 3 + INT(arg10)
!
!
    arg1 = recnum + maxrec - 2
    imx = MOD(arg1, maxrec) + 1
! Compute interpolation coefficients
! If data is located at the end of the time interval, then the
!  data value for the current record goes in slot 2
    arg1 = recnum - 1
    ixx = MOD(arg1, maxrec) + 1
!      ipx = mod(recnum,         maxrec) + 1
!
!
    recslot = 2
    ipx = 99
    CALL INTERP_COEFF(recnum, recslot, sec6hr)
! Read
!
    read6 = .false.
    IF (istep .EQ. 1 .OR. oldrecnum .NE. recnum) read6 = .true.
!
    CALL READ_DATA(read6, 0, fyear, imx, ixx, ipx, maxrec, tair_file, &
&             tair_data)
    CALL READ_DATA(read6, 0, fyear, imx, ixx, ipx, maxrec, uwind_file, &
&             uatm_data)
    CALL READ_DATA(read6, 0, fyear, imx, ixx, ipx, maxrec, vwind_file, &
&             vatm_data)
    CALL READ_DATA(read6, 0, fyear, imx, ixx, ipx, maxrec, rhoa_file, &
&             rhoa_data)
    CALL READ_DATA(read6, 0, fyear, imx, ixx, ipx, maxrec, humid_file, &
&             qa_data)
! Interpolate
!
    CALL INTERPOLATE_DATA(tair_data, tair)
    CALL INTERPOLATE_DATA(uatm_data, uatm)
    CALL INTERPOLATE_DATA(vatm_data, vatm)
    CALL INTERPOLATE_DATA(rhoa_data, rhoa)
    CALL INTERPOLATE_DATA(qa_data, qa)
! Save record number
!
    oldrecnum = recnum
  END SUBROUTINE NCAR_BULK_DAT
  SUBROUTINE READ_DATA(flag, recd, yr, imx, ixx, ipx, maxrec, data_file&
&    , field_data)
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    CHARACTER(len=char_len_long) :: data_file
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi, 2),INTENT(OUT) :: &
&    field_data
    LOGICAL(KIND=LOG_KIND),INTENT(IN) :: flag
    INTEGER(KIND=INT_KIND),INTENT(IN) :: imx
    INTEGER(KIND=INT_KIND),INTENT(IN) :: ipx
    INTEGER(KIND=INT_KIND),INTENT(IN) :: ixx
    INTEGER(KIND=INT_KIND),INTENT(IN) :: maxrec
    INTEGER(KIND=INT_KIND),INTENT(IN) :: recd
    INTEGER(KIND=INT_KIND),INTENT(IN) :: yr
    INTEGER(KIND=INT_KIND) :: arg1
    REAL(KIND=DBL_KIND) :: field_dummy(ilo:ihi, jlo:jhi)
    INTEGER(KIND=INT_KIND) :: arg, n2, n4, nbits, nrec
    LOGICAL(KIND=LOG_KIND) :: diag, scatter
!
! !DESCRIPTION:
!
! If data is at the beginning of a one-year record, get data from
!  the previous year.
! If data is at the end of a one-year record, get data from the 
!  following year.
! If no earlier data exists (beginning of fyear\_init), then \\
!  (1) For monthly data, get data from the end of fyear\_final. \\
!  (2) For more frequent data, let the imx value equal the 
!      first value of the year. \\
! If no later data exists (end of fyear\_final), then \\
!  (1) For monthly data, get data from the beginning of fyear\_init. \\
!  (2) For more frequent data, let the ipp and, if necessary, ipx 
!      values equal the last value of the year. \\
! In other words, we assume persistence when daily or 6-hourly
!   data is missing, and we assume periodicity when monthly data
!   is missing. \\
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !
!
! !INPUT/OUTPUT PARAMETERS:
!
!
! baseline record number
! year of forcing data
! record numbers of 3 data values
! relative to recd
! maximum record value
!
! 2 values needed for interpolation
!
!EOP
!
! data file to be read 
!
! = 32 for single precision, 64 for double
! record number to read
! like imx and ipx, but 
! adjusted at beginning and end of data
! value of time argument in field_data
!
!
! double precision data
    nbits = 64
!
! scatter data to processors
    scatter = .true.
! write diagnostic info
    diag = .false.
!
!! debugging
    IF (istep1 .GT. check_step) diag = .true.
! flag
!
!      if (my_task==master_task .and. (diag)) then
!         write(nu_diag,*) '  ',data_file
!      endif
!
    IF (flag) THEN
!-----------------------------------------------------------------
! Initialize record counters
! (n2, n4 will change only at the very beginning or end of
!  a forcing cycle.)
!-----------------------------------------------------------------
      n2 = imx
      n4 = ipx
!-----------------------------------------------------------------
! read data
!-----------------------------------------------------------------
      arg = 0
! imx ne 99
! always read ixx data from data file for current year
!
      IF (imx .NE. 99) THEN
! ixx <= 1
! currently in first half of data interval
        IF (ixx .LE. 1) THEN
! yr > fyear_init
          IF (yr .GT. fyear_init) THEN
! get data from previous year
            arg1 = yr - 1
            CALL FILE_YEAR(data_file, arg1)
          ELSE IF (maxrec .GT. 12) THEN
! yr = fyear_init, no prior data exists
! extrapolate from first record
            IF (ixx .EQ. 1) n2 = ixx
          ELSE
! go to end of fyear_final
            CALL FILE_YEAR(data_file, fyear_final)
          END IF
        END IF
        CALL ICE_OPEN(nu_forcing, data_file, nbits)
!
        arg = 1
        nrec = recd + n2
!jkim        call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                'rda8', scatter, diag)
        CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, &
&                diag)
        field_data(:, :, arg) = field_dummy(:, :)
!
        IF (ixx .EQ. 1 .AND. my_task .EQ. master_task) CLOSE(nu_forcing)&
&
      END IF
!
      CALL FILE_YEAR(data_file, yr)
      CALL ICE_OPEN(nu_forcing, data_file, nbits)
!
      arg = arg + 1
      nrec = recd + ixx
!jkim      call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                'rda8', scatter, diag)
      CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, diag&
&             )
      field_data(:, :, arg) = field_dummy(:, :)
! ipx ne 99
!
      IF (ipx .NE. 99) THEN
! ixx = maxrec
! currently in latter half of data interval
        IF (ixx .EQ. maxrec) THEN
! yr < fyear_final
          IF (yr .LT. fyear_final) THEN
! get data from following year
            IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
            arg1 = yr + 1
            CALL FILE_YEAR(data_file, arg1)
            CALL ICE_OPEN(nu_forcing, data_file, nbits)
          ELSE IF (maxrec .GT. 12) THEN
! yr = fyear_final, no more data exists
! extrapolate from ixx
            n4 = ixx
          ELSE
! go to beginning of fyear_init
            IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
            CALL FILE_YEAR(data_file, fyear_init)
            CALL ICE_OPEN(nu_forcing, data_file, nbits)
          END IF
        END IF
!
        arg = arg + 1
        nrec = recd + n4
!jkim        call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                'rda8', scatter, diag)
        CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, &
&                diag)
        field_data(:, :, arg) = field_dummy(:, :)
      END IF
!
      IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
    END IF
  END SUBROUTINE READ_DATA
  SUBROUTINE FILE_YEAR(data_file, yr)
    IMPLICIT NONE
    CHARACTER(len=char_len_long),INTENT(INOUT) :: data_file
    INTEGER(KIND=INT_KIND),INTENT(IN) :: yr
    INTEGER(KIND=INT_KIND) :: i
    CHARACTER(len=char_len_long) :: tmpname
    INTRINSIC INDEX
!
! !DESCRIPTION:
!
! Construct the correct name of the atmospheric data file
! to be read, given the year and assuming the naming convention
! that filenames end with 'yyyy.dat'.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!
!EOP
!
!
!
!
    i = INDEX(data_file, '.dat') - 5
    tmpname = data_file
    WRITE(data_file, '(a,i4.4,a)') tmpname(1:i), yr, '.dat'
  END SUBROUTINE FILE_YEAR
  SUBROUTINE INTERP_COEFF(recnum, recslot, secint)
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND),INTENT(IN) :: recnum
    INTEGER(KIND=INT_KIND),INTENT(IN) :: recslot
    REAL(KIND=DBL_KIND),INTENT(IN) :: secint
    REAL(KIND=DBL_KIND),PARAMETER :: secyr=c365*secday
    REAL(KIND=DBL_KIND) :: arg1
    REAL(KIND=DBL_KIND) :: t1, t2, tt
    INTRINSIC ABS, INT, DBLE
!
! !DESCRIPTION:
!
! Compute coefficients for interpolating data to current time step.
! Works for any data interval that divides evenly into a 365-day 
!  year (daily, 6-hourly, etc.)
! Use interp\_coef\_monthly for monthly data.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !USES:
!
! record number for current data value
! spline slot for current record 
!
! seconds in data interval
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! seconds in a 365-day year 
!
! seconds elapsed in current year
! seconds elapsed at data points
!
!
!jkim      tt = mod(ftime,secyr)
! Find neighboring times
    arg1 = ftime/secyr
    tt = ftime - INT(arg1)*secyr
! Compute coefficients
    IF (recslot .EQ. 2) THEN
! current record goes in slot 2 (NCEP)
      t2 = DBLE(recnum)*secint
!  - 1 interval
      t1 = t2 - secint
    ELSE
! recslot = 1
      t1 = DBLE(recnum)*secint
!  + 1 interval
      t2 = t1 + secint
    END IF
    IF ((t2-tt)/(t2-t1) .GE. 0.) THEN
      c1intp = (t2-tt)/(t2-t1)
    ELSE
      c1intp = -((t2-tt)/(t2-t1))
    END IF
    c2intp = c1 - c1intp
  END SUBROUTINE INTERP_COEFF
  SUBROUTINE PREPARE_FORCING()
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    REAL(KIND=DBL_KIND) :: arg1
    INTEGER(KIND=INT_KIND) :: i, j
    REAL(KIND=DBL_KIND) :: workx, worky
    REAL :: x1
    INTRINSIC COS, EXP, SIN, MAX, MIN, SQRT
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
    DO j=jlo,jhi
!-----------------------------------------------------------------
! Fix interpolation errors
!-----------------------------------------------------------------
      DO i=ilo,ihi
        IF (fsw(i, j) .LT. c0) THEN
          fsw(i, j) = c0
        ELSE
          fsw(i, j) = fsw(i, j)
        END IF
        IF (cldf(i, j) .GT. c1) THEN
          x1 = c1
        ELSE
          x1 = cldf(i, j)
        END IF
        IF (x1 .LT. c0) THEN
          cldf(i, j) = c0
        ELSE
          cldf(i, j) = x1
        END IF
        IF (fsnow(i, j) .LT. c0) THEN
          fsnow(i, j) = c0
        ELSE
          fsnow(i, j) = fsnow(i, j)
        END IF
        IF (rhoa(i, j) .LT. c0) THEN
          rhoa(i, j) = c0
        ELSE
          rhoa(i, j) = rhoa(i, j)
        END IF
        IF (qa(i, j) .LT. c0) THEN
          qa(i, j) = c0
        ELSE
          qa(i, j) = qa(i, j)
        END IF
!
! as in the dummy atm (latm)
        qa(i, j) = qa(i, j)*0.94_dbl_kind
! as in the dummy atm (latm)
!-----------------------------------------------------------------
! Compute other fields needed by model
!-----------------------------------------------------------------
        fsw(i, j) = fsw(i, j)*0.92_dbl_kind
!
!
        zlvl(i, j) = c10
! wind speed, m/s
        arg1 = uatm(i, j)**2 + vatm(i, j)**2
        wind(i, j) = SQRT(arg1)
! divide shortwave into spectral bands
        pott(i, j) = tair(i, j)
!
! visible direct
        swvdr(i, j) = fsw(i, j)*.28
! visible diffuse
        swvdf(i, j) = fsw(i, j)*.24
! near IR direct
        swidr(i, j) = fsw(i, j)*.31
! near IR diffuse
! as in the dummy atm (latm)
! longwave as in Parkinson and Washington (1979)
        swidf(i, j) = fsw(i, j)*.17
!
!!! downward longwave !!!
! determine whether precip is rain or snow
        arg1 = -(7.77e-4_dbl_kind*(tffresh-tair(i, j))**2)
        flw(i, j) = stefan_boltzmann*tair(i, j)**4*(c1-0.261_dbl_kind*&
&          EXP(arg1))*(c1+0.275_dbl_kind*cldf(i, j))
!
! mm/month -> kg/m^2 s
        fsnow(i, j) = fsnow(i, j)/2.592e+06_dbl_kind
        frain(i, j) = c0
!-----------------------------------------------------------------
! rotate zonal/meridional vectors to local coordinates
! Vector fields come in on T grid, but are oriented geographically
! need to rotate to pop-grid FIRST using ANGLET
! then interpolate to the U-cell centers  (otherwise we
! interpolate across the pole)
! use ANGLET which is on the T grid !
! atmo variables are needed in T cell centers in subroutine stability,
! and are interpolated to the U grid later as necessary
!-----------------------------------------------------------------
!        if (Tair(i,j) >= Tffresh) then
        IF (tair(i, j) .GE. tffresh) THEN
          frain(i, j) = fsnow(i, j)
          fsnow(i, j) = c0
        END IF
!
! wind velocity, m/s
        workx = uatm(i, j)
        worky = vatm(i, j)
! convert to POP grid
! note uatm, vatm, wind
        uatm(i, j) = workx*COS(anglet(i, j)) + worky*SIN(anglet(i, j))
!  are on the T-grid here
        vatm(i, j) = worky*COS(anglet(i, j)) - workx*SIN(anglet(i, j))
      END DO
    END DO
  END SUBROUTINE PREPARE_FORCING
  SUBROUTINE SSS_SST_RESTORE()
    IMPLICIT NONE
    INTEGER :: arg1
    INTEGER(KIND=INT_KIND) :: i, imx, ipx, j, maxrec, midmonth, recslot
    LOGICAL(KIND=LOG_KIND) :: readm
    REAL(KIND=DBL_KIND) :: sstdat(ilo:ihi, jlo:jhi)
    REAL(KIND=DBL_KIND) :: trestore
    INTRINSIC MOD
!
! !DESCRIPTION:
!
! Interpolate monthly sss, sst data to timestep.
! Restore sst computed by ice model to data.
! 
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! record numbers for neighboring months
! maximum record number
! spline slot for current record
! middle day of month
!
! data value toward which SST is restored
!
!
!
!jkim      sss_file = trim(ocn_data_dir)//'sss_Lev.mm'
!jkim      sst_file = trim(ocn_data_dir)//'sst_Lev.mm'
!
!      Trestore = dt       ! STRONG restoring 
!      Trestore = 1.296e6  ! 15 days
!      Trestore = 2.592e6  ! 30 days
! 90 days
    trestore = 7.776e6
!-------------------------------------------------------------------
! monthly data 
!
! Assume that monthly data values are located in the middle of the 
! month.
!-------------------------------------------------------------------
!
    IF (my_task .EQ. master_task .AND. istep .EQ. 1) THEN
      WRITE(nu_diag, *) ''
      WRITE(nu_diag, *) 'SSS and SST data interpolated to timestep:'
      WRITE(nu_diag, *) sss_file
      WRITE(nu_diag, *) sst_file
      WRITE(nu_diag, *) 'SST restoring timescale = ', trestore/secday, &
&      ' days'
    END IF
!
! data is given on 15th of every month
! Compute record numbers for surrounding months
    midmonth = 15
!      midmonth = fix(p5 * dble(daymo(month)))  ! exact middle
!
    maxrec = 12
    arg1 = month + maxrec - 2
    imx = MOD(arg1, maxrec) + 1
    ipx = MOD(month, maxrec) + 1
! other two points will be used
    IF (mday .GE. midmonth) imx = 99
! Determine whether interpolation will use values 1:2 or 2:3
! recslot = 2 means we use values 1:2, with the current value (2)
!  in the second slot
! recslot = 1 means we use values 2:3, with the current value (2)
!  in the first slot
    IF (mday .LT. midmonth) ipx = 99
!
! latter half of month
    recslot = 1
! first half of month
! Find interpolation coefficients
    IF (mday .LT. midmonth) recslot = 2
!
    CALL INTERP_COEFF_MONTHLY(recslot)
! Read 2 monthly values 
!
    readm = .false.
    IF (istep .EQ. 1 .OR. mday .EQ. midmonth .AND. sec .EQ. 0) readm = &
&        .true.
!
    CALL READ_CLIM_DATA(readm, 0, imx, month, ipx, sss_file, sss_data)
    CALL READ_CLIM_DATA(readm, 0, imx, month, ipx, sst_file, sst_data)
! Interpolate
!
    CALL INTERPOLATE_DATA(sss_data, sss)
    CALL INTERPOLATE_DATA(sst_data, sstdat)
!-------------------------------------------------------------------
! Restore sst to data
!-------------------------------------------------------------------
!
    DO j=jlo,jhi
      DO i=ilo,ihi
        IF (sss(i, j) .LT. c0) sss(i, j) = c0
        sst(i, j) = sst(i, j) + (sstdat(i, j)-sst(i, j))*dt/trestore
      END DO
    END DO
!
    CALL COMPLETE_GETFLUX_OCN()
  END SUBROUTINE SSS_SST_RESTORE
  SUBROUTINE INTERP_COEFF_MONTHLY(recslot)
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND),INTENT(IN) :: recslot
    REAL(KIND=DBL_KIND) :: arg1
    REAL(KIND=INT_KIND) :: daymid(0:13)
    REAL(KIND=DBL_KIND) :: t1, t2, tt
    INTRINSIC INT
!
! !DESCRIPTION:
!
! Compute coefficients for interpolating monthly data to current time step.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
! slot (1 or 2) for current record
!
!EOP
!
! seconds elapsed in current year
! seconds elapsed at month midpoint
!
! month mid-points
!
! time frame ends 0 sec into day 15
    daymid(1:13) = 14.
! Dec 15, 0 sec
! make time cyclic
    daymid(0) = -17.
!
!jkim      tt = mod(ftime/secday,c365)
! Find neighboring times
    arg1 = ftime/secday/c365
    tt = ftime/secday - INT(arg1)*c365
! Compute coefficients
    IF (recslot .EQ. 2) THEN
! first half of month
! midpoint, current month
      t2 = daycal(month) + daymid(month)
      IF (month .EQ. 1) THEN
! Dec 15 (0 sec)
        t1 = daymid(0)
      ELSE
! midpoint, previous month
        t1 = daycal(month-1) + daymid(month-1)
      END IF
    ELSE
! second half of month
! midpoint, current month
      t1 = daycal(month) + daymid(month)
! day 15 of next month (0 sec)
      t2 = daycal(month+1) + daymid(month+1)
    END IF
    c1intp = (t2-tt)/(t2-t1)
    c2intp = c1 - c1intp
  END SUBROUTINE INTERP_COEFF_MONTHLY
  SUBROUTINE READ_CLIM_DATA(readflag, recd, imx, ixx, ipx, data_file, &
&    field_data)
    IMPLICIT NONE
    CHARACTER(len=char_len_long),INTENT(IN) :: data_file
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi, 2),INTENT(OUT) :: &
&    field_data
    INTEGER(KIND=INT_KIND),INTENT(IN) :: imx
    INTEGER(KIND=INT_KIND),INTENT(IN) :: ipx
    INTEGER(KIND=INT_KIND),INTENT(IN) :: ixx
    LOGICAL(KIND=LOG_KIND),INTENT(IN) :: readflag
    INTEGER(KIND=INT_KIND),INTENT(IN) :: recd
    REAL(KIND=DBL_KIND) :: field_dummy(ilo:ihi, jlo:jhi)
    INTEGER(KIND=INT_KIND) :: arg, nbits, nrec
    LOGICAL(KIND=LOG_KIND) :: diag, scatter
!
! !DESCRIPTION:
!
! Read data needed for interpolation, as in read\_data.
! Assume a one-year cycle of climatological data, so that there is
!  no need to get data from other years or to extrapolate data beyond
!  the forcing time period.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!
! baseline record number
! record numbers of 3 data values
! relative to recd
!
!
! 2 values needed for interpolation
!
!EOP
!
! = 32 for single precision, 64 for double
! record number to read
! value of time argument in field_data
!
!
! double precision data
    nbits = 64
!
! scatter data to processors
    scatter = .true.
! write diagnostic info
    diag = .false.
!
!! debugging
    IF (istep1 .GT. check_step) diag = .true.
!
    IF (my_task .EQ. master_task .AND. diag) WRITE(nu_diag, *) '  ', &
&                                             data_file
! readflag
!
    IF (readflag) THEN
!-----------------------------------------------------------------
! read data
!-----------------------------------------------------------------
      CALL ICE_OPEN(nu_forcing, data_file, nbits)
!
      arg = 0
      IF (imx .NE. 99) THEN
        arg = 1
        nrec = recd + imx
!jkim        call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                 'rda8', scatter, diag)
        CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, &
&                diag)
        field_data(:, :, arg) = field_dummy(:, :)
      END IF
!
      arg = arg + 1
      nrec = recd + ixx
!jkim      call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                 'rda8', scatter, diag)
      CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, diag&
&             )
      field_data(:, :, arg) = field_dummy(:, :)
!
      IF (ipx .NE. 99) THEN
        arg = arg + 1
        nrec = recd + ipx
!jkim        call ice_read (nu_forcing, nrec, field_data(:,:,arg),
!jkim     &                 'rda8', scatter, diag)
        CALL ICE_READ(nu_forcing, nrec, field_dummy, 'rda8', scatter, &
&                diag)
        field_data(:, :, arg) = field_dummy(:, :)
      END IF
!
      IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
    END IF
  END SUBROUTINE READ_CLIM_DATA
  SUBROUTINE INTERPOLATE_DATA(field_data, field)
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: &
&    field
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi, 2),INTENT(IN) :: &
&    field_data
    INTEGER(KIND=INT_KIND) :: i, j
!
! !DESCRIPTION:
!
! Linear interpolation
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
! 2 values used for interpolation
!
! interpolated field
!
!EOP
!
!
    DO j=jlo,jhi
      DO i=ilo,ihi
        field(i, j) = c1intp*field_data(i, j, 1) + c2intp*field_data(i, &
&          j, 2)
      END DO
    END DO
  END SUBROUTINE INTERPOLATE_DATA
  SUBROUTINE SSS_CLIM()
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND) :: i, j, k, nbits
    LOGICAL(KIND=LOG_KIND) :: diag, scatter
!
! !DESCRIPTION:
!
! Creates annual mean climatology for Levitus sss from 12-month climatology.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
!
    scatter = .true.
    diag = .false.
!
! double precision data
    nbits = 64
!
!jkim      sss_file = trim(ocn_data_dir)//'sss_Lev.mm'
!
    IF (my_task .EQ. master_task) THEN
      WRITE(nu_diag, *) ''
      WRITE(nu_diag, *) 'SSS climatology computed from:'
      WRITE(nu_diag, *) sss_file
    END IF
!
    IF (my_task .EQ. master_task) CALL ICE_OPEN(nu_forcing, sss_file, &
&                                          nbits)
!-------------------------------------------------------------------
! create surface salinity climatology from monthly data
!-------------------------------------------------------------------
!
!
    DO j=jlo,jhi
      DO i=ilo,ihi
        worka(i, j) = c0
      END DO
    END DO
!
! loop over 12 months
    DO k=1,12
!
      CALL ICE_READ(nu_forcing, k, worka, 'rda8', scatter, diag)
      DO j=jlo,jhi
        DO i=ilo,ihi
          sss(i, j) = worka(i, j) + sss(i, j)
        END DO
      END DO
    END DO
!
! k
!
    DO j=jlo,jhi
      DO i=ilo,ihi
! annual average salinity
        sss(i, j) = sss(i, j)/c12
        IF (sss(i, j) .LT. c0) sss(i, j) = c0
      END DO
    END DO
! close file
!
    IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
!
    CALL COMPLETE_GETFLUX_OCN()
  END SUBROUTINE SSS_CLIM
  SUBROUTINE COMPLETE_GETFLUX_OCN()
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND) :: i, j
    INTRINSIC MAX
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!-------------------------------------------------------------------
! Compute freezing temperature based on salinity
! Make sure SST is not below Tf.
!-------------------------------------------------------------------
!
    DO j=jlo,jhi
      DO i=ilo,ihi
! freezing temp, deg C
        tf(i, j) = -(depresst*sss(i, j))
        IF (sst(i, j) .LT. tf(i, j)) THEN
          sst(i, j) = tf(i, j)
        ELSE
          sst(i, j) = sst(i, j)
        END IF
      END DO
    END DO
  END SUBROUTINE COMPLETE_GETFLUX_OCN
  SUBROUTINE SST_IC()
    IMPLICIT NONE
    INTEGER(KIND=INT_KIND) :: i, j, k, nbits
    LOGICAL(KIND=LOG_KIND) :: diag, scatter
    INTRINSIC MAX
!
! !DESCRIPTION:
!
! Reads sst data for current month, and adjusts sst based on freezing 
! temperature.  Does not interpolate.
!
! !REVISION HISTORY:
!
! authors Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
!
!
    scatter = .true.
    diag = .false.
!
! double precision data
!-----------------------------------------------------------------
! Shea, Trenberth and Reynolds SST
!-----------------------------------------------------------------
    nbits = 64
!
!      sst_file = '/n/filer3/climate/eclare/sst.sss/sst.mm.100x116.da'
!jkim      sst_file = trim(ocn_data_dir)//'sst_Lev.mm'
!
    IF (my_task .EQ. master_task) THEN
      WRITE(nu_diag, *) ''
      WRITE(nu_diag, *) 'SST initial condition:'
      WRITE(nu_diag, *) sst_file
    END IF
!
    IF (my_task .EQ. master_task) CALL ICE_OPEN(nu_forcing, sst_file, &
&                                          nbits)
!
    CALL ICE_READ(nu_forcing, month, sst, 'rda8', scatter, diag)
!
    IF (my_task .EQ. master_task) CLOSE(nu_forcing) 
!
    DO j=jlo,jhi
! Make sure sst is not less than Tf
      DO i=ilo,ihi
        IF (sst(i, j) .LT. tf(i, j)) THEN
          sst(i, j) = tf(i, j)
        ELSE
          sst(i, j) = sst(i, j)
        END IF
      END DO
    END DO
  END SUBROUTINE SST_IC
END MODULE ICE_FLUX_IN