model/src/external_forcing.F

C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.34 2005/12/15 17:47:54 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif

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
CEOP

      IF ( fluidIsAir ) THEN
       kSurface = 0
      ELSEIF ( usingPCoords ) THEN
       kSurface = Nr
      ELSE
       kSurface = 1
      ENDIF

C--   Forcing term
#ifdef ALLOW_AIM
      IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_AIM */

#ifdef ALLOW_FIZHI
      IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_FIZHI */

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

#if (defined (ALLOW_TAU_EDDY))
       CALL TAUEDDY_EXTERNAL_FORCING_U(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
#endif

#ifdef ALLOW_OBCS
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_U(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
      ENDIF
#endif

      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
#ifdef ALLOW_AIM
      IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_AIM */

#ifdef ALLOW_FIZHI
      IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_FIZHI */

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

#if (defined (ALLOW_TAU_EDDY))
       CALL TAUEDDY_EXTERNAL_FORCING_V(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
#endif

#ifdef ALLOW_OBCS
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_V(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
      ENDIF
#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
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

C--   Forcing term
#ifdef ALLOW_AIM
      IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_AIM */

#ifdef ALLOW_FIZHI
      IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_FIZHI */

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,
     I     myTime,myThid,
     U     swfracb)
       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

#ifdef ALLOW_RBCS
       if (useRBCS) then
          call RBCS_ADD_TENDENCY(bi,bj,klev, 1,
     &                            myTime, myThid )
       endif
#endif

#ifdef ALLOW_CLIMTEMP_RELAXATION
       IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN
        DO j=1,sNy
         DO i=1,sNx
          gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
     &     -1./tauThetaClimRelax3Dim
     &         *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj))
     &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)
         ENDDO
        ENDDO
       ENDIF
#endif

#ifdef ALLOW_OBCS
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_T(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
      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
#ifdef ALLOW_AIM
      IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_AIM */

#ifdef ALLOW_FIZHI
      IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
     &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
     &                      myTime, myThid )
#endif /* ALLOW_FIZHI */

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

#ifdef ALLOW_RBCS
       if (useRBCS) then
          call RBCS_ADD_TENDENCY(bi,bj,klev, 2,
     &                            myTime, myThid )
       endif
#endif

#ifdef ALLOW_CLIMSALT_RELAXATION
       IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN
        DO j=1,sNy
         DO i=1,sNx
          gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj)
     &     -1./tauSaltClimRelax3Dim
     &         *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj))
     &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)
         ENDDO
        ENDDO
       ENDIF
#endif

#ifdef ALLOW_OBCS
      IF (useOBCS) THEN
       CALL OBCS_SPONGE_S(
     I           iMin,iMax, jMin,jMax, bi,bj, kLev,
     I           myTime, myThid )
      ENDIF
#endif

      RETURN
      END
1	stephd	1.35	C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.34 2005/12/15 17:47:54 jmc Exp $
2	adcroft	1.12	C $Name: $
3	cnh	1.1
4	edhill	1.20	#include "PACKAGES_CONFIG.h"
5	cnh	1.1	#include "CPP_OPTIONS.h"
6	heimbach	1.33	#ifdef ALLOW_EXF
7	jmc	1.34	# include "EXF_OPTIONS.h"
8	heimbach	1.33	#endif
9	cnh	1.1
10	cnh	1.13	CBOP
11			C !ROUTINE: EXTERNAL_FORCING_U
12			C !INTERFACE:
13	cnh	1.1	SUBROUTINE EXTERNAL_FORCING_U(
14	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
15			I myTime, myThid )
16	cnh	1.13	C !DESCRIPTION: \bv
17			C ==========================================================
18	jmc	1.31	C \| S/R EXTERNAL_FORCING_U
19			C \| o Contains problem specific forcing for zonal velocity.
20	cnh	1.13	C ==========================================================
21	jmc	1.31	C \| Adds terms to gU for forcing by external sources
22			C \| e.g. wind stress, bottom friction etc ...
23	cnh	1.13	C ==========================================================
24			C \ev
25
26			C !USES:
27	cnh	1.2	IMPLICIT NONE
28	cnh	1.1	C == Global data ==
29			#include "SIZE.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32			#include "GRID.h"
33			#include "DYNVARS.h"
34	cnh	1.2	#include "FFIELDS.h"
35	cnh	1.13
36			C !INPUT/OUTPUT PARAMETERS:
37	cnh	1.1	C == Routine arguments ==
38	jmc	1.31	C iMin,iMax :: Working range of x-index for applying forcing.
39			C jMin,jMax :: Working range of y-index for applying forcing.
40			C bi,bj :: Current tile indices
41			C kLev :: Current vertical level index
42			C myTime :: Current time in simulation
43			C myThid :: Thread Id number
44	cnh	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
45	jmc	1.31	_RL myTime
46	adcroft	1.4	INTEGER myThid
47	cnh	1.1
48	cnh	1.13	C !LOCAL VARIABLES:
49	cnh	1.2	C == Local variables ==
50	jmc	1.31	C i,j :: Loop counters
51			C kSurface :: index of surface layer
52			INTEGER i, j
53	mlosch	1.17	INTEGER kSurface
54	cnh	1.13	CEOP
55	cnh	1.2
56	jmc	1.28	IF ( fluidIsAir ) THEN
57	jmc	1.21	kSurface = 0
58	jmc	1.28	ELSEIF ( usingPCoords ) THEN
59	mlosch	1.17	kSurface = Nr
60	jmc	1.28	ELSE
61	mlosch	1.17	kSurface = 1
62	jmc	1.28	ENDIF
63	mlosch	1.17
64	cnh	1.2	C-- Forcing term
65	jmc	1.21	#ifdef ALLOW_AIM
66			IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
67			& iMin,iMax, jMin,jMax, bi,bj, kLev,
68	jmc	1.31	& myTime, myThid )
69	jmc	1.21	#endif /* ALLOW_AIM */
70	jmc	1.31
71	molod	1.23	#ifdef ALLOW_FIZHI
72			IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
73			& iMin,iMax, jMin,jMax, bi,bj, kLev,
74	jmc	1.31	& myTime, myThid )
75	molod	1.23	#endif /* ALLOW_FIZHI */
76	jmc	1.21
77	cnh	1.2	C Add windstress momentum impulse into the top-layer
78	mlosch	1.17	IF ( kLev .EQ. kSurface ) THEN
79	jmc	1.32	c DO j=1,sNy
80			C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
81			DO j=0,sNy+1
82	jmc	1.31	DO i=1,sNx+1
83	cnh	1.2	gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
84	jmc	1.26	& +foFacMom*surfaceForcingU(i,j,bi,bj)
85			& recip_drF(kLev)recip_hFacW(i,j,kLev,bi,bj)
86	cnh	1.2	ENDDO
87			ENDDO
88			ENDIF
89
90	heimbach	1.30	#if (defined (ALLOW_TAU_EDDY))
91	heimbach	1.29	CALL TAUEDDY_EXTERNAL_FORCING_U(
92	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
93			I myTime, myThid )
94	heimbach	1.29	#endif
95
96	jmc	1.31	#ifdef ALLOW_OBCS
97	heimbach	1.16	IF (useOBCS) THEN
98			CALL OBCS_SPONGE_U(
99	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
100			I myTime, myThid )
101	heimbach	1.16	ENDIF
102	heimbach	1.14	#endif
103
104	cnh	1.1	RETURN
105			END
106	jmc	1.31
107			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
108	cnh	1.13	CBOP
109			C !ROUTINE: EXTERNAL_FORCING_V
110			C !INTERFACE:
111	cnh	1.1	SUBROUTINE EXTERNAL_FORCING_V(
112	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
113			I myTime, myThid )
114	cnh	1.13	C !DESCRIPTION: \bv
115			C ==========================================================
116	jmc	1.31	C \| S/R EXTERNAL_FORCING_V
117			C \| o Contains problem specific forcing for merid velocity.
118	cnh	1.13	C ==========================================================
119	jmc	1.31	C \| Adds terms to gV for forcing by external sources
120			C \| e.g. wind stress, bottom friction etc ...
121	cnh	1.13	C ==========================================================
122			C \ev
123
124			C !USES:
125	cnh	1.2	IMPLICIT NONE
126	cnh	1.1	C == Global data ==
127			#include "SIZE.h"
128			#include "EEPARAMS.h"
129			#include "PARAMS.h"
130			#include "GRID.h"
131			#include "DYNVARS.h"
132	cnh	1.2	#include "FFIELDS.h"
133
134	cnh	1.13	C !INPUT/OUTPUT PARAMETERS:
135	cnh	1.1	C == Routine arguments ==
136	jmc	1.31	C iMin,iMax :: Working range of x-index for applying forcing.
137			C jMin,jMax :: Working range of y-index for applying forcing.
138			C bi,bj :: Current tile indices
139			C kLev :: Current vertical level index
140			C myTime :: Current time in simulation
141			C myThid :: Thread Id number
142	cnh	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
143	jmc	1.31	_RL myTime
144	adcroft	1.4	INTEGER myThid
145	cnh	1.13
146			C !LOCAL VARIABLES:
147	cnh	1.2	C == Local variables ==
148	jmc	1.31	C i,j :: Loop counters
149			C kSurface :: index of surface layer
150			INTEGER i, j
151	mlosch	1.17	INTEGER kSurface
152	cnh	1.13	CEOP
153	cnh	1.2
154	jmc	1.28	IF ( fluidIsAir ) THEN
155	jmc	1.21	kSurface = 0
156	jmc	1.28	ELSEIF ( usingPCoords ) THEN
157	mlosch	1.17	kSurface = Nr
158	jmc	1.28	ELSE
159	mlosch	1.17	kSurface = 1
160	jmc	1.28	ENDIF
161	mlosch	1.17
162	cnh	1.2	C-- Forcing term
163	jmc	1.21	#ifdef ALLOW_AIM
164			IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
165			& iMin,iMax, jMin,jMax, bi,bj, kLev,
166	jmc	1.31	& myTime, myThid )
167	jmc	1.21	#endif /* ALLOW_AIM */
168
169	molod	1.23	#ifdef ALLOW_FIZHI
170			IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
171			& iMin,iMax, jMin,jMax, bi,bj, kLev,
172	jmc	1.31	& myTime, myThid )
173	molod	1.23	#endif /* ALLOW_FIZHI */
174	jmc	1.31
175	cnh	1.2	C Add windstress momentum impulse into the top-layer
176	mlosch	1.17	IF ( kLev .EQ. kSurface ) THEN
177	jmc	1.31	DO j=1,sNy+1
178	jmc	1.32	c DO i=1,sNx
179			C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
180			DO i=0,sNx+1
181	cnh	1.2	gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
182	jmc	1.26	& +foFacMom*surfaceForcingV(i,j,bi,bj)
183			& recip_drF(kLev)recip_hFacS(i,j,kLev,bi,bj)
184	cnh	1.2	ENDDO
185			ENDDO
186			ENDIF
187	cnh	1.1
188	heimbach	1.30	#if (defined (ALLOW_TAU_EDDY))
189	heimbach	1.29	CALL TAUEDDY_EXTERNAL_FORCING_V(
190	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
191			I myTime, myThid )
192	heimbach	1.29	#endif
193
194	jmc	1.31	#ifdef ALLOW_OBCS
195	heimbach	1.16	IF (useOBCS) THEN
196			CALL OBCS_SPONGE_V(
197	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
198			I myTime, myThid )
199	heimbach	1.16	ENDIF
200	heimbach	1.14	#endif
201
202	cnh	1.1	RETURN
203			END
204	jmc	1.31
205			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
206	cnh	1.13	CBOP
207			C !ROUTINE: EXTERNAL_FORCING_T
208			C !INTERFACE:
209	cnh	1.1	SUBROUTINE EXTERNAL_FORCING_T(
210	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
211			I myTime, myThid )
212	cnh	1.13	C !DESCRIPTION: \bv
213			C ==========================================================
214	jmc	1.31	C \| S/R EXTERNAL_FORCING_T
215			C \| o Contains problem specific forcing for temperature.
216	cnh	1.13	C ==========================================================
217	jmc	1.31	C \| Adds terms to gT for forcing by external sources
218			C \| e.g. heat flux, climatalogical relaxation, etc ...
219	cnh	1.13	C ==========================================================
220			C \ev
221
222			C !USES:
223	cnh	1.2	IMPLICIT NONE
224	cnh	1.1	C == Global data ==
225			#include "SIZE.h"
226			#include "EEPARAMS.h"
227			#include "PARAMS.h"
228			#include "GRID.h"
229			#include "DYNVARS.h"
230			#include "FFIELDS.h"
231	heimbach	1.7
232	cnh	1.13	C !INPUT/OUTPUT PARAMETERS:
233	cnh	1.1	C == Routine arguments ==
234	jmc	1.31	C iMin,iMax :: Working range of x-index for applying forcing.
235			C jMin,jMax :: Working range of y-index for applying forcing.
236			C bi,bj :: Current tile indices
237			C kLev :: Current vertical level index
238			C myTime :: Current time in simulation
239			C myThid :: Thread Id number
240	cnh	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
241	jmc	1.31	_RL myTime
242	adcroft	1.4	INTEGER myThid
243	cnh	1.1
244	cnh	1.13	C !LOCAL VARIABLES:
245	cnh	1.2	C == Local variables ==
246	jmc	1.31	C i,j :: Loop counters
247			C kSurface :: index of surface layer
248			INTEGER i, j
249	mlosch	1.17	INTEGER kSurface
250	jmc	1.31	CEOP
251	jmc	1.24	#ifdef SHORTWAVE_HEATING
252			integer two
253			_RL minusone
254			parameter (two=2,minusone=-1.)
255			_RL swfracb(two)
256			INTEGER kp1
257			#endif
258	cnh	1.2
259	jmc	1.28	IF ( fluidIsAir ) THEN
260	jmc	1.21	kSurface = 0
261	jmc	1.28	ELSEIF ( usingPCoords ) THEN
262	mlosch	1.17	kSurface = Nr
263	jmc	1.28	ELSE
264	mlosch	1.17	kSurface = 1
265	jmc	1.28	ENDIF
266	mlosch	1.17
267	cnh	1.2	C-- Forcing term
268	jmc	1.21	#ifdef ALLOW_AIM
269			IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
270			& iMin,iMax, jMin,jMax, bi,bj, kLev,
271	jmc	1.31	& myTime, myThid )
272	jmc	1.21	#endif /* ALLOW_AIM */
273
274	molod	1.23	#ifdef ALLOW_FIZHI
275			IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
276			& iMin,iMax, jMin,jMax, bi,bj, kLev,
277	jmc	1.31	& myTime, myThid )
278	molod	1.23	#endif /* ALLOW_FIZHI */
279	heimbach	1.25
280	cnh	1.2	C Add heat in top-layer
281	mlosch	1.17	IF ( kLev .EQ. kSurface ) THEN
282	jmc	1.31	DO j=1,sNy
283			DO i=1,sNx
284	cnh	1.2	gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
285	jmc	1.26	& +surfaceForcingT(i,j,bi,bj)
286			& recip_drF(kLev)recip_hFacC(i,j,kLev,bi,bj)
287	cnh	1.2	ENDDO
288			ENDDO
289			ENDIF
290	adcroft	1.5
291			#ifdef SHORTWAVE_HEATING
292			C Penetrating SW radiation
293	jmc	1.31	c IF ( usePenetratingSW ) THEN
294			swfracb(1)=abs(rF(klev))
295			swfracb(2)=abs(rF(klev+1))
296			CALL SWFRAC(
297	heimbach	1.8	I two,minusone,
298	jmc	1.31	I myTime,myThid,
299	dimitri	1.18	U swfracb)
300	jmc	1.31	kp1 = klev+1
301			IF (klev.EQ.Nr) THEN
302	jmc	1.24	kp1 = klev
303			swfracb(2)=0. _d 0
304	jmc	1.31	ENDIF
305			DO j=1,sNy
306			DO i=1,sNx
307			gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
308	jmc	1.24	& -Qsw(i,j,bi,bj)(swfracb(1)maskC(i,j,klev,bi,bj)
309			& -swfracb(2)*maskC(i,j,kp1, bi,bj))
310	jmc	1.27	& recip_Cprecip_rhoConst
311			& recip_drF(klev)recip_hFacC(i,j,kLev,bi,bj)
312	jmc	1.31	ENDDO
313	adcroft	1.5	ENDDO
314	jmc	1.31	c ENDIF
315	adcroft	1.5	#endif
316	heimbach	1.14
317	stephd	1.35	#ifdef ALLOW_RBCS
318			if (useRBCS) then
319			call RBCS_ADD_TENDENCY(bi,bj,klev, 1,
320			& myTime, myThid )
321			endif
322			#endif
323
324	heimbach	1.33	#ifdef ALLOW_CLIMTEMP_RELAXATION
325			IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN
326			DO j=1,sNy
327			DO i=1,sNx
328			gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
329			& -1./tauThetaClimRelax3Dim
330			& *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj))
331			& hFacC(i,j,klev,bi,bj)recip_hFacC(i,j,kLev,bi,bj)
332			ENDDO
333			ENDDO
334			ENDIF
335			#endif
336
337	jmc	1.31	#ifdef ALLOW_OBCS
338	heimbach	1.16	IF (useOBCS) THEN
339			CALL OBCS_SPONGE_T(
340	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
341			I myTime, myThid )
342	heimbach	1.16	ENDIF
343	heimbach	1.14	#endif
344
345	cnh	1.1	RETURN
346			END
347	jmc	1.31
348			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
349	cnh	1.13	CBOP
350			C !ROUTINE: EXTERNAL_FORCING_S
351			C !INTERFACE:
352	cnh	1.1	SUBROUTINE EXTERNAL_FORCING_S(
353	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
354			I myTime, myThid )
355	cnh	1.13
356			C !DESCRIPTION: \bv
357			C ==========================================================
358	jmc	1.31	C \| S/R EXTERNAL_FORCING_S
359			C \| o Contains problem specific forcing for merid velocity.
360	cnh	1.13	C ==========================================================
361	jmc	1.31	C \| Adds terms to gS for forcing by external sources
362			C \| e.g. fresh-water flux, climatalogical relaxation, etc ...
363	cnh	1.13	C ==========================================================
364			C \ev
365
366			C !USES:
367	cnh	1.2	IMPLICIT NONE
368	cnh	1.1	C == Global data ==
369			#include "SIZE.h"
370			#include "EEPARAMS.h"
371			#include "PARAMS.h"
372			#include "GRID.h"
373			#include "DYNVARS.h"
374	cnh	1.2	#include "FFIELDS.h"
375	cnh	1.1
376	cnh	1.13	C !INPUT/OUTPUT PARAMETERS:
377	cnh	1.1	C == Routine arguments ==
378	jmc	1.31	C iMin,iMax :: Working range of x-index for applying forcing.
379			C jMin,jMax :: Working range of y-index for applying forcing.
380			C bi,bj :: Current tile indices
381			C kLev :: Current vertical level index
382			C myTime :: Current time in simulation
383			C myThid :: Thread Id number
384	cnh	1.1	INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
385	jmc	1.31	_RL myTime
386	adcroft	1.4	INTEGER myThid
387	cnh	1.2
388	cnh	1.13	C !LOCAL VARIABLES:
389	cnh	1.2	C == Local variables ==
390	jmc	1.31	C i,j :: Loop counters
391			C kSurface :: index of surface layer
392			INTEGER i, j
393	mlosch	1.17	INTEGER kSurface
394	cnh	1.13	CEOP
395	cnh	1.2
396	jmc	1.28	IF ( fluidIsAir ) THEN
397	jmc	1.21	kSurface = 0
398	jmc	1.28	ELSEIF ( usingPCoords ) THEN
399	mlosch	1.17	kSurface = Nr
400	jmc	1.28	ELSE
401	mlosch	1.17	kSurface = 1
402	jmc	1.28	ENDIF
403	mlosch	1.17
404	cnh	1.2	C-- Forcing term
405	jmc	1.21	#ifdef ALLOW_AIM
406			IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
407			& iMin,iMax, jMin,jMax, bi,bj, kLev,
408	jmc	1.31	& myTime, myThid )
409	jmc	1.21	#endif /* ALLOW_AIM */
410
411	molod	1.23	#ifdef ALLOW_FIZHI
412			IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
413			& iMin,iMax, jMin,jMax, bi,bj, kLev,
414	jmc	1.31	& myTime, myThid )
415	molod	1.23	#endif /* ALLOW_FIZHI */
416	heimbach	1.25
417	cnh	1.2	C Add fresh-water in top-layer
418	mlosch	1.17	IF ( kLev .EQ. kSurface ) THEN
419	jmc	1.31	DO j=1,sNy
420			DO i=1,sNx
421	cnh	1.2	gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
422	jmc	1.26	& +surfaceForcingS(i,j,bi,bj)
423			& recip_drF(kLev)recip_hFacC(i,j,kLev,bi,bj)
424	cnh	1.2	ENDDO
425			ENDDO
426			ENDIF
427	heimbach	1.14
428	stephd	1.35	#ifdef ALLOW_RBCS
429			if (useRBCS) then
430			call RBCS_ADD_TENDENCY(bi,bj,klev, 2,
431			& myTime, myThid )
432			endif
433			#endif
434
435	heimbach	1.33	#ifdef ALLOW_CLIMSALT_RELAXATION
436			IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN
437			DO j=1,sNy
438			DO i=1,sNx
439			gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj)
440			& -1./tauSaltClimRelax3Dim
441			& *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj))
442			& hFacC(i,j,klev,bi,bj)recip_hFacC(i,j,kLev,bi,bj)
443			ENDDO
444			ENDDO
445			ENDIF
446			#endif
447
448	jmc	1.31	#ifdef ALLOW_OBCS
449	heimbach	1.16	IF (useOBCS) THEN
450			CALL OBCS_SPONGE_S(
451	jmc	1.31	I iMin,iMax, jMin,jMax, bi,bj, kLev,
452			I myTime, myThid )
453	heimbach	1.16	ENDIF
454	heimbach	1.14	#endif
455	cnh	1.1
456			RETURN
457			END