C $Header: /u/gcmpack/MITgcm/verification/tidal_basin_2d/code/external_forcing.F,v 1.8 2014/07/05 15:11:36 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: EXTERNAL_FORCING_U
C !INTERFACE:
SUBROUTINE EXTERNAL_FORCING_U(
I iMin,iMax, jMin,jMax, bi,bj, kLev,
I myTime, 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,iMax :: Working range of x-index for applying forcing.
C jMin,jMax :: Working range of y-index for applying forcing.
C bi,bj :: Current tile indices
C kLev :: Current vertical level index
C myTime :: Current time in simulation
C myThid :: Thread Id number
INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
_RL myTime
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C i,j :: Loop counters
C kSurface :: index of surface layer
INTEGER i, j
INTEGER kSurface
_RL tidal_freq,tidal_Hscale
_RL Coord2longitude,longitud1,longitud2
CEOP
IF ( fluidIsAir ) THEN
kSurface = 0
ELSEIF ( usingPCoords ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
C-- Forcing term
C Add windstress momentum impulse into the top-layer
IF ( kLev .EQ. kSurface ) THEN
c DO j=1,sNy
C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
DO j=0,sNy+1
DO i=1,sNx+1
gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
& +foFacMom*surfaceForcingU(i,j,bi,bj)
& *recip_drF(kLev)*recip_hFacW(i,j,kLev,bi,bj)
ENDDO
ENDDO
ENDIF
C-- Tidal body force: written as gradient of geopotential
C True M2 frequency is
c tidal_freq=2.*pi/(43200.+25.*60.)
C But for convenience we are using 12 hour period
tidal_freq=2.*pi/(43200.)
C Make the tide relatively strong (about 1 m)
tidal_Hscale=10.
IF ( usingCartesianGrid ) THEN
Coord2longitude=1./rSphere
ELSEIF ( usingSphericalPolarGrid ) THEN
Coord2longitude=pi/180.
ELSE
STOP 'Be careful about 2D!'
ENDIF
DO j=0,sNy+1
DO i=1,sNx+1
longitud1=XC(i-1,j,bi,bj)*Coord2longitude
longitud2=XC(i,j,bi,bj)*Coord2longitude
gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
& +gravity*tidal_Hscale*
& ( SIN( tidal_freq*myTime + 2.*longitud2 )
& -SIN( tidal_freq*myTime + 2.*longitud1 )
& )*recip_DXC(i,j,bi,bj)
& *_maskW(i,j,kLev,bi,bj)
c & *min( myTime/86400. , 1.)
ENDDO
ENDDO
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: EXTERNAL_FORCING_V
C !INTERFACE:
SUBROUTINE EXTERNAL_FORCING_V(
I iMin,iMax, jMin,jMax, bi,bj, kLev,
I myTime, 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,iMax :: Working range of x-index for applying forcing.
C jMin,jMax :: Working range of y-index for applying forcing.
C bi,bj :: Current tile indices
C kLev :: Current vertical level index
C myTime :: Current time in simulation
C myThid :: Thread Id number
INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
_RL myTime
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C i,j :: Loop counters
C kSurface :: index of surface layer
INTEGER i, j
INTEGER kSurface
CEOP
IF ( fluidIsAir ) THEN
kSurface = 0
ELSEIF ( usingPCoords ) 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=1,sNy+1
c DO i=1,sNx
C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
DO i=0,sNx+1
gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
& +foFacMom*surfaceForcingV(i,j,bi,bj)
& *recip_drF(kLev)*recip_hFacS(i,j,kLev,bi,bj)
ENDDO
ENDDO
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: EXTERNAL_FORCING_T
C !INTERFACE:
SUBROUTINE EXTERNAL_FORCING_T(
I iMin,iMax, jMin,jMax, bi,bj, kLev,
I myTime, 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, 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,iMax :: Working range of x-index for applying forcing.
C jMin,jMax :: Working range of y-index for applying forcing.
C bi,bj :: Current tile indices
C kLev :: Current vertical level index
C myTime :: Current time in simulation
C myThid :: Thread Id number
INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
_RL myTime
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C i,j :: Loop counters
C kSurface :: index of surface layer
INTEGER i, j
INTEGER kSurface
_RL recip_Cp
CEOP
#ifdef SHORTWAVE_HEATING
integer two
_RL minusone
parameter (two=2,minusone=-1.)
_RL swfracb(two)
INTEGER kp1
#endif
IF ( fluidIsAir ) THEN
kSurface = 0
ELSEIF ( usingPCoords ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
recip_Cp = 1. _d 0 / HeatCapacity_Cp
C-- Forcing term
C Add heat in top-layer
IF ( kLev .EQ. kSurface ) THEN
DO j=1,sNy
DO i=1,sNx
gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
& +surfaceForcingT(i,j,bi,bj)
& *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)
ENDDO
ENDDO
ENDIF
#ifdef SHORTWAVE_HEATING
C Penetrating SW radiation
c IF ( usePenetratingSW ) THEN
swfracb(1)=abs(rF(klev))
swfracb(2)=abs(rF(klev+1))
CALL SWFRAC(
I two, minusone,
U swfracb,
I myTime, 1, myThid )
kp1 = klev+1
IF (klev.EQ.Nr) THEN
kp1 = klev
swfracb(2)=0. _d 0
ENDIF
DO j=1,sNy
DO i=1,sNx
gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
& -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
& -swfracb(2)*maskC(i,j,kp1, bi,bj))
& *recip_Cp*recip_rhoConst
& *recip_drF(klev)*recip_hFacC(i,j,kLev,bi,bj)
ENDDO
ENDDO
c ENDIF
#endif
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: EXTERNAL_FORCING_S
C !INTERFACE:
SUBROUTINE EXTERNAL_FORCING_S(
I iMin,iMax, jMin,jMax, bi,bj, kLev,
I myTime, 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, 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,iMax :: Working range of x-index for applying forcing.
C jMin,jMax :: Working range of y-index for applying forcing.
C bi,bj :: Current tile indices
C kLev :: Current vertical level index
C myTime :: Current time in simulation
C myThid :: Thread Id number
INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
_RL myTime
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C i,j :: Loop counters
C kSurface :: index of surface layer
INTEGER i, j
INTEGER kSurface
CEOP
IF ( fluidIsAir ) THEN
kSurface = 0
ELSEIF ( usingPCoords ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
C-- Forcing term
C Add fresh-water in top-layer
IF ( kLev .EQ. kSurface ) THEN
DO j=1,sNy
DO i=1,sNx
gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
& +surfaceForcingS(i,j,bi,bj)
& *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)
ENDDO
ENDDO
ENDIF
RETURN
END