C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/natl_12/code/external_forcing.F,v 1.4 2003/08/07 13:03:31 cnh Exp $
C $Name:  $

#include "CPP_OPTIONS.h"
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: EXTERNAL_FORCING_U
C     !INTERFACE:
      SUBROUTINE EXTERNAL_FORCING_U(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R EXTERNAL_FORCING_U                                    
C     | o Contains problem specific forcing for zonal velocity.   
C     *==========================================================*
C     | Adds terms to gU for forcing by external sources          
C     | e.g. wind stress, bottom friction etc..................   
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     iMin - Working range of tile for applying forcing.
C     iMax
C     jMin
C     jMax
C     kLev
      INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
      _RL myCurrentTime
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     Loop counters
      INTEGER I, J
C     number of surface interface layer
      INTEGER kSurface
C     Cheap sponge layer
      _RS recip_tauSp(5)
      INTEGER spWidth 
      _RS curRecipTau 
      INTEGER jFromNBndy, jFromSBndy, 
     &        jNorthBndy, jSouthBndy, jG
CEOP

      if ( buoyancyRelation .eq. 'OCEANICP' ) then
       kSurface = Nr
      else
       kSurface = 1
      endif

C--   Forcing term
C     Add windstress momentum impulse into the top-layer
      IF ( kLev .EQ. kSurface ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
     &   +foFacMom*surfaceTendencyU(i,j,bi,bj)
     &   *_maskW(i,j,kLev,bi,bj)
        ENDDO
       ENDDO
      ENDIF

C--   Create a sponge layer where flow is linearly damped over entire water column
C     Damping time scale decreases away from boundary so that
C     tau = 1 day, 5days, 10days, 20days, 60days
      spWidth  = 5
      recip_tauSp(1)   = 1./(86400.*1. )
      recip_tauSp(2)   = 1./(86400.*5. )
      recip_tauSp(3)   = 1./(86400.*10.)
      recip_tauSp(4)   = 1./(86400.*20.)
      recip_tauSp(5)   = 1./(86400.*60.)
      jSouthBndy = 5
      jNorthBndy = ny-5+1
      DO j=1,sNy
       DO i=iMin,iMax
        jG = myYGlobalLo+(bj-1)*sNy+j-1
        jFromNBndy = jNorthBndy-jG
        jFromSBndy = jSouthBndy-jG
        curRecipTau=0.
        IF ( jFromNBndy .LE. 0 ) THEN
         curRecipTau = recip_tauSp(jFromNBndy+5)
        ENDIF
        IF ( jFromSBndy .GE. 0 ) THEN
         curRecipTau = recip_tauSp(-(jFromSBndy-5))
        ENDIF
        gu(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj) 
     &  -curRecipTau*uVel(i,j,Klev,bi,bj)
       ENDDO
      ENDDO
   
#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_U(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
      ENDIF
#endif

      RETURN
      END
CBOP
C     !ROUTINE: EXTERNAL_FORCING_V
C     !INTERFACE:
      SUBROUTINE EXTERNAL_FORCING_V(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R EXTERNAL_FORCING_V                                    
C     | o Contains problem specific forcing for merid velocity.   
C     *==========================================================*
C     | Adds terms to gV for forcing by external sources          
C     | e.g. wind stress, bottom friction etc..................   
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     iMin - Working range of tile for applying forcing.
C     iMax
C     jMin
C     jMax
C     kLev
      INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
      _RL myCurrentTime
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     Loop counters
      INTEGER I, J
C     number of surface interface layer
      INTEGER kSurface

C     == Cheap sponge layer ==
      _RS recip_tauSp(5)
      INTEGER spWidth
      _RS curRecipTau
      INTEGER jFromNBndy, jFromSBndy,
     &        jNorthBndy, jSouthBndy, jG


CEOP

      if ( buoyancyRelation .eq. 'OCEANICP' ) then
       kSurface = Nr
      else
       kSurface = 1
      endif

C--   Forcing term
C     Add windstress momentum impulse into the top-layer
      IF ( kLev .EQ. kSurface ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
     &   +foFacMom*surfaceTendencyV(i,j,bi,bj)
     &   *_maskS(i,j,kLev,bi,bj)
        ENDDO
       ENDDO
      ENDIF

#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_V(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
      ENDIF
#endif

C--   Create a sponge layer where flow is linearly damped over entire water column
C     Damping time scale decreases away from boundary so that
C     tau = 1 day, 5days, 10days, 20days, 60days
      spWidth  = 5
      recip_tauSp(1)   = 1./(86400.*1. )
      recip_tauSp(2)   = 1./(86400.*5. )
      recip_tauSp(3)   = 1./(86400.*10.)
      recip_tauSp(4)   = 1./(86400.*20.)
      recip_tauSp(5)   = 1./(86400.*60.)
      jSouthBndy = 5
      jNorthBndy = ny-5+1
      DO j=1,sNy
       DO i=iMin,iMax
        jG = myYGlobalLo+(bj-1)*sNy+j-1
        jFromNBndy = jNorthBndy-jG
        jFromSBndy = jSouthBndy-jG
        curRecipTau=0.
        IF ( jFromNBndy .LE. 0 ) THEN
         curRecipTau = recip_tauSp(jFromNBndy+5)
        ENDIF
        IF ( jFromSBndy .GE. 0 ) THEN
         curRecipTau = recip_tauSp(-(jFromSBndy-5))
        ENDIF
        gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
     &  -curRecipTau*vVel(i,j,Klev,bi,bj)
       ENDDO
      ENDDO

      RETURN
      END
CBOP
C     !ROUTINE: EXTERNAL_FORCING_T
C     !INTERFACE:
      SUBROUTINE EXTERNAL_FORCING_T(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R EXTERNAL_FORCING_T                                    
C     | o Contains problem specific forcing for temperature.      
C     *==========================================================*
C     | Adds terms to gT for forcing by external sources          
C     | e.g. heat flux, climatalogical relaxation..............   
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#ifdef SHORTWAVE_HEATING
      integer two
      _RL minusone
      parameter (two=2,minusone=-1.)
      _RL swfracb(two)
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     iMin - Working range of tile for applying forcing.
C     iMax
C     jMin
C     jMax
C     kLev
      INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
      _RL myCurrentTime
      INTEGER myThid
CEndOfInterface

C     !LOCAL VARIABLES:
C     == Local variables ==
C     Loop counters
      INTEGER I, J
C     number of surface interface layer
      INTEGER kSurface
C     Cheap sponge layer
      _RS recip_tauSp(5)
      INTEGER spWidth
      _RS curRecipTau
      INTEGER jFromNBndy, jFromSBndy,
     &        jNorthBndy, jSouthBndy, jG

CEOP

      if ( buoyancyRelation .eq. 'OCEANICP' ) then
       kSurface = Nr
      else
       kSurface = 1
      endif

C--   Forcing term
C     Add heat in top-layer
      IF ( kLev .EQ. kSurface ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
     &     +maskC(i,j,kLev,bi,bj)*surfaceTendencyT(i,j,bi,bj)
        ENDDO
       ENDDO
      ENDIF

C--   Forcing term
C     Add heat in top-layer ( 90 day climatalogical average relaxation )
      IF ( kLev .EQ. kSurface ) THEN
       curRecipTau=1./(86400.*90.)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
     &     +maskC(i,j,kLev,bi,bj)*(
     &      -curRecipTau*(theta(i,j,Klev,bi,bj)-thetaRef(i,j,kLev,bi,bj))
     &      )
        ENDDO
       ENDDO
      ENDIF

#ifdef SHORTWAVE_HEATING
C Penetrating SW radiation
      swfracb(1)=abs(rF(klev))
      swfracb(2)=abs(rF(klev+1))
      call SWFRAC(
     I     two,minusone,
     I     myCurrentTime,myThid,
     U     swfracb)
      DO j=jMin,jMax
       DO i=iMin,iMax
        gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj) 
     &   -maskC(i,j,klev,bi,bj)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))
     &    *recip_Cp*recip_rhoConst*recip_drF(klev)
       ENDDO
      ENDDO
#endif

#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_T(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
      ENDIF
#endif

C--   Create a sponge layer where flow is linearly damped over entire water column
C     Damping time scale decreases away from boundary so that
C     tau = 1 day, 5days, 10days, 20days, 60days
      spWidth  = 5
      recip_tauSp(1)   = 1./(86400.*1. )
      recip_tauSp(2)   = 1./(86400.*5. )
      recip_tauSp(3)   = 1./(86400.*10.)
      recip_tauSp(4)   = 1./(86400.*20.)
      recip_tauSp(5)   = 1./(86400.*60.)
      jSouthBndy = 5
      jNorthBndy = ny-5+1
      DO j=1,sNy
       DO i=iMin,iMax
        jG = myYGlobalLo+(bj-1)*sNy+j-1
        jFromNBndy = jNorthBndy-jG
        jFromSBndy = jSouthBndy-jG
        curRecipTau=0.
        IF ( jFromNBndy .LE. 0 ) THEN
         curRecipTau = recip_tauSp(jFromNBndy+5)
        ENDIF
        IF ( jFromSBndy .GE. 0 ) THEN
         curRecipTau = recip_tauSp(-(jFromSBndy-5))
        ENDIF
        gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
     &  -curRecipTau*(theta(i,j,Klev,bi,bj)-thetaRef(i,j,kLev,bi,bj))
C    &   *0.0000D0
       ENDDO
      ENDDO


      RETURN
      END
CBOP
C     !ROUTINE: EXTERNAL_FORCING_S
C     !INTERFACE:
      SUBROUTINE EXTERNAL_FORCING_S(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R EXTERNAL_FORCING_S                                    
C     | o Contains problem specific forcing for merid velocity.   
C     *==========================================================*
C     | Adds terms to gS for forcing by external sources          
C     | e.g. fresh-water flux, climatalogical relaxation.......   
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     iMin - Working range of tile for applying forcing.
C     iMax
C     jMin
C     jMax
C     kLev
      INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
      _RL myCurrentTime
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     Loop counters
      INTEGER I, J
C     number of surface interface layer
      INTEGER kSurface
C     Cheap sponge layer
      _RS recip_tauSp(5)
      INTEGER spWidth
      _RS curRecipTau
      INTEGER jFromNBndy, jFromSBndy,
     &        jNorthBndy, jSouthBndy, jG

CEOP

      if ( buoyancyRelation .eq. 'OCEANICP' ) then
       kSurface = Nr
      else
       kSurface = 1
      endif


C--   Forcing term
C     Add fresh-water in top-layer
      IF ( kLev .EQ. kSurface ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
     &   +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj)
        ENDDO
       ENDDO
      ENDIF

C--   Forcing term
C     Add freshening/salt in top-layer ( 90 day climatalogical average relaxation )
      IF ( kLev .EQ. kSurface ) THEN
       curRecipTau=1./(86400.*90.)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
     &     +maskC(i,j,kLev,bi,bj)*(
     &      -curRecipTau*(salt(i,j,Klev,bi,bj)-saltRef(i,j,kLev,bi,bj))
     &      )
        ENDDO
       ENDDO
      ENDIF

#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_S(
     I           iMin, iMax, jMin, jMax,bi,bj,kLev,
     I           myCurrentTime,myThid)
      ENDIF
#endif

C--   Create a sponge layer where flow is linearly damped over entire water column
C     Damping time scale decreases away from boundary so that
C     tau = 1 day, 5days, 10days, 20days, 60days
      spWidth  = 5
      recip_tauSp(1)   = 1./(86400.*1. )
      recip_tauSp(2)   = 1./(86400.*5. )
      recip_tauSp(3)   = 1./(86400.*10.)
      recip_tauSp(4)   = 1./(86400.*20.)
      recip_tauSp(5)   = 1./(86400.*60.)
      jSouthBndy = 5
      jNorthBndy = ny-5+1
      DO j=1,sNy
       DO i=iMin,iMax
        jG = myYGlobalLo+(bj-1)*sNy+j-1
        jFromNBndy = jNorthBndy-jG
        jFromSBndy = jSouthBndy-jG
        curRecipTau=0.
        IF ( jFromNBndy .LE. 0 ) THEN
         curRecipTau = recip_tauSp(jFromNBndy+5)
        ENDIF
        IF ( jFromSBndy .GE. 0 ) THEN
         curRecipTau = recip_tauSp(-(jFromSBndy-5))
        ENDIF
        gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
     &  -curRecipTau*(salt(i,j,Klev,bi,bj)-saltRef(i,j,kLev,bi,bj))
C    &   *0.0000D0
       ENDDO
      ENDDO

      RETURN
      END