/[MITgcm]/MITgcm/verification/tidal_basin_2d/code/external_forcing.F
ViewVC logotype

Annotation of /MITgcm/verification/tidal_basin_2d/code/external_forcing.F

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


Revision 1.7 - (hide annotations) (download)
Tue May 6 15:52:51 2014 UTC (10 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x
Changes since 1.6: +3 -2 lines
define "recip_Cp" as local variable (no longer in common block)

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

  ViewVC Help
Powered by ViewVC 1.1.22