/[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.44 by dimitri, Wed Nov 28 00:18:17 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    
38  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
39  C     === Routine arguments ===  C     === Routine arguments ===
40    C     bi,bj  :: tile indices
41    C     iMin,iMax, jMin,jMax :: Range of points for calculation
42    C     myTime :: Current time in simulation
43    C     myIter :: Current iteration number in simulation
44  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
       INTEGER myThid  
45        INTEGER bi,bj        INTEGER bi,bj
46        INTEGER iMin, iMax        INTEGER iMin, iMax
47        INTEGER jMin, jMax        INTEGER jMin, jMax
48          _RL myTime
49          INTEGER myIter
50          INTEGER myThid
51    
52  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
53  C     === Local variables ===  C     === Local variables ===
54    C     i,j    :: loop indices
55    C     ks     :: index of surface interface layer
56        INTEGER i,j        INTEGER i,j
57          INTEGER ks
58  CEOP  CEOP
59    #ifdef ALLOW_DIAGNOSTICS
60          _RL tmpFac
61    #endif /* ALLOW_DIAGNOSTICS */
62    
63          IF ( usingPCoords ) THEN
64           ks        = Nr
65          ELSE
66           ks        = 1
67          ENDIF
68    
69    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
70    
71          IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
72    C--   Start with surface restoring term :
73    
74    #ifdef ALLOW_SEAICE
75           IF ( useSEAICE ) THEN
76    C     Do not restore under sea-ice
77            DO j = jMin, jMax
78             DO i = iMin, iMax
79    C     Heat Flux (restoring term) :
80              surfaceForcingT(i,j,bi,bj) =
81         &      -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,1,bi,bj))
82         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
83         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
84    C     Salt Flux (restoring term) :
85              surfaceForcingS(i,j,bi,bj) =
86         &      -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,1,bi,bj))
87         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
88         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
89             ENDDO
90            ENDDO
91           ELSE
92    #endif /* ALLOW_SEAICE */
93            DO j = jMin, jMax
94             DO i = iMin, iMax
95    C     Heat Flux (restoring term) :
96              surfaceForcingT(i,j,bi,bj) =
97         &      -lambdaThetaClimRelax(i,j,bi,bj)
98         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
99         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
100    C     Salt Flux (restoring term) :
101              surfaceForcingS(i,j,bi,bj) =
102         &      -lambdaSaltClimRelax(i,j,bi,bj)
103         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
104         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
105             ENDDO
106            ENDDO
107    #ifdef ALLOW_SEAICE
108           ENDIF
109    #endif /* ALLOW_SEAICE */
110    
111    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
112    #ifdef NONLIN_FRSURF
113    C-    T,S surface forcing will be applied (thermodynamics) after the update
114    C     of surf.thickness (hFac): account for change in surf.thickness
115           IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
116            IF (select_rStar.GT.0) THEN
117    # ifndef DISABLE_RSTAR_CODE
118             DO j=jMin,jMax
119              DO i=iMin,iMax
120                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
121         &                                  * rStarExpC(i,j,bi,bj)
122                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
123         &                                  * rStarExpC(i,j,bi,bj)
124              ENDDO
125             ENDDO
126    # endif /* DISABLE_RSTAR_CODE */
127            ELSE
128             DO j=jMin,jMax
129              DO i=iMin,iMax
130               IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
131                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
132         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
133                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
134         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
135               ENDIF
136              ENDDO
137             ENDDO
138            ENDIF
139           ENDIF
140    #endif /* NONLIN_FRSURF */
141    
142    #ifdef ALLOW_DIAGNOSTICS
143           IF ( useDiagnostics ) THEN
144    
145    C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
146            tmpFac = HeatCapacity_Cp*rUnit2mass
147            CALL DIAGNOSTICS_SCALE_FILL(
148         &           surfaceForcingT(1-OLx,1-OLy,bi,bj),tmpFac,1,
149         &                             'TRELAX  ',0, 1,2,bi,bj,myThid)
150    
151    C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
152            tmpFac = rUnit2mass
153            CALL DIAGNOSTICS_SCALE_FILL(
154         &           surfaceForcingS(1-OLx,1-OLy,bi,bj),tmpFac,1,
155         &                             'SRELAX  ',0, 1,2,bi,bj,myThid)
156    
157           ENDIF
158    #endif /* ALLOW_DIAGNOSTICS */
159    
160          ELSE
161    C--   No restoring for T & S : set surfaceForcingT,S to zero :
162    
163           DO j = jMin, jMax
164            DO i = iMin, iMax
165              surfaceForcingT(i,j,bi,bj) = 0. _d 0
166              surfaceForcingS(i,j,bi,bj) = 0. _d 0
167            ENDDO
168           ENDDO
169    
170    C--   end restoring / no restoring block.
171          ENDIF
172    
173    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
174    
175    C--   Surface Fluxes :
176    
177        DO j = jMin, jMax        DO j = jMin, jMax
178           DO i = iMin, iMax           DO i = iMin, iMax
179    
180  c     Zonal wind stress fu:  C     Zonal wind stress fu:
181            surfaceTendencyU(i,j,bi,bj) = fu(i,j,bi,bj)            surfaceForcingU(i,j,bi,bj) =
182       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fu(i,j,bi,bj)*mass2rUnit
183  c     Meridional wind stress fv:  C     Meridional wind stress fv:
184            surfaceTendencyV(i,j,bi,bj) = fV(i,j,bi,bj)            surfaceForcingV(i,j,bi,bj) =
185       &         *horiVertRatio*recip_rhoNil*recip_dRf(1)       &      fv(i,j,bi,bj)*mass2rUnit
186  c     Net heat flux Qnet:  C     Net heat flux Qnet:
187            surfaceTendencyT(i,j,bi,bj) = -Qnet(i,j,bi,bj)            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
188       &        *recip_Cp*recip_rhoNil*recip_dRf(1)       &       - ( Qnet(i,j,bi,bj)
189       &         - lambdaThetaClimRelax*  #ifdef SHORTWAVE_HEATING
190       &           (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))       &          -Qsw(i,j,bi,bj)
191    #endif
192  #ifdef USE_NATURAL_BCS       &         ) *recip_Cp*mass2rUnit
193  c     Freshwater flux EmPmR:  C     Net Salt Flux :
194            surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
195       &         *recip_dRf(1)*salt(i,j,1,bi,bj)       &      -saltFlux(i,j,bi,bj)*mass2rUnit
      &         - 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            
196    
197           ENDDO           ENDDO
198        ENDDO        ENDDO
199    
200    #ifdef ALLOW_SALT_PLUME
201    C saltPlume is the amount of salt rejected by ice while freezing;
202    C it is here subtracted from surfaceForcingS and will be redistributed
203    C to multiple vertical levels later on as per Duffy et al. (GRL 1999)
204          IF ( useSaltPlume ) THEN
205             CALL SALT_PLUME_FORCING_SURF(
206         I        bi, bj, iMin, iMax, jMin, jMax,
207         I        myTime,myIter,myThid )
208          ENDIF
209    #endif /* ALLOW_SALT_PLUME */
210    
211    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
212    C--   Fresh-water flux
213    
214    C-    Apply mask on Fresh-Water flux
215    C      (needed for SSH forcing, whether or not exactConserv is used)
216          IF ( useRealFreshWaterFlux ) THEN
217            DO j=1-OLy,sNy+OLy
218             DO i=1-OLx,sNx+OLx
219               EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*maskH(i,j,bi,bj)
220             ENDDO
221            ENDDO
222          ENDIF
223    
224    #ifdef EXACT_CONSERV
225    C NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
226    C     to stay consitent with volume change (=d/dt etaH).
227          IF ( staggerTimeStep ) THEN
228            DO j=1-OLy,sNy+OLy
229             DO i=1-OLx,sNx+OLx
230               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
231             ENDDO
232            ENDDO
233          ENDIF
234    
235          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
236         &     .AND. useRealFreshWaterFlux ) THEN
237    
238    C--   NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
239    C     the water column height ; temp., salt, (tracer) flux associated
240    C     with this input/output of water is added here to the surface tendency.
241    
242           IF (temp_EvPrRn.NE.UNSET_RL) THEN
243            DO j = jMin, jMax
244             DO i = iMin, iMax
245              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
246         &      + PmEpR(i,j,bi,bj)
247         &          *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
248         &          *mass2rUnit
249             ENDDO
250            ENDDO
251           ENDIF
252    
253           IF (salt_EvPrRn.NE.UNSET_RL) THEN
254            DO j = jMin, jMax
255             DO i = iMin, iMax
256              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
257         &      + PmEpR(i,j,bi,bj)
258         &          *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
259         &          *mass2rUnit
260             ENDDO
261            ENDDO
262           ENDIF
263    
264    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
265          ELSE
266    #else /* EXACT_CONSERV */
267          IF (.TRUE.) THEN
268    #endif /* EXACT_CONSERV */
269    
270    C--   EmPmR does not really affect the water column height (for tracer budget)
271    C     and is converted to a salt tendency.
272    
273           IF (convertFW2Salt .EQ. -1.) THEN
274    C-    use local surface tracer field to calculate forcing term:
275    
276            IF (temp_EvPrRn.NE.UNSET_RL) THEN
277    C     account for Rain/Evap heat content (temp_EvPrRn) using local SST
278             DO j = jMin, jMax
279              DO i = iMin, iMax
280               surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
281         &       + EmPmR(i,j,bi,bj)
282         &           *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
283         &           *mass2rUnit
284              ENDDO
285             ENDDO
286            ENDIF
287            IF (salt_EvPrRn.NE.UNSET_RL) THEN
288    C     converts EmPmR to salinity tendency using surface local salinity
289             DO j = jMin, jMax
290              DO i = iMin, iMax
291               surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
292         &       + EmPmR(i,j,bi,bj)
293         &           *( salt(i,j,ks,bi,bj) - salt_EvPrRn )
294         &           *mass2rUnit
295              ENDDO
296             ENDDO
297            ENDIF
298    
299           ELSE
300    C-    use uniform tracer value to calculate forcing term:
301    
302            IF (temp_EvPrRn.NE.UNSET_RL) THEN
303    C     account for Rain/Evap heat content (temp_EvPrRn) assuming uniform SST (=tRef)
304             DO j = jMin, jMax
305              DO i = iMin, iMax
306               surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
307         &       + EmPmR(i,j,bi,bj)
308         &           *( tRef(ks) - temp_EvPrRn )
309         &           *mass2rUnit
310              ENDDO
311             ENDDO
312            ENDIF
313            IF (salt_EvPrRn.NE.UNSET_RL) THEN
314    C     converts EmPmR to virtual salt flux using uniform salinity (default=35)
315             DO j = jMin, jMax
316              DO i = iMin, iMax
317               surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
318         &       + EmPmR(i,j,bi,bj)
319         &           *( convertFW2Salt - salt_EvPrRn )
320         &           *mass2rUnit
321              ENDDO
322             ENDDO
323            ENDIF
324    
325    C-    end local-surface-tracer / uniform-value distinction
326           ENDIF
327    
328          ENDIF
329    
330    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
331    
332    #ifdef ALLOW_PTRACERS
333          IF ( usePTRACERS ) THEN
334             CALL PTRACERS_FORCING_SURF(
335         I        bi, bj, iMin, iMax, jMin, jMax,
336         I        myTime,myIter,myThid )
337          ENDIF
338    #endif /* ALLOW_PTRACERS */
339    
340    #ifdef ATMOSPHERIC_LOADING
341    
342    C-- Atmospheric surface Pressure loading :
343    
344          IF ( usingZCoords ) THEN
345           IF ( useRealFreshWaterFlux ) THEN
346            DO j = jMin, jMax
347             DO i = iMin, iMax
348              phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
349         &                        +sIceLoad(i,j,bi,bj)*gravity
350         &                          )*recip_rhoConst
351             ENDDO
352            ENDDO
353           ELSE
354            DO j = jMin, jMax
355             DO i = iMin, iMax
356              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
357             ENDDO
358            ENDDO
359           ENDIF
360          ELSEIF ( usingPCoords ) THEN
361    C-- This is a hack used to read phi0surf from a file (ploadFile)
362    C   instead of computing it from bathymetry & density ref. profile.
363    C   The true atmospheric P-loading is not yet implemented for P-coord
364    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
365            DO j = jMin, jMax
366             DO i = iMin, iMax
367              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
368             ENDDO
369            ENDDO
370          ENDIF
371    
372    #endif /* ATMOSPHERIC_LOADING */
373    
374    #ifdef ALLOW_SHELFICE
375          IF ( usingZCoords ) THEN
376           IF ( useSHELFICE) THEN
377            DO j = jMin, jMax
378             DO i = iMin, iMax
379              phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj)
380         &         + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst
381             ENDDO
382            ENDDO
383           ENDIF
384          ENDIF
385    #endif /* ALLOW_SHELFICE */
386    
387    #ifdef ALLOW_EBM
388    c--    Values for surfaceForcingT, surfaceForcingS
389    c      are overwritten by those produced by EBM
390    cph    AD recomputation problems if these IF useEBM are used
391    cph      IF ( useEBM ) THEN
392           CALL EBM_FORCING_SURF(
393         I        bi, bj, iMin, iMax, jMin, jMax,
394         I        myTime,myIter,myThid )
395    cph      ENDIF
396    #endif
397    
398        RETURN        RETURN
399        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22