/[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.37 by jmc, Tue May 22 02:35:58 2007 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
8  C     !ROUTINE: EXTERNAL_FORCING_SURF  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
16  C     | o Determines forcing terms based on external fields        C     | o Determines forcing terms based on external fields
17  C     |   relaxation terms etc.                                    C     |   relaxation terms etc.
18  C     *==========================================================*  C     *==========================================================*
19  C     \ev  C     \ev
20    
# 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    #ifdef ALLOW_SHELFICE
35    #include "SHELFICE.h"
36    #endif /* ALLOW_SHELFICE */
37  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
38  C     === Routine arguments ===  C     === Routine arguments ===
39    C     bi,bj  :: tile indices
40    C     iMin,iMax, jMin,jMax :: Range of points for calculation
41    C     myTime :: Current time in simulation
42    C     myIter :: Current iteration number in simulation
43  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
       INTEGER myThid  
44        INTEGER bi,bj        INTEGER bi,bj
45        INTEGER iMin, iMax        INTEGER iMin, iMax
46        INTEGER jMin, jMax        INTEGER jMin, jMax
47          _RL myTime
48          INTEGER myIter
49          INTEGER myThid
50    
51  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
52  C     === Local variables ===  C     === Local variables ===
53    C     i,j    :: loop indices
54    C     ks     :: index of surface interface layer
55        INTEGER i,j        INTEGER i,j
56          INTEGER ks
57  CEOP  CEOP
58    #ifdef ALLOW_DIAGNOSTICS
59          _RL tmpFac
60    #endif /* ALLOW_DIAGNOSTICS */
61    
62          IF ( usingPCoords ) THEN
63           ks        = Nr
64          ELSE
65           ks        = 1
66          ENDIF
67    
68    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
69    
70          IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
71    C--   Start with surface restoring term :
72    
73    #ifdef ALLOW_SEAICE
74           IF ( useSEAICE ) THEN
75    C     Do not restore under sea-ice
76            DO j = jMin, jMax
77             DO i = iMin, iMax
78    C     Heat Flux (restoring term) :
79              surfaceForcingT(i,j,bi,bj) =
80         &      -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,1,bi,bj))
81         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
82         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
83    C     Salt Flux (restoring term) :
84              surfaceForcingS(i,j,bi,bj) =
85         &      -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,1,bi,bj))
86         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
87         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
88             ENDDO
89            ENDDO
90           ELSE
91    #endif /* ALLOW_SEAICE */
92            DO j = jMin, jMax
93             DO i = iMin, iMax
94    C     Heat Flux (restoring term) :
95              surfaceForcingT(i,j,bi,bj) =
96         &      -lambdaThetaClimRelax(i,j,bi,bj)
97         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
98         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
99    C     Salt Flux (restoring term) :
100              surfaceForcingS(i,j,bi,bj) =
101         &      -lambdaSaltClimRelax(i,j,bi,bj)
102         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
103         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
104             ENDDO
105            ENDDO
106    #ifdef ALLOW_SEAICE
107           ENDIF
108    #endif /* ALLOW_SEAICE */
109    
110    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
111    #ifdef NONLIN_FRSURF
112    C-    T,S surface forcing will be applied (thermodynamics) after the update
113    C     of surf.thickness (hFac): account for change in surf.thickness
114           IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
115            IF (select_rStar.GT.0) THEN
116    # ifndef DISABLE_RSTAR_CODE
117             DO j=jMin,jMax
118              DO i=iMin,iMax
119                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
120         &                                  * rStarExpC(i,j,bi,bj)
121                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
122         &                                  * rStarExpC(i,j,bi,bj)
123              ENDDO
124             ENDDO
125    # endif /* DISABLE_RSTAR_CODE */
126            ELSE
127             DO j=jMin,jMax
128              DO i=iMin,iMax
129               IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
130                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
131         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
132                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
133         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
134               ENDIF
135              ENDDO
136             ENDDO
137            ENDIF
138           ENDIF
139    #endif /* NONLIN_FRSURF */
140    
141    #ifdef ALLOW_DIAGNOSTICS
142           IF ( useDiagnostics ) THEN
143    
144    C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
145            tmpFac = HeatCapacity_Cp*recip_horiVertRatio*rhoConst
146            CALL DIAGNOSTICS_SCALE_FILL(
147         &           surfaceForcingT(1-OLx,1-OLy,bi,bj),tmpFac,1,
148         &                             'TRELAX  ',0, 1,2,bi,bj,myThid)
149    
150    C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
151            tmpFac = recip_horiVertRatio*rhoConst
152            CALL DIAGNOSTICS_SCALE_FILL(
153         &           surfaceForcingS(1-OLx,1-OLy,bi,bj),tmpFac,1,
154         &                             'SRELAX  ',0, 1,2,bi,bj,myThid)
155    
156           ENDIF
157    #endif /* ALLOW_DIAGNOSTICS */
158    
159          ELSE
160    C--   No restoring for T & S : set surfaceForcingT,S to zero :
161    
162           DO j = jMin, jMax
163            DO i = iMin, iMax
164              surfaceForcingT(i,j,bi,bj) = 0. _d 0
165              surfaceForcingS(i,j,bi,bj) = 0. _d 0
166            ENDDO
167           ENDDO
168    
169    C--   end restoring / no restoring block.
170          ENDIF
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    C-    Apply mask on Fresh-Water flux
203    C      (needed for SSH forcing, whether or not exactConserv is used)
204          IF ( useRealFreshWaterFlux ) THEN
205            DO j=1-OLy,sNy+OLy
206             DO i=1-OLx,sNx+OLx
207               EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*maskH(i,j,bi,bj)
208             ENDDO
209            ENDDO
210          ENDIF
211    
212    #ifdef EXACT_CONSERV
213    c NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
214    c     to stay consitent with volume change (=d/dt etaH).
215          IF ( staggerTimeStep ) THEN
216            DO j=1-OLy,sNy+OLy
217             DO i=1-OLx,sNx+OLx
218               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
219             ENDDO
220            ENDDO
221          ENDIF
222    
223          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
224         &     .AND. useRealFreshWaterFlux ) THEN
225    
226    c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
227    c   the water column height ; temp., salt, (tracer) flux associated
228    c   with this input/output of water is added here to the surface tendency.
229    c
230    
231           IF (temp_EvPrRn.NE.UNSET_RL) THEN
232            DO j = jMin, jMax
233             DO i = iMin, iMax
234              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
235         &      + PmEpR(i,j,bi,bj)
236         &         *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
237         &         *convertEmP2rUnit
238             ENDDO
239            ENDDO
240           ENDIF
241    
242           IF (salt_EvPrRn.NE.UNSET_RL) THEN
243            DO j = jMin, jMax
244             DO i = iMin, iMax
245              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
246         &      + PmEpR(i,j,bi,bj)
247         &         *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
248         &         *convertEmP2rUnit
249             ENDDO
250            ENDDO
251           ENDIF
252    
253    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
254          ELSE
255    #else /* EXACT_CONSERV */
256          IF (.TRUE.) THEN
257    #endif /* EXACT_CONSERV */
258    
259    c- EmPmR does not really affect the water column height (for tracer budget)
260    c   and is converted to a salt tendency.
261    
262           IF (convertFW2Salt .EQ. -1.) THEN
263    c- converts EmPmR to salinity tendency using surface local salinity
264            DO j = jMin, jMax
265             DO i = iMin, iMax
266              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
267         &      + EmPmR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
268         &         *convertEmP2rUnit
269             ENDDO
270            ENDDO
271           ELSE
272    c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
273            DO j = jMin, jMax
274             DO i = iMin, iMax
275              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
276         &      + EmPmR(i,j,bi,bj)*convertFW2Salt
277         &         *convertEmP2rUnit
278             ENDDO
279            ENDDO
280           ENDIF
281    
282          ENDIF
283    
284    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
285    
286    #ifdef ALLOW_PTRACERS
287          IF ( usePTRACERS ) THEN
288             CALL PTRACERS_FORCING_SURF(
289         I        bi, bj, iMin, iMax, jMin, jMax,
290         I        myTime,myIter,myThid )
291          ENDIF
292    #endif /* ALLOW_PTRACERS */
293    
294    #ifdef ATMOSPHERIC_LOADING
295    
296    C-- Atmospheric surface Pressure loading :
297    
298          IF ( usingZCoords ) THEN
299           IF ( useRealFreshWaterFlux ) THEN
300            DO j = jMin, jMax
301             DO i = iMin, iMax
302              phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
303         &                        +sIceLoad(i,j,bi,bj)*gravity
304         &                          )*recip_rhoConst
305             ENDDO
306            ENDDO
307           ELSE
308            DO j = jMin, jMax
309             DO i = iMin, iMax
310              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
311             ENDDO
312            ENDDO
313           ENDIF
314          ELSEIF ( usingPCoords ) THEN
315    C-- This is a hack used to read phi0surf from a file (ploadFile)
316    C   instead of computing it from bathymetry & density ref. profile.
317    C   The true atmospheric P-loading is not yet implemented for P-coord
318    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
319            DO j = jMin, jMax
320             DO i = iMin, iMax
321              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
322             ENDDO
323            ENDDO
324          ENDIF
325    
326    #endif /* ATMOSPHERIC_LOADING */
327    
328    #ifdef ALLOW_SHELFICE
329          IF ( usingZCoords ) THEN
330           IF ( useSHELFICE) THEN
331            DO j = jMin, jMax
332             DO i = iMin, iMax
333              phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj)
334         &         + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst
335             ENDDO
336            ENDDO
337           ENDIF
338          ENDIF
339    #endif /* ALLOW_SHELFICE */
340    
341    #ifdef ALLOW_EBM
342    c--    Values for surfaceForcingT, surfaceForcingS
343    c      are overwritten by those produced by EBM
344    cph    AD recomputation problems if these IF useEBM are used
345    cph      IF ( useEBM ) THEN
346           CALL EBM_FORCING_SURF(
347         I        bi, bj, iMin, iMax, jMin, jMax,
348         I        myTime,myIter,myThid )
349    cph      ENDIF
350    #endif
351    
352        RETURN        RETURN
353        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22