model/src/external_forcing_surf.F

C $Header: /usr/local/gcmpack/MITgcm/model/src/external_forcing_surf.F,v 1.14 2003/10/16 23:46:47 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             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"
 
C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
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

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

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

C--   Surface Fluxes :

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j = 1-OLy, sNy+OLy
         DO i = 1-OLx, sNx+OLx

c     Zonal wind stress fu:
          surfaceTendencyU(i,j,bi,bj) = 
     &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
     &           *recip_drF(kSurface)*recip_hFacW(i,j,kSurface,bi,bj)
c     Meridional wind stress fv:
          surfaceTendencyV(i,j,bi,bj) = 
     &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
     &           *recip_drF(kSurface)*recip_hFacS(i,j,kSurface,bi,bj)
c     Net heat flux Qnet:
          surfaceTendencyT(i,j,bi,bj) = 
     &      -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst
     &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
C     Net Salt Flux : 
          surfaceTendencyS(i,j,bi,bj) = 
     &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst
     &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)

         ENDDO
        ENDDO
        
C--   Surface restoring term :

        IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
         DO j = 1-OLy, sNy+OLy
          DO i = 1-OLx, sNx+OLx
C     Heat Flux (restoring term) : 
           IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
              surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
     &             -lambdaThetaClimRelax
     &             *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))
C     Salt Flux (restoring term) : 
              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
     &             -lambdaSaltClimRelax
     &             *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))
           ENDIF
          ENDDO
         ENDDO
        ENDIF

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

#ifdef NONLIN_FRSURF
        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
c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay 
c     consitent with volume change (=d/dt etaN).

         IF (temp_EvPrRn.NE.UNSET_RL) THEN
          DO j = 1-OLy, sNy+OLy
           DO i = 1-OLx, sNx+OLx
              surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
     &             + PmEpR(i,j,bi,bj)
     &             *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )
     &             *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) 
     &             *convertEmP2rUnit
           ENDDO
          ENDDO
         ENDIF

         IF (salt_EvPrRn.NE.UNSET_RL) THEN
          DO j = 1-OLy, sNy+OLy
           DO i = 1-OLx, sNx+OLx
              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
     &             + PmEpR(i,j,bi,bj)
     &             *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) )
     &             *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
     &             *convertEmP2rUnit
           ENDDO
          ENDDO
         ENDIF

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

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 = 1-OLy, sNy+OLy
           DO i = 1-OLx, sNx+OLx
              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
     &             + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)
     &             *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
     &             *convertEmP2rUnit
           ENDDO
          ENDDO
         ELSE 
c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
          DO j = 1-OLy, sNy+OLy
           DO i = 1-OLx, sNx+OLx
              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
     &             + EmPmR(i,j,bi,bj)*convertFW2Salt
     &             *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) 
     &             *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, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
     I          myTime,myIter,myThid )
        ENDIF
#endif /* ALLOW_PTRACERS */

#ifdef ATMOSPHERIC_LOADING

C-- Atmospheric surface Pressure loading :

        IF (buoyancyRelation .eq. 'OCEANIC' ) THEN
         DO j = 1-OLy, sNy+OLy
          DO i = 1-OLx, sNx+OLx
             phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
          ENDDO
         ENDDO
        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 = 1-OLy, sNy+OLy
          DO i = 1-OLx, sNx+OLx
             phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF

#endif /* ATMOSPHERIC_LOADING */

       ENDDO
      ENDDO

      RETURN
      END
1	dimitri	1.15	C $Header: /usr/local/gcmpack/MITgcm/model/src/external_forcing_surf.F,v 1.14 2003/10/16 23:46:47 jmc Exp $
2	cnh	1.5	C $Name: $
3	heimbach	1.1
4	edhill	1.13	#include "PACKAGES_CONFIG.h"
5	heimbach	1.1	#include "CPP_OPTIONS.h"
6
7	cnh	1.5	CBOP
8			C !ROUTINE: EXTERNAL_FORCING_SURF
9			C !INTERFACE:
10	heimbach	1.2	SUBROUTINE EXTERNAL_FORCING_SURF(
11	dimitri	1.12	I myTime, myIter, myThid )
12	cnh	1.5	C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE EXTERNAL_FORCING_SURF
15			C \| o Determines forcing terms based on external fields
16			C \| relaxation terms etc.
17			C ==========================================================
18			C \ev
19
20			C !USES:
21	heimbach	1.1	IMPLICIT NONE
22			C === Global variables ===
23			#include "SIZE.h"
24			#include "EEPARAMS.h"
25			#include "PARAMS.h"
26			#include "FFIELDS.h"
27			#include "DYNVARS.h"
28			#include "GRID.h"
29	jmc	1.6	#include "SURFACE.h"
30	heimbach	1.1
31	cnh	1.5	C !INPUT/OUTPUT PARAMETERS:
32	heimbach	1.1	C === Routine arguments ===
33	dimitri	1.12	C myTime - Current time in simulation
34			C myIter - Current iteration number in simulation
35	cnh	1.5	C myThid :: Thread no. that called this routine.
36	dimitri	1.12	_RL myTime
37			INTEGER myIter
38	heimbach	1.1	INTEGER myThid
39
40	cnh	1.5	C !LOCAL VARIABLES:
41	heimbach	1.1	C === Local variables ===
42	dimitri	1.15	INTEGER bi,bj,i,j
43	mlosch	1.7	C number of surface interface layer
44			INTEGER kSurface
45	cnh	1.5	CEOP
46	heimbach	1.2
47	mlosch	1.7	if ( buoyancyRelation .eq. 'OCEANICP' ) then
48			kSurface = Nr
49			else
50			kSurface = 1
51			endif
52
53	jmc	1.14	C-- Surface Fluxes :
54
55	dimitri	1.15	DO bj = myByLo(myThid), myByHi(myThid)
56			DO bi = myBxLo(myThid), myBxHi(myThid)
57			DO j = 1-OLy, sNy+OLy
58			DO i = 1-OLx, sNx+OLx
59	heimbach	1.1
60			c Zonal wind stress fu:
61	jmc	1.6	surfaceTendencyU(i,j,bi,bj) =
62	mlosch	1.7	& fu(i,j,bi,bj)horiVertRatiorecip_rhoConst
63			& recip_drF(kSurface)recip_hFacW(i,j,kSurface,bi,bj)
64	heimbach	1.2	c Meridional wind stress fv:
65	jmc	1.6	surfaceTendencyV(i,j,bi,bj) =
66	mlosch	1.7	& fv(i,j,bi,bj)horiVertRatiorecip_rhoConst
67			& recip_drF(kSurface)recip_hFacS(i,j,kSurface,bi,bj)
68	heimbach	1.1	c Net heat flux Qnet:
69	jmc	1.6	surfaceTendencyT(i,j,bi,bj) =
70	mlosch	1.7	& -Qnet(i,j,bi,bj)recip_CphoriVertRatio*recip_rhoConst
71			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
72	jmc	1.14	C Net Salt Flux :
73			surfaceTendencyS(i,j,bi,bj) =
74			& -saltFlux(i,j,bi,bj)horiVertRatiorecip_rhoConst
75			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
76	jmc	1.6
77	heimbach	1.2	ENDDO
78	dimitri	1.15	ENDDO
79	jmc	1.14
80			C-- Surface restoring term :
81	jmc	1.6
82	dimitri	1.15	IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
83			DO j = 1-OLy, sNy+OLy
84			DO i = 1-OLx, sNx+OLx
85	jmc	1.14	C Heat Flux (restoring term) :
86	dimitri	1.15	IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
87			surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
88			& -lambdaThetaClimRelax
89			& *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))
90	jmc	1.14	C Salt Flux (restoring term) :
91	dimitri	1.15	surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
92			& -lambdaSaltClimRelax
93			& *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))
94			ENDIF
95			ENDDO
96			ENDDO
97			ENDIF
98	jmc	1.14
99			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
100			C-- Fresh-water flux
101	jmc	1.6
102			#ifdef NONLIN_FRSURF
103	dimitri	1.15	IF ( (nonlinFreeSurf.GT.0 .OR. buoyancyRelation.EQ.'OCEANICP')
104			& .AND. useRealFreshWaterFlux ) THEN
105	jmc	1.6
106			c- NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
107			c the water column height ; temp., salt, (tracer) flux associated
108			c with this input/output of water is added here to the surface tendency.
109			c
110			c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay
111			c consitent with volume change (=d/dt etaN).
112
113	dimitri	1.15	IF (temp_EvPrRn.NE.UNSET_RL) THEN
114			DO j = 1-OLy, sNy+OLy
115			DO i = 1-OLx, sNx+OLx
116			surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
117			& + PmEpR(i,j,bi,bj)
118			& *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )
119			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
120			& *convertEmP2rUnit
121			ENDDO
122			ENDDO
123			ENDIF
124	jmc	1.6
125	dimitri	1.15	IF (salt_EvPrRn.NE.UNSET_RL) THEN
126			DO j = 1-OLy, sNy+OLy
127			DO i = 1-OLx, sNx+OLx
128			surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
129			& + PmEpR(i,j,bi,bj)
130			& *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) )
131			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
132			& *convertEmP2rUnit
133			ENDDO
134			ENDDO
135			ENDIF
136	jmc	1.6
137	jmc	1.11	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
138	dimitri	1.15	ELSE
139	jmc	1.11	#else /* NONLIN_FRSURF */
140	dimitri	1.15	IF (.TRUE.) THEN
141	jmc	1.6	#endif /* NONLIN_FRSURF */
142	jmc	1.11
143			c- EmPmR does not really affect the water column height (for tracer budget)
144			c and is converted to a salt tendency.
145
146	dimitri	1.15	IF (convertFW2Salt .EQ. -1.) THEN
147	jmc	1.11	c- converts EmPmR to salinity tendency using surface local salinity
148	dimitri	1.15	DO j = 1-OLy, sNy+OLy
149			DO i = 1-OLx, sNx+OLx
150			surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
151			& + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)
152			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
153			& *convertEmP2rUnit
154			ENDDO
155			ENDDO
156			ELSE
157	jmc	1.11	c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
158	dimitri	1.15	DO j = 1-OLy, sNy+OLy
159			DO i = 1-OLx, sNx+OLx
160			surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
161			& + EmPmR(i,j,bi,bj)*convertFW2Salt
162			& recip_drF(kSurface)recip_hFacC(i,j,kSurface,bi,bj)
163			& *convertEmP2rUnit
164			ENDDO
165			ENDDO
166			ENDIF
167	jmc	1.11
168	dimitri	1.15	ENDIF
169	jmc	1.11	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
170	dimitri	1.12	#ifdef ALLOW_PTRACERS
171	dimitri	1.15	IF ( usePTRACERS ) THEN
172			CALL PTRACERS_FORCING_SURF(
173			I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
174			I myTime,myIter,myThid )
175			ENDIF
176	dimitri	1.12	#endif /* ALLOW_PTRACERS */
177	jmc	1.11
178			#ifdef ATMOSPHERIC_LOADING
179
180			C-- Atmospheric surface Pressure loading :
181
182	dimitri	1.15	IF (buoyancyRelation .eq. 'OCEANIC' ) THEN
183			DO j = 1-OLy, sNy+OLy
184			DO i = 1-OLx, sNx+OLx
185			phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
186			ENDDO
187	jmc	1.11	ENDDO
188	dimitri	1.15	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
189	jmc	1.11	C-- This is a hack used to read phi0surf from a file (ploadFile)
190			C instead of computing it from bathymetry & density ref. profile.
191			C The true atmospheric P-loading is not yet implemented for P-coord
192			C (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
193	dimitri	1.15	DO j = 1-OLy, sNy+OLy
194			DO i = 1-OLx, sNx+OLx
195			phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
196			ENDDO
197	jmc	1.11	ENDDO
198	dimitri	1.15	ENDIF
199	jmc	1.11
200			#endif /* ATMOSPHERIC_LOADING */
201	dimitri	1.15
202			ENDDO
203			ENDDO
204	heimbach	1.1
205			RETURN
206			END