      subroutine datefor(wday, iyr,mon,idy,ihr, yrflag)
      use mod_xc  ! HYCOM communication interface 
      implicit none
!
      real*8  wday
      integer iyr,mon,idy,ihr,yrflag
!
!**********
!*
!  1) convert date into 'model day', for yrflag=2,3,4 only.
!
!  2) for yrflag==4, all years are 365 days.
!     for example:
!      a) iyr=1,mon=1,idy=1 is the 1st day of the run,
!         so wday would be 1.0.
!      b) iyr=1,mon=1,idy=2 is the 2nd day of the run,
!         so wday would be 2.0.
!  
!  3) for yrflag==3, the 'model day' is the number of days since 
!     001/1901 (which is model day 1.0).
!     for example:
!      a) iyr=1901,mon=1,idy=1, represents 0000z hrs on 01/01/1901
!         so wday would be 1.0.
!      a) iyr=1901,mon=1,idy=2, represents 0000z hrs on 02/01/1901
!         so wday would be 2.0.
!     year must be no less than 1901.0, and no greater than 2099.0.
!     note that year 2000 is a leap year (but 1900 and 2100 are not).
!
!  4) for yrflag==2, all years are leap years.
!     for example:
!      a) iyr=1,mon=1,idy=1 is the 1st day of the run,
!         so wday would be 1.0.
!      a) iyr=1,mon=1,idy=2 is the 2nd day of the run,
!         so wday would be 2.0.
!*
!**********
!
      integer nleap
!
      integer month(13)
      data    month / 0,  31,  59,  90, 120, 151, 181, &
                         212, 243, 273, 304, 334, 365 /
!
      if     (yrflag.eq.4) then
        wday = 365.d0*(iyr-1) + month(mon) + idy-1 + ihr/24.0d0
      elseif (yrflag.eq.3) then
        nleap = (iyr-1901)/4
        wday  = 365.0d0*(iyr-1901) + nleap + month(mon) + idy + ihr/24.0d0
        if     (mod(iyr,4).eq.0 .and. mon.gt.2) then
          wday  = wday + 1.0d0
        endif
      elseif (yrflag.eq.2) then
        wday = 366.d0*(iyr-1) + month(mon) + idy-1 + ihr/24.0d0
        if     (mon.gt.2) then
          wday  = wday + 1.0d0
        endif
      else 
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in datefor - unsupported yrflag value'
        write(lp,*)
        endif !1st tile
        call xcstop('(datefor)')
               stop '(datefor)'
      endif !yrflag
      return
!     end of datefor
      end
!
!
      subroutine forday(dtime,yrflag, iyear,iday,ihour)
      use mod_xc  ! HYCOM communication interface
      implicit none
!
      real*8  dtime
      integer yrflag, iyear,iday,ihour
!
! --- converts model day to "calendar" date (year,ordinal-day,hour).
!
      real*8  dtim1,day
      integer iyr,nleap
!
      if     (yrflag.eq.0) then
! ---   360 days per model year, starting Jan 16
        iyear =  int((dtime+15.001d0)/360.d0) + 1
        iday  =  mod( dtime+15.001d0 ,360.d0) + 1
        ihour = (mod( dtime+15.001d0 ,360.d0) + 1.d0 - iday)*24.d0
!
      elseif (yrflag.eq.1) then
! ---   366 days per model year, starting Jan 16
        iyear =  int((dtime+15.001d0)/366.d0) + 1
        iday  =  mod( dtime+15.001d0 ,366.d0) + 1
        ihour = (mod( dtime+15.001d0 ,366.d0) + 1.d0 - iday)*24.d0
!
      elseif (yrflag.eq.2) then
! ---   366 days per model year, starting Jan 01
        iyear =  int((dtime+ 0.001d0)/366.d0) + 1
        iday  =  mod( dtime+ 0.001d0 ,366.d0) + 1
        ihour = (mod( dtime+ 0.001d0 ,366.d0) + 1.d0 - iday)*24.d0
!
      elseif (yrflag.eq.3) then
! ---   model day is calendar days since 01/01/1901
        iyr   = (dtime-1.d0)/365.25d0
        nleap = iyr/4
        dtim1 = 365.d0*iyr + nleap + 1.d0
        day   = dtime - dtim1 + 1.d0
        if     (dtim1.gt.dtime) then
          iyr = iyr - 1
        elseif (day.ge.367.d0) then
          iyr = iyr + 1
        elseif (day.ge.366.d0 .and. mod(iyr,4).ne.3) then
          iyr = iyr + 1
        endif
        nleap = iyr/4
        dtim1 = 365.d0*iyr + nleap + 1.d0
!
        iyear =  1901 + iyr
        iday  =  dtime - dtim1 + 1.001d0
        ihour = (dtime - dtim1 + 1.001d0 - iday)*24.d0
!
      elseif (yrflag.eq.4) then
! ---   365 days per model year, starting Jan 01 -No Leap year-
        iyear =  int((dtime+ 0.001d0)/365.d0) + 1
        iday  =  mod( dtime+ 0.001d0 ,365.d0) + 1
        ihour = (mod( dtime+ 0.001d0 ,365.d0) + 1.d0 - iday)*24.d0
!
      else
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in forday - unsupported yrflag value'
        write(lp,*)
        endif !1st tile
        call xcstop('(forday)')
               stop '(forday)'
      endif
      return
      end
!
!
      subroutine fordate(dtime,yrflag, iyear,month,iday,ihour)
      implicit none
!
      real*8  dtime
      integer yrflag, iyear,month,iday,ihour
!
! --- converts model day to "calendar" date (year,month,day,hour).
!
      integer jday,k,m
!
      integer month0(13,3)
      data    month0 / 1,  31,  61,  91, 121, 151, 181, &
                          211, 241, 271, 301, 331, 361, &
                       1,  32,  60,  91, 121, 152, 182, &
                          213, 244, 274, 305, 335, 366, &
                       1,  32,  61,  92, 122, 153, 183, &
                          214, 245, 275, 306, 336, 367 /
!
      call forday(dtime,yrflag, iyear,jday,ihour)
!
      if (yrflag.eq.3) then
        if     (mod(iyear,4).eq.0) then
          k = 3  !leap year
        else
          k = 2  !standard year
        endif
      elseif (yrflag.eq.4) then
        k = 2  !365-day year
      elseif (yrflag.eq.0) then
        k = 1  !360-day year
      else
        k = 3  !leap year
      endif
      do m= 1,12
        if     (jday.ge.month0(m,  k) .and. &
                jday.lt.month0(m+1,k)      ) then
          month = m
          iday  = jday - month0(m,k) + 1
        endif
      enddo
      return
      end
!
!
      subroutine forfuna
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize input of atmospheric forcing fields
!
! --- units of tau_x  are N/m^2  (positive eastwards  w.r.t. the grid)
! --- units of tau_y  are N/m^2  (positive northwards w.r.t. the grid)
! --- units of wnd_x  are m/s    (positive eastwards  w.r.t. the grid)
! --- units of wnd_y  are m/s    (positive northwards w.r.t. the grid)
! --- units of wndspd are m/s
! --- units of ustar  are m/s
! --- units of airtmp are degC
! --- units of surtmp are degC
! --- units of seatmp are degC
! --- units of vapmix are kg/k 
! --- units of spchum are kg/k
! --- units of mslprs are Pa     (anomaly, offset from total by prsbas)
! --- units of precip are m/s    (positive into ocean)
! --- units of radflx are w/m^2  (positive into ocean)
! --- units of swflx  are w/m^2  (positive into ocean)
! --- units of offlux are w/m^2  (positive into ocean)
! --- units of oftaux are N/m^2  (positive eastwards  w.r.t. the grid)
! --- units of oftauy are N/m^2  (positive northwards w.r.t. the grid)
!
! --- tau_x and tau_y are either on u&v grids or both on the p grid,
! --- depending on the value of blkdat input parameter "wndflg".
! --- in any case, they are always oriented along the local grid
! --- which need not be east-west and north-south.
! --- all other fields, including wnd* and ustar, are always on the p grid.
!
! --- I/O and array I/O units 899-910 are reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer   lgth,mo
      character preambl(5)*79
      real      pcmax,one,oneps
      integer   i,j
!
      mreca=1
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening forcing fields ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      if     (mslprf .or. flxflg.eq.6) then
      call zaiopf(flnmfor(1:lgth)//'forcing.mslprs.a', 'old', 899)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+899,file=flnmfor(1:lgth)//'forcing.mslprs.b', &
         status='old', action='read')
      read (uoff+899,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 899)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1., &
!diag      'mslprs (Pa)')
      endif !mslprf
!
      if (windf) then
!
      if     (wndflg.lt.4) then
        call zaiopf(flnmfor(1:lgth)//'forcing.tauewd.a', 'old', 901)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+901,file=flnmfor(1:lgth)//'forcing.tauewd.b', &
              status='old', action='read')
        read (uoff+901,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 901)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1000., &
!diag        'tau_x (x 1000 N/m^2 ) ')
!
        call zaiopf(flnmfor(1:lgth)//'forcing.taunwd.a', 'old', 902)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+902,file=flnmfor(1:lgth)//'forcing.taunwd.b', &
              status='old', action='read')
        read (uoff+902,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 902)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1000., &
!diag        'tau_y (x 1000 N/m^2 ) ')
      else !read 10m wind components
        call zaiopf(flnmfor(1:lgth)//'forcing.wndewd.a', 'old', 901)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+901,file=flnmfor(1:lgth)//'forcing.wndewd.b', &
              status='old', action='read')
        read (uoff+901,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 901)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1., &
!diag        'wnd_x (m/s ) ')
!
        call zaiopf(flnmfor(1:lgth)//'forcing.wndnwd.a', 'old', 902)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+902,file=flnmfor(1:lgth)//'forcing.wndnwd.b', &
              status='old', action='read')
        read (uoff+902,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 902)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1., &
!diag        'wnd_y (m/s) ')
        endif !stress:wind
!
        if     (stroff) then
          call zaiopf(flnmfor(1:lgth)//'forcing.ofstrs.a', 'old', 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+916,file=flnmfor(1:lgth)//'forcing.ofstrs.b', &
             status='old', action='read')
          read (uoff+916,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
          call rdmonth(oftaux, 916)
          call rdmonth(oftauy, 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close( unit=uoff+916)
          endif
          call zaiocl(916)
!diag     call prtmsk(ip,oftaux,util2,idm,idm,jdm,  0.,1.0, &
!diag          'taux offset (N)  ')
!diag     call prtmsk(ip,oftauy,util2,idm,idm,jdm,  0.,1.0, &
!diag          'tauy offset (N)  ')
          else
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
              oftaux(i,j) = 0.0
              oftauy(i,j) = 0.0
            enddo !i
          enddo !j
        endif !stroff:else
!
      endif				!  windf = .true.
!
      if (thermo) then
!
      if     (ustflg.eq.3) then
      call zaiopf(flnmfor(1:lgth)//'forcing.ustar.a', 'old', 900)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+900,file=flnmfor(1:lgth)//'forcing.ustar.b', &
         status='old', action='read')
      read (uoff+900,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 900)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1000., &
!diag      'ustar (x 1000 ???)')
      endif !ustflg.eq.3
!
      if (wndflg.lt.3) then
      call zaiopf(flnmfor(1:lgth)//'forcing.wndspd.a', 'old', 903)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+903,file=flnmfor(1:lgth)//'forcing.wndspd.b', &
         status='old', action='read')
      read (uoff+903,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 903)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,10., &
!diag      'wind speed  (x 10 m/s)')
      endif !wndflg.lt.3
!
      if     (flxflg.ne.3) then
      call zaiopf(flnmfor(1:lgth)//'forcing.airtmp.a', 'old', 904)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+904,file=flnmfor(1:lgth)//'forcing.airtmp.b', &
         status='old', action='read')
      read (uoff+904,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 904)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,10., &
!diag      'air temperature  (0.1 c)')
!
      if     (flxflg.ne.5) then
        call zaiopf(flnmfor(1:lgth)//'forcing.vapmix.a', 'old', 905)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+905,file=flnmfor(1:lgth)//'forcing.vapmix.b', &
           status='old', action='read')
        read (uoff+905,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 905)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,10000., &
!diag        'mixing ratio  (0.1 g/kg)')
      else !flgflg==5
        call zaiopf(flnmfor(1:lgth)//'forcing.spchum.a', 'old', 905)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+905,file=flnmfor(1:lgth)//'forcing.spchum.b', &
           status='old', action='read')
        read (uoff+905,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 905)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,10000., &
!diag        'specific humidity  (0.1 g/kg)')
      endif !flxflg 4,5
      endif !flxflg.ne.3
!
      if     (pcipf) then
        call zaiopf(flnmfor(1:lgth)//'forcing.precip.a', 'old', 906)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+906,file=flnmfor(1:lgth)//'forcing.precip.b', &
           status='old', action='read')
        read (uoff+906,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        pcmax = -huge(pcmax)
        call rdmonth(util1, 906)
        do j=1,jj
          do i=1,ii
            pcmax = max(pcmax,util1(i,j))
          enddo
        enddo
        call xcmaxr(pcmax)
!diag   call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,86400.*36000., &
!diag        'precipitation (cm/year) ')
!
! ---   zero fields implies no surface salinity flux
        pcipf = pcmax.ne.0.0
        if     (.not.pcipf) then
          if     (mnproc.eq.1) then
          write (lp,*)
          write (lp,*) '***** no surface salinity flux *****'
          write (lp,*)
          endif !1st tile
          call xcsync(flush_lp)
        endif  !pcipf actually .false.
      endif  !pcipf initially .true.
!
      call zaiopf(flnmfor(1:lgth)//'forcing.radflx.a', 'old', 907)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+907,file=flnmfor(1:lgth)//'forcing.radflx.b', &
         status='old', action='read')
      read (uoff+907,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 907)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1.0, &
!diag      'net radiation (w/m^2 )  ')
!
      call zaiopf(flnmfor(1:lgth)//'forcing.shwflx.a', 'old', 908)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+908,file=flnmfor(1:lgth)//'forcing.shwflx.b', &
         status='old', action='read')
      read (uoff+908,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 908)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1.0, &
!diag      'sw radiation (w/m^2 )  ')
!
      endif                    !  thermo = .true.
!
      if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
              icmflg.eq.2 .or. ticegr.eq.0.0     ) then
      call zaiopf(flnmfor(1:lgth)//'forcing.surtmp.a', 'old', 909)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+909,file=flnmfor(1:lgth)//'forcing.surtmp.b', &
         status='old', action='read')
      read (uoff+909,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 909)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1.0, &
!diag      'SST from lw rad (degC) ')
      endif !surtmp
!
      if     (sstflg.eq.3) then
      call zaiopf(flnmfor(1:lgth)//'forcing.seatmp.a', 'old', 910)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+910,file=flnmfor(1:lgth)//'forcing.seatmp.b', &
         status='old', action='read')
      read (uoff+910,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      call rdmonth(util1, 910)
!diag call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,1.0, &
!diag      'SST from obs. (degC)   ')
      endif !sstflg.eq.3
!
      if     (flxoff) then
        call zaiopf(flnmfor(1:lgth)//'forcing.offlux.a', 'old', 916)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+916,file=flnmfor(1:lgth)//'forcing.offlux.b', &
           status='old', action='read')
        read (uoff+916,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(offlux, 916)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close( unit=uoff+916)
        endif
        call zaiocl(916)
!diag   call prtmsk(ip,offlux,util2,idm,idm,jdm,  0.,1.0, &
!diag        'heat flux offset (w/m^2 )  ')
        else
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-nbdy,jj+nbdy
          do i=1-nbdy,ii+nbdy
            offlux(i,j) = 0.0
          enddo !i
        enddo !j
      endif !flxoff:else
!
! --- jerlov used in call to swfrac_ij
      if     (jerlv0.gt.0) then
! ---   calculate jerlov water type,
! ---   which governs the penetration depth of shortwave radiation.
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-nbdy,jj+nbdy
          do i=1-nbdy,ii+nbdy
! ---       map shallow depths to high jerlov numbers
            jerlov(i,j)=6-max(1,min(5,int(depths(i,j)/15.0)))
            jerlov(i,j)=max(jerlv0,jerlov(i,j))
          enddo
        enddo
      else
! ---   jerlv0= 0 uses an input annual/monthly kpar field
! ---   jerlv0=-1 uses an input annual/monthly chl  field
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-nbdy,jj+nbdy
          do i=1-nbdy,ii+nbdy
            jerlov(i,j)=jerlv0
          enddo
        enddo
      endif
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening forcing fields'
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfund(tiddrg)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer tiddrg
!
! --- initialize tidal drag tensor
!
! --- units of dragrh and drgten are m/s on the p-grid.
! ---          dragrh and drgten should be positive.
!
! --- I/O and array I/O unit 925 used here, but not reserved.
!
! --- all input fields must be defined at all grid points
!
      integer   i,j,k,l,lgth
!
      if     (mnproc.eq.1) then
      if     (tiddrg.eq.1) then
        write (lp,*) ' now opening tidal drag rh field  ...'
      else
        write (lp,*) ' now opening tidal drag tensor fields  ...'
      endif !tiddrg
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      if     (tiddrg.eq.1) then
! ---   original filename, 1 field (rh)
        call zaiopf(flnmfor(1:lgth)//'tidal.rh.a', 'old', 925)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+925,file=flnmfor(1:lgth)//'tidal.rh.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util1, 925)
        call xctilr( util1,1,1, nbdy,nbdy, halo_ps)
! ---   cast scalar to tensor drag.
        drgten(1,1,:,:) = drgscl*util1(:,:) !uu
        drgten(1,2,:,:) = 0.0  !uv
        drgten(2,1,:,:) = 0.0  !vu
        drgten(2,2,:,:) = drgscl*util1(:,:) !vv
      else
! ---   drag tensor filename, 4 fields uu, uv, vv, vu
        call zaiopf(flnmfor(1:lgth)//'tidal.tensor.a', 'old', 925)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+925,file=flnmfor(1:lgth)//'tidal.tensor.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util1, 925)
        call xctilr( util1,1,1, nbdy,nbdy, halo_ps)
        drgten(1,1,:,:) = drgscl*util1(:,:) !uu
        call rdmonth(util1, 925)
        call xctilr( util1,1,1, nbdy,nbdy, halo_ps)
        drgten(1,2,:,:) = drgscl*util1(:,:) !uv
        call rdmonth(util1, 925)
        call xctilr( util1,1,1, nbdy,nbdy, halo_ps)
        drgten(2,1,:,:) = drgscl*util1(:,:) !vu
        call rdmonth(util1, 925)
        call xctilr( util1,1,1, nbdy,nbdy, halo_ps)
        drgten(2,2,:,:) = drgscl*util1(:,:) !vv
      endif !tiddrg
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close (unit=uoff+925)
      endif
      call zaiocl(925)
!
      if     (mnproc.eq.1) then
      if     (tiddrg.eq.1) then
        write (lp,*) ' ...finished opening tidal drag rh field'
      else
        write (lp,*) ' ...finished opening tidal drag tensor fields'
      endif !tiddrg
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunh(dtime)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
!
! --- high frequency atmospheric forcing field processing.
!
! --- units of tau_x  are N/m^2  (positive eastwards  w.r.t. the grid)
! --- units of tau_y  are N/m^2  (positive northwards w.r.t. the grid)
! --- units of wnd_x  are m/s    (positive eastwards  w.r.t. the grid)
! --- units of wnd_y  are m/s    (positive northwards w.r.t. the grid)
! --- units of wndspd are m/s
! --- units of ustar  are m/s
! --- units of airtmp are degC
! --- units of surtmp are degC
! --- units of seatmp are degC
! --- units of stoc_t are degC/hr
! --- units of stoc_s are psu/hr
! --- units of stoc_u are m/s/hr on the u-grid
! --- units of stoc_v are m/s/hr on the v-grid
! --- units of vapmix are kg/kg
! --- units of spchum are kg/k
! --- units of mslprs are Pa     (anomaly, offset from total by prsbas)
! --- units of precip are m/s    (positive into ocean)
! --- units of radflx are w/m^2  (positive into ocean)
! --- units of swflx  are w/m^2  (positive into ocean)
! --- units of offlux are w/m^2  (positive into ocean)
!
! --- tau_x and tau_y are either on u&v grids or both on the p grid,
! --- depending on the value of blkdat input parameter "wndflg".
! --- in any case, they are always oriented along the local grid
! --- which need not be east-west and north-south.
! --- all other fields, including wnd* and ustar, are always on the p grid.
!
! --- I/O and array I/O units 895-910 are reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      real*8    dtime0,dtime1
      save      dtime0,dtime1
!
      character preambl(5)*79,cline*80
      real      pcmax
      integer   i,ios,iunit,j,lgth,nrec
! ESPC --- add
#if ! defined(ESPC_ATM) && ! defined(ESPC_NAVGEM) && ! defined(ESPC_DATA_ATM)

!
! --- w0 negative on first call only.
      if     (w0.lt.-1.0) then
!
! ---   initialize forcing fields
!
        if      (.not.windf) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in forfunh - windf must be .true.'
          write(lp,*)
          endif !1st tile
          call xcstop('(forfunh)')
                 stop '(forfunh)'
        elseif (.not.thermo) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in forfunh - thermo must be .true.'
          write(lp,*)
          endif !1st tile
          call xcstop('(forfunh)')
                 stop '(forfunh)'
        endif
!
        if     (natm.eq.4) then
! ---     linear interpolation in time, so slots 3 and 4 are zero.
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
                taux(i,j,natm-1) = 0.0
                taux(i,j,natm)   = 0.0
                tauy(i,j,natm-1) = 0.0
                tauy(i,j,natm)   = 0.0
              wndspd(i,j,natm-1) = 0.0
              wndspd(i,j,natm)   = 0.0
              ustara(i,j,natm-1) = 0.0
              ustara(i,j,natm)   = 0.0
              airtmp(i,j,natm-1) = 0.0
              airtmp(i,j,natm)   = 0.0
              vapmix(i,j,natm-1) = 0.0
              vapmix(i,j,natm)   = 0.0
              mslprs(i,j,natm-1) = 0.0
              mslprs(i,j,natm)   = 0.0
              precip(i,j,natm-1) = 0.0
              precip(i,j,natm)   = 0.0
              radflx(i,j,natm-1) = 0.0
              radflx(i,j,natm)   = 0.0
               swflx(i,j,natm-1) = 0.0
               swflx(i,j,natm)   = 0.0
              surtmp(i,j,natm-1) = 0.0
              surtmp(i,j,natm)   = 0.0
              seatmp(i,j,natm-1) = 0.0
              seatmp(i,j,natm)   = 0.0
              stoc_t(i,j,natm-1) = 0.0
              stoc_t(i,j,natm)   = 0.0
              stoc_s(i,j,natm-1) = 0.0
              stoc_s(i,j,natm)   = 0.0
              stoc_u(i,j,natm-1) = 0.0
              stoc_u(i,j,natm)   = 0.0
              stoc_v(i,j,natm-1) = 0.0
              stoc_v(i,j,natm)   = 0.0
            enddo
          enddo
        endif !natm.eq.4
!
! ---   open all forcing files.
        if     (mnproc.eq.1) then
        write (lp,*) ' now initializing forcing fields ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        if     (stfflg.gt.0) then
        call zaiopf(flnmfor(1:lgth)//'forcing.stoc_t.a', 'old', 896)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+896,file=flnmfor(1:lgth)//'forcing.stoc_t.b', &
           status='old', action='read')
        read (uoff+896,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call zaiopf(flnmfor(1:lgth)//'forcing.stoc_s.a', 'old', 897)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+897,file=flnmfor(1:lgth)//'forcing.stoc_s.b', &
           status='old', action='read')
        read (uoff+897,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !stfflg>0
!
        if     (stfflg.eq.2) then
        call zaiopf(flnmfor(1:lgth)//'forcing.stoc_u.a', 'old', 895)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+895,file=flnmfor(1:lgth)//'forcing.stoc_u.b', &
           status='old', action='read')
        read (uoff+895,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call zaiopf(flnmfor(1:lgth)//'forcing.stoc_v.a', 'old', 898)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+898,file=flnmfor(1:lgth)//'forcing.stoc_v.b', &
           status='old', action='read')
        read (uoff+898,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !stfflg==2
!
        if     (mslprf .or. flxflg.eq.6) then
        call zaiopf(flnmfor(1:lgth)//'forcing.mslprs.a', 'old', 899)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+899,file=flnmfor(1:lgth)//'forcing.mslprs.b', &
           status='old', action='read')
        read (uoff+899,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !mslprf
!
        if     (wndflg.lt.4) then
          call zaiopf(flnmfor(1:lgth)//'forcing.tauewd.a', 'old', 901)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+901,file=flnmfor(1:lgth)//'forcing.tauewd.b', &
                status='old', action='read')
          read (uoff+901,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
!
          call zaiopf(flnmfor(1:lgth)//'forcing.taunwd.a', 'old', 902)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+902,file=flnmfor(1:lgth)//'forcing.taunwd.b', &
             status='old', action='read')
          read (uoff+902,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        else !read 10m wind components
          call zaiopf(flnmfor(1:lgth)//'forcing.wndewd.a', 'old', 901)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+901,file=flnmfor(1:lgth)//'forcing.wndewd.b', &
                status='old', action='read')
          read (uoff+901,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
!             
          call zaiopf(flnmfor(1:lgth)//'forcing.wndnwd.a', 'old', 902)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+902,file=flnmfor(1:lgth)//'forcing.wndnwd.b', &
                status='old', action='read')
          read (uoff+902,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        endif !stress:wind
!
        if     (stroff) then
          call zaiopf(flnmfor(1:lgth)//'forcing.ofstrs.a', 'old', 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+916,file=flnmfor(1:lgth)//'forcing.ofstrs.b', &
             status='old', action='read')
          read (uoff+916,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
          call rdmonth(oftaux, 916)
          call rdmonth(oftauy, 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close( unit=uoff+916)
          endif
          call zaiocl(916)
!diag     call prtmsk(ip,oftaux,util2,idm,idm,jdm,  0.,1.0, &
!diag          'taux offset (N)  ')
!diag     call prtmsk(ip,oftauy,util2,idm,idm,jdm,  0.,1.0, &
!diag          'tauy offset (N)  ')
          else
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
              oftaux(i,j) = 0.0
              oftauy(i,j) = 0.0
            enddo !i
          enddo !j
        endif !stroff:else
!
        if     (ustflg.eq.3) then
        call zaiopf(flnmfor(1:lgth)//'forcing.ustar.a', 'old', 900)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+900,file=flnmfor(1:lgth)//'forcing.ustar.b', &
           status='old', action='read')
        read (uoff+900,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !ustflg.eq.3
!
        if (wndflg.lt.3) then
        call zaiopf(flnmfor(1:lgth)//'forcing.wndspd.a', 'old', 903)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+903,file=flnmfor(1:lgth)//'forcing.wndspd.b', &
           status='old', action='read')
        read (uoff+903,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !wndflg.lt.3
!
        if     (flxflg.ne.3) then
        call zaiopf(flnmfor(1:lgth)//'forcing.airtmp.a', 'old', 904)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+904,file=flnmfor(1:lgth)//'forcing.airtmp.b', &
           status='old', action='read')
        read (uoff+904,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
!
        if     (flxflg.ne.5) then
          call zaiopf(flnmfor(1:lgth)//'forcing.vapmix.a', 'old', 905)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+905,file=flnmfor(1:lgth)//'forcing.vapmix.b', &
             status='old', action='read')
          read (uoff+905,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        else !flgflg==5
          call zaiopf(flnmfor(1:lgth)//'forcing.spchum.a', 'old', 905)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+905,file=flnmfor(1:lgth)//'forcing.spchum.b', &
             status='old', action='read')
          read (uoff+905,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        endif !flxflg 4,5
        endif !flxflg.ne.3
!
        if     (pcipf) then
        call zaiopf(flnmfor(1:lgth)//'forcing.precip.a', 'old', 906)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+906,file=flnmfor(1:lgth)//'forcing.precip.b', &
           status='old', action='read')
        read (uoff+906,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        endif !pcipf
!
        call zaiopf(flnmfor(1:lgth)//'forcing.radflx.a', 'old', 907)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+907,file=flnmfor(1:lgth)//'forcing.radflx.b', &
           status='old', action='read')
        read (uoff+907,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
!
        call zaiopf(flnmfor(1:lgth)//'forcing.shwflx.a', 'old', 908)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+908,file=flnmfor(1:lgth)//'forcing.shwflx.b', &
           status='old', action='read')
        read (uoff+908,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
!
        if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
                icmflg.eq.2 .or. ticegr.eq.0.0     ) then
          call zaiopf(flnmfor(1:lgth)//'forcing.surtmp.a', 'old', 909)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+909,file=flnmfor(1:lgth)//'forcing.surtmp.b', &
             status='old', action='read')
          read (uoff+909,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        endif !surtmp
!
        if     (sstflg.eq.3) then
          call zaiopf(flnmfor(1:lgth)//'forcing.seatmp.a', 'old', 910)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+910,file=flnmfor(1:lgth)//'forcing.seatmp.b', &
             status='old', action='read')
          read (uoff+910,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
        endif
!
! ESPC --- add
# endif

        if     (flxoff) then
          call zaiopf(flnmfor(1:lgth)//'forcing.offlux.a', 'old', 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          open (unit=uoff+916,file=flnmfor(1:lgth)//'forcing.offlux.b', &
             status='old', action='read')
          read (uoff+916,'(a79)') preambl
          endif !1st tile
          call preambl_print(preambl)
          call rdmonth(offlux, 916)
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close( unit=uoff+916)
          endif
          call zaiocl(916)
        else
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
              offlux(i,j) = 0.0
            enddo !i
          enddo !j
        endif !flxoff:else

! ESPC ---add
#if ! defined(ESPC_ATM) && ! defined(ESPC_NAVGEM) && ! defined(ESPC_DATA_ATM)

!
! ---   skip ahead to the start time.
        nrec   = 0
        dtime1 = huge(dtime1)
        do  ! infinate loop, with exit at end
          dtime0 = dtime1
          nrec   = nrec + 1
          call zagetc(cline,ios, uoff+901)
          if     (ios.ne.0) then
            if     (mnproc.eq.1) then
              write(lp,*)
              write(lp,*) 'error in forfunh - hit end of input'
              write(lp,*) 'dtime0,dtime1 = ',dtime0,dtime1
              write(lp,*) 'dtime = ',dtime
              write(lp,*)
            endif !1st tile
            call xcstop('(forfunh)')
                   stop '(forfunh)'
          endif
          i = index(cline,'=')
          read (cline(i+1:),*) dtime1
          if     (yrflag.eq.2) then
            if     (nrec.eq.1 .and. abs(dtime1-1096.0d0).gt.0.01) then
!
! ---         climatology must start on wind day 1096.0, 01/01/1904.
              if     (mnproc.eq.1) then
              write(lp,'(a)')  cline
              write(lp,'(/ a,a / a,g15.6 /)') &
                'error in forfunh - forcing climatology', &
                ' must start on wind day 1096', &
                'dtime1 = ',dtime1
              endif !1st tile
              call xcstop('(forfunh)')
                     stop '(forfunh)'
            endif
            dtime1 = (dtime1 - 1096.0d0) +  &
                     wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
            if     (nrec.ne.1 .and. dtime1.lt.dtime0) then
              dtime1 = dtime1 + wndrep
            endif
          elseif (yrflag.eq.4) then 
            if     (nrec.eq.1 .and. abs(dtime1-731.0d0).gt.0.01) then
!
! ---         climatology must start on wind day 731.0, 01/01/1903.
              if     (mnproc.eq.1) then
              write(lp,'(a)')  cline
              write(lp,'(/ a,a / a,g15.6 /)') &
                'error in forfunh - forcing climatology', &
                ' must start on wind day 731', &
                'dtime1 = ',dtime1
              endif !1st tile
              call xcstop('(forfunh)')
                     stop '(forfunh)'
            endif
            dtime1 = (dtime1 - 731.0d0) +  &
                     wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
            if     (nrec.ne.1 .and. dtime1.lt.dtime0) then
              dtime1 = dtime1 + wndrep
            endif
          elseif (nrec.eq.1 .and. dtime1.lt.1462.0d0) then
!
! ---       otherwise, must start after wind day 1462.0, 01/01/1905.
            if     (mnproc.eq.1) then
            write(lp,'(a)')  cline
            write(lp,'(/ a,a / a,g15.6 /)') &
              'error in forfunh - actual forcing', &
              ' must start after wind day 1462', &
              'dtime1 = ',dtime1
            endif !1st tile
            call xcstop('(forfunh)')
                   stop '(forfunh)'
          endif
          if     (dtime0.le.dtime .and. dtime1.gt.dtime) then
            exit
          endif
        enddo   ! infinate loop, with exit above
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          rewind(unit=uoff+901)
          read (uoff+901,'(a79)') preambl
        endif
!
        do iunit= 895,908
          if     (iunit.eq.896 .and. stfflg.eq.0) then
            cycle  !no stoc_t
          endif
          if     (iunit.eq.897 .and. stfflg.eq.0) then
            cycle  !no stoc_s
          endif
          if     (iunit.eq.895 .and. stfflg.lt.2) then
            cycle  !no stoc_u
          endif
          if     (iunit.eq.898 .and. stfflg.lt.2) then
            cycle  !no stoc_v
          endif
          if     (iunit.eq.899 .and.  &
                  .not.(mslprf .or. flxflg.eq.6)) then
            cycle  !no mslprs
          endif
          if     (iunit.eq.900 .and. ustflg.ne.3) then
            cycle  !no ustar 
          endif
          if     (iunit.eq.903 .and. wndflg.ge.3) then
            cycle  !no wndspd
          endif
          do i= 1,nrec-2
            call skmonth(iunit)
          enddo
        enddo
        if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
                icmflg.eq.2 .or. ticegr.eq.0.0     ) then
          do i= 1,nrec-2
            call skmonth(909)
          enddo
        endif !surtmp
        if     (sstflg.eq.3) then
          do i= 1,nrec-2
            call skmonth(910)
          enddo
        endif
        dtime1 = huge(dtime1)
        call rdpall(dtime0,dtime1)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)
        endif
        call rdpall(dtime0,dtime1)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        endif
!
! ---   zero precip field implies no surface salinity flux
        if     (pcipf) then
          pcmax = -huge(pcmax)
          do j=1,jj
            do i=1,ii
              pcmax = max(pcmax,precip(i,j,1),precip(i,j,2))
            enddo
          enddo
          call xcmaxr(pcmax)
          pcipf = pcmax.ne.0.0
          if     (.not.pcipf) then
            if     (mnproc.eq.1) then
            write (lp,*)
            write (lp,*) '***** no surface salinity flux *****'
            write (lp,*)
            endif !1st tile
            call xcsync(flush_lp)
          endif  !pcipf actually .false.
        endif  !pcipf initially .true.
!
! ---   jerlov used in call to swfrac_ij
        if     (jerlv0.gt.0) then
! ---     calculate jerlov water type,
! ---     which governs the penetration depth of shortwave radiation.
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
! ---         map shallow depths to high jerlov numbers
              jerlov(i,j)=6-max(1,min(5,int(depths(i,j)/15.0)))
              jerlov(i,j)=max(jerlv0,jerlov(i,j))
            enddo
          enddo
        else
! ---     jerlv0= 0 uses an input annual/monthly kpar field
! ---     jerlv0=-1 uses an input annual/monthly chl  field
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
              jerlov(i,j)=jerlv0
            enddo
          enddo
        endif
        if     (mnproc.eq.1) then
        write (lp,*) 
        write (lp,*) ' dtime,dtime0,dtime1 = ',dtime,dtime0,dtime1
        write (lp,*) 
        write (lp,*) ' ...finished initializing forcing fields'
        endif !1st tile
        call xcsync(flush_lp)
      endif  ! initialization
!
      if     (dtime.gt.dtime1) then
!
! ---   get the next set of fields.
!           if     (mnproc.eq.1) then
!           write(lp,*) 'enter rdpall - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
        call rdpall(dtime0,dtime1)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        endif
!           if     (mnproc.eq.1) then
!           write(lp,*) ' exit rdpall - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
      endif
!
! --- linear interpolation in time.
      w0 = (dtime1-dtime)/(dtime1-dtime0)
      w1 = 1.0 - w0
!           if     (mnproc.eq.1) then
!           write(lp,*) 'rdpall - dtime,w0,w1 = ',dtime,w0,w1
!           endif !1st tile
!           call xcsync(flush_lp)
      return
! ESPC --- add
# endif
      end
!
!
      subroutine forfunhp(dtime)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
!
! --- high frequency atmospheric pressure forcing field processing.
! --- call either forfunh or forfunhp, not both.
!
! --- units of mslprs are Pa (anomaly, offset from total by prsbas)
!
! --- mslprs is always on the p grid.
! --- mslprs must be defined at all grid points
!
! --- I/O and array I/O unit 899 is reserved for the entire run.
!
!
      real*8    dtime0,dtime1
      save      dtime0,dtime1
!
      character preambl(5)*79,cline*80
      integer   i,ios,iunit,j,lgth,nrec
!
! --- w0 negative on first call only.
      if     (w0.lt.-1.0) then
!
! ---   initialize forcing fields
!
        if      (.not.mslprf) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in forfunhp - mslprf must be .true.'
          write(lp,*)
          endif !1st tile
          call xcstop('(forfunhp)')
                 stop '(forfunhp)'
        elseif  (windf) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in forfunhp - windf must be .false.'
          write(lp,*)
          endif !1st tile
          call xcstop('(forfunhp)')
                 stop '(forfunhp)'
        elseif (thermo) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in forfunhp - thermo must be .false.'
          write(lp,*)
          endif !1st tile
          call xcstop('(forfunhp)')
                 stop '(forfunhp)'
        endif
!
        if     (natm.eq.4) then
! ---     linear interpolation in time, so slots 3 and 4 are zero.
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j=1-nbdy,jj+nbdy
            do i=1-nbdy,ii+nbdy
              mslprs(i,j,natm-1) = 0.0
              mslprs(i,j,natm)   = 0.0
            enddo
          enddo
        endif !natm.eq.4
!
! ---   open pressure forcing file.
        if     (mnproc.eq.1) then
        write (lp,*) ' now initializing forcing fields ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'forcing.mslprs.a', 'old', 899)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+899,file=flnmfor(1:lgth)//'forcing.mslprs.b', &
           status='old', action='read')
        read (uoff+899,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
!
! ---   skip ahead to the start time.
        nrec   = 0
        dtime1 = huge(dtime1)
        do  ! infinate loop, with exit at end
          dtime0 = dtime1
          nrec   = nrec + 1
          call zagetc(cline,ios, uoff+899)
          if     (ios.ne.0) then
            if     (mnproc.eq.1) then
              write(lp,*)
              write(lp,*) 'error in forfunhp - hit end of input'
              write(lp,*) 'dtime0,dtime1 = ',dtime0,dtime1
              write(lp,*) 'dtime = ',dtime
              write(lp,*)
            endif !1st tile
            call xcstop('(forfunhp)')
                   stop '(forfunhp)'
          endif
          i = index(cline,'=')
          read (cline(i+1:),*) dtime1
          if     (yrflag.eq.2) then
            if     (nrec.eq.1 .and. abs(dtime1-1096.0d0).gt.0.01) then
!
! ---         climatology must start on wind day 1096.0, 01/01/1904.
              if     (mnproc.eq.1) then
              write(lp,'(a)')  cline
              write(lp,'(/ a,a / a,g15.6 /)') &
                'error in forfunhp - forcing climatology', &
                ' must start on wind day 1096', &
                'dtime1 = ',dtime1
              endif !1st tile
              call xcstop('(forfunhp)')
                     stop '(forfunhp)'
            endif
            dtime1 = (dtime1 - 1096.0d0) +  &
                     wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
            if     (nrec.ne.1 .and. dtime1.lt.dtime0) then
              dtime1 = dtime1 + wndrep
            endif
          elseif (yrflag.eq.4) then
            if     (nrec.eq.1 .and. abs(dtime1-731.0d0).gt.0.01) then
!
! ---         climatology must start on wind day 731.0, 01/01/1903.
              if     (mnproc.eq.1) then
              write(lp,'(a)')  cline
              write(lp,'(/ a,a / a,g15.6 /)') &
                'error in forfunhp - forcing climatology', &
                ' must start on wind day 731', &
                'dtime1 = ',dtime1
              endif !1st tile
              call xcstop('(forfunhp)')
                     stop '(forfunhp)'
            endif
            dtime1 = (dtime1 - 731.0d0) +  &
                     wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
            if     (nrec.ne.1 .and. dtime1.lt.dtime0) then
              dtime1 = dtime1 + wndrep
            endif
          elseif (nrec.eq.1 .and. dtime1.lt.1462.0d0) then
!
! ---       otherwise, must start after wind day 1462.0, 01/01/1905.
            if     (mnproc.eq.1) then
            write(lp,'(a)')  cline
            write(lp,'(/ a,a / a,g15.6 /)') &
              'error in forfunhp - actual forcing', &
              ' must start after wind day 1462', &
              'dtime1 = ',dtime1
            endif !1st tile
            call xcstop('(forfunhp)')
                   stop '(forfunhp)'
          endif
          if     (dtime0.le.dtime .and. dtime1.gt.dtime) then
            exit
          endif
        enddo   ! infinate loop, with exit above
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          rewind(unit=uoff+899)
          read (uoff+899,'(a79)') preambl
        endif
!
        do i= 1,nrec-2
          call skmonth(899)
        enddo
        dtime0 = huge(dtime1)
        call rdpall1(mslprs,dtime1,899,.true.)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)
        endif
        dtime0 = dtime1
        call rdpall1(mslprs,dtime1,899,.true.)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        endif
        if     (mnproc.eq.1) then
        write (lp,*) 
        write (lp,*) ' dtime,dtime0,dtime1 = ',dtime,dtime0,dtime1
        write (lp,*) 
        write (lp,*) ' ...finished initializing forcing fields'
        endif !1st tile
        call xcsync(flush_lp)
      endif  ! initialization
!
      if     (dtime.gt.dtime1) then
!
! ---   get the next set of fields.
!           if     (mnproc.eq.1) then
!           write(lp,*) 'enter rdpall - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
        dtime0 = dtime1
        call rdpall1(mslprs,dtime1,899,.true.)
        if     (yrflag.eq.2) then
          dtime1 = (dtime1 - 1096.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=366 or 732
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        elseif (yrflag.eq.4) then
          dtime1 = (dtime1 - 731.0d0) +  &
                   wndrep*int((dtime+0.00001d0)/wndrep)  !wndrep=365 or 731
          if     (dtime1.lt.dtime0) then
            dtime1 = dtime1 + wndrep
          endif
        endif
!           if     (mnproc.eq.1) then
!           write(lp,*) ' exit rdpall - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
      endif
!
! --- linear interpolation in time.
      w0 = (dtime1-dtime)/(dtime1-dtime0)
      w1 = 1.0 - w0
!           if     (mnproc.eq.1) then
!           write(lp,*) 'rdpall - dtime,w0,w1 = ',dtime,w0,w1
!           endif !1st tile
!           call xcsync(flush_lp)
      return
      end
!
!
      subroutine forfunhz
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- high frequency atmospheric forcing field processing.
! --- set all fields to zero.
!
      integer   i,j,l
!
      if     (mnproc.eq.1) then
      write (lp,*) ' now zeroing forcing fields ...'
      endif !1st tile
      call xcsync(flush_lp)
!
!$OMP   PARALLEL DO PRIVATE(l,j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
      do l=1,natm
        do j=1-nbdy,jj+nbdy
          do i=1-nbdy,ii+nbdy
              taux(i,j,l) = 0.0
              tauy(i,j,l) = 0.0
            wndspd(i,j,l) = 0.0
            ustara(i,j,l) = 0.0
            airtmp(i,j,l) = 0.0
            vapmix(i,j,l) = 0.0
            mslprs(i,j,l) = 0.0
            precip(i,j,l) = 0.0
            radflx(i,j,l) = 0.0
             swflx(i,j,l) = 0.0
            surtmp(i,j,l) = 0.0
            seatmp(i,j,l) = 0.0
            stoc_t(i,j,l) = 0.0
            stoc_s(i,j,l) = 0.0
            stoc_u(i,j,l) = 0.0
            stoc_v(i,j,l) = 0.0
          enddo !i
        enddo !j
      enddo !l
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished zeroing forcing fields'
      endif !1st tile
      call xcsync(flush_lp)
      return
      end
!
!
      subroutine forfunc
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize input of chl forcing field (into akpar array)
! --- call either forfunc or forfunk, not both
!
! --- units of chl are mg/m^3
! --- chl    is always on the p grid.
!
! --- I/O and array I/O unit 919 is reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer   l,lgth
      character preambl(5)*79
!
      mreck=1
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening chl field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      if (thermo) then
        call zaiopf(flnmfor(1:lgth)//'forcing.chl.a', 'old', 919)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+919,file=flnmfor(1:lgth)//'forcing.chl.b', &
           status='old', action='read')
        read (uoff+919,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 919)
        if     (kparan) then
          if     (mnproc.eq.1) then
          write (lp,*)
          write (lp,*) '***** annual chl *****'
          write (lp,*)
          endif !1st tile
          call xcsync(flush_lp)
          do l= 1,4
            akpar(:,:,l) = util1(:,:)
          enddo
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close (unit=uoff+919)
          endif
          call zaiocl(919)
        endif !kparan
!
      else  ! .not.thermo
        kparan = .true.  
      endif                    !  thermo
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening chl field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunk
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize input of kpar forcing field
! --- call either forfunc or forfunk, not both
!
! --- units of akpar are 1/m 
! --- akpar  is always on the p grid.
!
! --- I/O and array I/O unit 919 is reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer   l,lgth
      character preambl(5)*79
!
      mreck=1
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening kpar field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      if (thermo) then
        call zaiopf(flnmfor(1:lgth)//'forcing.kpar.a', 'old', 919)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+919,file=flnmfor(1:lgth)//'forcing.kpar.b', &
           status='old', action='read')
        read (uoff+919,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 919)
        if     (kparan) then
          if     (mnproc.eq.1) then
          write (lp,*)
          write (lp,*) '***** annual kpar *****'
          write (lp,*)
          endif !1st tile
          call xcsync(flush_lp)
          do l= 1,4
            akpar(:,:,l) = util1(:,:)
          enddo
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close (unit=uoff+919)
          endif
          call zaiocl(919)
!diag     call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,86400.*36000., &
!diag          'river precip. (cm/year) ')
        endif !kparan
!
      else  ! .not.thermo
        kparan = .true.  
      endif                    !  thermo
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening kpar field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunp
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize input of river (precip bogas) forcing field
!
! --- units of rivers are m/s    (positive into ocean)
! --- rivers is always on the p grid.
!
! --- I/O and array I/O unit 918 is reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      logical   lexist
      integer   l,lgth
      character preambl(5)*79
!
      mrecr=1
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening rivers  field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      if (thermo) then
!
      if     (.not.priver) then
        if     (mnproc.eq.1) then
        write (lp,*)
        write (lp,*) '***** no river precipitation *****'
        write (lp,*)
        endif !1st tile
        call xcsync(flush_lp)
        rivers(:,:,:) = 0.0
        rivera = .true.  
      else
        call zaiopf(flnmfor(1:lgth)//'forcing.rivers.a', 'old', 918)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+918,file=flnmfor(1:lgth)//'forcing.rivers.b', &
           status='old', action='read')
        read (uoff+918,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(util1, 918)
        if     (rivera) then
          if     (mnproc.eq.1) then
          write (lp,*)
          write (lp,*) '***** annual river precipitation *****'
          write (lp,*)
          endif !1st tile
          call xcsync(flush_lp)
          do l= 1,4
            rivers(:,:,l) = util1(:,:)
          enddo
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
          close (unit=uoff+918)
          endif
          call zaiocl(918)
!diag     call prtmsk(ip,util1,util2,idm,idm,jdm,  0.,86400.*36000., &
!diag          'river precip. (cm/year) ')
        endif
      endif
!
      else  ! .not.thermo
        rivers(:,:,:) = 0.0
        priver = .false.
        rivera = .true.  
      endif                    !  thermo
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening river   field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunr
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize input of thermal/tracer relaxation forcing fields
!
! --- thmean is long term mean vertically averaged potential density
! --- sshgmn is long term mean sea surface height (m)
!
! --- sefold is e-folding time of SSS relaxation (days)
! ---         with 0.0 indicating no relaxation at that point,
! ---         converted here to rmus (1/s)
! --- sssrmx is maximum SSS difference for SSS relaxation (psu)
!
! --- rmu    is a single field specifying 1/e-folding time (1/s)
! ---         set to zero where there is no thermal boundary relaxation
! --- twall  is temperature climatology for all layers
! --- swall  is salinity    climatology for all layers
! --- pwall  is interface   climatology for all layers (pressure units, Pa)
!
! --- rmutr  is a single field specifying 1/e-folding time (1/s)
! ---         set to zero where there is no tracer boundary relaxation
! --- trwall is tracer climatology for all layers and tracers
!
! --- I/O and array I/O units 911-914 are reserved for the entire run.
! --- I/O and array I/O unit  915 is used but not reserved.
!
! --- all input fields must be defined at all grid points
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer   lgth,mo
      character preambl(5)*79
      character cline*80
      real      one,oneps,sefsec
      integer   i,incmon,ios,j,k,ktr,nrec
!
! --- initialize all rmu fields to zero
!
      rmu(   :,:) = 0.0  !needed for thermf
      rmutra(:,:) = 0.0
      rmunp( :,:) = 0.0
      rmunv( :,:) = 0.0  !needed for thermf
!
      lgth = len_trim(flnmforw)
!
      if     (sshflg.eq.1) then
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening mean SSH   fields ...'
        endif !1st tile
        call xcsync(flush_lp)
        call zaiopf(flnmforw(1:lgth)//'relax.ssh.a', 'old', 915)
        call rdmonth(thmean, -915)  !no .b file (mean depth averaged density)
        call rdmonth(sshgmn, -915)  !no .b file
        call zaiocl(915)
        do j= 1,jj
          do i= 1,ii
            if     (ip(i,j).eq.1) then
              sshgmn(i,j) = sshgmn(i,j)*g  !input is mean ssh in m
            else
              thmean(i,j) = 0.0
              sshgmn(i,j) = 0.0
            endif
          enddo !i
        enddo !j
        call xctilr(thmean,1,1, nbdy,nbdy, halo_ps)
        call xctilr(sshgmn,1,1, nbdy,nbdy, halo_ps)
      elseif (sshflg.eq.2) then
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening mean SSH & Montg. Pot. fields ...'
        endif !1st tile
        call xcsync(flush_lp)
        call zaiopf(flnmforw(1:lgth)//'relax.montg.a', 'old', 915)
        call rdmonth(thmean, -915)  !no .b file (Montg. Pot. correction)
        call rdmonth(sshgmn, -915)  !no .b file
        call zaiocl(915)
        do j= 1,jj
          do i= 1,ii
            if     (ip(i,j).eq.1) then
              montg_c(i,j) = (sshgmn(i,j)-thmean(i,j))*g  !input is mean ssh/montg in m
            else
              montg_c(i,j) = 0.0
            endif
          enddo !i
        enddo !j
        call xctilr(montg_c,1,1, nbdy,nbdy, halo_ps)

      endif !sshflg
!
      if (.not.relaxf) then
        return
      endif
!
! --- read fields needed for boundary and surface relaxation
!
      mrecc=1
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening relaxation fields ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      if     (sefold.gt.0.0) then  !constant rmus
        rmus(:,:) = 1.0/(sefold*86400.0)   !from days to 1/s
      else
        call xcsync(flush_lp)
        lgth = len_trim(flnmfor)
        call zaiopf(flnmfor(1:lgth)//'relax.sefold.a', 'old', 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+915,file=flnmfor(1:lgth)//'relax.sefold.b', &
           status='old', action='read')
        read (uoff+915,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(rmus, 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+915)
        endif !1st tile
        call zaiocl(915)
!$OMP   PARALLEL DO PRIVATE(j,i,sefsec) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1,jj
          do i= 1,ii
            sefsec = rmus(i,j)*86400.0  !from days to s
            if     (sefsec.ge.epsil) then
              rmus(i,j)=1.0/sefsec  !from s to 1/s
            else
              rmus(i,j)=0.0         !sefold==0 turns off relaxation
            endif !sefold>0:else
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening sefold field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !cbar
!
      if     (sssflg.eq.-1) then  ! sss relaxation limiter
        call zaiopf(flnmforw(1:lgth)//'relax.sssrmx.a', 'old', 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+915,file=flnmforw(1:lgth)//'relax.sssrmx.b', &
              status='old', action='read')
        read (uoff+915,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(sssrmx, 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+915)
        endif
        call zaiocl(915)
      else
        if     (mnproc.eq.1) then
        write (lp,*) 'No sss relaxation limiter.'
        endif !1st tile
        sssrmx(:,:) = 99.9  !needed for thermf, set to no limit
      endif
!
      if     (relax) then  ! boundary thermal relaxation
        call zaiopf(flnmforw(1:lgth)//'relax.rmu.a', 'old', 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+915,file=flnmforw(1:lgth)//'relax.rmu.b', &
              status='old', action='read')
        read (uoff+915,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(rmu, 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+915)
        endif
        call zaiocl(915)
      else
        if     (mnproc.eq.1) then
        write (lp,*) 'No thermal relaxation mask.'
        endif !1st tile
      endif
!
      if     (trcrlx) then  ! boundary tracer relaxation
        call zaiopf(flnmforw(1:lgth)//'relax.rmutr.a', 'old', 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+915,file=flnmforw(1:lgth)//'relax.rmutr.b', &
              status='old', action='read')
        endif !1st tile
        call zagetc(cline,ios, uoff+915)  !1st line of the header on all tiles
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        rewind uoff+915
        read  (uoff+915,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        call rdmonth(rmutra, 915)
        do j= 1,jj
          do i= 1,ii
            rmutr(i,j,1) = rmutra(i,j)
          enddo !i
        enddo !j
        if     (cline.eq.'Relaxation Masks') then  !multiple masks
          do ktr= 2,ntracr
            call rdmonth(rmutr(1-nbdy,1-nbdy,ktr), 915)
            do j= 1,jj
              do i= 1,ii
                rmutra(i,j) = max(rmutra(i,j), rmutr(i,j,ktr) )
              enddo !i
            enddo !j
          enddo !ktr
        else !one mask (repeated)
          do ktr= 2,ntracr
            do j= 1,jj
              do i= 1,ii
                rmutr(i,j,ktr) = rmutra(i,j)
              enddo !i
            enddo !j
          enddo !ktr
        endif !multiple or one mask(s)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+915)
        endif
        call zaiocl(915)
      elseif (ntracr.gt.0) then
        if     (mnproc.eq.1) then
        write (lp,*) 'No tracer relaxation mask.'
        endif !1st tile
      endif
!
      call zaiopf(flnmforw(1:lgth)//'relax.temp.a', 'old', 911)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+911,file=flnmforw(1:lgth)//'relax.temp.b', &
            status='old', action='read')
      read (uoff+911,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      if     (relaxs) then  ! surface only
        call rdmonth(util1, 911)
        do k=2,kk
          call skmonth(       911)
        enddo
      else
        do k=1,kk
          call rdmonth(util1, 911)
        enddo
      endif
!
      call zaiopf(flnmforw(1:lgth)//'relax.saln.a', 'old', 912)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+912,file=flnmforw(1:lgth)//'relax.saln.b', &
            status='old', action='read')
      read (uoff+912,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      if     (relaxs) then  ! surface only
        call rdmonth(util1, 912)
        do k=2,kk
          call skmonth(       912)
        enddo
      else
        do k=1,kk
          call rdmonth(util1, 912)
        enddo
      endif
!
      if     (relaxt .or. .not.relaxs) then
        call zaiopf(flnmforw(1:lgth)//'relax.intf.a', 'old', 913)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+913,file=flnmforw(1:lgth)//'relax.intf.b', &
              status='old', action='read')
        read (uoff+913,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        do k=1,kk
          call rdmonth(util1, 913)
        enddo
      endif
!
      if     (relaxt) then
        call zaiopf(flnmforw(1:lgth)//'relax.trcr.a', 'old', 914)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+914,file=flnmforw(1:lgth)//'relax.trcr.b', &
              status='old', action='read')
        read (uoff+914,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        do ktr= 1,ntracr
          do k=1,kk
            call rdmonth(util1, 914)
          enddo
        enddo
      endif
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening relaxation fields'
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfuns
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize SAL factor field
!
! --- salfac is unitless and on the p-grid.
!
! --- I/O and array I/O unit 925 used here, but not reserved.
!
! --- all input fields must be defined at all grid points
!
      integer   i,j,k,l,lgth
!
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening salfac field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      call zaiopf(flnmfor(1:lgth)//'tidal.sal.a', 'old', 925)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+925,file=flnmfor(1:lgth)//'tidal.sal.b', &
         status='old', action='read')
      endif !1st tile
      call rdmonth(salfac, 925)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close (unit=uoff+925)
      endif
      call zaiocl(925)
!
      call xctilr(salfac,1,1, nbdy,nbdy, halo_ps)
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening salfac field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunt
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize spacially varying minimum depth for isopycnal layers
!
! --- units of isotop are m on the p-grid.
!
! --- I/O and array I/O unit 924 used here, but not reserved.
!
! --- all input fields must be defined at all grid points
!
      integer   k,l,lgth
!
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening isotop field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      call zaiopf(flnmfor(1:lgth)//'iso.top.a', 'old', 924)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+924,file=flnmfor(1:lgth)//'iso.top.b', &
         status='old', action='read')
      endif !1st tile
      call rdmonth(topiso, 924)
      vland = -isotop  !should be the typical shallowest depth (m)
      call xctilr( topiso,1,1, nbdy,nbdy, halo_ps)
      vland = 0.0
! --- convert to pressure units (Pa).
      topiso(:,:) = onem*topiso(:,:)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close (unit=uoff+924)
      endif
      call zaiocl(924)
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening isotop field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfunv
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize spacially varying isopycnal target densities
!
! --- units of sigma are sigma-theta or sigma-2.
! --- sigma  is always on the p grid.
!
! --- I/O and array I/O unit 922 used here, but not reserved.
!
! --- all input fields must be defined at all grid points
!
      integer   k,l,lgth
!
      if     (mnproc.eq.1) then
      write (lp,*) ' now opening sigma field  ...'
      endif !1st tile
      call xcsync(flush_lp)
!
      lgth = len_trim(flnmfor)
!
      call zaiopf(flnmfor(1:lgth)//'iso.sigma.a', 'old', 922)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      open (unit=uoff+922,file=flnmfor(1:lgth)//'iso.sigma.b', &
         status='old', action='read')
      endif !1st tile
      do k= 1,kk
        call rdmonth(util1, 922)
        vland = sigma(k)
        call xctilr( util1, 1,1, nbdy,nbdy, halo_ps)
        vland = 0.0
! ---   subtract constant 'thbase' from theta to reduce roundoff errors
        theta(:,:,k) = util1(:,:) - thbase
      enddo
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close (unit=uoff+922)
      endif
      call zaiocl(922)
!
      if     (mnproc.eq.1) then
      write (lp,*) ' ...finished opening sigma field '
      endif !1st tile
      call xcsync(flush_lp)
!
      return
      end
!
!
      subroutine forfundf
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
! --- initialize spacially varying veldf2, veldf4, thkdf[24], cbar and cb,
! ---  if necessary,
! ---  and calculate laplacian and biharmonic diffusion coefficients.
! --- initialize spacially varying diws, diwm, diwbot and diwqh0,
! ---   if necessary.
! --- alll are on the p-grid.
!
! --- veldf2 is diffusion velocity for laplacian  background diffusion
! --- veldf4 is diffusion velocity for biharmonic background diffusion
! --- at most one of thkdf2 and thkdf4 are active
! --- thkdf2 is diffusion velocity for laplacian  thickness  diffusion
! --- thkdf4 is diffusion velocity for biharmonic thickness  diffusion
! --- cbar   is rms flow speed     for linear bottom friction
! --- all units are (m/s)
!
! --- cb     is coefficient of quadratic bottom friction (unitless)
!
! --- diws   is background/internal wave diffusivity (m**2/s)
! --- diwm   is background/internal wave viscosity   (m**2/s)
! --- diwbot is background diffusivity at the bottom (m**2/s)
! --- diwqh0 is vertical scale for background diffusivity (m)
! ---           but diwqh0 converted here to 1/pressure-units
!
! --- I/O and array I/O unit 923 used here, but not reserved
!
! --- all input fields must be defined at all grid points
!
      logical   lthkdf4
      integer   i,j,l,lgth
!
      if     (veldf2.ge.0.0) then
        util2(:,:) = veldf2
      else
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening veldf2 field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'veldf2.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'veldf2.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util2, 923)
        call xctilr( util2, 1,1,    1,   1, halo_ps)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening veldf2 field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !veldf2
!
      if     (veldf4.ge.0.0) then
        util3(:,:) = veldf4
      else
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening veldf4 field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'veldf4.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'veldf4.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util3, 923)
        call xctilr( util3, 1,1,    1,   1, halo_ps)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening veldf4 field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !veldf4
!
      if     (thkdf4.gt.0.0) then
        util1(:,:) = thkdf4
        lthkdf4    = .true.
      elseif (thkdf2.gt.0.0) then
        util1(:,:) = thkdf2
        lthkdf4    = .false.
      elseif (thkdf2.eq.thkdf4) then  !both are zero
        util1(:,:) = 0.0
        lthkdf4    = .false.
      elseif (thkdf2.lt.0.0) then
        lthkdf4 = .false.
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening thkdf2 field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'thkdf2.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'thkdf2.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util1, 923)
        call xctilr( util1, 1,1,    1,   1, halo_ps)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening thkdf2 field '
        endif !1st tile
        call xcsync(flush_lp)
      else   !thkdf4.lt.0.0
        lthkdf4 = .true.
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening thkdf4 field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'thkdf4.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'thkdf4.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(util1, 923)
        call xctilr( util1, 1,1,    1,   1, halo_ps)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening thkdf4 field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !thkdf[24]
!
      if     (cbar.ge.0.0) then
        cbarp(:,:) = cbar
      else
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening cbar field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'cbar.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'cbar.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(cbarp, 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        call xctilr( cbarp, 1,1, nbdy,nbdy, halo_ps)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening cbar field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !cbar
!
      if     (cb.ge.0.0) then
        cbp(:,:) = cb
      else
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening cb field  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'cb.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'cb.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(cbp, 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        call xctilr( cbp, 1,1, nbdy,nbdy, halo_ps)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening cb field '
        endif !1st tile
        call xcsync(flush_lp)
      endif !cbar
!
      if     (difsiw.ge.0.0) then !difmiw also +ve
        diws(:,:) = difsiw
        diwm(:,:) = difmiw
      else
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening diws and diwm fields  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'diwsm.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'diwsm.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(diws, 923)
        call xctilr( diws, 1,1,    1,   1, halo_ps)
        call rdmonth(diwm, 923)
        call xctilr( diwm, 1,1,    1,   1, halo_ps)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening diws and diwm fields '
        endif !1st tile
        call xcsync(flush_lp)
      endif !difsiw
!
      if     (botdiw) then
        if     (mnproc.eq.1) then
        write (lp,*) ' now opening diwbot fields  ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        lgth = len_trim(flnmfor)
!
        call zaiopf(flnmfor(1:lgth)//'diwbot.a', 'old', 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+923,file=flnmfor(1:lgth)//'diwbot.b', &
           status='old', action='read')
        endif !1st tile
        call rdmonth(diwbot, 923)
        call rdmonth(diwqh0, 923)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+923)
        endif !1st tile
        call zaiocl(923)
!
        call xctilr( diwbot, 1,1, nbdy,nbdy, halo_ps)
        call xctilr( diwqh0, 1,1, nbdy,nbdy, halo_ps)
        diwqh0(:,:) = 1.0 / (diwqh0(:,:)*onem) !1/h0 in pressure units
!
        if     (mnproc.eq.1) then
        write (lp,*) ' ...finished opening diwbot fields '
        endif !1st tile
        call xcsync(flush_lp)
      endif !botdiw
!
! --- diffusion coeficients
!
      if     (momtyp.eq.2) then
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1,jj
          do i= 1,ii
            if     (iu(i,j).eq.1) then
              veldf2u(i,j) = 0.5*(util2(i,j)+util2(i-1,j))* &
                                  aspux(i,j)
              veldf4u(i,j) = 0.5*(util3(i,j)+util3(i-1,j))* &
                                  aspux(i,j)
              if     (lthkdf4) then
              thkdf4u(i,j) = 0.5*(util1(i,j)+util1(i-1,j))* &
                                 (aspux(i,j)**3)*scuy(i,j)
              else  !use thkdf4u array for thkdf2u
              thkdf4u(i,j) = 0.5*(util1(i,j)+util1(i-1,j))* &
                                  aspux(i,j)    *scuy(i,j)
              endif !lthkdf4
            else
              veldf2u(i,j) = 0.0
              veldf4u(i,j) = 0.0
              thkdf4u(i,j) = 0.0
            endif
            if     (iv(i,j).eq.1) then
              veldf2v(i,j) = 0.5*(util2(i,j)+util2(i,j-1))* &
                                  aspvy(i,j)
              veldf4v(i,j) = 0.5*(util3(i,j)+util3(i,j-1))* &
                                  aspvy(i,j)
              if     (lthkdf4) then
              thkdf4v(i,j) = 0.5*(util1(i,j)+util1(i,j-1))* &
                                 (aspvy(i,j)**3)*scvx(i,j)
              else  !use thkdf4v array for thkdf2v
              thkdf4v(i,j) = 0.5*(util1(i,j)+util1(i,j-1))* &
                                  aspvy(i,j)    *scvx(i,j)
              endif !lthkdf4
            else
              veldf2v(i,j) = 0.0
              veldf4v(i,j) = 0.0
              thkdf4v(i,j) = 0.0
            endif
          enddo !i
        enddo !j
        call xctilr(veldf2u, 1,1, nbdy,nbdy, halo_us)
        call xctilr(veldf4u, 1,1, nbdy,nbdy, halo_us)
        call xctilr(thkdf4u, 1,1, nbdy,nbdy, halo_us)
        call xctilr(veldf2v, 1,1, nbdy,nbdy, halo_vs)
        call xctilr(veldf4v, 1,1, nbdy,nbdy, halo_vs)
        call xctilr(thkdf4v, 1,1, nbdy,nbdy, halo_vs)
      else !momtyp==3,4
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1,jj
          do i= 1,ii
! ---       veldf[24]u are veldf[24]p, and veldf[24]v not used
            if     (ip(i,j).eq.1) then
              veldf2u(i,j) = util2(i,j) !veldf2p
              veldf4u(i,j) = util3(i,j) !veldf4p
            else
              veldf2u(i,j) = 0.0
              veldf4u(i,j) = 0.0
            endif
            if     (iu(i,j).eq.1) then
              if     (lthkdf4) then
              thkdf4u(i,j) = 0.5*(util1(i,j)+util1(i-1,j))* &
                                 (aspux(i,j)**3)*scuy(i,j)
              else  !use thkdf4u array for thkdf2u
              thkdf4u(i,j) = 0.5*(util1(i,j)+util1(i-1,j))* &
                                  aspux(i,j)    *scuy(i,j)
              endif !lthkdf4
            else
              thkdf4u(i,j) = 0.0
            endif
            if     (iv(i,j).eq.1) then
              if     (lthkdf4) then
              thkdf4v(i,j) = 0.5*(util1(i,j)+util1(i,j-1))* &
                                 (aspvy(i,j)**3)*scvx(i,j)
              else  !use thkdf4v array for thkdf2v
              thkdf4v(i,j) = 0.5*(util1(i,j)+util1(i,j-1))* &
                                  aspvy(i,j)    *scvx(i,j)
              endif !lthkdf4
            else
              thkdf4v(i,j) = 0.0
            endif
          enddo !i
        enddo !j
        call xctilr(veldf2u, 1,1, nbdy,nbdy, halo_ps)  !veldf2p
        call xctilr(veldf4u, 1,1, nbdy,nbdy, halo_ps)  !veldf4p
        call xctilr(thkdf4u, 1,1, nbdy,nbdy, halo_us)
        call xctilr(thkdf4v, 1,1, nbdy,nbdy, halo_vs)
      endif !momtyp
!
      return
      end
!
!
      subroutine preambl_print(preambl)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      character preambl(5)*79
!
! --- print non-blank lines of preambl
!
      integer i
!
      if     (mnproc.eq.1) then
        write(lp,*)
        do i= 1,5
          if     (len_trim(preambl(i)).ne.0) then
            write(lp,'(a)') trim(preambl(i))
          endif
        enddo
      endif !1st tile
      call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdmonth(field, iunit)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer   iunit
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
                field
!
! --- read a single array field from unit iunit (see also rdmonthck).
! --- ignore the input month, so long as it is between 1 and 12.
!
! --- iunit=900-910; atmospheric forcing field
! --- iunit=911-914; relaxation  forcing field
! --- iunit=915;     relaxation  time scale field
! --- iunit=918;     river       forcing field
! --- iunit=919;     kpar or chl forcing field
! --- iunit=922;     isopycnal target density field
! --- iunit=923;     laplacian or biharmonic diffusion velocity field
! --- iunit=924;     minimum depth for isopycnal layers
! --- iunit=925;     tidal drag (dragrh or drgten) or SAL (salfac)
!
! --- most of work now done by rdmonthck
!
      call rdmonthck(field, iunit, 0)
      return
      end
!
!
      subroutine rdmonthck(field, iunit, mnthck)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer   iunit,mnthck
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
                field
!
! --- read a single array field from unit iunit (see also rdmonth).
! --- check the input month against mnthck, if mnthck>0
!
! --- iunit=900-910; atmospheric forcing field
! --- iunit=911-914; relaxation  forcing field
! --- iunit=915;     relaxation  time scale field
! --- iunit=916;     offset      forcing field
! --- iunit=918;     river       forcing field
! --- iunit=919;     kpar or chl forcing field
! --- iunit=922;     isopycnal target density field
! --- iunit=923;     laplacian or biharmonic diffusion velocity field
! --- iunit=924;     minimum depth for isopycnal layers
! --- iunit=925;     tidal drag (dragrh or drgten) or SAL (salfac)
!
      integer   i,ios,layer,mnth
      real      denlay,hmina,hminb,hmaxa,hmaxb
      character cline*80
!
      if     (iunit.lt.0) then
!
! ---   special case, no .b file.
!
        call zaiord(field,ip,.false., hmina,hmaxa, &
                    -iunit)
        return
      endif
!
      call zagetc(cline,ios, uoff+iunit)
      if     (ios.ne.0) then
        if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in rdmonth - hit end of input'
          write(lp,*) 'iunit,ios = ',iunit,ios
          write(lp,*)
        endif !1st tile
        call xcstop('(rdmonth)')
               stop '(rdmonth)'
      endif
      if     (mnproc.eq.1) then
      write (lp,'(a)')  cline  !print input array info
      endif !1st tile
      i = index(cline,'=')
      if     (iunit.ge.900 .and. iunit.le.910) then
! ---   atmospheric forcing
        read (cline(i+1:),*) mnth,hminb,hmaxb
        if     (mnth.lt.1 .or. mnth.gt.12) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,a /)') &
            'error on unit',iunit,' - not monthly atmospheric data'
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
        if     (mnthck.gt.0 .and. mnth.ne.mnthck) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,a,a,2i4,a /)') &
            'error on unit',iunit,' - wrong atmospheric month', &
            ' (expected,input =',mnthck,mnth,')'
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
      elseif (iunit.eq.916) then
! ---   time-invarient heat flux correction
        read (cline(i+1:),*) hminb,hmaxb
      elseif (iunit.ge.911 .and. iunit.le.914) then
! ---   relaxation forcing
        read (cline(i+1:),*) mnth,layer,denlay,hminb,hmaxb
        if     (mnth.lt.1 .or. mnth.gt.12) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,a /)')  &
            'error on unit',iunit,' - not monthly relaxation data'
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
        if     (mnthck.gt.0 .and. mnth.ne.mnthck) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,a,a,2i4,a /)') &
            'error on unit',iunit,' - wrong relaxation month', &
            ' (expected,input =',mnthck,mnth,')'
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
      elseif (iunit.eq.919) then
! ---   kpar or chl forcing
        kparan = cline(i-8:i) .eq. ': range ='
        if     (kparan) then
! ---     annual
          read (cline(i+1:),*) hminb,hmaxb
        else
! ---     monthly
          read (cline(i+1:),*) mnth,hminb,hmaxb
          if     (mnth.lt.1 .or. mnth.gt.12) then
            if     (mnproc.eq.1) then
            write(lp,'(/ a,i4,a /)')  &
              'error on unit',iunit,' - not monthly kpar or chl data'
            endif !1st tile
            call xcstop('(rdmonth)')
                   stop '(rdmonth)'
          endif
          if     (mnthck.gt.0 .and. mnth.ne.mnthck) then
            if     (mnproc.eq.1) then
            write(lp,'(/ a,i4,a,a,2i4,a /)') &
             'error on unit',iunit,' - wrong kpar month', &
             ' (expected,input =',mnthck,mnth,')'
           endif !1st tile
           call xcstop('(rdmonth)')
                  stop '(rdmonth)'
         endif
        endif
      elseif (iunit.eq.918) then
! ---   river forcing
        rivera = cline(i-8:i) .eq. ': range ='
        if     (rivera) then
! ---     annual
          read (cline(i+1:),*) hminb,hmaxb
        else
! ---     monthly
          read (cline(i+1:),*) mnth,hminb,hmaxb
          if     (mnth.lt.1 .or. mnth.gt.12) then
            if     (mnproc.eq.1) then
            write(lp,'(/ a,i4,a /)')  &
              'error on unit',iunit,' - not monthly river data'
            endif !1st tile
            call xcstop('(rdmonth)')
                   stop '(rdmonth)'
          endif
          if     (mnthck.gt.0 .and. mnth.ne.mnthck) then
            if     (mnproc.eq.1) then
            write(lp,'(/ a,i4,a,a,2i4,a /)') &
             'error on unit',iunit,' - wrong river month', &
             ' (expected,input =',mnthck,mnth,')'
           endif !1st tile
           call xcstop('(rdmonth)')
                  stop '(rdmonth)'
         endif
        endif
      elseif (iunit.eq.915) then
! ---   relaxation time scale
        read (cline(i+1:),*) hminb,hmaxb
      elseif (iunit.eq.922) then
! ---   target density field.
        read (cline(i+1:),*) layer,hminb,hmaxb
        if     (hminb.gt.sigma(layer)+0.005 .or. &
                hmaxb.lt.sigma(layer)-0.005     ) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,a /)')  &
            'error on unit',iunit,' - not consistent with sigma(k)'
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
      elseif (iunit.eq.923) then
! ---   laplacian or biharmonic diffusion velocity field
        read (cline(i+1:),*) hminb,hmaxb
      elseif (iunit.eq.924) then
! ---   minimum depth for isopycnal layers
        read (cline(i+1:),*) hminb,hmaxb
      elseif (iunit.eq.925) then
! ---   tidal drag roughness or SAL
        read (cline(i+1:),*) hminb,hmaxb
      else
        if     (mnproc.eq.1) then
        write(lp,'(a,a / a,i5)') &
          'error - iunit must be 900-910 or 911-916', &
                             'or 918-919 or 922-925', &
          'iunit =',iunit
        endif !1st tile
        call xcstop('(rdmonth)')
               stop '(rdmonth)'
      endif
!
      if     (hminb.eq.hmaxb) then  !constant field
        field(:,:) = hminb
        call zaiosk(iunit)
      else
        call zaiord(field,ip,.false., hmina,hmaxa, &
                    iunit)
!
        if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or. &
                abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,i3 / a / a,1p3e14.6 / a,1p3e14.6 /)') &
            'error - .a and .b files not consistent:', &
            'iunit = ',iunit, &
            cline, &
            '.a,.b min = ',hmina,hminb,hmina-hminb, &
            '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
          endif !1st tile
          call xcstop('(rdmonth)')
                 stop '(rdmonth)'
        endif
      endif
!
      return
      end
!
!
      subroutine skmonth(iunit)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer   iunit
!
! --- skip a single array field from unit iunit.
!
! --- iunit=900-910; atmospheric forcing field
! --- iunit=911-914; relaxation  forcing field
! --- iunit=915;     relaxation strength field
! --- iunit=918;     river       forcing field
! --- iunit=919;     kpar or chl forcing field
!
      character cline*80
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        read (uoff+iunit,'(a)')  cline
!       write(lp,   '(a)')  cline
      endif
!
      call zaiosk(iunit)
!
      return
      end
!
!
      subroutine rdpall(dtime0,dtime1)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      implicit none
!
      real*8  dtime0,dtime1
!
! --- copy slot 2 into slot 1, and 
! --- read a set of high frequency forcing fields into slot 2.
! --- on exit, dtime0 and dtime1 are the associated times (wind days).
!
      integer i,j,k
      real    albw,degtorad
      real*8  dtime(895:910)
!
      integer, save :: icall = -1
!
      real, parameter :: sstmin = -1.8
      real, parameter :: sstmax = 35.0
!
      icall = icall + 1
!
      call rdpall1(  taux,dtime(901),901,mod(icall,3).eq.0)
      call rdpall1(  tauy,dtime(902),902,mod(icall,3).eq.0)
      if     (ustflg.eq.3) then
        call rdpall1(ustara,dtime(900),900,mod(icall,3).eq.0)
      else
        dtime(900) = dtime(901)
      endif
      if     (wndflg.lt.3) then
        call rdpall1(wndspd,dtime(903),903,mod(icall,3).eq.0)
      elseif (wndflg.ge.4) then
! ---   taux,tauy contains wndx,wndy
        dtime(903) = dtime(902)
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            wndspd(i,j,1) = wndspd(i,j,2)
            wndspd(i,j,2) = sqrt( taux(i,j,2)**2 + tauy(i,j,2)**2 )
          enddo
        enddo
      else
        dtime(903) = dtime(902)
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            wndspd(i,j,1) = wndspd(i,j,2)
          enddo
        enddo
        call str2spd(wndspd(1-nbdy,1-nbdy,2), &
                       taux(1-nbdy,1-nbdy,2), &
                       tauy(1-nbdy,1-nbdy,2) )
      endif !wndspd
      if     (flxflg.ne.3) then
        call rdpall1(airtmp,dtime(904),904,mod(icall,3).eq.1)
        call rdpall1(vapmix,dtime(905),905,mod(icall,3).eq.1)
      else
        dtime(904) = dtime(900)
        dtime(905) = dtime(900)
      endif
      if     (mslprf .or. flxflg.eq.6) then
        call rdpall1(mslprs,dtime(899),899,mod(icall,3).eq.2)
      else
        dtime(899) = dtime(905)
      endif
      if     (pcipf) then
        call rdpall1(precip,dtime(906),906,mod(icall,3).eq.1)
      else
        dtime(906) = dtime(905)
      endif
      call rdpall1(radflx,dtime(907),907,mod(icall,3).eq.2)
      call rdpall1( swflx,dtime(908),908,mod(icall,3).eq.2)
      if     (albflg.ne.0) then  !swflx is Qswdn
! ---   convert swflx to net shortwave into the ocean
! ---   shortwave through sea ice is handled separately
        if     (albflg.eq.1) then
          do j= 1-nbdy,jj+nbdy
            do i= 1-nbdy,ii+nbdy
              swflx(i,j,2) = swflx(i,j,2)*(1.0-0.09)  !NAVGEM albedo
            enddo
          enddo
        else   !albflg.eq.2
          degtorad = 4.d0*atan(1.d0)/180.d0
          do j= 1-nbdy,jj+nbdy
            do i= 1-nbdy,ii+nbdy
! ---         latitudinally-varying ocean albedo (Large and Yeager, 2009)
! ---         5.8% at the equator and 8% at the poles
              albw = ( 0.069 - 0.011*cos(2.0*degtorad*plat(i,j) ) )
              swflx(i,j,2) = swflx(i,j,2)*(1.0-albw)
            enddo
          enddo
        endif !albflg
      endif !Qswdn
      if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
              icmflg.eq.2 .or. ticegr.eq.0.0     ) then
        call rdpall1(surtmp,dtime(909),909,mod(icall,3).eq.2)
        if     (sstflg.ne.3) then  !use atmos. sst as "truth"
          do j= 1-nbdy,jj+nbdy
            do i= 1-nbdy,ii+nbdy
              seatmp(i,j,1) = max( sstmin, min(surtmp(i,j,1), sstmax ) )
              seatmp(i,j,2) = max( sstmin, min(surtmp(i,j,2), sstmax ) )
            enddo
          enddo
        endif
      else
        dtime(909) = dtime(905)
      endif !surtmp:else
      if     (sstflg.eq.3) then
        call rdpall1(seatmp,dtime(910),910,mod(icall,3).eq.2)
      else
        dtime(910) = dtime(905)
      endif
      if     (stfflg.gt.0) then
        call rdpall1(stoc_t,dtime(896),896,mod(icall,3).eq.0)
        call rdpall1(stoc_s,dtime(897),897,mod(icall,3).eq.0)
      else
        dtime(896) = dtime(905)
        dtime(897) = dtime(905)
      endif
      if     (stfflg.eq.2) then
        call rdpall1(stoc_u,dtime(895),895,mod(icall,3).eq.0)
        call rdpall1(stoc_v,dtime(898),898,mod(icall,3).eq.0)
      else
        dtime(895) = dtime(905)
        dtime(898) = dtime(905)
      endif
!
      dtime0 = dtime1
      dtime1 = dtime(901)
!
! --- check the input times.
      do k= 895,910
        if     (dtime(k).ne.dtime1) then
          if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in rdpall - inconsistent forcing times'
          write(lp,*) 'dtime0,dtime1 = ',dtime0,dtime1
          write(lp,*) 'dtime = ',dtime
          write(lp,*)
          endif !1st tile
          call xcstop('(rdpall)')
                 stop '(rdpall)'
        endif
      enddo
      return
      end
!
!
      subroutine rdpall1(field,dtime,iunit,lprint)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2) :: &
              field
      real*8  dtime
      integer iunit
      logical lprint
!
! --- copy field(:,:,2) into field(:,:,1), and
! --- read a high frequency forcing field into field(:,:,2).
! --- on exit, dtime is the time (wind day) of the forcing field.
!
      integer   i,ios,j
      character cline*80
      real      hmina,hminb,hmaxa,hmaxb,span
!
      call zagetc(cline,ios, uoff+iunit)
      if     (ios.lt.0) then  ! e-o-f
        if     (yrflag.eq.2 .or. yrflag.eq.4) then
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
!           write(lp,*) 'rdpall1 - rewind unit ',iunit
!           call flush(lp)
            rewind uoff+iunit
            read (uoff+iunit,*)
            read (uoff+iunit,*)
            read (uoff+iunit,*)
            read (uoff+iunit,*)
            read (uoff+iunit,*)
          endif
          call zaiorw(iunit)
          call zagetc(cline,ios, uoff+iunit)
        else
          if     (mnproc.eq.1) then
          write(lp,'(/ a,i4 /)') &
            'end of file from zagetc, iunit = ',iunit
          endif !1st tile
          call xcstop('(rdpall1)')
                 stop '(rdpall1)'
        endif
      endif  ! e-o-f
      if     (ios.gt.0) then
        if     (mnproc.eq.1) then
          write(lp,'(/ a,i4,i9 /)') &
            'I/O error from zagetc, iunit,ios = ',iunit,ios
        endif !1st tile
        call xcstop('(rdpall1)')
               stop '(rdpall1)'
      endif
      if     (lprint) then
        if     (mnproc.eq.1) then
        write (lp,'(a)') cline
        endif !1st tile
      endif !lprint
!
      i = index(cline,'=')
      read (cline(i+1:),*) dtime,span,hminb,hmaxb
      if     (span.gt.0.008d0) then
! ---   correct wind day to nearest 15 minutes
        dtime = nint(dtime*96.d0)/96.d0
      endif
!
      if     (hminb.eq.hmaxb) then  !constant field
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            field(i,j,1) = field(i,j,2)
            field(i,j,2) = hminb
          enddo
        enddo
        call zaiosk(iunit)
      else  !input field
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            field(i,j,1) = field(i,j,2)
          enddo
        enddo
        call zaiord(field(1-nbdy,1-nbdy,2),ip,.false., hmina,hmaxa, &
                    iunit)
!
        if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or. &
                abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,i3 / a / a,1p3e14.6 / a,1p3e14.6 /)') &
            'error - .a and .b files not consistent:', &
            'iunit = ',iunit, &
            cline, &
            '.a,.b min = ',hmina,hminb,hmina-hminb, &
            '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
          endif !1st tile
          call xcstop('(rdpall1)')
                 stop '(rdpall1)'
        endif
      endif
!
! --- msl pressure uses the the halo.
!
      if     (iunit.eq.899) then
        call xctilr(field(1-nbdy,1-nbdy,2),1,1, nbdy,nbdy, halo_ps)
      endif
!
! --- wind stress uses the the halo.
!
      if     (iunit.eq.901 .and. wndflg.eq.1) then  ! taux on u-grid
        call xctilr(field(1-nbdy,1-nbdy,2),1,1, nbdy,nbdy, halo_uv)
      elseif (iunit.eq.901) then                    ! taux on p-grid
        call xctilr(field(1-nbdy,1-nbdy,2),1,1, nbdy,nbdy, halo_pv)
      elseif (iunit.eq.902 .and. wndflg.eq.1) then  ! tauy on v-grid
        call xctilr(field(1-nbdy,1-nbdy,2),1,1, nbdy,nbdy, halo_vv)
      elseif (iunit.eq.902) then                    ! tauy on p-grid
        call xctilr(field(1-nbdy,1-nbdy,2),1,1, nbdy,nbdy, halo_pv)
      endif
      return
      end
!
!
      subroutine rdforf(mnth,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer lslot,mnth
!
! --- read forcing functions for one month.
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer i,irec,iunit,j
      real    albw,degtorad
!
      real, parameter :: sstmin = -1.8
      real, parameter :: sstmax = 35.0
!
      if     (mnth.le.mreca) then
!
! ---   rewind all units
!
        if     (mslprf .or. flxflg.eq.6) then
          iunit = 899
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              rewind uoff+iunit
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
            endif
            call zaiorw(iunit)
        endif
        if     (ustflg.eq.3) then
          iunit = 900
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              rewind uoff+iunit
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
            endif
            call zaiorw(iunit)
        endif
        if     (windf) then
          do iunit= 901,902
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              rewind uoff+iunit
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
            endif
            call zaiorw(iunit)
          enddo
        endif
        if     (thermo) then
          do iunit= max(903,901+min(wndflg,3)),908  !904,908 when wndflg>=3
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              rewind uoff+iunit
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
            endif
            call zaiorw(iunit)
          enddo
        endif
        if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
                icmflg.eq.2 .or. ticegr.eq.0.0     ) then
          iunit = 909
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            rewind uoff+iunit
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
          endif
          call zaiorw(iunit)
        endif !surtmp
        if     (sstflg.eq.3) then
          iunit = 910
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            rewind uoff+iunit
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
            read  (uoff+iunit,*)
          endif
          call zaiorw(iunit)
        endif
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdforf: mreca,mnth = ',mreca,mnth,'  (rewind)'
!       endif !1st tile
!       call xcsync(flush_lp)
        mreca = 0
      endif
!
! --- skip forward to desired month
!
      do irec= mreca+1,mnth-1
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdforf: mreca,mnth = ',mreca,mnth,
!    &              '  (skipping ',irec,')'
!       endif !1st tile
!       call xcsync(flush_lp)
        if     (mslprf .or. flxflg.eq.6) then
          call skmonth(899)
        endif
        if     (ustflg.eq.3) then
          call skmonth(900)
        endif
        if     (windf) then
          do iunit= 901,902
            call skmonth(iunit)
          enddo
        endif
        if     (thermo) then
          do iunit= max(903,901+min(wndflg,3)),908  !904,908 when wndflg>=3
            call skmonth(iunit)
          enddo
        endif
        if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
                icmflg.eq.2 .or. ticegr.eq.0.0     ) then
          call skmonth(909)
        endif !surtmp
        if     (sstflg.eq.3) then
          call skmonth(910)
        endif
      enddo
!
! --- read desired month
!
      if     (mslprf .or. flxflg.eq.6) then
        call rdmonthck(mslprs(1-nbdy,1-nbdy,lslot),899,mnth)
        call xctilr(mslprs(1-nbdy,1-nbdy,lslot),1,1, nbdy,nbdy, halo_ps)
      endif
      if     (windf) then
        call rdmonthck(taux(1-nbdy,1-nbdy,lslot),901,mnth)
        call rdmonthck(tauy(1-nbdy,1-nbdy,lslot),902,mnth)
        if     (wndflg.eq.1) then  !on uv-grids
          call xctilr(taux(1-nbdy,1-nbdy,lslot),1,1, nbdy,nbdy, halo_uv)
          call xctilr(tauy(1-nbdy,1-nbdy,lslot),1,1, nbdy,nbdy, halo_vv)
        else                       !on p-grid
          call xctilr(taux(1-nbdy,1-nbdy,lslot),1,1, nbdy,nbdy, halo_pv)
          call xctilr(tauy(1-nbdy,1-nbdy,lslot),1,1, nbdy,nbdy, halo_pv)
        endif
      else
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            taux(i,j,lslot) = 0.0
            tauy(i,j,lslot) = 0.0
          enddo
        enddo
      endif
      if     (thermo) then
        if     (ustflg.eq.3) then
          call rdmonthck(ustara(1-nbdy,1-nbdy,lslot),900,mnth)
        endif
        if (wndflg.lt.3) then
          call rdmonthck(wndspd(1-nbdy,1-nbdy,lslot),903,mnth)
        elseif (wndflg.ge.4) then
! ---     taux,tauy contains wndx,wndy
!$OMP     PARALLEL DO PRIVATE(j,i) &
!$OMP              SCHEDULE(STATIC,jblk)
          do j= 1-nbdy,jj+nbdy
            do i= 1-nbdy,ii+nbdy
              wndspd(i,j,lslot) = sqrt( taux(i,j,lslot)**2 + &
                                        tauy(i,j,lslot)**2  )
            enddo
          enddo
        else
          call str2spd(wndspd(1-nbdy,1-nbdy,lslot), &
                         taux(1-nbdy,1-nbdy,lslot), &
                         tauy(1-nbdy,1-nbdy,lslot) )
        endif !wndspd
        if     (flxflg.ne.3) then
          call rdmonthck(airtmp(1-nbdy,1-nbdy,lslot),904,mnth)
          call rdmonthck(vapmix(1-nbdy,1-nbdy,lslot),905,mnth)
        endif
        if     (pcipf) then
          call rdmonthck(precip(1-nbdy,1-nbdy,lslot),906,mnth)
        endif
        call rdmonthck(radflx(1-nbdy,1-nbdy,lslot),907,mnth)
        call rdmonthck( swflx(1-nbdy,1-nbdy,lslot),908,mnth)
        if     (albflg.ne.0) then  !swflx is Qswdn
! ---     convert swflx to net shortwave into the ocean
! ---     shortwave through sea ice is handled separately
          if     (albflg.eq.1) then
            do j= 1-nbdy,jj+nbdy
              do i= 1-nbdy,ii+nbdy
                swflx(i,j,lslot) = swflx(i,j,lslot)*(1.0-0.09)  !NAVGEM albedo
              enddo
            enddo
          else   !albflg.eq.2
            degtorad = 4.d0*atan(1.d0)/180.d0
            do j= 1-nbdy,jj+nbdy
              do i= 1-nbdy,ii+nbdy
! ---           latitudinally-varying ocean albedo (Large and Yeager, 2009)
! ---           5.8% at the equator and 8% at the poles
                albw = ( 0.069 - 0.011*cos(2.0*degtorad*plat(i,j) ) )
                swflx(i,j,lslot) = swflx(i,j,lslot)*(1.0-albw)
              enddo
            enddo
          endif !albflg
        endif !Qswdn
      else
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            wndspd(i,j,lslot) = 0.0
            airtmp(i,j,lslot) = 0.0
            vapmix(i,j,lslot) = 0.0
            precip(i,j,lslot) = 0.0
            radflx(i,j,lslot) = 0.0
             swflx(i,j,lslot) = 0.0
          enddo
        enddo
      endif
!
      if     (lwflag.eq.2 .or. sstflg.eq.2   .or. &
              icmflg.eq.2 .or. ticegr.eq.0.0     ) then
        call rdmonthck(surtmp(1-nbdy,1-nbdy,lslot),909,mnth)
        if     (sstflg.ne.3) then
          do j= 1-nbdy,jj+nbdy
            do i= 1-nbdy,ii+nbdy
              seatmp(i,j,lslot) = max( sstmin, &
                                       min(sstmax, surtmp(i,j,lslot) ) )
            enddo
          enddo
        endif
      endif !surtmp
!
      if     (sstflg.eq.3) then
        call rdmonthck(seatmp(1-nbdy,1-nbdy,lslot),910,mnth)
      endif
!
      mreca = mnth
!
      if     (mnproc.eq.1) then
      write (lp,'(2(a,i3))') ' forcing functions for month',mnth, &
         ' written into slot',lslot
      endif !1st tile
      call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdkpar(mnth,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer lslot,mnth
!
! --- read kpar or chl forcing for one month.
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer i,irec,iunit,j
!
      if     (kparan) then
        return   ! annual (constant) kpar forcing
      endif
!
      if     (mnth.le.mreck) then
!
! ---   rewind
!
        iunit=919
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          rewind uoff+iunit
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
        endif
        call zaiorw(iunit)
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdkpar: mreck,mnth = ',mreck,mnth,'  (rewind)'
!       endif !1st tile
!       call xcsync(flush_lp)
        mreck = 0
      endif
!
! --- skip forward to desired month
!
      do irec= mreck+1,mnth-1
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdkpar: mreck,mnth = ',mreck,mnth,
!    &              '  (skipping ',irec,')'
!       endif !1st tile
!       call xcsync(flush_lp)
        call skmonth(919)
      enddo
!
! --- read desired month
!
      call rdmonthck(akpar(1-nbdy,1-nbdy,lslot),919,mnth)
!
      mreck = mnth
!
      if     (mnproc.eq.1) then
      write (lp,'(2(a,i3))') ' kpar    forcing   for month',mnth, &
         ' written into slot',lslot
      endif !1st tile
      call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdrivr(mnth,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer lslot,mnth
!
! --- read river forcing for one month.
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      integer i,irec,iunit,j
!
      if     (rivera) then
        return   ! annual (constant) river forcing
      endif
!
      if     (mnth.le.mrecr) then
!
! ---   rewind
!
        iunit=918
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          rewind uoff+iunit
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
          read  (uoff+iunit,*)
        endif
        call zaiorw(iunit)
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdrivr: mreca,mnth = ',mreca,mnth,'  (rewind)'
!       endif !1st tile
!       call xcsync(flush_lp)
        mrecr = 0
      endif
!
! --- skip forward to desired month
!
      do irec= mrecr+1,mnth-1
!       if     (mnproc.eq.1) then
!       write(lp,*) 'rdforf: mrecr,mnth = ',mrecr,mnth,
!    &              '  (skipping ',irec,')'
!       endif !1st tile
!       call xcsync(flush_lp)
        call skmonth(918)
      enddo
!
! --- read desired month
!
      call rdmonthck(rivers(1-nbdy,1-nbdy,lslot),918,mnth)
!
      mrecr = mnth
!
      if     (mnproc.eq.1) then
      write (lp,'(2(a,i3))') ' rivers  forcing   for month',mnth, &
         ' written into slot',lslot
      endif !1st tile
      call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdrlax(month,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      integer lslot,month
!
! --- read relaxation fields for one month,
! --- monthly (clmflg==12) or bi-monthly (clmflg==6) data.
!
      integer         mreca,mrecc,mreck,mrecr
      common/rdforfi/ mreca,mrecc,mreck,mrecr
      save  /rdforfi/
!
      logical lexist,lfatal
      integer i,irec,iunit,j,k,ktr,mrec,mnth,mxunit
      real    p23min(2),shallow
!
      mnth=mod(month-1,12)+1
      if     (mnproc.eq.1) then
      write(lp,*) 'rdrlax - month = ',month,mnth
      endif !1st tile
      call xcsync(flush_lp)
!
      if     (relaxf) then
        if     (clmflg.eq.12) then
          mrec = mnth
        else
          mrec = (mnth+1)/2
        endif
        if     (relaxt) then
          mxunit = 914  ! tracers
        elseif (relaxs) then
          mxunit = 912  ! T&S only
        else
          mxunit = 913
        endif
!
        if     (mrec.le.mrecc) then
!
! ---     rewind all units
!
          do iunit= 911,mxunit
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              rewind uoff+iunit
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
              read  (uoff+iunit,*)
            endif
            call zaiorw(iunit)
          enddo
!         if     (mnproc.eq.1) then
!         write(lp,*) 'rdrlax: mrecc,mrec = ',mrecc,mrec,'  (rewind)'
!         endif !1st tile
!         call xcsync(flush_lp)
          mrecc = 0
        endif
!
! ---   skip forward to desired month
!
        do irec= mrecc+1,mrec-1
!         if     (mnproc.eq.1) then
!         write(lp,*) 'rdrlax: mrecc,mrec = ',mrecc,mrec,
!    &                '  (skipping ',irec,')'
!         endif !1st tile
!         call xcsync(flush_lp)
          do iunit= 911,mxunit
            if     (iunit.lt.914) then
              do k= 1,kk
                call skmonth(iunit)
              enddo
            else
              do ktr= 1,ntracr
                do k= 1,kk
                  call skmonth(iunit)
                enddo
              enddo
            endif
          enddo
        enddo
!
! --- read desired month
!
        if     (relaxs) then  ! surface only
          k=1
            call rdmonthck(twall(1-nbdy,1-nbdy,k,lslot),911,mnth)
          do k= 2,kk
            call skmonth(                               911)
          enddo
          k=1
            call rdmonthck(swall(1-nbdy,1-nbdy,k,lslot),912,mnth)
          do k= 2,kk
            call skmonth(                               912)
          enddo
          if (relaxt) then  !need pwall for tracers
            do k= 1,kk
              call rdmonthck(pwall(1-nbdy,1-nbdy,k,lslot),913,mnth)
            enddo
          endif
        else
          do k= 1,kk
            call rdmonthck(twall(1-nbdy,1-nbdy,k,lslot),911,mnth)
          enddo
          do k= 1,kk
            call rdmonthck(swall(1-nbdy,1-nbdy,k,lslot),912,mnth)
          enddo
          do k= 1,kk
            call rdmonthck(pwall(1-nbdy,1-nbdy,k,lslot),913,mnth)
          enddo
        endif
!
        if (relaxt) then
          do ktr= 1,ntracr
            do k= 1,kk
              call rdmonthck(trwall(1-nbdy,1-nbdy,k,lslot,ktr),914,mnth)
            enddo
          enddo
        endif
!
        mrecc = mrec
!
! ---   sanity check.
!
        if     (lslot.eq.1 .and. (relaxt .or. .not.relaxs)) then
          if     (isopyc) then
            shallow = dp0k(1)*5.0*qonem
          else
            shallow = sum(dp0k(1:min(max(5,nsigma+2),kk)))*qonem
          endif
          p23min(1) = huge(p23min(1))
          p23min(2) = huge(p23min(2))
          do j= 1,jj
            do i= 1,ii
              if     (depths(i,j).gt.shallow) then
                p23min(1) = min( p23min(1),  &
                                 pwall(i,j,min(2,kkwall),lslot) )
                p23min(2) = min( p23min(2),  &
                                 pwall(i,j,min(3,kkwall),lslot) )
              endif
            enddo
          enddo
          call xcminr(p23min)
          if     (p23min(1).eq.huge(p23min(1))) then
            if     (mnproc.eq.1) then
            write (lp,'(2a,f7.1)') &
              'rdrlax: could not check pwall.2, all depths below', &
              shallow
            endif !1st tile
          else
            shallow = dp0k(1)+dp0k(min(2,kk))
            lfatal = .false.
            if    (abs(p23min(1)-dp0k(1)).le.dp0k(1)*0.01) then
              if     (mnproc.eq.1) then
              write (lp,'(a,2f7.2)') &
                'rdrlax: pwall.2 ok; expected,input min depth =', &
                dp0k(1)*qonem,p23min(1)*qonem
              endif !1st tile
            else
              lfatal = .true.
              if     (mnproc.eq.1) then
              write (lp,'(a,2f7.2,a)') &
                'rdrlax: pwall.2 NOT ok; expected,input min depth =', &
                dp0k(1)*qonem,p23min(1)*qonem, &
                ' (bad climatology?)'
              endif !1st tile
            endif
            if     (abs(p23min(2)-shallow).le.shallow*0.01) then
              if     (mnproc.eq.1) then
              write (lp,'(a,2f7.2)') &
                'rdrlax: pwall.3 ok; expected,input min depth =', &
                shallow*qonem,p23min(2)*qonem
              endif !1st tile
            elseif (.not.isopyc .and. kk.gt.2) then
              lfatal = .true.
              if     (mnproc.eq.1) then
              write (lp,'(a,2f7.2,a)') &
                'rdrlax: pwall.3 NOT ok; expected,input min depth =', &
                shallow*qonem,p23min(2)*qonem, &
                ' (bad climatology?)'
              endif !1st tile
            endif
            if     (lfatal) then
              inquire( &
                file=flnmforw(1:len_trim(flnmforw))//'relax.weird', &
                exist=lexist)
              if     (lexist) then
                if     (mnproc.eq.1) then
                write (lp,'(3a)') &
                  'rdrlax: continuing because file ', &
                  flnmforw(1:len_trim(flnmforw))//'relax.weird', &
                  ' exists (ignore the "bad" climatology)'
                endif !1st tile
                call xcsync(flush_lp)
              else
                call xcsync(flush_lp)
                if     (mnproc.eq.1) then
                write (lp,'(3a)') &
                  'rdrlax: create an empty file ', &
                  flnmforw(1:len_trim(flnmforw))//'relax.weird', &
                  ' to ignore the "bad" climatology'
                endif !1st tile
                call xcstop('(rdrlax)')
                       stop '(rdrlax)'
              endif
            endif
          endif
        endif  ! sanity check
!
        if     (mnproc.eq.1) then
        write (lp,'(2(a,i3))') ' relaxation fields for month',mnth, &
           ' written into slot',lslot
        endif !1st tile
        call xcsync(flush_lp)
      endif
      return
      end
!
!
      subroutine rdbaro(dtime)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
!
! --- baroclinic velocity nesting archive input processing
!
! --- filenames  nest/arch[vm].????_???_??.[ab]
!
! --- I/O and array I/O unit 921 is reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      logical   larchm,lmonth
      save      larchm,lmonth
      integer   iarch
      save      iarch
      real*8    dbnstf,dbnsti,dtimei,dtime0,dtime1
      save      dbnstf,dbnsti,dtimei,dtime0,dtime1
!
      integer   iyr,mon,idy,ihr
      integer   i,j,k
      real*8    wdaymn(-1:2)
!
! --- wb0 negative on first call only.
      if     (wb0.lt.-1.0) then
!
! ---   initialize nesting fields
!
        if     (mnproc.eq.1) then
        write (lp,*) ' now initializing baro nesting fields ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        dbnstf = abs(bnstfq)
        if     (dbnstf.lt.1.0) then
          dbnstf = (baclin/86400.0d0)* &
                   max(1,nint((86400.0d0*dbnstf)/baclin))
        endif
!
        larchm = bnstfq.lt.0.0  !mean archives spaning -bnstfq days
        lmonth = yrflag.ge.2 .and. abs(bnstfq+30.5).lt.0.1  !monthly mean archives
        if     (lmonth) then
          call fordate(dtime,yrflag, iyr,mon,idy,ihr)
          if     (mon.eq.1) then
            call datefor(wdaymn(-1), iyr-1,   12,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,      1,1,0, yrflag)
            call datefor(wdaymn( 1), iyr,      2,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,      3,1,0, yrflag)
          elseif (mon.eq.11) then
            call datefor(wdaymn(-1), iyr,     10,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,     11,1,0, yrflag)
            call datefor(wdaymn( 1), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    1,1,0, yrflag)
          elseif (mon.eq.12) then
            call datefor(wdaymn(-1), iyr,     11,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 1), iyr+1,    1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    2,1,0, yrflag)
          else
            call datefor(wdaymn(-1), iyr,  mon-1,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,  mon,  1,0, yrflag)
            call datefor(wdaymn( 1), iyr,  mon+1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,  mon+2,1,0, yrflag)
          endif
          if     (dtime.lt.0.5d0*(wdaymn(0)+wdaymn(1))) then
            dtime0 = 0.5d0*(wdaymn(-1)+wdaymn( 0))
            dtime1 = 0.5d0*(wdaymn( 0)+wdaymn( 1))
          else
            dtime0 = 0.5d0*(wdaymn( 0)+wdaymn( 1))
            dtime1 = 0.5d0*(wdaymn( 1)+wdaymn( 2))
          endif
          lb0    = 1
          call rdbaro_in(dtime0,larchm,1)
          lb1    = 2
          call rdbaro_in(dtime1,larchm,2)
        else
          if     (larchm) then
            dbnsti = 0.5*dbnstf
          else
            dbnsti = 0.0
          endif
          dtimei = int((dtime-dbnsti)/dbnstf)*dbnstf + dbnsti
!
          dtime0 = dtimei
          lb0    = 1
          call rdbaro_in(dtime0,larchm,1)
!
          iarch  = 1
          dtime1 = dtimei + dbnstf
          lb1    = 2
          call rdbaro_in(dtime1,larchm,2)
        endif
!
        if     (mnproc.eq.1) then
        write (lp,*)
        write (lp,*) ' dtime,dtime0,dtime1 = ',dtime,dtime0,dtime1
        write (lp,*)
        write (lp,*) ' ...finished initializing baro nesting fields'
        endif !1st tile
        call xcsync(flush_lp)
      endif  ! initialization
!
      if     (dtime.gt.dtime1) then
!
! ---   get the next set of fields.
        do j= 1-nbdy,jj+nbdy
          do i= 1-nbdy,ii+nbdy
            ubnest(i,j,1) = ubnest(i,j,2)
            vbnest(i,j,1) = vbnest(i,j,2)
            ubpnst(i,j,1) = ubpnst(i,j,2)
            vbpnst(i,j,1) = vbpnst(i,j,2)
            pbnest(i,j,1) = pbnest(i,j,2)
          enddo
        enddo
        if     (lmonth) then
          dtime0 = dtime1
          call fordate(dtime1,yrflag, iyr,mon,idy,ihr)
          if     (mon.eq.11) then
            call datefor(wdaymn( 1), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    1,1,0, yrflag)
          elseif (mon.eq.12) then
            call datefor(wdaymn( 1), iyr+1,    1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    2,1,0, yrflag)
          else
            call datefor(wdaymn( 1), iyr,  mon+1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,  mon+2,1,0, yrflag)
          endif
          dtime1 = 0.5d0*(wdaymn( 1)+wdaymn( 2))
          call rdbaro_in(dtime1,larchm,2)
        else
          dtime0 = dtime1
          iarch  = iarch + 1
          dtime1 = dtimei + dbnstf*iarch
          call rdbaro_in(dtime1,larchm,2)
        endif
!
!           if     (mnproc.eq.1) then
!           write(lp,*) ' exit rdbaro_in - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
      endif  ! next set of fields.
!
! --- linear interpolation in time.
      wb0 = (dtime1-dtime)/(dtime1-dtime0)
      wb1 = 1.0 - wb0
!           if     (mnproc.eq.1) then
!           write(lp,*) 'rdbaro - dtime,wb0,wb1 = ',dtime,wb0,wb1
!           endif !1st tile
!           call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdbaro_in(dtime,larchm,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
      integer   lslot
      logical   larchm
!
! --- input barotropic fields from archive on model day dtime.
! --- filenames  nest/arch[vm].????_???_??.[ab]
! --- I/O and array I/O unit 921 is reserved for the entire run.
!
      character flnm*22, cline*80, cvarin*6, cfield*8
      integer   i,idmtst,ios,j,jdmtst,k,layer
      integer   iyear,iday,ihour
      logical   nodens
      real      hqpbot
!
      call forday(dtime, yrflag, iyear,iday,ihour)
!
      if     (larchm) then
        write(flnm,'("nest/archm.",i4.4,"_",i3.3,"_",i2.2)') &
                                   iyear,iday,ihour
      else
        write(flnm,'("nest/archv.",i4.4,"_",i3.3,"_",i2.2)') &
                                   iyear,iday,ihour
      endif
!
      if     (mnproc.eq.1) then
      write (lp,"(a,a,a,f12.5,a)") 'rdbaro_in: ',flnm," (",dtime,")"
      endif !1st tile
      call xcsync(flush_lp)
!
      call zaiopf(flnm//'.a','old', 921)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+921,file=flnm//'.b',form='formatted', &
              status='old',action='read')
!
        read(uoff+921,'(a)') cline
        read(uoff+921,'(a)') cline
        read(uoff+921,'(a)') cline
        read(uoff+921,'(a)') cline
!
        read(uoff+921,'(a)') cline
        read(uoff+921,'(a)') cline
        read(uoff+921,'(a)') cline
      endif !1st tile
!
      call zagetc(cline,ios, uoff+921)
      read(cline,*) idmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',idmtst
!     endif !1st tile
      if (cvarin.ne.'idm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdbaro_in - input ',cvarin, &
                              ' but should be idm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      call zagetc(cline,ios, uoff+921)
      read(cline,*) jdmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',jdmtst
!     endif !1st tile
      if (cvarin.ne.'jdm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdbaro_in - input ',cvarin, &
                              ' but should be jdm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
!
      if (idmtst.ne.itdm .or. jdmtst.ne.jtdm) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdbaro_in - input idm,jdm', &
                              ' not consistent with parameters'
        write(lp,*) 'idm,jdm = ',itdm,  jtdm,  '  (dimensions.h)'
        write(lp,*) 'idm,jdm = ',idmtst,jdmtst,'  (input)'
        write(lp,*)
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
!
      call zagetc(cline,ios, uoff+921)
!
! --- skip some surface fields.
!
      call rd_archive(util1, cfield,layer, 921)  ! montg1
      if     (cfield.ne.'montg1  ') then
        if     (mnproc.eq.1) then
        write(lp,'(/ a / a,a /)') cfield, &
               'error in rdbaro_in - expected ','montg1  '
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      nodens = layer.ne.0  !new or original archive type
!
      call rd_archive(util2, cfield,layer, 921)  ! srfhgt=montg1+svref*pbnest
      if     (cfield.ne.'srfhgt  ') then
        if     (mnproc.eq.1) then
        write(lp,'(/ a / a,a /)') cfield, &
               'error in rdbaro_in - expected ','srfhgt  '
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      call zagetc(cline,ios, uoff+921)  !steric or surflx
      call zaiosk(921)
      if     (cline(1:8).eq.'steric  ') then  !surflx
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+921,*)
        endif
        call zaiosk(921)
      endif
      if     (nodens) then
        call zagetc(cline,ios, uoff+921)  !wtrflx or salflx
        call zaiosk(921)
        if     (cline(1:8).eq.'wtrflx  ') then
          do i= 1,3 !salflx,dpbl,dpmixl
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              read (uoff+921,*)
            endif
            call zaiosk(921)
          enddo
        else  !.eq.salflx
          do i= 1,2 !dpbl,dpmixl
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              read (uoff+921,*)
            endif
            call zaiosk(921)
          enddo
        endif !wtrflx/salflx
      else
        do i= 1,8 !salflx,dpbl,dpmixl,tmix,smix,thmix,umix,vmix
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+921,*)
          endif
          call zaiosk(921)
        enddo
      endif !nodens:else
      call rd_archive(ubnest(1-nbdy,1-nbdy,lslot), cfield,layer, 921)
      if     (cfield.eq.'kemix   ') then
        call rd_archive(ubnest(1-nbdy,1-nbdy,lslot), cfield,layer, 921)
      endif
      if     (cfield.eq.'covice  ') then
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+921,*)
        endif
        call zaiosk(921)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+921,*)
        endif
        call zaiosk(921)
        call rd_archive(ubnest(1-nbdy,1-nbdy,lslot), cfield,layer, 921)
      endif
      if     (cfield.ne.'u_btrop ') then
        if     (mnproc.eq.1) then
        write(lp,'(/ a / a,a /)') cfield, &
               'error in rdbaro_in - expected ','u_btrop '
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      call rd_archive(vbnest(1-nbdy,1-nbdy,lslot), cfield,layer, 921)
      if     (cfield.ne.'v_btrop ') then
        if     (mnproc.eq.1) then
        write(lp,'(/ a / a,a /)') cfield, &
               'error in rdbaro_in - expected ','v_btrop '
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close( unit=uoff+921)
      endif
      call zaiocl(921)
!
      call xctilr(ubnest(1-nbdy,1-nbdy,lslot),1,1, 1,1, halo_uv)
      call xctilr(vbnest(1-nbdy,1-nbdy,lslot),1,1, 1,1, halo_vv)
!
!$OMP PARALLEL DO PRIVATE(j,i,hqpbot) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
        do i=1,ii
          if     (ip(i,j).eq.1) then
            hqpbot = 0.5/pbot(i,j)
            if     (sshflg.eq.2 .and. .not.lbmont) then
!             apply montgomery correction to input pb
              pbnest(i,j,lslot) = (util2(i,j)-util1(i,j)-montg_c(i,j)) &
                                  *rhoref
            else
!             montgomery correction either not needed (sshflg<2) 
!                                  or already applied (lbmont)
              pbnest(i,j,lslot) = (util2(i,j)-util1(i,j))*rhoref
            endif
            ubpnst(i,j,lslot) = (ubnest(i,  j,lslot)*depthu(i,  j)+ &
                                 ubnest(i+1,j,lslot)*depthu(i+1,j) ) &
                                *hqpbot
            vbpnst(i,j,lslot) = (vbnest(i,j,  lslot)*depthv(i,j  )+ &
                                 vbnest(i,j+1,lslot)*depthv(i,j+1) ) &
                                *hqpbot
          else
            pbnest(i,j,lslot) = 0.0
            ubpnst(i,j,lslot) = 0.0
            vbpnst(i,j,lslot) = 0.0
          endif
        enddo
      enddo
!
      if     (.false. .and. ittest.ne.-1 .and. jttest.ne.-1) then
        call xcsync(flush_lp)
        if     (i0.lt.ittest .and. i0+ii.ge.ittest .and. &
                j0.lt.jttest .and. j0+jj.ge.jttest      ) then
          write(lp,'(i5,i4,a,1p5e13.5)') &
             itest+i0,jtest+j0,' rdbaro: ub,vb,pb,ubp,vbp = ', &
             ubnest(itest,jtest,lslot), &
             vbnest(itest,jtest,lslot), &
             pbnest(itest,jtest,lslot), &
             ubpnst(itest,jtest,lslot), &
             vbpnst(itest,jtest,lslot)
          write(lp,'(i5,i4,a,1p5e13.5)') &
             itest+1+i0,jtest+j0,' rdbaro: ub,vb,pb,ubp,vbp = ', &
             ubnest(itest+1,jtest,lslot), &
             vbnest(itest+1,jtest,lslot), &
             pbnest(itest+1,jtest,lslot), &
             ubpnst(itest+1,jtest,lslot), &
             vbpnst(itest+1,jtest,lslot)
          write(lp,'(i5,i4,a,1p5e13.5)') &
             itest+i0,jtest+1+j0,' rdbaro: ub,vb,pb,ubp,vbp = ', &
             ubnest(itest,jtest+1,lslot), &
             vbnest(itest,jtest+1,lslot), &
             pbnest(itest,jtest+1,lslot), &
             ubpnst(itest,jtest+1,lslot), &
             vbpnst(itest,jtest+1,lslot)
        endif
        call xcsync(flush_lp)
      endif
      return
      end
!
!
      subroutine rdnest(dtime)
      use mod_xc           ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za           ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
!
! --- 3-d nesting archive input processing
!
! --- filenames  ./nest/arch[vm].????_???_??.[ab]
! ---            ./nest/rmu.[ab]
!
! --- I/O and array I/O unit 915 is used for rmun[pv] only (not reserved).
! --- I/O and array I/O unit 920 is reserved for the entire run.
!
! --- all input fields must be defined at all grid points
!
      logical   larchm,lmonth
      save      larchm,lmonth
      integer   iarch
      save      iarch
      real*8    dnestf,dnesti,dtimei,dtime0,dtime1
      save      dnestf,dnesti,dtimei,dtime0,dtime1
!
      integer   iyr,mon,idy,ihr
      integer   i,ios,j,k
      character preambl(5)*79,cline*80
      real*8    wdaymn(-1:2)
!
! --- wn0 negative on first call only.
      if     (wn0.lt.-1.0) then
!
! ---   initialize nesting fields
!
        if     (mnproc.eq.1) then
        write (lp,*) ' now initializing 3-d nesting fields ...'
        endif !1st tile
        call xcsync(flush_lp)
!
        call zaiopf('nest/rmu.a', 'old', 915)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+915,file='nest/rmu.b', &
              status='old',action='read')
        endif !1st tile
        call zagetc(cline,ios, uoff+915)  !1st line of the header on all tiles
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        rewind uoff+915
        read  (uoff+915,'(a79)') preambl
        endif !1st tile
        call preambl_print(preambl)
        if     (cline.eq.'Relaxation Masks') then  !two masks
          call rdmonth(rmunp, 915)
          call rdmonth(rmunv, 915)
          call xctilr(rmunp,1,1, nbdy,nbdy, halo_ps)
          call xctilr(rmunv,1,1, nbdy,nbdy, halo_ps)
        else !one mask
          call rdmonth(rmunp, 915)
          call xctilr(rmunp,1,1, nbdy,nbdy, halo_ps)
          rmunv(:,:) = rmunp(:,:)
        endif !1 or 2 masks
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        close (unit=uoff+915)
        endif !1st tile
        call zaiocl(915)
!
! ---   rmunv[uv] may get masked in rdnest
!
        do j= 1,jj
          do i= 1,ii
            rmunvu(i,j) = max(rmunv(i,j),rmunv(i-1,j))
            rmunvv(i,j) = max(rmunv(i,j),rmunv(i,j-1))
          enddo
        enddo
!
        dnestf = abs(nestfq)
        if     (dnestf.lt.1.0) then
          dnestf = (baclin/86400.0d0)* &
                   max(1,nint((86400.0d0*dnestf)/baclin))
        endif
!
        larchm = nestfq.lt.0.0  !mean archives spaning -nestfq days
        lmonth = yrflag.ge.2 .and. abs(nestfq+30.5).lt.0.1  !monthly mean archives
        if     (lmonth) then
          call fordate(dtime,yrflag, iyr,mon,idy,ihr)
          if     (mon.eq.1) then
            call datefor(wdaymn(-1), iyr-1,   12,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,      1,1,0, yrflag)
            call datefor(wdaymn( 1), iyr,      2,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,      3,1,0, yrflag)
          elseif (mon.eq.11) then
            call datefor(wdaymn(-1), iyr,     10,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,     11,1,0, yrflag)
            call datefor(wdaymn( 1), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    1,1,0, yrflag)
          elseif (mon.eq.12) then
            call datefor(wdaymn(-1), iyr,     11,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 1), iyr+1,    1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    2,1,0, yrflag)
          else
            call datefor(wdaymn(-1), iyr,  mon-1,1,0, yrflag)
            call datefor(wdaymn( 0), iyr,  mon,  1,0, yrflag)
            call datefor(wdaymn( 1), iyr,  mon+1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,  mon+2,1,0, yrflag)
          endif
          if     (dtime.lt.0.5d0*(wdaymn(0)+wdaymn(1))) then
            dtime0 = 0.5d0*(wdaymn(-1)+wdaymn( 0))
            dtime1 = 0.5d0*(wdaymn( 0)+wdaymn( 1))
          else
            dtime0 = 0.5d0*(wdaymn( 0)+wdaymn( 1))
            dtime1 = 0.5d0*(wdaymn( 1)+wdaymn( 2))
          endif
          ln0    = 1
          call rdnest_in(dtime0,larchm,1)
          ln1    = 2
          call rdnest_in(dtime1,larchm,2)
        else
          if     (larchm) then
            dnesti = 0.5d0*dnestf
          else
            dnesti = 0.0
          endif
          dtimei = int((dtime-dnesti)/dnestf)*dnestf + dnesti
!
          dtime0 = dtimei
          ln0    = 1
          call rdnest_in(dtime0,larchm,1)
!
          iarch  = 1
          dtime1 = dtimei + dnestf
          ln1    = 2
          call rdnest_in(dtime1,larchm,2)
        endif
!
        if     (mnproc.eq.1) then
        write (lp,*)
        write (lp,*) ' dtime,dtime0,dtime1 = ',dtime,dtime0,dtime1
        write (lp,*)
        write (lp,*) ' ...finished initializing 3-d nesting fields'
        endif !1st tile
        call xcsync(flush_lp)
      endif  ! initialization
!
      if     (dtime.gt.dtime1) then
!
! ---   get the next set of fields.
        do k= 1,kk
          do j= 1,jj
            do i= 1,ii
              tnest(i,j,k,1) = tnest(i,j,k,2)
              snest(i,j,k,1) = snest(i,j,k,2)
              pnest(i,j,k,1) = pnest(i,j,k,2)
              unest(i,j,k,1) = unest(i,j,k,2)
              vnest(i,j,k,1) = vnest(i,j,k,2)
            enddo
          enddo
        enddo
        un1min(1) = un1min(2)
        un1max(1) = un1max(2)
        vn1min(1) = vn1min(2)
        vn1max(1) = vn1max(2)
        if     (lmonth) then
          dtime0 = dtime1
          call fordate(dtime1,yrflag, iyr,mon,idy,ihr)
          if     (mon.eq.11) then
            call datefor(wdaymn( 1), iyr,     12,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    1,1,0, yrflag)
          elseif (mon.eq.12) then
            call datefor(wdaymn( 1), iyr+1,    1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr+1,    2,1,0, yrflag)
          else
            call datefor(wdaymn( 1), iyr,  mon+1,1,0, yrflag)
            call datefor(wdaymn( 2), iyr,  mon+2,1,0, yrflag)
          endif
          dtime1 = 0.5d0*(wdaymn( 1)+wdaymn( 2))
          call rdnest_in(dtime1,larchm,2)
        else
          dtime0 = dtime1
          iarch  = iarch + 1
          dtime1 = dtimei + dnestf*iarch
          call rdnest_in(dtime1,larchm,2)
        endif
!
!           if     (mnproc.eq.1) then
!           write(lp,*) ' exit rdnest_in - ',dtime,dtime0,dtime1
!           endif !1st tile
!           call xcsync(flush_lp)
      endif  ! next set of fields.
!
! --- linear interpolation in time.
      wn0 = (dtime1-dtime)/(dtime1-dtime0)
      wn1 = 1.0 - wn0
!           if     (mnproc.eq.1) then
!           write(lp,*) 'rdnest - dtime,wn0,wn1 = ',dtime,wn0,wn1
!           endif !1st tile
!           call xcsync(flush_lp)
      return
      end
!
!
      subroutine rdnest_in(dtime,larchm,lslot)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real*8    dtime
      integer   lslot
      logical   larchm
!
! --- input 3-d nesting fields from archive on model day dtime.
! --- filenames  nest/arch[vm].????_???_??.[ab]
! --- I/O and array I/O unit 920 is reserved for the entire run.
!
      logical, parameter :: ldebug_rdnest = .false.  !usually .false.
#if defined(RDNEST_MASK)
      logical, parameter :: lmask_rdnest  = .true.   !set by a CPP macro
                                                     !mask velocity outliers
#else
      logical, parameter :: lmask_rdnest  = .false.  !usually .false.
                                                     !mask velocity outliers
#endif
!
      character flnm*22, cline*80, cvarin*6, cfield*8
      integer   i,idmtst,ios,j,jdmtst,k,layer
      integer   iyear,iday,ihour,nucnt,nvcnt
      logical   meanar,nodens
      real      sumn(2)
!
      call forday(dtime, yrflag, iyear,iday,ihour)
!
      if     (larchm) then
        write(flnm,'("nest/archm.",i4.4,"_",i3.3,"_",i2.2)') &
                                   iyear,iday,ihour
      else
        write(flnm,'("nest/archv.",i4.4,"_",i3.3,"_",i2.2)') &
                                   iyear,iday,ihour
      endif
!
      if     (mnproc.eq.1) then
      write (lp,"(a,a,a,f12.5,a)") 'rdnest_in: ',flnm," (",dtime,")"
      endif !1st tile
      call xcsync(flush_lp)
!
      call zaiopf(flnm//'.a','old', 920)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+920,file=flnm//'.b',form='formatted', &
              status='old',action='read')
!
        read(uoff+920,'(a)') cline
        read(uoff+920,'(a)') cline
        read(uoff+920,'(a)') cline
        read(uoff+920,'(a)') cline
!
        read(uoff+920,'(a)') cline
        read(uoff+920,'(a)') cline
        read(uoff+920,'(a)') cline
      endif !1st tile
!
      call zagetc(cline,ios, uoff+920)
      read(cline,*) idmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',idmtst
!     endif !1st tile
      if (cvarin.ne.'idm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdnest_in - input ',cvarin, &
                              ' but should be idm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(rdnest_in)')
               stop '(rdnest_in)'
      endif
      call zagetc(cline,ios, uoff+920)
      read(cline,*) jdmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',jdmtst
!     endif !1st tile
      if (cvarin.ne.'jdm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdnest_in - input ',cvarin, &
                              ' but should be jdm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(rdnest_in)')
               stop '(rdnest_in)'
      endif
!
      if (idmtst.ne.itdm .or. jdmtst.ne.jtdm) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in rdnest_in - input idm,jdm', &
                              ' not consistent with parameters'
        write(lp,*) 'idm,jdm = ',itdm,  jtdm,  '  (dimensions.h)'
        write(lp,*) 'idm,jdm = ',idmtst,jdmtst,'  (input)'
        write(lp,*)
        endif !1st tile
        call xcstop('(rdnest_in)')
               stop '(rdnest_in)'
      endif
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        read (uoff+920,*)
      endif
!
! --- skip surface fields.
!
      call rd_archive(util1, cfield,layer, 920)  ! montg1 (discarded)
      if     (cfield.ne.'montg1  ') then
        if     (mnproc.eq.1) then
        write(lp,'(/ a / a,a /)') cfield, &
               'error in rdnest_in - expected ','montg1  '
        endif !1st tile
        call xcstop('(rdbaro_in)')
               stop '(rdbaro_in)'
      endif
      nodens = layer.ne.0  !new or original archive type
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        read (uoff+920,*)
      endif
      call zaiosk(920)                  !srfhgt
      call zagetc(cline,ios, uoff+920)  !steric or surflx
      call zaiosk(920)
      if     (cline(1:8).eq.'steric  ') then  !surflx
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+920,*)
        endif
        call zaiosk(920)
      endif
      if     (nodens) then
        call zagetc(cline,ios, uoff+920)  !wtrflx or salflx
        call zaiosk(920)
        if     (cline(1:8).eq.'wtrflx  ') then
          do i= 1,3 !salflx,dpbl,dpmixl
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              read (uoff+920,*)
            endif
            call zaiosk(920)
          enddo
        else  !.eq.salflx
          do i= 1,2 !dpbl,dpmixl
            if     (mnproc.eq.1) then  ! .b file from 1st tile only
              read (uoff+920,*)
            endif
            call zaiosk(920)
          enddo
        endif !wtrflx/salflx
      else
        do i= 1,8 !salflx,dpbl,dpmixl,tmix,smix,thmix,umix,vmix
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920)
        enddo
      endif !nodens:else
      call zagetc(cline,ios, uoff+920)  !kemix or covice or u_btrop
      meanar = cline(1:8).eq.'kemix   '
      if     (meanar) then
        call zaiosk(920)  !skip kemix
        call rd_archive(util1, cfield,layer, 920) !covice or u_btrop
        if     (cfield.eq.'covice  ') then
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920) !skip thkice
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920) !skip temice
          call rd_archive(util1, cfield,layer, 920) !u_btrop
        endif
        call rd_archive(util2, cfield,layer, 920) !v_btrop
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+920,*)
        endif
        call zaiosk(920)  !skip kebtrop
      else !standard archive file
        if     (cline(1:8).eq.'covice  ') then
          call zaiosk(920) !skip covice
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920) !skip thkice
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920) !skip temice
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
        endif
        call zaiosk(920)  !skip u_btrop
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+920,*)
        endif
        call zaiosk(920)  !skip v_btrop
      endif !meanar:else
!
! --- 3-d fields.
!
      do k=1,kk
        if     (k.eq.1) then
          call rd_archive_rng(unest(1-nbdy,1-nbdy,k,lslot), &
                              un1min(lslot),un1max(lslot), &
                              cfield,layer, 920)
          un1min(lslot) = un1min(lslot) - 0.1
          un1max(lslot) = un1max(lslot) + 0.1
        else
          call rd_archive(    unest(1-nbdy,1-nbdy,k,lslot), &
                                           cfield,layer, 920)
        endif !k==1:else
        if     (cfield.ne.'u-vel.  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in rdnest_in - expected ','u-vel.  '
          endif !1st tile
          call xcstop('(rdnest_in)')
                 stop '(rdnest_in)'
        endif
        if     (k.eq.1) then
          call rd_archive_rng(vnest(1-nbdy,1-nbdy,k,lslot), &
                              vn1min(lslot),vn1max(lslot), &
                              cfield,layer, 920)
          vn1min(lslot) = vn1min(lslot) - 0.1
          vn1max(lslot) = vn1max(lslot) + 0.1
        else
          call rd_archive(    vnest(1-nbdy,1-nbdy,k,lslot), &
                                           cfield,layer, 920)
        endif
        if     (cfield.ne.'v-vel.  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in rdnest_in - expected ','v-vel.  '
          endif !1st tile
          call xcstop('(rdnest_in)')
                 stop '(rdnest_in)'
        endif
        if     (meanar) then
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920)  !skip k.e.
        endif
        if     (k.ne.kk) then
          call rd_archive(pnest(1-nbdy,1-nbdy,k+1,lslot), &
                          cfield,layer, 920)
          if     (cfield.ne.'thknss  ') then
            if     (mnproc.eq.1) then
            write(lp,'(/ a / a,a /)') cfield, &
                   'error in rdnest_in - expected ','thknss  '
            endif !1st tile
            call xcstop('(rdnest_in)')
                   stop '(rdnest_in)'
          endif
        else
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920)
        endif
        call rd_archive(tnest(1-nbdy,1-nbdy,k,lslot), cfield,layer, 920)
        if     (cfield.ne.'temp    ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in rdnest_in - expected ','temp    '
          endif !1st tile
          call xcstop('(rdnest_in)')
                 stop '(rdnest_in)'
        endif
        call rd_archive(snest(1-nbdy,1-nbdy,k,lslot), cfield,layer, 920)
        if     (cfield.ne.'salin   ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in rdnest_in - expected ','salin   '
          endif !1st tile
          call xcstop('(rdnest_in)')
                 stop '(rdnest_in)'
        endif
        if     (.not. nodens) then
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+920,*)
          endif
          call zaiosk(920)  !skip density
        endif !.not.nodens
      enddo
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close( unit=uoff+920)
      endif
      call zaiocl(920)

      if     (meanar) then
        call xctilr(pnest(1-nbdy,1-nbdy,1,lslot),1,kk, 1,1, halo_ps)
      endif
!
      nucnt = 0
      nvcnt = 0
!$OMP PARALLEL DO PRIVATE(j,i,k) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
        if     (lmask_rdnest) then
! ---     mask rmunv[uv] where velocities exceed near-surface range
          do i=1,ii
            if     (iu(i,j).eq.1) then
              rmunvu(i,j) = max(rmunv(i,j),rmunv(i-1,j))
              if    (rmunvu(i,j).gt.0.0) then
                do k= 2,kk
                  if     (unest(i,j,k,1).lt.un1min(1) .or. &
                          unest(i,j,k,1).gt.un1max(1) .or. &
                          unest(i,j,k,2).lt.un1min(2) .or. &
                          unest(i,j,k,2).gt.un1max(2)     ) then
                    rmunvu(i,j) = 0.0  !mask
                    nucnt = nucnt + 1
                    exit !k
                  endif
                enddo !k
              endif !rmunvu
            endif !iu
            if     (iv(i,j).eq.1) then
              rmunvv(i,j) = max(rmunv(i,j),rmunv(i,j-1))
              if    (rmunvv(i,j).gt.0.0) then
                do k= 2,kk
                  if     (vnest(i,j,k,1).lt.vn1min(1) .or. &
                          vnest(i,j,k,1).gt.vn1max(1) .or. &
                          vnest(i,j,k,2).lt.vn1min(2) .or. &
                          vnest(i,j,k,2).gt.vn1max(2)     ) then
                    rmunvv(i,j) = 0.0  !mask
                    nvcnt = nvcnt + 1
                    exit !k
                  endif
                enddo !k
              endif !rmunvv
            endif !iv
          enddo !i
        endif !lmask_rdnest
!
        if     (meanar) then  !mean archive
!         for thin layers, take baroclinic velocity from above
!         otherwise, convert from total to baroclinic velocity
          do k= 1,kk
            do i=1,ii
              if     (iu(i,j).eq.1) then
                if     (min(pnest(i,  j,k,lslot), &
                            pnest(i-1,j,k,lslot) ).lt.tencm) then
                  unest(i,j,k,lslot) = unest(i,j,max(1,k-1),lslot)
                else
                  unest(i,j,k,lslot) = unest(i,j,k,lslot) - util1(i,j)
                endif !thin layer:else
              endif !iu
              if     (iv(i,j).eq.1) then
                if     (min(pnest(i,j,  k,lslot), &
                            pnest(i,j-1,k,lslot) ).lt.tencm) then
                  vnest(i,j,k,lslot) = vnest(i,j,max(1,k-1),lslot)
                else
                  vnest(i,j,k,lslot) = vnest(i,j,k,lslot) - util2(i,j)
                endif !thin layer:else
              endif !iv
            enddo !i
          enddo !k
        endif !meanar
!       convert from layer thickness to interface depth (pressure)
        do i=1,ii
          pnest(i,j,1,lslot) = 0.0
          do k= 3,kk
            pnest(i,j,k,lslot) = pnest(i,j,k,  lslot) + &
                                 pnest(i,j,k-1,lslot)
          enddo !k
        enddo !i
      enddo !j
!
      if     (lmask_rdnest) then
        sumn(1) = nucnt
        sumn(2) = nvcnt
        call xcsumr(sumn, 1)
        if     (sumn(1)+sumn(2).gt.0.0 .and. mnproc.eq.1) then
          write(lp,'(f12.5,a,2i9)') &
            dtime,'  nest mask counts =',int(sumn(1)),int(sumn(2))
        endif 
      endif !lmask_rdnest
!
      if     (ldebug_rdnest .and. ittest.ne.-1 .and. jttest.ne.-1) then
        call xcsync(flush_lp)
        if     (i0.lt.ittest .and. i0+ii.ge.ittest .and. &
                j0.lt.jttest .and. j0+jj.ge.jttest      ) then
 103      format(i8,i5,i4,1x,a,a/ &
                 (i8,5x,i4,1x,a,a,2f7.3,2f7.3,f8.4,f9.3,f9.2))
          write(lp,103) &
             nstep,itest+i0,jtest+j0,'rdnest', &
             ':   utot   vtot   temp   saln    dens    thkns     dpth', &
            (nstep,k,                'rdnest',':', &
             unest(itest,jtest,k,lslot)+ubnest(itest,jtest,lslot), &
             vnest(itest,jtest,k,lslot)+vbnest(itest,jtest,lslot), &
             tnest(itest,jtest,k,lslot), &
             snest(itest,jtest,k,lslot), &
             0.0, &
             (pnest(itest,jtest,k+1,lslot)- &
              pnest(itest,jtest,k,  lslot) )*qonem, &
             pnest(itest,jtest,k+1,lslot)*qonem, &
             k=1,kk-1), &
            (nstep,k,                'rdnest',':', &
             unest(itest,jtest,k,lslot)+ubnest(itest,jtest,lslot), &
             vnest(itest,jtest,k,lslot)+vbnest(itest,jtest,lslot), &
             tnest(itest,jtest,k,lslot), &
             snest(itest,jtest,k,lslot), &
             0.0, &
             depths(i,j)-pnest(itest,jtest,k,lslot)*qonem, &
             depths(i,j), &
             k=kk,kk)
        endif
        call xcsync(flush_lp)
      endif
!
      return
      end
!
!
      subroutine rd_archive(field, cfield,layer, iunit)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      character cfield*8
      integer   layer,iunit
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
                field
!
! --- read a single archive array field from unit iunit.
!
      integer   i,ios,nnstep
      real      hmina,hminb,hmaxa,hmaxb,timein,thet
      character cline*80
!
      call zagetc(cline,ios, uoff+iunit)
      if     (ios.ne.0) then
        if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in rd_archive - hit end of input'
          write(lp,*) 'iunit,ios = ',iunit,ios
          write(lp,*)
        endif !1st tile
        call xcstop('(rd_archive)')
               stop '(rd_archive)'
      endif
!     if     (mnproc.eq.1) then
!     write(lp,'(a)') cline
!     endif !1st tile
!
      cfield = cline(1:8)
!
      i = index(cline,'=')
      read(cline(i+1:),*) nnstep,timein,layer,thet,hminb,hmaxb
!
      if     (hminb.eq.hmaxb) then  !constant field
        field(:,:) = hminb
        call zaiosk(iunit)
      else
        call zaiord(field,ip,.false., hmina,hmaxa, &
                    iunit)
!
        if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or. &
                abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') &
            'error - .a and .b files not consistent:', &
            '.a,.b min = ',hmina,hminb,hmina-hminb, &
            '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
          endif !1st tile
!nostop   call xcstop('(rd_archive)')
!nostop          stop '(rd_archive)'
        endif
      endif
      return
      end
!
!
      subroutine rd_archive_rng(field,fmin,fmax, cfield,layer, iunit)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      character cfield*8
      integer   layer,iunit
      real      fmin,fmax
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
                field
!
! --- read a single archive array field from unit iunit.
! --- return its range in fmin,fmax
!
      integer   i,ios,nnstep
      real      hmina,hminb,hmaxa,hmaxb,timein,thet
      character cline*80
!
      call zagetc(cline,ios, uoff+iunit)
      if     (ios.ne.0) then
        if     (mnproc.eq.1) then
          write(lp,*)
          write(lp,*) 'error in rd_archive - hit end of input'
          write(lp,*) 'iunit,ios = ',iunit,ios
          write(lp,*)
        endif !1st tile
        call xcstop('(rd_archive)')
               stop '(rd_archive)'
      endif
!     if     (mnproc.eq.1) then
!     write(lp,'(a)') cline
!     endif !1st tile
!
      cfield = cline(1:8)
!
      i = index(cline,'=')
      read(cline(i+1:),*) nnstep,timein,layer,thet,hminb,hmaxb
!
      if     (hminb.eq.hmaxb) then  !constant field
        field(:,:) = hminb
        call zaiosk(iunit)
      else
        call zaiord(field,ip,.false., hmina,hmaxa, &
                    iunit)
!
        if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or. &
                abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') &
            'error - .a and .b files not consistent:', &
            '.a,.b min = ',hmina,hminb,hmina-hminb, &
            '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
          endif !1st tile
!nostop   call xcstop('(rd_archive)')
!nostop          stop '(rd_archive)'
        endif
      endif
      fmin = hminb
      fmax = hmaxb
      return
      end
!
!
      subroutine str2spd(wspd, tx,ty)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      implicit none
!
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
              wspd, tx,ty
!
! --- calculate wind speed from wind stress
!
! --- speed-dependent scale factor from stress to speed is based
! --- on the Kara (neutral) wind-speed dependent drag coefficient 
!
      integer   i,j
      real      strspd,wndstr
!
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
      do j= 1-nbdy,jj+nbdy
        do i= 1-nbdy,ii+nbdy
          wndstr = sqrt( tx(i,j)**2 + ty(i,j)**2 )
          if     (wndstr.LE.0.7711) THEN
            strspd = 1.0/(1.22*(((3.236E-3 *wndstr - &
                                  5.230E-3)*wndstr + &
                                  3.218E-3)*wndstr + &
                                  0.926E-3)       )
          else
            strspd = 1.0/(1.22*(((0.007E-3 *wndstr - &
                                  0.092E-3)*wndstr + &
                                  0.485E-3)*wndstr + &
                                  1.461E-3)       )
          endif
          wspd(I,J) = sqrt( strspd*wndstr )
        enddo
      enddo
      return
      end
!
!
!> Revision history:
!>
!> Mar. 1995 - added logical variable 'windf'
!> Oct. 1997 - made necessary changes to reduce time dimension from 12 to 4
!> Oct. 1999 - added code to read and store shortwave heat flux used for
!>             penetrating shortwave radiation
!> Jan. 2000 - removed all conversion factors (apply before input)
!> Jan. 2000 - removed biasrd and biaspc      (apply before input)
!> May. 2000 - conversion to SI units, positive flux into ocean
!> Aug. 2000 - added option for high frequency atmospheric forcing
!> Jan. 2001 - converted from pakk to array input file type
!> July 2001 - added skmonth and support for relaxs (surface only relax)
!> July 2001 - added rdopen and rdbaro
!> Aug. 2001 - added constant field logic (skip the array input)
!> Mar. 2003 - added surtmp and seatmp
!> Mar. 2005 - added wndflg==3 and str2spd
!> Aug. 2005 - added tracer climatology
!> Nov. 2006 - [uv]pnst now from interpolation of transports
!> Mar. 2010 - support for daily mean nesting archive files
!> Mar. 2010 - added diwbot
!> Apr. 2010 - added sssrmx
!> Apr  2010 - added diwqh0 and removed diwlat
!> Nov. 2010 - added wndrep (yrflag==2) for one or two year forcing repeat
!> Apr. 2011 - added cbarp
!> July 2011 - fixed a cbarp input bug when cbar<0
!> July 2011 - added forfuns for salfac
!> Sep. 2011 - added cbp
!> Nov. 2012 - added wndflg=4 for reading 10m wind components
!> Nov. 2012 - added iftaux,oftauy, primarily for wndflg=4
!> Jan. 2013 - replaced dragrh with drgten
!> July 2013 - added diws and diwm
!> Oct. 2013 - added jerlv0=-1 and forfunc
!> Nov. 2013 - added wndflg=5 (also) for reading 10m wind components
!> Nov. 2013 - added lwflag.eq.-1 for radflx=Qlwdn, swflx=Qswdn
!> Jan. 2014 - added mslprf and mslprs and forfunhp
!> Jan. 2014 - added natm logic
!> Jan. 2014 - modified natm and pwall logic to avoid gfortran warnings
!> May  2014 - added forfunhz
!> Oct. 2014 - flxflg==6 requires mslprs input
!> Aug. 2015 - flxflg==3 does not read airtmp and vapmix
!> June 2016 - nestfq and bnstfq -30.5 now for actual monthly means
!> July 2017 - always define rivers array
!> July 2017 - datefor and lmonth now also for yrflag=2.
!> July 2017 - for montyp=3,4 veldf[24]u is actually veldf[24]p
!> Sep. 2017 - bugfix for mslprs skip to 1st day
!> Sep. 2017 - added stoc_[tsuv]
!> Jan. 2018 - overlay thkdf2[uv] on thkdf4[uv]
!> Apr. 2018 - read in spchum when flxflg and wndflg are 5.
!> Nov. 2018 - added sefold (rmus)
!> Nov. 2018 - input archive files can now contain wtrflx
!> Dec. 2018 - cycle when no ustar
!> Dec  2018 - add yrflag=4 for 365 days no-leap calendar (CESM)
!> Feb. 2019 - add sshflg=2 for steric Montg. Potential and montg_c in rdbaro_in
!> May  2019 - add yrflag=4 logic to datefor and fordate
!> Oct  2019 - added lbmont
!> Oct  2019 - bugfix in datefor for yrflag=2,4
!> Oct  2019 - optionally mask nest velocities if outside near-surface range
!> Oct  2019 - added a CPP macro to set lmask_rdnest
