C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_fields_load.F,v 1.10 2005/04/06 18:36:22 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
 
C     !ROUTINE: CHEAPAML_FIELDS_LOAD
C     !INTERFACE:
      SUBROUTINE CHEAPAML_FIELDS_LOAD( myTime, myIter, myThid )
C     *==========================================================*
C     | SUBROUTINE CHEAPAML_FIELDS_LOAD                           
C     | o Control reading of fields from external source.         
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
c #include "GRID.h"
c #include "DYNVARS.h"
C #include "BULKF.h"
#ifdef ALLOW_THSICE
#include "THSICE_VARS.h"
#endif
#include "CHEAPAML.h"
 
C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
C     myTime - Simulation time
C     myIter - Simulation timestep number
      INTEGER myThid
      _RL     myTime
      _RL     local
      INTEGER myIter
	INTEGER jg
 
C     !LOCAL VARIABLES:
C     === Local arrays ===
C     trair[01]  :: Relaxation temp. profile for air temperature
C     solar[01]  :: short wave flux
C     uwind[01]  :: zonal wind
C     vwind[01]  :: meridional wind

C     aWght, bWght :: Interpolation weights
      COMMON /BULKFFIELDS/
     &                 trair0,
     &                 trair1
 
      _RS  trair0    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  trair1    (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  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)

      INTEGER bi,bj,i,j,intime0,intime1
      _RL aWght,bWght,rdt
      INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm


      IF ( periodicExternalForcing ) THEN

C First call requires that we initialize everything to zero for safety
      IF ( myIter .EQ. nIter0 ) THEN
       CALL LEF_ZERO( trair0 ,myThid )
       CALL LEF_ZERO( trair1 ,myThid )
       CALL LEF_ZERO( solar0 ,myThid )
       CALL LEF_ZERO( solar1 ,myThid )
       CALL LEF_ZERO( uwind0 ,myThid )
       CALL LEF_ZERO( uwind1 ,myThid )
       CALL LEF_ZERO( vwind0 ,myThid )
       CALL LEF_ZERO( vwind1 ,myThid )
      ENDIF

C Now calculate whether it is time to update the forcing arrays
      rdt=1. _d 0 / deltaTclock
      nForcingPeriods=
     &  int(externForcingCycle/externForcingPeriod+0.5)
      Imytm=int(myTime*rdt+0.5)
      Ifprd=int(externForcingPeriod*rdt+0.5)
      Ifcyc=int(externForcingCycle*rdt+0.5)
      Iftm=mod( Imytm+Ifcyc-Ifprd/2 ,Ifcyc)

      intime0=int(Iftm/Ifprd)
      intime1=mod(intime0+1,nForcingPeriods)
c     aWght=float( Iftm-Ifprd*intime0 )/float( Ifprd )
      aWght=dfloat( Iftm-Ifprd*intime0 )/dfloat( Ifprd )
      bWght=1.-aWght

      intime0=intime0+1
      intime1=intime1+1

      IF (
     &  Iftm-Ifprd*(intime0-1) .EQ. 0
     &  .OR. myIter .EQ. nIter0
     & ) THEN

       _BEGIN_MASTER(myThid)

C      If the above condition is met then we need to read in
C      data for the period ahead and the period behind myTime.
       WRITE(*,*)
     &  'S/R CHEAPAML_FIELDS_LOAD'
      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 ( TRFile .NE. ' '  ) THEN
       CALL READ_REC_XY_RS( TRFile,trair0,intime0,
     &       myIter,myThid )
       CALL READ_REC_XY_RS( TRFile,trair1,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

       _END_MASTER(myThid)
C
       _EXCH_XY_R4(trair0 , myThid )
       _EXCH_XY_R4(solar0 , myThid )
       _EXCH_XY_R4(uwind0 , myThid )
       _EXCH_XY_R4(vwind0 , myThid )
       _EXCH_XY_R4(trair1 , myThid )
       _EXCH_XY_R4(solar1 , myThid )
       _EXCH_XY_R4(uwind1 , myThid )
       _EXCH_XY_R4(vwind1 , myThid )
C
      ENDIF

C--   Interpolate SOLAR, WINDS and LAND TEMPERATURE
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          TR(i,j,bi,bj)   = bWght*trair0(i,j,bi,bj)  
     &                     +aWght*trair1(i,j,bi,bj)   !+273.15
          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)   
          solar(i,j,bi,bj)   = bWght*solar0(i,j,bi,bj)  
     &                     +aWght*solar1(i,j,bi,bj)   
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      ELSE
       IF ( myIter .EQ. nIter0 ) THEN
        IF ( SolarFile .NE. ' '  ) THEN
         CALL READ_FLD_XY_RL( SolarFile,solar,' ',0,myThid )
	ELSE
      DO bj=1,nSy
         DO bi=1,nSx
            DO j=1,sNy
               DO i=1,sNx
                  jG = myYGlobalLo-1+(bj-1)*sNy+j
	local=250.d0-float((jg-1))/float((ny-1))*75.d0
                  Solar(i,j,bi,bj) = local
               ENDDO
            ENDDO
         ENDDO
      ENDDO
	_EXCH_XY_RL(solar, mythid)
        ENDIF
        IF ( TRFile .NE. ' '  ) THEN
         CALL READ_FLD_XY_RL( TRFile,tr,' ',0,myThid )
	ELSE
      DO bj=1,nSy
         DO bi=1,nSx
            DO j=1,sNy
               DO i=1,sNx
c                  jG = myYGlobalLo-1+(bj-1)*sNy+j
	local=solar(i,j,bi,bj)
	local=(2.d0*local/stefan)**(0.25d0)-273.16
                  TR(i,j,bi,bj) = local
               ENDDO
            ENDDO
         ENDDO
      ENDDO
	_EXCH_XY_RL(TR, mythid)
        ENDIF
        IF ( UWindFile .NE. ' '  ) THEN
         CALL READ_FLD_XY_RL( UWindFile,uwind,' ',0,myThid )
	ELSE
      DO bj=1,nSy
         DO bi=1,nSx
            DO j=1,sNy
               DO i=1,sNx
                  jG = myYGlobalLo-1+(bj-1)*sNy+j
        	local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
c                  local=-5.d0
                  uwind(i,j,bi,bj) = local
               ENDDO
            ENDDO
         ENDDO
      ENDDO
	_EXCH_XY_RL(uwind, mythid)
        ENDIF
        IF ( VWindFile .NE. ' '  ) THEN
         CALL READ_FLD_XY_RL( VWindFile,vwind,' ',0,myThid )
	ELSE
      DO bj=1,nSy
         DO bi=1,nSx
            DO j=1,sNy
               DO i=1,sNx
                  jG = myYGlobalLo-1+(bj-1)*sNy+j
                  vwind(i,j,bi,bj) =-0.d0
               ENDDO
            ENDDO
         ENDDO
      ENDDO
	_EXCH_XY_RL(vwind, mythid)
        ENDIF
       ENDIF
C endif for periodicForcing
      ENDIF

      RETURN
      END
