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

Contents 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.8 - (show annotations) (download)
Sat Jul 5 15:11:36 2014 UTC (9 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64z, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.7: +1 -89 lines
remove PKG S/R calls not used in this experiment: this is closer
 to the original version (and easier to maintained, specially
 since this set-up is not tested).

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

  ViewVC Help
Powered by ViewVC 1.1.22