#include "BULK_FORCE_OPTIONS.h"

C     !ROUTINE: BULKF_FIELDS_LOAD
C     !INTERFACE:
      SUBROUTINE BULKF_FIELDS_LOAD( myTime, myIter, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE BULKF_FIELDS_LOAD
C     | o Control reading of fields from external source.
C     *==========================================================*
C     | Bulk formula External source field loading routine.
C     | This routine is called every time we want to
C     | load a a set of external fields. The routine decides
C     | which fields to load and then reads them in.
C     | This routine needs to be customised for particular
C     | experiments.
C     | Notes
C     | =====
C     | Two-dimensional and three-dimensional I/O are handled in
C     | the following way under MITgcmUV. A master thread
C     | performs I/O using system calls. This threads reads data
C     | into a temporary buffer. At present the buffer is loaded
C     | with the entire model domain. This is probably OK for now
C     | Each thread then copies data from the buffer to the
C     | region of the proper array it is responsible for.
C     | =====
C     | Conversion of flux fields are described in FFIELDS.h
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"
#include "BULKF.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: Thread no. that called this routine.
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_BULK_FORCE

C     !LOCAL VARIABLES:
C     === Local arrays in common block ===
C     bulkfRec   :: time-record currently loaded (in temp arrays *[1])
C     tair[01]   :: Temp. for air temperature
C     qair[01]   :: Temp. for air specific humidity
C     solar[01]  :: Temp. for incoming solar radition
C     flwdwn[01] :: Temp. for downward longwave flux
C     cloud[01]  :: Temp for cloud
C     wspeed[01] :: Temp. for wind speed
C     uwind[01]  :: Temp. for zonal wind speed
C     vwind[01]  :: Temp. for meridional wind speed
C     rain[01]   :: Temp. for rain
C     runoff[01] :: Temp. for runoff
C     snow[01]   :: Temp  for snow
C     thAir[01]  :: Temp for Air potential temp. in the BL [K]
c     qnetch[01] :: Temp for qnet (cheating)
C     empch[01]  :: Temp for empmr (cheating)
C     [01]       :: End points for interpolation

      COMMON /BULKF_FIELDS_LOAD_I/ bulkfRec
      COMMON /BULKF_FIELDS_LOAD_RS/
     &                 tair0, tair1, qair0, qair1,
     &                 solar0, solar1, flwdwn0, flwdwn1,
     &                 cloud0, cloud1, wspeed0, wspeed1,
     &                 uwind0, uwind1, vwind0, vwind1,
     &                 rain0, rain1, runoff0, runoff1,
     &                 snow0, snow1,
#ifdef ALLOW_FORMULA_AIM
     &                 thAir0, thAir1,
#endif
     &                 qnetch0, qnetch1, empch0, empch1

      INTEGER bulkfRec(nSx,nSy)
      _RS  tair0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  tair1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  qair0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  qair1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  solar0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  solar1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  flwdwn0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  flwdwn1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  cloud0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  cloud1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  wspeed0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  wspeed1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  uwind0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  uwind1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  vwind0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  vwind1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  rain0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  rain1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  runoff0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  runoff1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  snow0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  snow1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef ALLOW_FORMULA_AIM
      _RS  thAir0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  thAir1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#endif
      _RS  qnetch0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  qnetch1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  empch0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  empch1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C     === Local variables ===
C     aWght, bWght :: Interpolation weights
      _RL  aWght,bWght
      INTEGER intimeP, intime0, intime1
      INTEGER bi, bj, i, j
CEOP

      IF ( periodicExternalForcing ) THEN

C First call requires that we initialize everything to zero for safety
      IF ( myIter .EQ. nIter0 ) THEN
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
           bulkfRec(bi,bj) = 0
           DO j = 1-Oly,sNy+Oly
            DO i = 1-Olx,sNx+Olx
              tair0(i,j,bi,bj)   = 0. _d 0
              tair1(i,j,bi,bj)   = 0. _d 0
              qair0(i,j,bi,bj)   = 0. _d 0
              qair1(i,j,bi,bj)   = 0. _d 0
              rain0(i,j,bi,bj)   = 0. _d 0
              rain1(i,j,bi,bj)   = 0. _d 0
              solar0(i,j,bi,bj)  = 0. _d 0
              solar1(i,j,bi,bj)  = 0. _d 0
              flwdwn0(i,j,bi,bj) = 0. _d 0
              flwdwn1(i,j,bi,bj) = 0. _d 0
              cloud0(i,j,bi,bj)  = 0. _d 0
              cloud1(i,j,bi,bj)  = 0. _d 0
              wspeed0(i,j,bi,bj) = 0. _d 0
              wspeed1(i,j,bi,bj) = 0. _d 0
              uwind0(i,j,bi,bj)  = 0. _d 0
              uwind1(i,j,bi,bj)  = 0. _d 0
              vwind0(i,j,bi,bj)  = 0. _d 0
              vwind1(i,j,bi,bj)  = 0. _d 0
              runoff0(i,j,bi,bj) = 0. _d 0
              runoff1(i,j,bi,bj) = 0. _d 0
              snow0(i,j,bi,bj)   = 0. _d 0
              snow1(i,j,bi,bj)   = 0. _d 0
#ifdef ALLOW_FORMULA_AIM
              thAir0(i,j,bi,bj)  = 0. _d 0
              thAir1(i,j,bi,bj)  = 0. _d 0
#endif
              qnetch0(i,j,bi,bj) = 0. _d 0
              qnetch1(i,j,bi,bj) = 0. _d 0
              empch0(i,j,bi,bj)  = 0. _d 0
              empch1(i,j,bi,bj)  = 0. _d 0
            ENDDO
           ENDDO
         ENDDO
        ENDDO
      ENDIF

C--   Now calculate whether it is time to update the forcing arrays
      CALL GET_PERIODIC_INTERVAL(
     O                  intimeP, intime0, intime1, bWght, aWght,
     I                  externForcingCycle, externForcingPeriod,
     I                  deltaTclock, myTime, myThid )

      bi = myBxLo(myThid)
      bj = myByLo(myThid)
#ifdef ALLOW_DEBUG
      IF ( debugLevel.GE.debLevB ) THEN
        _BEGIN_MASTER(myThid)
        WRITE(standardMessageUnit,'(A,I10,A,4I5,A,2F14.10)')
     &   ' BULKF_FIELDS_LOAD,', myIter,
     &   ' : iP,iLd,i0,i1=', intimeP,bulkfRec(bi,bj), intime0,intime1,
     &   ' ; Wght=', bWght, aWght
        _END_MASTER(myThid)
      ENDIF
#endif /* ALLOW_DEBUG */

C-    Make no assumption on sequence of calls to BULKF_FIELDS_LOAD ;
C     This is the correct formulation (works in Adjoint run).
      IF ( intime1.NE.bulkfRec(bi,bj) ) THEN

       _BARRIER

C--   If the above condition is met then we need to read in
C     data for the period ahead and the period behind myTime.
       IF ( debugLevel.GE.debLevZero ) THEN
         _BEGIN_MASTER(myThid)
        WRITE(standardMessageUnit,'(A,I10,A,2(2I5,A))')
     &   ' BULKF_FIELDS_LOAD, it=', myIter,
     &   ' : Reading new data, i0,i1=', intime0, intime1,
     &    ' (prev=', intimeP, bulkfRec(bi,bj), ' )'
        _END_MASTER(myThid)
       ENDIF

      IF ( AirTempFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( AirTempFile,tair0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( AirTempFile,tair1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( AirHumidityFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( AirhumidityFile,qair0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( AirhumidityFile,qair1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( SolarFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( SolarFile,solar0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( SolarFile,solar1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( LongwaveFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( LongwaveFile,flwdwn0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( LongwaveFile,flwdwn1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( CloudFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( CloudFile,cloud0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( CloudFile,cloud1,intime1,
     &       myIter,myThid )
      ENDIF

      IF ( WSpeedFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( WSpeedFile,wspeed0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( WSpeedFile,wspeed1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( UwindFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( UWindFile,uwind0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( UWindFile,uwind1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( VwindFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( VWindFile,vwind0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( VWindFile,vwind1,intime1,
     &       myIter,myThid )
      ENDIF

      IF ( RainFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( RainFile,rain0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( RainFile,rain1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( RunoffFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( RunoffFile,runoff0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( RunoffFile,runoff1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( SnowFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( SnowFile,snow0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( SnowFile,snow1,intime1,
     &       myIter,myThid )
      ENDIF

#ifdef ALLOW_FORMULA_AIM
      IF ( airPotTempFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( airPotTempFile, thAir0, intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( airPotTempFile, thAir1, intime1,
     &       myIter,myThid )
      ENDIF
#endif

      IF ( QnetFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( QnetFile,qnetch0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( QnetFile,qnetch1,intime1,
     &       myIter,myThid )
      ENDIF
      IF ( EmPFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( EmpFile,empch0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( EmpFile,empch1,intime1,
     &       myIter,myThid )
c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
        _BARRIER
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            empch0(i,j,bi,bj) = empch0(i,j,bi,bj)*rhoConstFresh
            empch1(i,j,bi,bj) = empch1(i,j,bi,bj)*rhoConstFresh
           ENDDO
          ENDDO
         ENDDO
        ENDDO
c      ENDIF
      ENDIF

C note: no need for extra barrier here since EXCH contains threads synchronization
c      _BARRIER

       _EXCH_XY_RS(tair0 , myThid )
       _EXCH_XY_RS(tair1 , myThid )
       _EXCH_XY_RS(qair0 , myThid )
       _EXCH_XY_RS(qair1 , myThid )
       _EXCH_XY_RS(solar0, myThid )
       _EXCH_XY_RS(solar1, myThid )
       _EXCH_XY_RS(flwdwn0, myThid )
       _EXCH_XY_RS(flwdwn1, myThid )
       _EXCH_XY_RS(cloud0, myThid )
       _EXCH_XY_RS(cloud1, myThid )
       _EXCH_XY_RS(wspeed0, myThid )
       _EXCH_XY_RS(wspeed1, myThid )
c      _EXCH_XY_RS(uwind0, myThid )
c      _EXCH_XY_RS(uwind1, myThid )
c      _EXCH_XY_RS(vwind0, myThid )
c      _EXCH_XY_RS(vwind1, myThid )
       CALL EXCH_UV_AGRID_3D_RS( uwind0, vwind0, .TRUE., 1, myThid )
       CALL EXCH_UV_AGRID_3D_RS( uwind1, vwind1, .TRUE., 1, myThid )
       _EXCH_XY_RS(rain0, myThid )
       _EXCH_XY_RS(rain1, myThid )
       _EXCH_XY_RS(runoff0, myThid )
       _EXCH_XY_RS(runoff1, myThid )
       _EXCH_XY_RS(snow0 , myThid )
       _EXCH_XY_RS(snow1 , myThid )
#ifdef ALLOW_FORMULA_AIM
       IF ( useFluxFormula_AIM ) THEN
         _EXCH_XY_RS( thAir0, myThid )
         _EXCH_XY_RS( thAir1, myThid )
       ENDIF
#endif
       _EXCH_XY_RS(qnetch0, myThid )
       _EXCH_XY_RS(qnetch1, myThid )
       _EXCH_XY_RS(empch0, myThid )
       _EXCH_XY_RS(empch1, myThid )

C-    save newly loaded time-record
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
           bulkfRec(bi,bj) = intime1
         ENDDO
        ENDDO

C--   end if-block for loading new time-records
      ENDIF

C--   Interpolate TAIR, QAIR, RAIN, SOLAR
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
Cswdblkf -- QQQQQ check if tair is K or C ------
C        -- dasilva data in C, ncep data in K ---
          TAIR(i,j,bi,bj)   = bWght*tair0(i,j,bi,bj)
     &                       +aWght*tair1(i,j,bi,bj)
c                            +273.15 _d 0
Cswdblkf -- QQQQQ set to kg.kg??? ---
C       -- dasilva data in g, ncep in kg ---
          QAIR(i,j,bi,bj)    =(bWght*qair0(i,j,bi,bj)
     &                        +aWght*qair1(i,j,bi,bj) )
          SOLAR(i,j,bi,bj)   = bWght*solar0(i,j,bi,bj)
     &                        +aWght*solar1(i,j,bi,bj)
          FLWDWN(i,j,bi,bj)  = bWght*flwdwn0(i,j,bi,bj)
     &                        +aWght*flwdwn1(i,j,bi,bj)
          CLOUD(i,j,bi,bj)   = bWght*cloud0(i,j,bi,bj)
     &                        +aWght*cloud1(i,j,bi,bj)
          WSPEED(i,j,bi,bj)  = bWght*wspeed0(i,j,bi,bj)
     &                        +aWght*wspeed1(i,j,bi,bj)
          UWIND(i,j,bi,bj)   = bWght*uwind0(i,j,bi,bj)
     &                        +aWght*uwind1(i,j,bi,bj)
          VWIND(i,j,bi,bj)   = bWght*vwind0(i,j,bi,bj)
     &                        +aWght*vwind1(i,j,bi,bj)
          RAIN(i,j,bi,bj)    = bWght*rain0(i,j,bi,bj)
     &                        +aWght*rain1(i,j,bi,bj)
          RUNOFF(i,j,bi,bj)  = bWght*runoff0(i,j,bi,bj)
     &                        +aWght*runoff1(i,j,bi,bj)
c#ifdef ALLOW_THSICE
c         SNOW(i,j,bi,bj)    = bWght*snow0(i,j,bi,bj)
c    &                        +aWght*snow1(i,j,bi,bj)
c#endif
          Qnetch(i,j,bi,bj)  = bWght*qnetch0(i,j,bi,bj)
     &                        +aWght*qnetch1(i,j,bi,bj)
          EmPch(i,j,bi,bj)   = bWght*empch0(i,j,bi,bj)
     &                        +aWght*empch1(i,j,bi,bj)
         ENDDO
        ENDDO
#ifdef ALLOW_FORMULA_AIM
        IF ( useFluxFormula_AIM ) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            thAir(i,j,bi,bj)  = bWght*thAir0(i,j,bi,bj)
     &                        + aWght*thAir1(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF
#endif
       ENDDO
      ENDDO

C- jmc: not needed since local fields have been exchanged.
c      _EXCH_XY_RL(tair , myThid )

C-- Diagnostics
c      IF (myThid.EQ.1 .AND. myTime.LT.62208000.) THEN
c        write(*,'(a,1p5e12.4,2i6,2e12.4)')
c     &   'time,TAIR,QAIR,RAIN,SOLAR,i0,i1,a,b = ',
c     &   myTime,
c     &   TAIR(1,sNy,1,1),QAIR(1,sNy,1,1),
c     &   RAIN(1,sNy,1,1),SOLAR(1,sNy,1,1),
c     &   intime0,intime1,aWght,bWght
c        write(*,'(a,1p4e12.4,2e12.4)')
c     &   'time,tair0,tair1,TAIR = ',
c     &   myTime,
c     &   tair0(1,sNy,1,1),tair1(1,sNy,1,1),TAIR(1,sNy,1,1),
c     &   aWght,bWght
c      ENDIF

C endif for periodicForcing
      ENDIF

#endif /*ALLOW_BULK_FORCE*/

      RETURN
      END
