/[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.1 by heimbach, Mon Sep 11 20:45:57 2000 UTC revision 1.30 by heimbach, Wed Apr 6 20:18:19 2005 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8        SUBROUTINE EXTERNAL_FORCING_SURF_U(  C     !ROUTINE: EXTERNAL_FORCING_SURF
9       I           iMin, iMax, jMin, jMax,bi,bj,myThid )  C     !INTERFACE:
10  C     /==========================================================\        SUBROUTINE EXTERNAL_FORCING_SURF(
11  C     | SUBROUTINE EXTERNAL_FORCING_SURF_U                       |       I             bi, bj, iMin, iMax, jMin, jMax,
12  C     | o Determines forcing terms based on external fields      |       I             myTime, myIter, myThid )
13  C     |   relaxation terms etc.                                  |  C     !DESCRIPTION: \bv
14  C     \==========================================================/  C     *==========================================================*
15    C     | SUBROUTINE EXTERNAL_FORCING_SURF                          
16    C     | o Determines forcing terms based on external fields      
17    C     |   relaxation terms etc.                                  
18    C     *==========================================================*
19    C     \ev
20    
21    C     !USES:
22        IMPLICIT NONE        IMPLICIT NONE
   
23  C     === Global variables ===  C     === Global variables ===
24  #include "SIZE.h"  #include "SIZE.h"
25  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 18  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:
36  C     === Routine arguments ===  C     === Routine arguments ===
37  C     myThid - Thread no. that called this routine.  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.
42          _RL myTime
43          INTEGER myIter
44        INTEGER myThid        INTEGER myThid
45  C     iMin - Working range of tile for applying forcing.        INTEGER bi,bj
46  C     iMax        INTEGER iMin, iMax
47  C     jMin        INTEGER jMin, jMax
 C     jMax  
       INTEGER iMin, iMax, jMin, jMax, bi, bj  
 CEndOfInterface  
   
 C     === Functions ===  
   
 C     === Local arrays ===  
48    
49    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  c     Zonal wind stress fu:        INTEGER ks
54         DO j=jMin,jMax  #ifdef ALLOW_DIAGNOSTICS
55          DO i=iMin,iMax        LOGICAL  DIAGNOSTICS_IS_ON
56            surfaceTendencyU(i,j,bi,bj) = fu(i,j,bi,bj)        EXTERNAL DIAGNOSTICS_IS_ON
57          _RL tmp1k(1:sNx,1:sNy)
58    #endif /* ALLOW_DIAGNOSTICS */
59    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(i,j,bi,bj) * (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(i,j,bi,bj) * (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              surfaceForcingT(i,j,bi,bj) =
89         &      -lambdaThetaClimRelax(i,j,bi,bj)
90         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
91         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
92    C     Salt Flux (restoring term) :
93              surfaceForcingS(i,j,bi,bj) =
94         &      -lambdaSaltClimRelax(i,j,bi,bj)
95         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
96         &         *drF(ks)*hFacC(i,j,ks,bi,bj)
97    #endif /* ALLOW_SEAICE */
98          ENDDO          ENDDO
99         ENDDO         ENDDO
100    
101        RETURN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
102        END  #ifdef NONLIN_FRSURF
103    C-    T,S surface forcing will be applied (thermodynamics) after the update
104  CStartOfInterface  C     of surf.thickness (hFac): account for change in surf.thickness
105        SUBROUTINE EXTERNAL_FORCING_SURF_V(         IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
106       I           iMin, iMax, jMin, jMax,bi,bj,myThid )          IF (select_rStar.GT.0) THEN
107  C     /==========================================================\           DO j=jMin,jMax
108  C     | SUBROUTINE EXTERNAL_FORCING_SURF_V                       |            DO i=iMin,iMax
109  C     | o Determines forcing terms based on external fields      |              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
110  C     |   relaxation terms etc.                                  |       &                                  * rStarExpC(i,j,bi,bj)
111  C     \==========================================================/              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
112        IMPLICIT NONE       &                                  * rStarExpC(i,j,bi,bj)
113              ENDDO
114  C     === Global variables ===           ENDDO
115  #include "SIZE.h"          ELSE
116  #include "EEPARAMS.h"           DO j=jMin,jMax
117  #include "PARAMS.h"            DO i=iMin,iMax
118  #include "FFIELDS.h"             IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
119  #include "DYNVARS.h"              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
120  #include "GRID.h"       &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
121                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
122  C     === Routine arguments ===       &             *recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
123  C     myThid - Thread no. that called this routine.             ENDIF
124        INTEGER myThid            ENDDO
125  C     iMin - Working range of tile for applying forcing.           ENDDO
126  C     iMax          ENDIF
127  C     jMin         ENDIF
128  C     jMax  #endif /* NONLIN_FRSURF */
129        INTEGER iMin, iMax, jMin, jMax, bi, bj  
130  CEndOfInterface        ELSE
131    C--   No restoring for T & S : set surfaceForcingT,S to zero :
132  C     === Functions ===  
133           DO j = jMin, jMax
134  C     === Local arrays ===          DO i = iMin, iMax
135              surfaceForcingT(i,j,bi,bj) = 0. _d 0
136  C     === Local variables ===            surfaceForcingS(i,j,bi,bj) = 0. _d 0
       INTEGER i,j  
   
 c     Zonal wind stress fv:  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
           surfaceTendencyV(i,j,bi,bj) = fV(i,j,bi,bj)  
137          ENDDO          ENDDO
138         ENDDO         ENDDO
139    
140        RETURN  C--   end restoring / no restoring block.
141        END        ENDIF
142    
143  CStartOfInterface  #ifdef ALLOW_DIAGNOSTICS
144        SUBROUTINE EXTERNAL_FORCING_SURF_T(  C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
145       I           iMin, iMax, jMin, jMax,bi,bj,myThid )        IF ( DIAGNOSTICS_IS_ON('TRELAX  ',myThid) ) THEN
146  C     /==========================================================\           DO j = 1,sNy
147  C     | SUBROUTINE EXTERNAL_FORCING_SURF_T                       |              DO i = 1,sNx
148  C     | o Determines forcing terms based on external fields      |                 tmp1k(i,j) = surfaceForcingT(i,j,bi,bj)
149  C     |   relaxation terms etc.                                  |       &              *HeatCapacity_Cp*recip_horiVertRatio*rhoConst
150  C     \==========================================================/              ENDDO
151        IMPLICIT NONE           ENDDO
152             CALL DIAGNOSTICS_FILL(tmp1k,'TRELAX  ',0,1,3,bi,bj,myThid)
153  C     === Global variables ===        ENDIF
154  #include "SIZE.h"  
155  #include "EEPARAMS.h"  C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
156  #include "PARAMS.h"        IF ( DIAGNOSTICS_IS_ON('SRELAX  ',myThid) ) THEN
157  #include "FFIELDS.h"           DO j = 1,sNy
158  #include "DYNVARS.h"              DO i = 1,sNx
159                   tmp1k(i,j) = surfaceForcingS(i,j,bi,bj)*
160  C     === Routine arguments ===       &              recip_horiVertRatio*rhoConst
161  C     myThid - Thread no. that called this routine.              ENDDO
162        INTEGER myThid           ENDDO
163  C     iMin - Working range of tile for applying forcing.           CALL DIAGNOSTICS_FILL(tmp1k,'SRELAX  ',0,1,3,bi,bj,myThid)
164  C     iMax        ENDIF
165  C     jMin  #endif /* ALLOW_DIAGNOSTICS */
166  C     jMax  
167        INTEGER iMin, iMax, jMin, jMax, bi, bj  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168  CEndOfInterface  
169    C--   Surface Fluxes :
170  C     === Functions ===  
171          DO j = jMin, jMax
172  C     === Local arrays ===           DO i = iMin, iMax
173    
174    C     Zonal wind stress fu:
175              surfaceForcingU(i,j,bi,bj) =
176         &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
177    C     Meridional wind stress fv:
178              surfaceForcingV(i,j,bi,bj) =
179         &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
180    C     Net heat flux Qnet:
181              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
182         &       - ( Qnet(i,j,bi,bj)
183    #ifdef SHORTWAVE_HEATING
184         &          -Qsw(i,j,bi,bj)
185    #endif
186         &         ) *recip_Cp*horiVertRatio*recip_rhoConst
187    C     Net Salt Flux :
188              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
189         &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst
190    
191             ENDDO
192          ENDDO
193            
194    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
195    C--   Fresh-water flux
196    
197    #ifdef EXACT_CONSERV
198    c NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
199    c     to stay consitent with volume change (=d/dt etaH).
200          IF ( staggerTimeStep ) THEN
201            DO j=1,sNy
202             DO i=1,sNx
203               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
204             ENDDO
205            ENDDO
206          ENDIF
207    
208  C     === Local variables ===        IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
209        INTEGER i,j       &     .AND. useRealFreshWaterFlux ) THEN
210    
211  c     Net heat flux Qnet:  c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
212         DO j=jMin,jMax  c   the water column height ; temp., salt, (tracer) flux associated
213          DO i=iMin,iMax  c   with this input/output of water is added here to the surface tendency.
214            surfaceTendencyT(i,j,bi,bj) = Qnet(i,j,bi,bj)  c
215       &         - lambdaThetaClimRelax*  
216       &           (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))         IF (temp_EvPrRn.NE.UNSET_RL) THEN
217            DO j = jMin, jMax
218             DO i = iMin, iMax
219              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
220         &      + PmEpR(i,j,bi,bj)
221         &         *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
222         &         *convertEmP2rUnit
223             ENDDO
224          ENDDO          ENDDO
225         ENDDO         ENDIF
226    
227        RETURN         IF (salt_EvPrRn.NE.UNSET_RL) THEN
228        END          DO j = jMin, jMax
229             DO i = iMin, iMax
230              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
231         &      + PmEpR(i,j,bi,bj)
232         &         *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
233         &         *convertEmP2rUnit
234             ENDDO
235            ENDDO
236           ENDIF
237    
238  CStartOfInterface  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
239        SUBROUTINE EXTERNAL_FORCING_SURF_S(        ELSE
240       I           iMin, iMax, jMin, jMax,bi,bj,myThid )  #else /* EXACT_CONSERV */
241  C     /==========================================================\        IF (.TRUE.) THEN
242  C     | SUBROUTINE EXTERNAL_FORCING_SURF_S                       |  #endif /* EXACT_CONSERV */
243  C     | o Determines forcing terms based on external fields      |  
244  C     |   relaxation terms etc.                                  |  c- EmPmR does not really affect the water column height (for tracer budget)
245  C     \==========================================================/  c   and is converted to a salt tendency.
246        IMPLICIT NONE  
247           IF (convertFW2Salt .EQ. -1.) THEN
248  C     === Global variables ===  c- converts EmPmR to salinity tendency using surface local salinity
249  #include "SIZE.h"          DO j = jMin, jMax
250  #include "EEPARAMS.h"           DO i = iMin, iMax
251  #include "PARAMS.h"            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
252  #include "FFIELDS.h"       &      + EmPmR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
253  #include "DYNVARS.h"       &         *convertEmP2rUnit
254             ENDDO
255  C     === Routine arguments ===          ENDDO
256  C     myThid - Thread no. that called this routine.         ELSE
257        INTEGER myThid  c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
258  C     iMin - Working range of tile for applying forcing.          DO j = jMin, jMax
259  C     iMax           DO i = iMin, iMax
260  C     jMin            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
261  C     jMax       &      + EmPmR(i,j,bi,bj)*convertFW2Salt
262        INTEGER iMin, iMax, jMin, jMax, bi, bj       &         *convertEmP2rUnit
263  CEndOfInterface           ENDDO
264            ENDDO
265  C     === Functions ===         ENDIF
266    
267  C     === Local arrays ===        ENDIF
268    
269  C     === Local variables ===  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       INTEGER i,j  
270    
271  c     Net heat flux Qnet:  #ifdef ALLOW_PTRACERS
272         DO j=jMin,jMax        IF ( usePTRACERS ) THEN
273          DO i=iMin,iMax           CALL PTRACERS_FORCING_SURF(
274  #ifdef USE_NATURAL_BCS       I        bi, bj, iMin, iMax, jMin, jMax,
275  c     Freshwater flux EmPmR:       I        myTime,myIter,myThid )
276            surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)        ENDIF
277       &         - lambdaSaltClimRelax*  #endif /* ALLOW_PTRACERS */
278       &           (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))  
279  #else  #ifdef ATMOSPHERIC_LOADING
280  c     Freshwater flux EmPmR:  
281            surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)  C-- Atmospheric surface Pressure loading :
282       &         - lambdaSaltClimRelax*  
283       &           (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))        IF ( usingZCoords ) THEN
284  #endif                   IF ( useRealFreshWaterFlux ) THEN
285            DO j = jMin, jMax
286             DO i = iMin, iMax
287              phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
288         &                        +sIceLoad(i,j,bi,bj)*gravity
289         &                          )*recip_rhoConst
290             ENDDO
291          ENDDO          ENDDO
292         ENDDO         ELSE
293            DO j = jMin, jMax
294             DO i = iMin, iMax
295              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
296             ENDDO
297            ENDDO
298           ENDIF
299          ELSEIF ( usingPCoords ) THEN
300    C-- This is a hack used to read phi0surf from a file (ploadFile)
301    C   instead of computing it from bathymetry & density ref. profile.
302    C   The true atmospheric P-loading is not yet implemented for P-coord
303    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
304            DO j = jMin, jMax
305             DO i = iMin, iMax
306              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
307             ENDDO
308            ENDDO
309          ENDIF
310    
311    #endif /* ATMOSPHERIC_LOADING */
312    
313    #ifdef ALLOW_EBM
314    c--    Values for surfaceForcingT, surfaceForcingS
315    c      are overwritten by those produced by EBM
316    cph    AD recomputation problems if these IF useEBM are used
317    cph      IF ( useEBM ) THEN
318           CALL EBM_FORCING_SURF(
319         I        bi, bj, iMin, iMax, jMin, jMax,
320         I        myTime,myIter,myThid )
321    cph      ENDIF
322    #endif
323    
324        RETURN        RETURN
325        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22