C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/external_forcing_surf.F,v 1.26 2004/07/18 01:04:23 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
 
CBOP
C     !ROUTINE: EXTERNAL_FORCING_SURF
C     !INTERFACE:
      SUBROUTINE EXTERNAL_FORCING_SURF( 
     I             bi, bj, iMin, iMax, jMin, jMax,
     I             myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE EXTERNAL_FORCING_SURF                          
C     | o Determines forcing terms based on external fields       
C     |   relaxation terms etc.                                   
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#ifdef ALLOW_SEAICE
#include "SEAICE.h"
#endif /* ALLOW_SEAICE */
 
C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     bi,bj  :: tile indices
C     iMin,iMax, jMin,jMax :: Range of points for calculation
C     myTime :: Current time in simulation
C     myIter :: Current iteration number in simulation
C     myThid :: Thread no. that called this routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax

C     !LOCAL VARIABLES:
C     === Local variables ===
      INTEGER i,j
C     number of surface interface layer
      INTEGER ks
CEOP

      IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
       ks        = Nr 
      ELSE
       ks        = 1
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
C--   Start with surface restoring term :

       DO j = jMin, jMax
        DO i = iMin, iMax
#ifdef ALLOW_SEAICE
C     Don't restore under sea-ice
C     Heat Flux (restoring term) : 
          surfaceForcingT(i,j,bi,bj) = 
     &      -lambdaThetaClimRelax * (1-AREA(i,j,1,bi,bj))
     &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
     &         *drF(ks)*hFacC(i,j,ks,bi,bj)
C     Salt Flux (restoring term) : 
          surfaceForcingS(i,j,bi,bj) =
     &      -lambdaSaltClimRelax * (1-AREA(i,j,1,bi,bj))
     &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
     &         *drF(ks)*hFacC(i,j,ks,bi,bj)
#else /* ifndef ALLOW_SEAICE */
C     Heat Flux (restoring term) : 
         IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
          surfaceForcingT(i,j,bi,bj) = 
     &      -lambdaThetaClimRelax
     &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
     &         *drF(ks)*hFacC(i,j,ks,bi,bj)
C     Salt Flux (restoring term) : 
          surfaceForcingS(i,j,bi,bj) =
     &      -lambdaSaltClimRelax
     &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
     &         *drF(ks)*hFacC(i,j,ks,bi,bj)
         ELSE
          surfaceForcingT(i,j,bi,bj) = 0. _d 0
          surfaceForcingS(i,j,bi,bj) = 0. _d 0
         ENDIF
#endif /* ALLOW_SEAICE */
        ENDDO
       ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef NONLIN_FRSURF
C-    T,S surface forcing will be applied (thermodynamics) after the update
C     of surf.thickness (hFac): account for change in surf.thickness
       IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
        IF (select_rStar.GT.0) THEN
         DO j=jMin,jMax
          DO i=iMin,iMax
            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
     &                                  * rStarExpC(i,j,bi,bj)
            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &                                  * rStarExpC(i,j,bi,bj)
          ENDDO
         ENDDO
        ELSE
         DO j=jMin,jMax
          DO i=iMin,iMax
           IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
     &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
           ENDIF
          ENDDO
         ENDDO
        ENDIF
       ENDIF
#endif /* NONLIN_FRSURF */

      ELSE
C--   No restoring for T & S : set surfaceForcingT,S to zero :

       DO j = jMin, jMax
        DO i = iMin, iMax
          surfaceForcingT(i,j,bi,bj) = 0. _d 0
          surfaceForcingS(i,j,bi,bj) = 0. _d 0
        ENDDO
       ENDDO

C--   end restoring / no restoring block.
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Surface Fluxes :

      DO j = jMin, jMax
         DO i = iMin, iMax

C     Zonal wind stress fu:
          surfaceForcingU(i,j,bi,bj) = 
     &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
C     Meridional wind stress fv:
          surfaceForcingV(i,j,bi,bj) = 
     &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
C     Net heat flux Qnet:
          surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
     &       - ( Qnet(i,j,bi,bj)
#ifdef SHORTWAVE_HEATING
     &          -Qsw(i,j,bi,bj)
#endif
     &         ) *recip_Cp*horiVertRatio*recip_rhoConst
C     Net Salt Flux : 
          surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst

         ENDDO
      ENDDO
        
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Fresh-water flux 

#ifdef EXACT_CONSERV
c NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
c     to stay consitent with volume change (=d/dt etaH).
      IF ( staggerTimeStep ) THEN
        DO j=1,sNy
         DO i=1,sNx
           PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
         ENDDO
        ENDDO
      ENDIF

      IF ( (nonlinFreeSurf.GT.0 .OR. buoyancyRelation.EQ.'OCEANICP')
     &     .AND. useRealFreshWaterFlux ) THEN

c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes 
c   the water column height ; temp., salt, (tracer) flux associated 
c   with this input/output of water is added here to the surface tendency.
c

       IF (temp_EvPrRn.NE.UNSET_RL) THEN
        DO j = jMin, jMax
         DO i = iMin, iMax
          surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
     &      + PmEpR(i,j,bi,bj)
     &         *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
     &         *convertEmP2rUnit
         ENDDO
        ENDDO
       ENDIF

       IF (salt_EvPrRn.NE.UNSET_RL) THEN
        DO j = jMin, jMax
         DO i = iMin, iMax
          surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &      + PmEpR(i,j,bi,bj)
     &         *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
     &         *convertEmP2rUnit
         ENDDO
        ENDDO
       ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      ELSE
#else /* EXACT_CONSERV */
      IF (.TRUE.) THEN
#endif /* EXACT_CONSERV */

c- EmPmR does not really affect the water column height (for tracer budget)
c   and is converted to a salt tendency.

       IF (convertFW2Salt .EQ. -1.) THEN
c- converts EmPmR to salinity tendency using surface local salinity
        DO j = jMin, jMax
         DO i = iMin, iMax
          surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &      + EmPmR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
     &         *convertEmP2rUnit
         ENDDO
        ENDDO
       ELSE 
c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
        DO j = jMin, jMax
         DO i = iMin, iMax
          surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
     &      + EmPmR(i,j,bi,bj)*convertFW2Salt
     &         *convertEmP2rUnit
         ENDDO
        ENDDO
       ENDIF

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_PTRACERS
      IF ( usePTRACERS ) THEN
         CALL PTRACERS_FORCING_SURF(
     I        bi, bj, iMin, iMax, jMin, jMax,
     I        myTime,myIter,myThid )
      ENDIF
#endif /* ALLOW_PTRACERS */

#ifdef ATMOSPHERIC_LOADING

C-- Atmospheric surface Pressure loading :

      IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
       IF ( useRealFreshWaterFlux ) THEN
        DO j = jMin, jMax
         DO i = iMin, iMax
          phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
     &                        +sIceLoad(i,j,bi,bj)*gravity
     &                          )*recip_rhoConst
         ENDDO
        ENDDO
       ELSE
        DO j = jMin, jMax
         DO i = iMin, iMax
          phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
         ENDDO
        ENDDO
       ENDIF
      ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN 
C-- This is a hack used to read phi0surf from a file (ploadFile)
C   instead of computing it from bathymetry & density ref. profile.
C   The true atmospheric P-loading is not yet implemented for P-coord
C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
        DO j = jMin, jMax
         DO i = iMin, iMax
          phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
         ENDDO
        ENDDO
      ENDIF

#endif /* ATMOSPHERIC_LOADING */

#ifdef ALLOW_EBM
c--    Values for surfaceForcingT, surfaceForcingS
c      are overwritten by those produced by EBM
cph    AD recomputation problems if these IF useEBM are used
cph      IF ( useEBM ) THEN
       CALL EBM_FORCING_SURF(
     I        bi, bj, iMin, iMax, jMin, jMax,
     I        myTime,myIter,myThid )
cph      ENDIF
#endif

#ifdef ALLOW_TIMEAVE
c     IF ( taveFreq .NE. 0. _d 0 ) THEN
c       CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
c     ENDIF
#endif /* ALLOW_TIMEAVE */

      RETURN
      END