/[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.29 by dimitri, Fri Dec 17 19:17:56 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     bi,bj  :: tile indices
38    C     iMin,iMax, jMin,jMax :: Range of points for calculation
39    C     myTime :: Current time in simulation
40    C     myIter :: Current iteration number in simulation
41  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
42          _RL myTime
43          INTEGER myIter
44        INTEGER myThid        INTEGER myThid
45        INTEGER bi,bj        INTEGER bi,bj
46        INTEGER iMin, iMax        INTEGER iMin, iMax
# Line 38  C     myThid :: Thread no. that called t Line 49  C     myThid :: Thread no. that called t
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     === Local variables ===  C     === Local variables ===
51        INTEGER i,j        INTEGER i,j
52    C     number of surface interface layer
53          INTEGER ks
54    #ifdef ALLOW_DIAGNOSTICS
55          LOGICAL  DIAGNOSTICS_IS_ON
56          EXTERNAL DIAGNOSTICS_IS_ON
57          _RL tmp1k(1:sNx,1:sNy)
58    #endif /* ALLOW_DIAGNOSTICS */
59  CEOP  CEOP
60    
61          IF ( usingPCoords ) THEN
62           ks        = Nr
63          ELSE
64           ks        = 1
65          ENDIF
66    
67    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
68    
69          IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
70    C--   Start with surface restoring term :
71    
72           DO j = jMin, jMax
73            DO i = iMin, iMax
74    #ifdef ALLOW_SEAICE
75    C     Do not restore under sea-ice
76    C     Heat Flux (restoring term) :
77              surfaceForcingT(i,j,bi,bj) =
78         &      -lambdaThetaClimRelax * (1-AREA(i,j,1,bi,bj))
79         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
80         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
81    C     Salt Flux (restoring term) :
82              surfaceForcingS(i,j,bi,bj) =
83         &      -lambdaSaltClimRelax * (1-AREA(i,j,1,bi,bj))
84         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
85         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
86    #else /* ifndef ALLOW_SEAICE */
87    C     Heat Flux (restoring term) :
88             IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
89              surfaceForcingT(i,j,bi,bj) =
90         &      -lambdaThetaClimRelax
91         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
92         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
93    C     Salt Flux (restoring term) :
94              surfaceForcingS(i,j,bi,bj) =
95         &      -lambdaSaltClimRelax
96         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
97         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
98             ELSE
99              surfaceForcingT(i,j,bi,bj) = 0. _d 0
100              surfaceForcingS(i,j,bi,bj) = 0. _d 0
101             ENDIF
102    #endif /* ALLOW_SEAICE */
103            ENDDO
104           ENDDO
105    
106    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
107    #ifdef NONLIN_FRSURF
108    C-    T,S surface forcing will be applied (thermodynamics) after the update
109    C     of surf.thickness (hFac): account for change in surf.thickness
110           IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
111            IF (select_rStar.GT.0) THEN
112             DO j=jMin,jMax
113              DO i=iMin,iMax
114                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
115         &                                  * rStarExpC(i,j,bi,bj)
116                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
117         &                                  * rStarExpC(i,j,bi,bj)
118              ENDDO
119             ENDDO
120            ELSE
121             DO j=jMin,jMax
122              DO i=iMin,iMax
123               IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
124                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
125         &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
126                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
127         &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
128               ENDIF
129              ENDDO
130             ENDDO
131            ENDIF
132           ENDIF
133    #endif /* NONLIN_FRSURF */
134    
135          ELSE
136    C--   No restoring for T & S : set surfaceForcingT,S to zero :
137    
138           DO j = jMin, jMax
139            DO i = iMin, iMax
140              surfaceForcingT(i,j,bi,bj) = 0. _d 0
141              surfaceForcingS(i,j,bi,bj) = 0. _d 0
142            ENDDO
143           ENDDO
144    
145    C--   end restoring / no restoring block.
146          ENDIF
147    
148    #ifdef ALLOW_DIAGNOSTICS
149    C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
150          IF ( DIAGNOSTICS_IS_ON('TRELAX  ',myThid) ) THEN
151             DO j = 1,sNy
152                DO i = 1,sNx
153                   tmp1k(i,j) = surfaceForcingT(i,j,bi,bj)
154         &              *HeatCapacity_Cp*recip_horiVertRatio*rhoConst
155                ENDDO
156             ENDDO
157             CALL DIAGNOSTICS_FILL(tmp1k,'TRELAX  ',0,1,3,bi,bj,myThid)
158          ENDIF
159    
160    C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
161          IF ( DIAGNOSTICS_IS_ON('SRELAX  ',myThid) ) THEN
162             DO j = 1,sNy
163                DO i = 1,sNx
164                   tmp1k(i,j) = surfaceForcingS(i,j,bi,bj)*
165         &              recip_horiVertRatio*rhoConst
166                ENDDO
167             ENDDO
168             CALL DIAGNOSTICS_FILL(tmp1k,'SRELAX  ',0,1,3,bi,bj,myThid)
169          ENDIF
170    #endif /* ALLOW_DIAGNOSTICS */
171    
172    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
173    
174    C--   Surface Fluxes :
175    
176        DO j = jMin, jMax        DO j = jMin, jMax
177           DO i = iMin, iMax           DO i = iMin, iMax
178    
179  c     Zonal wind stress fu:  C     Zonal wind stress fu:
180            surfaceTendencyU(i,j,bi,bj) = fu(i,j,bi,bj)            surfaceForcingU(i,j,bi,bj) =
181       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
182  c     Meridional wind stress fv:  C     Meridional wind stress fv:
183            surfaceTendencyV(i,j,bi,bj) = fV(i,j,bi,bj)            surfaceForcingV(i,j,bi,bj) =
184       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
185  c     Net heat flux Qnet:  C     Net heat flux Qnet:
186            surfaceTendencyT(i,j,bi,bj) = -Qnet(i,j,bi,bj)            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
187       &        *recip_Cp*recip_rhoNil*recip_dRf(1)       &       - ( Qnet(i,j,bi,bj)
188       &         - lambdaThetaClimRelax*  #ifdef SHORTWAVE_HEATING
189       &           (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))       &          -Qsw(i,j,bi,bj)
190    #endif
191  #ifdef USE_NATURAL_BCS       &         ) *recip_Cp*horiVertRatio*recip_rhoConst
192  c     Freshwater flux EmPmR:  C     Net Salt Flux :
193            surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
194       &         *recip_dRf(1)*salt(i,j,1,bi,bj)       &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst
      &         - lambdaSaltClimRelax*  
      &           (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))  
 #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            
195    
196           ENDDO           ENDDO
197        ENDDO        ENDDO
198            
199    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
200    C--   Fresh-water flux
201    
202    #ifdef EXACT_CONSERV
203    c NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
204    c     to stay consitent with volume change (=d/dt etaH).
205          IF ( staggerTimeStep ) THEN
206            DO j=1,sNy
207             DO i=1,sNx
208               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
209             ENDDO
210            ENDDO
211          ENDIF
212    
213          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
214         &     .AND. useRealFreshWaterFlux ) THEN
215    
216    c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
217    c   the water column height ; temp., salt, (tracer) flux associated
218    c   with this input/output of water is added here to the surface tendency.
219    c
220    
221           IF (temp_EvPrRn.NE.UNSET_RL) THEN
222            DO j = jMin, jMax
223             DO i = iMin, iMax
224              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
225         &      + PmEpR(i,j,bi,bj)
226         &         *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
227         &         *convertEmP2rUnit
228             ENDDO
229            ENDDO
230           ENDIF
231    
232           IF (salt_EvPrRn.NE.UNSET_RL) THEN
233            DO j = jMin, jMax
234             DO i = iMin, iMax
235              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
236         &      + PmEpR(i,j,bi,bj)
237         &         *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
238         &         *convertEmP2rUnit
239             ENDDO
240            ENDDO
241           ENDIF
242    
243    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
244          ELSE
245    #else /* EXACT_CONSERV */
246          IF (.TRUE.) THEN
247    #endif /* EXACT_CONSERV */
248    
249    c- EmPmR does not really affect the water column height (for tracer budget)
250    c   and is converted to a salt tendency.
251    
252           IF (convertFW2Salt .EQ. -1.) THEN
253    c- converts EmPmR to salinity tendency using surface local salinity
254            DO j = jMin, jMax
255             DO i = iMin, iMax
256              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
257         &      + EmPmR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
258         &         *convertEmP2rUnit
259             ENDDO
260            ENDDO
261           ELSE
262    c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
263            DO j = jMin, jMax
264             DO i = iMin, iMax
265              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
266         &      + EmPmR(i,j,bi,bj)*convertFW2Salt
267         &         *convertEmP2rUnit
268             ENDDO
269            ENDDO
270           ENDIF
271    
272          ENDIF
273    
274    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
275    
276    #ifdef ALLOW_PTRACERS
277          IF ( usePTRACERS ) THEN
278             CALL PTRACERS_FORCING_SURF(
279         I        bi, bj, iMin, iMax, jMin, jMax,
280         I        myTime,myIter,myThid )
281          ENDIF
282    #endif /* ALLOW_PTRACERS */
283    
284    #ifdef ATMOSPHERIC_LOADING
285    
286    C-- Atmospheric surface Pressure loading :
287    
288          IF ( usingZCoords ) THEN
289           IF ( useRealFreshWaterFlux ) THEN
290            DO j = jMin, jMax
291             DO i = iMin, iMax
292              phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
293         &                        +sIceLoad(i,j,bi,bj)*gravity
294         &                          )*recip_rhoConst
295             ENDDO
296            ENDDO
297           ELSE
298            DO j = jMin, jMax
299             DO i = iMin, iMax
300              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
301             ENDDO
302            ENDDO
303           ENDIF
304          ELSEIF ( usingPCoords ) THEN
305    C-- This is a hack used to read phi0surf from a file (ploadFile)
306    C   instead of computing it from bathymetry & density ref. profile.
307    C   The true atmospheric P-loading is not yet implemented for P-coord
308    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
309            DO j = jMin, jMax
310             DO i = iMin, iMax
311              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
312             ENDDO
313            ENDDO
314          ENDIF
315    
316    #endif /* ATMOSPHERIC_LOADING */
317    
318    #ifdef ALLOW_EBM
319    c--    Values for surfaceForcingT, surfaceForcingS
320    c      are overwritten by those produced by EBM
321    cph    AD recomputation problems if these IF useEBM are used
322    cph      IF ( useEBM ) THEN
323           CALL EBM_FORCING_SURF(
324         I        bi, bj, iMin, iMax, jMin, jMax,
325         I        myTime,myIter,myThid )
326    cph      ENDIF
327    #endif
328    
329        RETURN        RETURN
330        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22