/[MITgcm]/MITgcm/model/src/external_forcing.F
ViewVC logotype

Diff of /MITgcm/model/src/external_forcing.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.13 by cnh, Wed Sep 26 18:09:14 2001 UTC revision 1.34 by jmc, Thu Dec 15 17:47:54 2005 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    #ifdef ALLOW_EXF
7    # include "EXF_OPTIONS.h"
8    #endif
9    
10  CBOP  CBOP
11  C     !ROUTINE: EXTERNAL_FORCING_U  C     !ROUTINE: EXTERNAL_FORCING_U
12  C     !INTERFACE:  C     !INTERFACE:
13        SUBROUTINE EXTERNAL_FORCING_U(        SUBROUTINE EXTERNAL_FORCING_U(
14       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
15       I           myCurrentTime,myThid)       I           myTime, myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | S/R EXTERNAL_FORCING_U                                      C     | S/R EXTERNAL_FORCING_U
19  C     | o Contains problem specific forcing for zonal velocity.    C     | o Contains problem specific forcing for zonal velocity.
20  C     *==========================================================*  C     *==========================================================*
21  C     | Adds terms to gU for forcing by external sources            C     | Adds terms to gU for forcing by external sources
22  C     | e.g. wind stress, bottom friction etc..................    C     | e.g. wind stress, bottom friction etc ...
23  C     *==========================================================*  C     *==========================================================*
24  C     \ev  C     \ev
25    
# Line 31  C     == Global data == Line 35  C     == Global data ==
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
39  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
40  C     jMin  C     bi,bj     :: Current tile indices
41  C     jMax  C     kLev      :: Current vertical level index
42  C     kLev  C     myTime    :: Current time in simulation
43    C     myThid    :: Thread Id number
44        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
45        _RL myCurrentTime        _RL myTime
46        INTEGER myThid        INTEGER myThid
47    
48  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
49  C     == Local variables ==  C     == Local variables ==
50  C     Loop counters  C     i,j       :: Loop counters
51        INTEGER I, J  C     kSurface  :: index of surface layer
52          INTEGER i, j
53          INTEGER kSurface
54  CEOP  CEOP
55    
56          IF ( fluidIsAir ) THEN
57           kSurface = 0
58          ELSEIF ( usingPCoords ) THEN
59           kSurface = Nr
60          ELSE
61           kSurface = 1
62          ENDIF
63    
64  C--   Forcing term  C--   Forcing term
65    #ifdef ALLOW_AIM
66          IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
67         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
68         &                      myTime, myThid )
69    #endif /* ALLOW_AIM */
70    
71    #ifdef ALLOW_FIZHI
72          IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
73         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
74         &                      myTime, myThid )
75    #endif /* ALLOW_FIZHI */
76    
77  C     Add windstress momentum impulse into the top-layer  C     Add windstress momentum impulse into the top-layer
78        IF ( kLev .EQ. 1 ) THEN        IF ( kLev .EQ. kSurface ) THEN
79         DO j=jMin,jMax  c      DO j=1,sNy
80          DO i=iMin,iMax  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
81           DO j=0,sNy+1
82            DO i=1,sNx+1
83           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
84       &   +foFacMom*surfaceTendencyU(i,j,bi,bj)       &   +foFacMom*surfaceForcingU(i,j,bi,bj)
85       &   *_maskW(i,j,kLev,bi,bj)       &   *recip_drF(kLev)*recip_hFacW(i,j,kLev,bi,bj)
86          ENDDO          ENDDO
87         ENDDO         ENDDO
88        ENDIF        ENDIF
89    
90    #if (defined (ALLOW_TAU_EDDY))
91           CALL TAUEDDY_EXTERNAL_FORCING_U(
92         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
93         I           myTime, myThid )
94    #endif
95    
96    #ifdef ALLOW_OBCS
97          IF (useOBCS) THEN
98           CALL OBCS_SPONGE_U(
99         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
100         I           myTime, myThid )
101          ENDIF
102    #endif
103    
104        RETURN        RETURN
105        END        END
106    
107    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108  CBOP  CBOP
109  C     !ROUTINE: EXTERNAL_FORCING_V  C     !ROUTINE: EXTERNAL_FORCING_V
110  C     !INTERFACE:  C     !INTERFACE:
111        SUBROUTINE EXTERNAL_FORCING_V(        SUBROUTINE EXTERNAL_FORCING_V(
112       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
113       I           myCurrentTime,myThid)       I           myTime, myThid )
114  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
115  C     *==========================================================*  C     *==========================================================*
116  C     | S/R EXTERNAL_FORCING_V                                      C     | S/R EXTERNAL_FORCING_V
117  C     | o Contains problem specific forcing for merid velocity.    C     | o Contains problem specific forcing for merid velocity.
118  C     *==========================================================*  C     *==========================================================*
119  C     | Adds terms to gV for forcing by external sources            C     | Adds terms to gV for forcing by external sources
120  C     | e.g. wind stress, bottom friction etc..................    C     | e.g. wind stress, bottom friction etc ...
121  C     *==========================================================*  C     *==========================================================*
122  C     \ev  C     \ev
123    
# Line 88  C     == Global data == Line 133  C     == Global data ==
133    
134  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
135  C     == Routine arguments ==  C     == Routine arguments ==
136  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
137  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
138  C     jMin  C     bi,bj     :: Current tile indices
139  C     jMax  C     kLev      :: Current vertical level index
140  C     kLev  C     myTime    :: Current time in simulation
141    C     myThid    :: Thread Id number
142        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
143        _RL myCurrentTime        _RL myTime
144        INTEGER myThid        INTEGER myThid
145    
146  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
147  C     == Local variables ==  C     == Local variables ==
148  C     Loop counters  C     i,j       :: Loop counters
149        INTEGER I, J  C     kSurface  :: index of surface layer
150          INTEGER i, j
151          INTEGER kSurface
152  CEOP  CEOP
153    
154          IF ( fluidIsAir ) THEN
155           kSurface = 0
156          ELSEIF ( usingPCoords ) THEN
157           kSurface = Nr
158          ELSE
159           kSurface = 1
160          ENDIF
161    
162  C--   Forcing term  C--   Forcing term
163    #ifdef ALLOW_AIM
164          IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
165         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
166         &                      myTime, myThid )
167    #endif /* ALLOW_AIM */
168    
169    #ifdef ALLOW_FIZHI
170          IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
171         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
172         &                      myTime, myThid )
173    #endif /* ALLOW_FIZHI */
174    
175  C     Add windstress momentum impulse into the top-layer  C     Add windstress momentum impulse into the top-layer
176        IF ( kLev .EQ. 1 ) THEN        IF ( kLev .EQ. kSurface ) THEN
177         DO j=jMin,jMax         DO j=1,sNy+1
178          DO i=iMin,iMax  c       DO i=1,sNx
179    C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
180            DO i=0,sNx+1
181           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
182       &   +foFacMom*surfaceTendencyV(i,j,bi,bj)       &   +foFacMom*surfaceForcingV(i,j,bi,bj)
183       &   *_maskS(i,j,kLev,bi,bj)       &   *recip_drF(kLev)*recip_hFacS(i,j,kLev,bi,bj)
184          ENDDO          ENDDO
185         ENDDO         ENDDO
186        ENDIF        ENDIF
187    
188    #if (defined (ALLOW_TAU_EDDY))
189           CALL TAUEDDY_EXTERNAL_FORCING_V(
190         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
191         I           myTime, myThid )
192    #endif
193    
194    #ifdef ALLOW_OBCS
195          IF (useOBCS) THEN
196           CALL OBCS_SPONGE_V(
197         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
198         I           myTime, myThid )
199          ENDIF
200    #endif
201    
202        RETURN        RETURN
203        END        END
204    
205    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
206  CBOP  CBOP
207  C     !ROUTINE: EXTERNAL_FORCING_T  C     !ROUTINE: EXTERNAL_FORCING_T
208  C     !INTERFACE:  C     !INTERFACE:
209        SUBROUTINE EXTERNAL_FORCING_T(        SUBROUTINE EXTERNAL_FORCING_T(
210       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
211       I           myCurrentTime,myThid)       I           myTime, myThid )
212  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
213  C     *==========================================================*  C     *==========================================================*
214  C     | S/R EXTERNAL_FORCING_T                                      C     | S/R EXTERNAL_FORCING_T
215  C     | o Contains problem specific forcing for temperature.        C     | o Contains problem specific forcing for temperature.
216  C     *==========================================================*  C     *==========================================================*
217  C     | Adds terms to gT for forcing by external sources            C     | Adds terms to gT for forcing by external sources
218  C     | e.g. heat flux, climatalogical relaxation..............    C     | e.g. heat flux, climatalogical relaxation, etc ...
219  C     *==========================================================*  C     *==========================================================*
220  C     \ev  C     \ev
221    
# Line 142  C     == Global data == Line 228  C     == Global data ==
228  #include "GRID.h"  #include "GRID.h"
229  #include "DYNVARS.h"  #include "DYNVARS.h"
230  #include "FFIELDS.h"  #include "FFIELDS.h"
 #ifdef SHORTWAVE_HEATING  
       integer two  
       _RL minusone  
       parameter (two=2,minusone=-1.)  
       _RL swfracb(two)  
 #endif  
231    
232  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
233  C     == Routine arguments ==  C     == Routine arguments ==
234  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
235  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
236  C     jMin  C     bi,bj     :: Current tile indices
237  C     jMax  C     kLev      :: Current vertical level index
238  C     kLev  C     myTime    :: Current time in simulation
239    C     myThid    :: Thread Id number
240        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
241        _RL myCurrentTime        _RL myTime
242        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
243    
244  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
245  C     == Local variables ==  C     == Local variables ==
246  C     Loop counters  C     i,j       :: Loop counters
247        INTEGER I, J  C     kSurface  :: index of surface layer
248          INTEGER i, j
249          INTEGER kSurface
250  CEOP  CEOP
251    #ifdef SHORTWAVE_HEATING
252          integer two
253          _RL minusone
254          parameter (two=2,minusone=-1.)
255          _RL swfracb(two)
256          INTEGER kp1
257    #endif
258    
259          IF ( fluidIsAir ) THEN
260           kSurface = 0
261          ELSEIF ( usingPCoords ) THEN
262           kSurface = Nr
263          ELSE
264           kSurface = 1
265          ENDIF
266    
267  C--   Forcing term  C--   Forcing term
268    #ifdef ALLOW_AIM
269          IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
270         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
271         &                      myTime, myThid )
272    #endif /* ALLOW_AIM */
273    
274    #ifdef ALLOW_FIZHI
275          IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
276         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
277         &                      myTime, myThid )
278    #endif /* ALLOW_FIZHI */
279    
280  C     Add heat in top-layer  C     Add heat in top-layer
281        IF ( kLev .EQ. 1 ) THEN        IF ( kLev .EQ. kSurface ) THEN
282         DO j=jMin,jMax         DO j=1,sNy
283          DO i=iMin,iMax          DO i=1,sNx
284           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
285       &     +maskC(i,j,kLev,bi,bj)*surfaceTendencyT(i,j,bi,bj)       &     +surfaceForcingT(i,j,bi,bj)
286         &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)
287          ENDDO          ENDDO
288         ENDDO         ENDDO
289        ENDIF        ENDIF
290    
291  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
292  C Penetrating SW radiation  C Penetrating SW radiation
293        swfracb(1)=abs(rF(klev))  c     IF ( usePenetratingSW ) THEN
294        swfracb(2)=abs(rF(klev+1))         swfracb(1)=abs(rF(klev))
295        call SWFRAC(         swfracb(2)=abs(rF(klev+1))
296           CALL SWFRAC(
297       I     two,minusone,       I     two,minusone,
298       I     myCurrentTime,myThid,       I     myTime,myThid,
299       O     swfracb)       U     swfracb)
300        DO j=jMin,jMax         kp1 = klev+1
301         DO i=iMin,iMax         IF (klev.EQ.Nr) THEN
302          gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)          kp1 = klev
303       &   -maskC(i,j,klev,bi,bj)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))          swfracb(2)=0. _d 0
304       &    *recip_Cp*recip_rhoNil*recip_drF(klev)         ENDIF
305           DO j=1,sNy
306            DO i=1,sNx
307             gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
308         &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
309         &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))
310         &    *recip_Cp*recip_rhoConst
311         &    *recip_drF(klev)*recip_hFacC(i,j,kLev,bi,bj)
312            ENDDO
313         ENDDO         ENDDO
314        ENDDO  c     ENDIF
315    #endif
316    
317    #ifdef ALLOW_CLIMTEMP_RELAXATION
318           IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN
319            DO j=1,sNy
320             DO i=1,sNx
321              gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
322         &     -1./tauThetaClimRelax3Dim
323         &         *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj))
324         &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)
325             ENDDO
326            ENDDO
327           ENDIF
328    #endif
329    
330    #ifdef ALLOW_OBCS
331          IF (useOBCS) THEN
332           CALL OBCS_SPONGE_T(
333         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
334         I           myTime, myThid )
335          ENDIF
336  #endif  #endif
337    
338        RETURN        RETURN
339        END        END
340    
341    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
342  CBOP  CBOP
343  C     !ROUTINE: EXTERNAL_FORCING_S  C     !ROUTINE: EXTERNAL_FORCING_S
344  C     !INTERFACE:  C     !INTERFACE:
345        SUBROUTINE EXTERNAL_FORCING_S(        SUBROUTINE EXTERNAL_FORCING_S(
346       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
347       I           myCurrentTime,myThid)       I           myTime, myThid )
348    
349  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
350  C     *==========================================================*  C     *==========================================================*
351  C     | S/R EXTERNAL_FORCING_S                                      C     | S/R EXTERNAL_FORCING_S
352  C     | o Contains problem specific forcing for merid velocity.    C     | o Contains problem specific forcing for merid velocity.
353  C     *==========================================================*  C     *==========================================================*
354  C     | Adds terms to gS for forcing by external sources            C     | Adds terms to gS for forcing by external sources
355  C     | e.g. fresh-water flux, climatalogical relaxation.......    C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
356  C     *==========================================================*  C     *==========================================================*
357  C     \ev  C     \ev
358    
# Line 225  C     == Global data == Line 368  C     == Global data ==
368    
369  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
370  C     == Routine arguments ==  C     == Routine arguments ==
371  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
372  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
373  C     jMin  C     bi,bj     :: Current tile indices
374  C     jMax  C     kLev      :: Current vertical level index
375  C     kLev  C     myTime    :: Current time in simulation
376    C     myThid    :: Thread Id number
377        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
378        _RL myCurrentTime        _RL myTime
379        INTEGER myThid        INTEGER myThid
380    
381  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
382  C     == Local variables ==  C     == Local variables ==
383  C     Loop counters  C     i,j       :: Loop counters
384        INTEGER I, J  C     kSurface  :: index of surface layer
385          INTEGER i, j
386          INTEGER kSurface
387  CEOP  CEOP
388    
389          IF ( fluidIsAir ) THEN
390           kSurface = 0
391          ELSEIF ( usingPCoords ) THEN
392           kSurface = Nr
393          ELSE
394           kSurface = 1
395          ENDIF
396    
397  C--   Forcing term  C--   Forcing term
398    #ifdef ALLOW_AIM
399          IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
400         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
401         &                      myTime, myThid )
402    #endif /* ALLOW_AIM */
403    
404    #ifdef ALLOW_FIZHI
405          IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
406         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
407         &                      myTime, myThid )
408    #endif /* ALLOW_FIZHI */
409    
410  C     Add fresh-water in top-layer  C     Add fresh-water in top-layer
411        IF ( kLev .EQ. 1 ) THEN        IF ( kLev .EQ. kSurface ) THEN
412         DO j=jMin,jMax         DO j=1,sNy
413          DO i=iMin,iMax          DO i=1,sNx
414           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
415       &   +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj)       &     +surfaceForcingS(i,j,bi,bj)
416         &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)
417          ENDDO          ENDDO
418         ENDDO         ENDDO
419        ENDIF        ENDIF
420    
421    #ifdef ALLOW_CLIMSALT_RELAXATION
422           IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN
423            DO j=1,sNy
424             DO i=1,sNx
425              gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj)
426         &     -1./tauSaltClimRelax3Dim
427         &         *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj))
428         &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)
429             ENDDO
430            ENDDO
431           ENDIF
432    #endif
433    
434    #ifdef ALLOW_OBCS
435          IF (useOBCS) THEN
436           CALL OBCS_SPONGE_S(
437         I           iMin,iMax, jMin,jMax, bi,bj, kLev,
438         I           myTime, myThid )
439          ENDIF
440    #endif
441    
442        RETURN        RETURN
443        END        END
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22