/[MITgcm]/MITgcm_contrib/natl_12/code/external_forcing.F
ViewVC logotype

Contents of /MITgcm_contrib/natl_12/code/external_forcing.F

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


Revision 1.4 - (show annotations) (download)
Thu Aug 7 13:03:31 2003 UTC (21 years, 11 months ago) by cnh
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +6 -6 lines
Mods to allow laplacian as well as biharmonic viscosity with varying resolution

1 C $Header: /u/u0/gcmpack/MITgcm/model/src/external_forcing.F,v 1.19 2003/06/19 15:00:45 heimbach Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #ifdef ALLOW_OBCS
6 # include "OBCS_OPTIONS.h"
7 #endif
8
9 CBOP
10 C !ROUTINE: EXTERNAL_FORCING_U
11 C !INTERFACE:
12 SUBROUTINE EXTERNAL_FORCING_U(
13 I iMin, iMax, jMin, jMax,bi,bj,kLev,
14 I myCurrentTime,myThid)
15 C !DESCRIPTION: \bv
16 C *==========================================================*
17 C | S/R EXTERNAL_FORCING_U
18 C | o Contains problem specific forcing for zonal velocity.
19 C *==========================================================*
20 C | Adds terms to gU for forcing by external sources
21 C | e.g. wind stress, bottom friction etc..................
22 C *==========================================================*
23 C \ev
24
25 C !USES:
26 IMPLICIT NONE
27 C == Global data ==
28 #include "SIZE.h"
29 #include "EEPARAMS.h"
30 #include "PARAMS.h"
31 #include "GRID.h"
32 #include "DYNVARS.h"
33 #include "FFIELDS.h"
34
35 C !INPUT/OUTPUT PARAMETERS:
36 C == Routine arguments ==
37 C iMin - Working range of tile for applying forcing.
38 C iMax
39 C jMin
40 C jMax
41 C kLev
42 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
43 _RL myCurrentTime
44 INTEGER myThid
45
46 C !LOCAL VARIABLES:
47 C == Local variables ==
48 C Loop counters
49 INTEGER I, J
50 C number of surface interface layer
51 INTEGER kSurface
52 C Cheap sponge layer
53 _RS recip_tauSp(5)
54 INTEGER spWidth
55 _RS curRecipTau
56 INTEGER jFromNBndy, jFromSBndy,
57 & jNorthBndy, jSouthBndy, jG
58 CEOP
59
60 if ( buoyancyRelation .eq. 'OCEANICP' ) then
61 kSurface = Nr
62 else
63 kSurface = 1
64 endif
65
66 C-- Forcing term
67 C Add windstress momentum impulse into the top-layer
68 IF ( kLev .EQ. kSurface ) THEN
69 DO j=jMin,jMax
70 DO i=iMin,iMax
71 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
72 & +foFacMom*surfaceTendencyU(i,j,bi,bj)
73 & *_maskW(i,j,kLev,bi,bj)
74 ENDDO
75 ENDDO
76 ENDIF
77
78 C-- Create a sponge layer where flow is linearly damped over entire water column
79 C Damping time scale decreases away from boundary so that
80 C tau = 1 day, 5days, 10days, 20days, 60days
81 spWidth = 5
82 recip_tauSp(1) = 1./(86400.*1. )
83 recip_tauSp(2) = 1./(86400.*5. )
84 recip_tauSp(3) = 1./(86400.*10.)
85 recip_tauSp(4) = 1./(86400.*20.)
86 recip_tauSp(5) = 1./(86400.*60.)
87 jSouthBndy = 5
88 jNorthBndy = ny-5+1
89 DO j=1,sNy
90 DO i=iMin,iMax
91 jG = myYGlobalLo+(bj-1)*sNy+j-1
92 jFromNBndy = jNorthBndy-jG
93 jFromSBndy = jSouthBndy-jG
94 curRecipTau=0.
95 IF ( jFromNBndy .LE. 0 ) THEN
96 curRecipTau = recip_tauSp(jFromNBndy+5)
97 ENDIF
98 IF ( jFromSBndy .GE. 0 ) THEN
99 curRecipTau = recip_tauSp(-(jFromSBndy-5))
100 ENDIF
101 gu(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
102 & -curRecipTau*uVel(i,j,Klev,bi,bj)
103 ENDDO
104 ENDDO
105
106 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
107 IF (useOBCS) THEN
108 CALL OBCS_SPONGE_U(
109 I iMin, iMax, jMin, jMax,bi,bj,kLev,
110 I myCurrentTime,myThid)
111 ENDIF
112 #endif
113
114 RETURN
115 END
116 CBOP
117 C !ROUTINE: EXTERNAL_FORCING_V
118 C !INTERFACE:
119 SUBROUTINE EXTERNAL_FORCING_V(
120 I iMin, iMax, jMin, jMax,bi,bj,kLev,
121 I myCurrentTime,myThid)
122 C !DESCRIPTION: \bv
123 C *==========================================================*
124 C | S/R EXTERNAL_FORCING_V
125 C | o Contains problem specific forcing for merid velocity.
126 C *==========================================================*
127 C | Adds terms to gV for forcing by external sources
128 C | e.g. wind stress, bottom friction etc..................
129 C *==========================================================*
130 C \ev
131
132 C !USES:
133 IMPLICIT NONE
134 C == Global data ==
135 #include "SIZE.h"
136 #include "EEPARAMS.h"
137 #include "PARAMS.h"
138 #include "GRID.h"
139 #include "DYNVARS.h"
140 #include "FFIELDS.h"
141
142 C !INPUT/OUTPUT PARAMETERS:
143 C == Routine arguments ==
144 C iMin - Working range of tile for applying forcing.
145 C iMax
146 C jMin
147 C jMax
148 C kLev
149 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
150 _RL myCurrentTime
151 INTEGER myThid
152
153 C !LOCAL VARIABLES:
154 C == Local variables ==
155 C Loop counters
156 INTEGER I, J
157 C number of surface interface layer
158 INTEGER kSurface
159
160 C == Cheap sponge layer ==
161 _RS recip_tauSp(5)
162 INTEGER spWidth
163 _RS curRecipTau
164 INTEGER jFromNBndy, jFromSBndy,
165 & jNorthBndy, jSouthBndy, jG
166
167
168 CEOP
169
170 if ( buoyancyRelation .eq. 'OCEANICP' ) then
171 kSurface = Nr
172 else
173 kSurface = 1
174 endif
175
176 C-- Forcing term
177 C Add windstress momentum impulse into the top-layer
178 IF ( kLev .EQ. kSurface ) THEN
179 DO j=jMin,jMax
180 DO i=iMin,iMax
181 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
182 & +foFacMom*surfaceTendencyV(i,j,bi,bj)
183 & *_maskS(i,j,kLev,bi,bj)
184 ENDDO
185 ENDDO
186 ENDIF
187
188 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
189 IF (useOBCS) THEN
190 CALL OBCS_SPONGE_V(
191 I iMin, iMax, jMin, jMax,bi,bj,kLev,
192 I myCurrentTime,myThid)
193 ENDIF
194 #endif
195
196 C-- Create a sponge layer where flow is linearly damped over entire water column
197 C Damping time scale decreases away from boundary so that
198 C tau = 1 day, 5days, 10days, 20days, 60days
199 spWidth = 5
200 recip_tauSp(1) = 1./(86400.*1. )
201 recip_tauSp(2) = 1./(86400.*5. )
202 recip_tauSp(3) = 1./(86400.*10.)
203 recip_tauSp(4) = 1./(86400.*20.)
204 recip_tauSp(5) = 1./(86400.*60.)
205 jSouthBndy = 5
206 jNorthBndy = ny-5+1
207 DO j=1,sNy
208 DO i=iMin,iMax
209 jG = myYGlobalLo+(bj-1)*sNy+j-1
210 jFromNBndy = jNorthBndy-jG
211 jFromSBndy = jSouthBndy-jG
212 curRecipTau=0.
213 IF ( jFromNBndy .LE. 0 ) THEN
214 curRecipTau = recip_tauSp(jFromNBndy+5)
215 ENDIF
216 IF ( jFromSBndy .GE. 0 ) THEN
217 curRecipTau = recip_tauSp(-(jFromSBndy-5))
218 ENDIF
219 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
220 & -curRecipTau*vVel(i,j,Klev,bi,bj)
221 ENDDO
222 ENDDO
223
224 RETURN
225 END
226 CBOP
227 C !ROUTINE: EXTERNAL_FORCING_T
228 C !INTERFACE:
229 SUBROUTINE EXTERNAL_FORCING_T(
230 I iMin, iMax, jMin, jMax,bi,bj,kLev,
231 I myCurrentTime,myThid)
232 C !DESCRIPTION: \bv
233 C *==========================================================*
234 C | S/R EXTERNAL_FORCING_T
235 C | o Contains problem specific forcing for temperature.
236 C *==========================================================*
237 C | Adds terms to gT for forcing by external sources
238 C | e.g. heat flux, climatalogical relaxation..............
239 C *==========================================================*
240 C \ev
241
242 C !USES:
243 IMPLICIT NONE
244 C == Global data ==
245 #include "SIZE.h"
246 #include "EEPARAMS.h"
247 #include "PARAMS.h"
248 #include "GRID.h"
249 #include "DYNVARS.h"
250 #include "FFIELDS.h"
251 #ifdef SHORTWAVE_HEATING
252 integer two
253 _RL minusone
254 parameter (two=2,minusone=-1.)
255 _RL swfracb(two)
256 #endif
257
258 C !INPUT/OUTPUT PARAMETERS:
259 C == Routine arguments ==
260 C iMin - Working range of tile for applying forcing.
261 C iMax
262 C jMin
263 C jMax
264 C kLev
265 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
266 _RL myCurrentTime
267 INTEGER myThid
268 CEndOfInterface
269
270 C !LOCAL VARIABLES:
271 C == Local variables ==
272 C Loop counters
273 INTEGER I, J
274 C number of surface interface layer
275 INTEGER kSurface
276 C Cheap sponge layer
277 _RS recip_tauSp(5)
278 INTEGER spWidth
279 _RS curRecipTau
280 INTEGER jFromNBndy, jFromSBndy,
281 & jNorthBndy, jSouthBndy, jG
282
283 CEOP
284
285 if ( buoyancyRelation .eq. 'OCEANICP' ) then
286 kSurface = Nr
287 else
288 kSurface = 1
289 endif
290
291 C-- Forcing term
292 C Add heat in top-layer
293 IF ( kLev .EQ. kSurface ) THEN
294 DO j=jMin,jMax
295 DO i=iMin,iMax
296 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
297 & +maskC(i,j,kLev,bi,bj)*surfaceTendencyT(i,j,bi,bj)
298 ENDDO
299 ENDDO
300 ENDIF
301
302 C-- Forcing term
303 C Add heat in top-layer ( 90 day climatalogical average relaxation )
304 IF ( kLev .EQ. kSurface ) THEN
305 curRecipTau=1./(86400.*90.)
306 DO j=jMin,jMax
307 DO i=iMin,iMax
308 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
309 & +maskC(i,j,kLev,bi,bj)*(
310 & -curRecipTau*(theta(i,j,Klev,bi,bj)-thetaRef(i,j,kLev,bi,bj))
311 & )
312 ENDDO
313 ENDDO
314 ENDIF
315
316 #ifdef SHORTWAVE_HEATING
317 C Penetrating SW radiation
318 swfracb(1)=abs(rF(klev))
319 swfracb(2)=abs(rF(klev+1))
320 call SWFRAC(
321 I two,minusone,
322 I myCurrentTime,myThid,
323 U swfracb)
324 DO j=jMin,jMax
325 DO i=iMin,iMax
326 gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
327 & -maskC(i,j,klev,bi,bj)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))
328 & *recip_Cp*recip_rhoConst*recip_drF(klev)
329 ENDDO
330 ENDDO
331 #endif
332
333 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
334 IF (useOBCS) THEN
335 CALL OBCS_SPONGE_T(
336 I iMin, iMax, jMin, jMax,bi,bj,kLev,
337 I myCurrentTime,myThid)
338 ENDIF
339 #endif
340
341 C-- Create a sponge layer where flow is linearly damped over entire water column
342 C Damping time scale decreases away from boundary so that
343 C tau = 1 day, 5days, 10days, 20days, 60days
344 spWidth = 5
345 recip_tauSp(1) = 1./(86400.*1. )
346 recip_tauSp(2) = 1./(86400.*5. )
347 recip_tauSp(3) = 1./(86400.*10.)
348 recip_tauSp(4) = 1./(86400.*20.)
349 recip_tauSp(5) = 1./(86400.*60.)
350 jSouthBndy = 5
351 jNorthBndy = ny-5+1
352 DO j=1,sNy
353 DO i=iMin,iMax
354 jG = myYGlobalLo+(bj-1)*sNy+j-1
355 jFromNBndy = jNorthBndy-jG
356 jFromSBndy = jSouthBndy-jG
357 curRecipTau=0.
358 IF ( jFromNBndy .LE. 0 ) THEN
359 curRecipTau = recip_tauSp(jFromNBndy+5)
360 ENDIF
361 IF ( jFromSBndy .GE. 0 ) THEN
362 curRecipTau = recip_tauSp(-(jFromSBndy-5))
363 ENDIF
364 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
365 & -curRecipTau*(theta(i,j,Klev,bi,bj)-thetaRef(i,j,kLev,bi,bj))
366 C & *0.0000D0
367 ENDDO
368 ENDDO
369
370
371 RETURN
372 END
373 CBOP
374 C !ROUTINE: EXTERNAL_FORCING_S
375 C !INTERFACE:
376 SUBROUTINE EXTERNAL_FORCING_S(
377 I iMin, iMax, jMin, jMax,bi,bj,kLev,
378 I myCurrentTime,myThid)
379
380 C !DESCRIPTION: \bv
381 C *==========================================================*
382 C | S/R EXTERNAL_FORCING_S
383 C | o Contains problem specific forcing for merid velocity.
384 C *==========================================================*
385 C | Adds terms to gS for forcing by external sources
386 C | e.g. fresh-water flux, climatalogical relaxation.......
387 C *==========================================================*
388 C \ev
389
390 C !USES:
391 IMPLICIT NONE
392 C == Global data ==
393 #include "SIZE.h"
394 #include "EEPARAMS.h"
395 #include "PARAMS.h"
396 #include "GRID.h"
397 #include "DYNVARS.h"
398 #include "FFIELDS.h"
399
400 C !INPUT/OUTPUT PARAMETERS:
401 C == Routine arguments ==
402 C iMin - Working range of tile for applying forcing.
403 C iMax
404 C jMin
405 C jMax
406 C kLev
407 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
408 _RL myCurrentTime
409 INTEGER myThid
410
411 C !LOCAL VARIABLES:
412 C == Local variables ==
413 C Loop counters
414 INTEGER I, J
415 C number of surface interface layer
416 INTEGER kSurface
417 C Cheap sponge layer
418 _RS recip_tauSp(5)
419 INTEGER spWidth
420 _RS curRecipTau
421 INTEGER jFromNBndy, jFromSBndy,
422 & jNorthBndy, jSouthBndy, jG
423
424 CEOP
425
426 if ( buoyancyRelation .eq. 'OCEANICP' ) then
427 kSurface = Nr
428 else
429 kSurface = 1
430 endif
431
432
433 C-- Forcing term
434 C Add fresh-water in top-layer
435 IF ( kLev .EQ. kSurface ) THEN
436 DO j=jMin,jMax
437 DO i=iMin,iMax
438 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
439 & +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj)
440 ENDDO
441 ENDDO
442 ENDIF
443
444 C-- Forcing term
445 C Add freshening/salt in top-layer ( 90 day climatalogical average relaxation )
446 IF ( kLev .EQ. kSurface ) THEN
447 curRecipTau=1./(86400.*90.)
448 DO j=jMin,jMax
449 DO i=iMin,iMax
450 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
451 & +maskC(i,j,kLev,bi,bj)*(
452 & -curRecipTau*(salt(i,j,Klev,bi,bj)-saltRef(i,j,kLev,bi,bj))
453 & )
454 ENDDO
455 ENDDO
456 ENDIF
457
458 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
459 IF (useOBCS) THEN
460 CALL OBCS_SPONGE_S(
461 I iMin, iMax, jMin, jMax,bi,bj,kLev,
462 I myCurrentTime,myThid)
463 ENDIF
464 #endif
465
466 C-- Create a sponge layer where flow is linearly damped over entire water column
467 C Damping time scale decreases away from boundary so that
468 C tau = 1 day, 5days, 10days, 20days, 60days
469 spWidth = 5
470 recip_tauSp(1) = 1./(86400.*1. )
471 recip_tauSp(2) = 1./(86400.*5. )
472 recip_tauSp(3) = 1./(86400.*10.)
473 recip_tauSp(4) = 1./(86400.*20.)
474 recip_tauSp(5) = 1./(86400.*60.)
475 jSouthBndy = 5
476 jNorthBndy = ny-5+1
477 DO j=1,sNy
478 DO i=iMin,iMax
479 jG = myYGlobalLo+(bj-1)*sNy+j-1
480 jFromNBndy = jNorthBndy-jG
481 jFromSBndy = jSouthBndy-jG
482 curRecipTau=0.
483 IF ( jFromNBndy .LE. 0 ) THEN
484 curRecipTau = recip_tauSp(jFromNBndy+5)
485 ENDIF
486 IF ( jFromSBndy .GE. 0 ) THEN
487 curRecipTau = recip_tauSp(-(jFromSBndy-5))
488 ENDIF
489 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
490 & -curRecipTau*(salt(i,j,Klev,bi,bj)-saltRef(i,j,kLev,bi,bj))
491 C & *0.0000D0
492 ENDDO
493 ENDDO
494
495 RETURN
496 END

  ViewVC Help
Powered by ViewVC 1.1.22