tidal_basin_2d/code/external_forcing.F

C $Header: /u/gcmpack/MITgcm/verification/tidal_basin_2d/code/external_forcing.F,v 1.7 2014/05/06 15:52:51 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
1	jmc	1.8	C $Header: /u/gcmpack/MITgcm/verification/tidal_basin_2d/code/external_forcing.F,v 1.7 2014/05/06 15:52:51 jmc Exp $
2	adcroft	1.1	C $Name: $
3
4	jmc	1.5	#include "PACKAGES_CONFIG.h"
5	adcroft	1.1	#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: EXTERNAL_FORCING_U
9			C !INTERFACE:
10			SUBROUTINE EXTERNAL_FORCING_U(
11	jmc	1.5	I iMin,iMax, jMin,jMax, bi,bj, kLev,
12			I myTime, myThid )
13	adcroft	1.1	C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.5	C \| S/R EXTERNAL_FORCING_U
16			C \| o Contains problem specific forcing for zonal velocity.
17	adcroft	1.1	C ==========================================================
18	jmc	1.5	C \| Adds terms to gU for forcing by external sources
19			C \| e.g. wind stress, bottom friction etc ...
20	adcroft	1.1	C ==========================================================
21			C \ev
22
23			C !USES:
24			IMPLICIT NONE
25			C == Global data ==
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "DYNVARS.h"
31			#include "FFIELDS.h"
32
33			C !INPUT/OUTPUT PARAMETERS:
34			C == Routine arguments ==
35	jmc	1.5	C iMin,iMax :: Working range of x-index for applying forcing.
36			C jMin,jMax :: Working range of y-index for applying forcing.
37			C bi,bj :: Current tile indices
38			C kLev :: Current vertical level index
39			C myTime :: Current time in simulation
40			C myThid :: Thread Id number
41	adcroft	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
42	jmc	1.5	_RL myTime
43	adcroft	1.1	INTEGER myThid
44
45			C !LOCAL VARIABLES:
46			C == Local variables ==
47	jmc	1.5	C i,j :: Loop counters
48			C kSurface :: index of surface layer
49			INTEGER i, j
50	adcroft	1.1	INTEGER kSurface
51			_RL tidal_freq,tidal_Hscale
52			_RL Coord2longitude,longitud1,longitud2
53			CEOP
54
55	jmc	1.5	IF ( fluidIsAir ) THEN
56			kSurface = 0
57			ELSEIF ( usingPCoords ) THEN
58	adcroft	1.1	kSurface = Nr
59	jmc	1.5	ELSE
60	adcroft	1.1	kSurface = 1
61	jmc	1.5	ENDIF
62	adcroft	1.1
63			C-- Forcing term
64	jmc	1.5
65	adcroft	1.1	C Add windstress momentum impulse into the top-layer
66			IF ( kLev .EQ. kSurface ) THEN
67	jmc	1.5	c DO j=1,sNy
68			C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
69			DO j=0,sNy+1
70			DO i=1,sNx+1
71	adcroft	1.1	gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
72	jmc	1.4	& +foFacMom*surfaceForcingU(i,j,bi,bj)
73			& recip_drF(kLev)recip_hFacW(i,j,kLev,bi,bj)
74	adcroft	1.1	ENDDO
75			ENDDO
76			ENDIF
77
78			C-- Tidal body force: written as gradient of geopotential
79			C True M2 frequency is
80			c tidal_freq=2.pi/(43200.+25.60.)
81			C But for convenience we are using 12 hour period
82			tidal_freq=2.*pi/(43200.)
83			C Make the tide relatively strong (about 1 m)
84	adcroft	1.2	tidal_Hscale=10.
85	adcroft	1.1	IF ( usingCartesianGrid ) THEN
86			Coord2longitude=1./rSphere
87			ELSEIF ( usingSphericalPolarGrid ) THEN
88			Coord2longitude=pi/180.
89			ELSE
90			STOP 'Be careful about 2D!'
91			ENDIF
92	jmc	1.5	DO j=0,sNy+1
93			DO i=1,sNx+1
94	adcroft	1.1	longitud1=XC(i-1,j,bi,bj)*Coord2longitude
95			longitud2=XC(i,j,bi,bj)*Coord2longitude
96			gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
97			& +gravitytidal_Hscale
98	jmc	1.5	& ( SIN( tidal_freqmyTime + 2.longitud2 )
99			& -SIN( tidal_freqmyTime + 2.longitud1 )
100	adcroft	1.1	& )*recip_DXC(i,j,bi,bj)
101			& *_maskW(i,j,kLev,bi,bj)
102	jmc	1.5	c & *min( myTime/86400. , 1.)
103	adcroft	1.1	ENDDO
104			ENDDO
105
106			RETURN
107			END
108	jmc	1.5
109			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
110	adcroft	1.1	CBOP
111			C !ROUTINE: EXTERNAL_FORCING_V
112			C !INTERFACE:
113			SUBROUTINE EXTERNAL_FORCING_V(
114	jmc	1.5	I iMin,iMax, jMin,jMax, bi,bj, kLev,
115			I myTime, myThid )
116	adcroft	1.1	C !DESCRIPTION: \bv
117			C ==========================================================
118	jmc	1.5	C \| S/R EXTERNAL_FORCING_V
119			C \| o Contains problem specific forcing for merid velocity.
120	adcroft	1.1	C ==========================================================
121	jmc	1.5	C \| Adds terms to gV for forcing by external sources
122			C \| e.g. wind stress, bottom friction etc ...
123	adcroft	1.1	C ==========================================================
124			C \ev
125
126			C !USES:
127			IMPLICIT NONE
128			C == Global data ==
129			#include "SIZE.h"
130			#include "EEPARAMS.h"
131			#include "PARAMS.h"
132			#include "GRID.h"
133			#include "DYNVARS.h"
134			#include "FFIELDS.h"
135
136			C !INPUT/OUTPUT PARAMETERS:
137			C == Routine arguments ==
138	jmc	1.5	C iMin,iMax :: Working range of x-index for applying forcing.
139			C jMin,jMax :: Working range of y-index for applying forcing.
140			C bi,bj :: Current tile indices
141			C kLev :: Current vertical level index
142			C myTime :: Current time in simulation
143			C myThid :: Thread Id number
144	adcroft	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
145	jmc	1.5	_RL myTime
146	adcroft	1.1	INTEGER myThid
147
148			C !LOCAL VARIABLES:
149			C == Local variables ==
150	jmc	1.5	C i,j :: Loop counters
151			C kSurface :: index of surface layer
152			INTEGER i, j
153	adcroft	1.1	INTEGER kSurface
154			CEOP
155
156	jmc	1.5	IF ( fluidIsAir ) THEN
157			kSurface = 0
158			ELSEIF ( usingPCoords ) THEN
159	adcroft	1.1	kSurface = Nr
160	jmc	1.5	ELSE
161	adcroft	1.1	kSurface = 1
162	jmc	1.5	ENDIF
163	adcroft	1.1
164			C-- Forcing term
165	jmc	1.5
166	adcroft	1.1	C Add windstress momentum impulse into the top-layer
167			IF ( kLev .EQ. kSurface ) THEN
168	jmc	1.5	DO j=1,sNy+1
169			c DO i=1,sNx
170			C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
171			DO i=0,sNx+1
172	adcroft	1.1	gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
173	jmc	1.4	& +foFacMom*surfaceForcingV(i,j,bi,bj)
174			& recip_drF(kLev)recip_hFacS(i,j,kLev,bi,bj)
175	adcroft	1.1	ENDDO
176			ENDDO
177			ENDIF
178
179			RETURN
180			END
181	jmc	1.5
182			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
183	adcroft	1.1	CBOP
184			C !ROUTINE: EXTERNAL_FORCING_T
185			C !INTERFACE:
186			SUBROUTINE EXTERNAL_FORCING_T(
187	jmc	1.5	I iMin,iMax, jMin,jMax, bi,bj, kLev,
188			I myTime, myThid )
189	adcroft	1.1	C !DESCRIPTION: \bv
190			C ==========================================================
191	jmc	1.5	C \| S/R EXTERNAL_FORCING_T
192			C \| o Contains problem specific forcing for temperature.
193	adcroft	1.1	C ==========================================================
194	jmc	1.5	C \| Adds terms to gT for forcing by external sources
195			C \| e.g. heat flux, climatalogical relaxation, etc ...
196	adcroft	1.1	C ==========================================================
197			C \ev
198
199			C !USES:
200			IMPLICIT NONE
201			C == Global data ==
202			#include "SIZE.h"
203			#include "EEPARAMS.h"
204			#include "PARAMS.h"
205			#include "GRID.h"
206			#include "DYNVARS.h"
207			#include "FFIELDS.h"
208
209			C !INPUT/OUTPUT PARAMETERS:
210			C == Routine arguments ==
211	jmc	1.5	C iMin,iMax :: Working range of x-index for applying forcing.
212			C jMin,jMax :: Working range of y-index for applying forcing.
213			C bi,bj :: Current tile indices
214			C kLev :: Current vertical level index
215			C myTime :: Current time in simulation
216			C myThid :: Thread Id number
217	adcroft	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
218	jmc	1.5	_RL myTime
219	adcroft	1.1	INTEGER myThid
220
221			C !LOCAL VARIABLES:
222			C == Local variables ==
223	jmc	1.5	C i,j :: Loop counters
224			C kSurface :: index of surface layer
225			INTEGER i, j
226	adcroft	1.1	INTEGER kSurface
227	jmc	1.7	_RL recip_Cp
228	adcroft	1.1	CEOP
229	jmc	1.5	#ifdef SHORTWAVE_HEATING
230			integer two
231			_RL minusone
232			parameter (two=2,minusone=-1.)
233			_RL swfracb(two)
234			INTEGER kp1
235			#endif
236	adcroft	1.1
237	jmc	1.5	IF ( fluidIsAir ) THEN
238			kSurface = 0
239			ELSEIF ( usingPCoords ) THEN
240	adcroft	1.1	kSurface = Nr
241	jmc	1.5	ELSE
242	adcroft	1.1	kSurface = 1
243	jmc	1.5	ENDIF
244	jmc	1.7	recip_Cp = 1. _d 0 / HeatCapacity_Cp
245	adcroft	1.1
246			C-- Forcing term
247	jmc	1.5
248	adcroft	1.1	C Add heat in top-layer
249			IF ( kLev .EQ. kSurface ) THEN
250	jmc	1.5	DO j=1,sNy
251			DO i=1,sNx
252	adcroft	1.1	gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
253	jmc	1.4	& +surfaceForcingT(i,j,bi,bj)
254			& recip_drF(kLev)recip_hFacC(i,j,kLev,bi,bj)
255	adcroft	1.1	ENDDO
256			ENDDO
257			ENDIF
258
259			#ifdef SHORTWAVE_HEATING
260			C Penetrating SW radiation
261	jmc	1.5	c IF ( usePenetratingSW ) THEN
262			swfracb(1)=abs(rF(klev))
263			swfracb(2)=abs(rF(klev+1))
264			CALL SWFRAC(
265	jmc	1.6	I two, minusone,
266			U swfracb,
267			I myTime, 1, myThid )
268	jmc	1.5	kp1 = klev+1
269			IF (klev.EQ.Nr) THEN
270			kp1 = klev
271			swfracb(2)=0. _d 0
272			ENDIF
273			DO j=1,sNy
274			DO i=1,sNx
275			gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
276			& -Qsw(i,j,bi,bj)(swfracb(1)maskC(i,j,klev,bi,bj)
277			& -swfracb(2)*maskC(i,j,kp1, bi,bj))
278			& recip_Cprecip_rhoConst
279			& recip_drF(klev)recip_hFacC(i,j,kLev,bi,bj)
280			ENDDO
281	adcroft	1.1	ENDDO
282	jmc	1.5	c ENDIF
283	adcroft	1.1	#endif
284
285			RETURN
286			END
287	jmc	1.5
288			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
289	adcroft	1.1	CBOP
290			C !ROUTINE: EXTERNAL_FORCING_S
291			C !INTERFACE:
292			SUBROUTINE EXTERNAL_FORCING_S(
293	jmc	1.5	I iMin,iMax, jMin,jMax, bi,bj, kLev,
294			I myTime, myThid )
295	adcroft	1.1
296			C !DESCRIPTION: \bv
297			C ==========================================================
298	jmc	1.5	C \| S/R EXTERNAL_FORCING_S
299			C \| o Contains problem specific forcing for merid velocity.
300	adcroft	1.1	C ==========================================================
301	jmc	1.5	C \| Adds terms to gS for forcing by external sources
302			C \| e.g. fresh-water flux, climatalogical relaxation, etc ...
303	adcroft	1.1	C ==========================================================
304			C \ev
305
306			C !USES:
307			IMPLICIT NONE
308			C == Global data ==
309			#include "SIZE.h"
310			#include "EEPARAMS.h"
311			#include "PARAMS.h"
312			#include "GRID.h"
313			#include "DYNVARS.h"
314			#include "FFIELDS.h"
315
316			C !INPUT/OUTPUT PARAMETERS:
317			C == Routine arguments ==
318	jmc	1.5	C iMin,iMax :: Working range of x-index for applying forcing.
319			C jMin,jMax :: Working range of y-index for applying forcing.
320			C bi,bj :: Current tile indices
321			C kLev :: Current vertical level index
322			C myTime :: Current time in simulation
323			C myThid :: Thread Id number
324	adcroft	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
325	jmc	1.5	_RL myTime
326	adcroft	1.1	INTEGER myThid
327
328			C !LOCAL VARIABLES:
329			C == Local variables ==
330	jmc	1.5	C i,j :: Loop counters
331			C kSurface :: index of surface layer
332			INTEGER i, j
333	adcroft	1.1	INTEGER kSurface
334			CEOP
335
336	jmc	1.5	IF ( fluidIsAir ) THEN
337			kSurface = 0
338			ELSEIF ( usingPCoords ) THEN
339	adcroft	1.1	kSurface = Nr
340	jmc	1.5	ELSE
341	adcroft	1.1	kSurface = 1
342	jmc	1.5	ENDIF
343	adcroft	1.1
344	jmc	1.5	C-- Forcing term
345	adcroft	1.1
346			C Add fresh-water in top-layer
347			IF ( kLev .EQ. kSurface ) THEN
348	jmc	1.5	DO j=1,sNy
349			DO i=1,sNx
350	adcroft	1.1	gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
351	jmc	1.4	& +surfaceForcingS(i,j,bi,bj)
352			& recip_drF(kLev)recip_hFacC(i,j,kLev,bi,bj)
353	adcroft	1.1	ENDDO
354			ENDDO
355			ENDIF
356
357			RETURN
358			END