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