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