/[MITgcm]/MITgcm/model/src/external_forcing_surf.F
ViewVC logotype

Diff of /MITgcm/model/src/external_forcing_surf.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.5 by cnh, Wed Sep 26 18:09:14 2001 UTC revision 1.20 by heimbach, Fri May 14 21:08:28 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 8  C     !ROUTINE: EXTERNAL_FORCING_SURF Line 9  C     !ROUTINE: EXTERNAL_FORCING_SURF
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE EXTERNAL_FORCING_SURF(        SUBROUTINE EXTERNAL_FORCING_SURF(
11       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
12       I             myThid )       I             myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | SUBROUTINE EXTERNAL_FORCING_SURF                            C     | SUBROUTINE EXTERNAL_FORCING_SURF                          
# Line 26  C     === Global variables === Line 27  C     === Global variables ===
27  #include "FFIELDS.h"  #include "FFIELDS.h"
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "GRID.h"  #include "GRID.h"
30    #include "SURFACE.h"
31    #ifdef ALLOW_SEAICE
32    #include "SEAICE.h"
33    #endif /* ALLOW_SEAICE */
34    
35  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
36  C     === Routine arguments ===  C     === Routine arguments ===
37    C     myTime - Current time in simulation
38    C     myIter - Current iteration number in simulation
39  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
40          _RL myTime
41          INTEGER myIter
42        INTEGER myThid        INTEGER myThid
43        INTEGER bi,bj        INTEGER bi,bj
44        INTEGER iMin, iMax        INTEGER iMin, iMax
# Line 38  C     myThid :: Thread no. that called t Line 47  C     myThid :: Thread no. that called t
47  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
48  C     === Local variables ===  C     === Local variables ===
49        INTEGER i,j        INTEGER i,j
50    C     number of surface interface layer
51          INTEGER kSurface
52  CEOP  CEOP
53    
54          if ( buoyancyRelation .eq. 'OCEANICP' ) then
55           kSurface        = Nr
56          else
57           kSurface        = 1
58          endif
59    
60    C--   Surface Fluxes :
61    
62        DO j = jMin, jMax        DO j = jMin, jMax
63           DO i = iMin, iMax           DO i = iMin, iMax
64    
65  c     Zonal wind stress fu:  c     Zonal wind stress fu:
66            surfaceTendencyU(i,j,bi,bj) = fu(i,j,bi,bj)            surfaceTendencyU(i,j,bi,bj) =
67       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
68         &           *recip_drF(kSurface)*recip_hFacW(i,j,kSurface,bi,bj)
69  c     Meridional wind stress fv:  c     Meridional wind stress fv:
70            surfaceTendencyV(i,j,bi,bj) = fV(i,j,bi,bj)            surfaceTendencyV(i,j,bi,bj) =
71       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
72         &           *recip_drF(kSurface)*recip_hFacS(i,j,kSurface,bi,bj)
73  c     Net heat flux Qnet:  c     Net heat flux Qnet:
74            surfaceTendencyT(i,j,bi,bj) = -Qnet(i,j,bi,bj)            surfaceTendencyT(i,j,bi,bj) =
75       &        *recip_Cp*recip_rhoNil*recip_dRf(1)       &      -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst
76       &         - lambdaThetaClimRelax*       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
77       &           (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))  C     Net Salt Flux :
78              surfaceTendencyS(i,j,bi,bj) =
79  #ifdef USE_NATURAL_BCS       &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst
80  c     Freshwater flux EmPmR:       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
81            surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)  
82       &         *recip_dRf(1)*salt(i,j,1,bi,bj)  #ifdef ALLOW_PASSIVE_TRACER
83       &         - lambdaSaltClimRelax*  c ***  define the tracer surface tendency here ***
84       &           (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))  #endif /* ALLOW_PASSIVE_TRACER */
 #else  
 c     Freshwater flux EmPmR:  
           surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)  
      &         *recip_dRf(1)*35.  
      &         - lambdaSaltClimRelax*  
      &           (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))  
 #endif            
85    
86           ENDDO           ENDDO
87        ENDDO        ENDDO
88            
89    C--   Surface restoring term :
90    
91          IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
92           DO j = jMin, jMax
93            DO i = iMin, iMax
94    #ifdef ALLOW_SEAICE
95    C     Don't restore under sea-ice
96    C     Heat Flux (restoring term) :
97              surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
98         &      -lambdaThetaClimRelax * (1-AREA(i,j,1,bi,bj))
99         &         *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))
100    C     Salt Flux (restoring term) :
101              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
102         &      -lambdaSaltClimRelax * (1-AREA(i,j,1,bi,bj))
103         &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))
104    #else /* ifndef ALLOW_SEAICE */
105    C     Heat Flux (restoring term) :
106             IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
107              surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
108         &      -lambdaThetaClimRelax
109         &         *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))
110    C     Salt Flux (restoring term) :
111              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
112         &      -lambdaSaltClimRelax
113         &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))
114             ENDIF
115    #endif /* ALLOW_SEAICE */
116            ENDDO
117           ENDDO
118          ENDIF
119    
120    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
121    C--   Fresh-water flux
122    
123    #ifdef NONLIN_FRSURF
124          IF ( (nonlinFreeSurf.GT.0 .OR. buoyancyRelation.EQ.'OCEANICP')
125         &     .AND. useRealFreshWaterFlux ) THEN
126    
127    c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
128    c   the water column height ; temp., salt, (tracer) flux associated
129    c   with this input/output of water is added here to the surface tendency.
130    c
131    c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay
132    c     consitent with volume change (=d/dt etaN).
133    
134           IF (temp_EvPrRn.NE.UNSET_RL) THEN
135            DO j = jMin, jMax
136             DO i = iMin, iMax
137              surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)
138         &      + PmEpR(i,j,bi,bj)
139         &         *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )
140         &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
141         &         *convertEmP2rUnit
142             ENDDO
143            ENDDO
144           ENDIF
145    
146           IF (salt_EvPrRn.NE.UNSET_RL) THEN
147            DO j = jMin, jMax
148             DO i = iMin, iMax
149              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
150         &      + PmEpR(i,j,bi,bj)
151         &         *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) )
152         &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
153         &         *convertEmP2rUnit
154             ENDDO
155            ENDDO
156           ENDIF
157    
158    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
159          ELSE
160    #else /* NONLIN_FRSURF */
161          IF (.TRUE.) THEN
162    #endif /* NONLIN_FRSURF */
163    
164    c- EmPmR does not really affect the water column height (for tracer budget)
165    c   and is converted to a salt tendency.
166    
167           IF (convertFW2Salt .EQ. -1.) THEN
168    c- converts EmPmR to salinity tendency using surface local salinity
169            DO j = jMin, jMax
170             DO i = iMin, iMax
171              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
172         &      + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)
173         &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
174         &         *convertEmP2rUnit
175             ENDDO
176            ENDDO
177           ELSE
178    c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
179            DO j = jMin, jMax
180             DO i = iMin, iMax
181              surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)
182         &      + EmPmR(i,j,bi,bj)*convertFW2Salt
183         &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
184         &         *convertEmP2rUnit
185             ENDDO
186            ENDDO
187           ENDIF
188    
189          ENDIF
190    
191    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192    
193    #ifdef ALLOW_PTRACERS
194          IF ( usePTRACERS ) THEN
195             CALL PTRACERS_FORCING_SURF(
196         I        bi, bj, iMin, iMax, jMin, jMax,
197         I        myTime,myIter,myThid )
198          ENDIF
199    #endif /* ALLOW_PTRACERS */
200    
201    #ifdef ATMOSPHERIC_LOADING
202    
203    C-- Atmospheric surface Pressure loading :
204    
205          IF (buoyancyRelation .eq. 'OCEANIC' ) THEN
206            DO j = jMin, jMax
207             DO i = iMin, iMax
208              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
209             ENDDO
210            ENDDO
211          ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
212    C-- This is a hack used to read phi0surf from a file (ploadFile)
213    C   instead of computing it from bathymetry & density ref. profile.
214    C   The true atmospheric P-loading is not yet implemented for P-coord
215    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
216            DO j = jMin, jMax
217             DO i = iMin, iMax
218              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
219             ENDDO
220            ENDDO
221          ENDIF
222    
223    #endif /* ATMOSPHERIC_LOADING */
224    
225    #ifdef ALLOW_EBM
226    c--    Values for surfaceTendencyT, surfaceTendencyS
227    c      are overwritten by those produced by EBM
228          IF ( useEBM ) THEN
229           CALL EBM_FORCING_SURF(
230         I        bi, bj, iMin, iMax, jMin, jMax,
231         I        myTime,myIter,myThid )
232          ENDIF
233    #endif
234    
235    #ifdef ALLOW_TIMEAVE
236          IF ( taveFreq .NE. 0. _d 0 ) THEN
237            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
238          ENDIF
239    #endif /* ALLOW_TIMEAVE */
240    
241        RETURN        RETURN
242        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22