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

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

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


Revision 1.43 - (show annotations) (download)
Sun Jul 22 23:51:16 2007 UTC (16 years, 9 months ago) by dimitri
Branch: MAIN
Changes since 1.42: +35 -3 lines
o starting implementation of Duffy et al. (GRL 1999) salt plume scheme
  (#ifdef ALLOW_SALT_PLUME)

1 C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.42 2007/05/03 21:41:35 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 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
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 #ifdef ALLOW_MYPACKAGE
75 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_U(
76 & iMin,iMax, jMin,jMax, bi,bj, kLev,
77 & myTime, myThid )
78 #endif /* ALLOW_MYPACKAGE */
79
80 C Add windstress momentum impulse into the top-layer
81 IF ( kLev .EQ. kSurface ) THEN
82 c DO j=1,sNy
83 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
84 DO j=0,sNy+1
85 DO i=1,sNx+1
86 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
87 & +foFacMom*surfaceForcingU(i,j,bi,bj)
88 & *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
89 ENDDO
90 ENDDO
91 ENDIF
92
93 #if (defined (ALLOW_TAU_EDDY))
94 CALL TAUEDDY_EXTERNAL_FORCING_U(
95 I iMin,iMax, jMin,jMax, bi,bj, kLev,
96 I myTime, myThid )
97 #endif
98
99 #ifdef ALLOW_OBCS
100 IF (useOBCS) THEN
101 CALL OBCS_SPONGE_U(
102 I iMin,iMax, jMin,jMax, bi,bj, kLev,
103 I myTime, myThid )
104 ENDIF
105 #endif
106
107 RETURN
108 END
109
110 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
111 CBOP
112 C !ROUTINE: EXTERNAL_FORCING_V
113 C !INTERFACE:
114 SUBROUTINE EXTERNAL_FORCING_V(
115 I iMin,iMax, jMin,jMax, bi,bj, kLev,
116 I myTime, myThid )
117 C !DESCRIPTION: \bv
118 C *==========================================================*
119 C | S/R EXTERNAL_FORCING_V
120 C | o Contains problem specific forcing for merid velocity.
121 C *==========================================================*
122 C | Adds terms to gV for forcing by external sources
123 C | e.g. wind stress, bottom friction etc ...
124 C *==========================================================*
125 C \ev
126
127 C !USES:
128 IMPLICIT NONE
129 C == Global data ==
130 #include "SIZE.h"
131 #include "EEPARAMS.h"
132 #include "PARAMS.h"
133 #include "GRID.h"
134 #include "DYNVARS.h"
135 #include "FFIELDS.h"
136
137 C !INPUT/OUTPUT PARAMETERS:
138 C == Routine arguments ==
139 C iMin,iMax :: Working range of x-index for applying forcing.
140 C jMin,jMax :: Working range of y-index for applying forcing.
141 C bi,bj :: Current tile indices
142 C kLev :: Current vertical level index
143 C myTime :: Current time in simulation
144 C myThid :: Thread Id number
145 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
146 _RL myTime
147 INTEGER myThid
148
149 C !LOCAL VARIABLES:
150 C == Local variables ==
151 C i,j :: Loop counters
152 C kSurface :: index of surface layer
153 INTEGER i, j
154 INTEGER kSurface
155 CEOP
156
157 IF ( fluidIsAir ) THEN
158 kSurface = 0
159 ELSEIF ( usingPCoords ) THEN
160 kSurface = Nr
161 ELSE
162 kSurface = 1
163 ENDIF
164
165 C-- Forcing term
166 #ifdef ALLOW_AIM
167 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
168 & iMin,iMax, jMin,jMax, bi,bj, kLev,
169 & myTime, myThid )
170 #endif /* ALLOW_AIM */
171
172 #ifdef ALLOW_FIZHI
173 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
174 & iMin,iMax, jMin,jMax, bi,bj, kLev,
175 & myTime, myThid )
176 #endif /* ALLOW_FIZHI */
177
178 #ifdef ALLOW_MYPACKAGE
179 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_V(
180 & iMin,iMax, jMin,jMax, bi,bj, kLev,
181 & myTime, myThid )
182 #endif /* ALLOW_MYPACKAGE */
183
184 C Add windstress momentum impulse into the top-layer
185 IF ( kLev .EQ. kSurface ) THEN
186 DO j=1,sNy+1
187 c DO i=1,sNx
188 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
189 DO i=0,sNx+1
190 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
191 & +foFacMom*surfaceForcingV(i,j,bi,bj)
192 & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
193 ENDDO
194 ENDDO
195 ENDIF
196
197 #if (defined (ALLOW_TAU_EDDY))
198 CALL TAUEDDY_EXTERNAL_FORCING_V(
199 I iMin,iMax, jMin,jMax, bi,bj, kLev,
200 I myTime, myThid )
201 #endif
202
203 #ifdef ALLOW_OBCS
204 IF (useOBCS) THEN
205 CALL OBCS_SPONGE_V(
206 I iMin,iMax, jMin,jMax, bi,bj, kLev,
207 I myTime, myThid )
208 ENDIF
209 #endif
210
211 RETURN
212 END
213
214 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
215 CBOP
216 C !ROUTINE: EXTERNAL_FORCING_T
217 C !INTERFACE:
218 SUBROUTINE EXTERNAL_FORCING_T(
219 I iMin,iMax, jMin,jMax, bi,bj, kLev,
220 I myTime, myThid )
221 C !DESCRIPTION: \bv
222 C *==========================================================*
223 C | S/R EXTERNAL_FORCING_T
224 C | o Contains problem specific forcing for temperature.
225 C *==========================================================*
226 C | Adds terms to gT for forcing by external sources
227 C | e.g. heat flux, climatalogical relaxation, etc ...
228 C *==========================================================*
229 C \ev
230
231 C !USES:
232 IMPLICIT NONE
233 C == Global data ==
234 #include "SIZE.h"
235 #include "EEPARAMS.h"
236 #include "PARAMS.h"
237 #include "GRID.h"
238 #include "DYNVARS.h"
239 #include "FFIELDS.h"
240 #include "SURFACE.h"
241
242 C !INPUT/OUTPUT PARAMETERS:
243 C == Routine arguments ==
244 C iMin,iMax :: Working range of x-index for applying forcing.
245 C jMin,jMax :: Working range of y-index for applying forcing.
246 C bi,bj :: Current tile indices
247 C kLev :: Current vertical level index
248 C myTime :: Current time in simulation
249 C myThid :: Thread Id number
250 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
251 _RL myTime
252 INTEGER myThid
253
254 C !LOCAL VARIABLES:
255 C == Local variables ==
256 C i,j :: Loop counters
257 C kSurface :: index of surface layer
258 INTEGER i, j
259 INTEGER kSurface
260 CEOP
261 #ifdef SHORTWAVE_HEATING
262 integer two
263 _RL minusone
264 parameter (two=2,minusone=-1.)
265 _RL swfracb(two)
266 INTEGER kp1
267 #endif
268
269 IF ( fluidIsAir ) THEN
270 kSurface = 0
271 ELSEIF ( usingPCoords ) THEN
272 kSurface = Nr
273 ELSE
274 kSurface = 1
275 ENDIF
276
277 C-- Forcing term
278 #ifdef ALLOW_AIM
279 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
280 & iMin,iMax, jMin,jMax, bi,bj, kLev,
281 & myTime, myThid )
282 #endif /* ALLOW_AIM */
283
284 #ifdef ALLOW_FIZHI
285 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
286 & iMin,iMax, jMin,jMax, bi,bj, kLev,
287 & myTime, myThid )
288 #endif /* ALLOW_FIZHI */
289
290 #ifdef ALLOW_MYPACKAGE
291 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
292 & iMin,iMax, jMin,jMax, bi,bj, kLev,
293 & myTime, myThid )
294 #endif /* ALLOW_MYPACKAGE */
295
296 C Add heat in top-layer
297 IF ( kLev .EQ. kSurface ) THEN
298 DO j=1,sNy
299 DO i=1,sNx
300 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
301 & +surfaceForcingT(i,j,bi,bj)
302 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
303 ENDDO
304 ENDDO
305 ENDIF
306
307 #ifndef ALLOW_AUTODIFF_TAMC
308 IF (linFSConserveTr) THEN
309 DO j=1,sNy
310 DO i=1,sNx
311 IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
312 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
313 & +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
314 ENDIF
315 ENDDO
316 ENDDO
317 ENDIF
318 #endif /* ndfef ALLOW_AUTODIFF_TAMC */
319
320 #ifdef ALLOW_SHELFICE
321 IF ( useShelfIce )
322 & CALL SHELFICE_FORCING_T(
323 I iMin,iMax, jMin,jMax, bi,bj, kLev,
324 I myTime, myThid )
325 #endif /* ALLOW_SHELFICE */
326
327 #ifdef SHORTWAVE_HEATING
328 C Penetrating SW radiation
329 c IF ( usePenetratingSW ) THEN
330 swfracb(1)=abs(rF(klev))
331 swfracb(2)=abs(rF(klev+1))
332 CALL SWFRAC(
333 I two, minusone,
334 U swfracb,
335 I myTime, 1, myThid )
336 kp1 = klev+1
337 IF (klev.EQ.Nr) THEN
338 kp1 = klev
339 swfracb(2)=0. _d 0
340 ENDIF
341 DO j=1,sNy
342 DO i=1,sNx
343 gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
344 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
345 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
346 & *recip_Cp*recip_rhoConst
347 & *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
348 ENDDO
349 ENDDO
350 c ENDIF
351 #endif
352
353 #ifdef ALLOW_RBCS
354 if (useRBCS) then
355 call RBCS_ADD_TENDENCY(bi,bj,klev, 1,
356 & myTime, myThid )
357 endif
358 #endif
359
360 #ifdef ALLOW_OBCS
361 IF (useOBCS) THEN
362 CALL OBCS_SPONGE_T(
363 I iMin,iMax, jMin,jMax, bi,bj, kLev,
364 I myTime, myThid )
365 ENDIF
366 #endif
367
368 RETURN
369 END
370
371 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
372 CBOP
373 C !ROUTINE: EXTERNAL_FORCING_S
374 C !INTERFACE:
375 SUBROUTINE EXTERNAL_FORCING_S(
376 I iMin,iMax, jMin,jMax, bi,bj, kLev,
377 I myTime, myThid )
378
379 C !DESCRIPTION: \bv
380 C *==========================================================*
381 C | S/R EXTERNAL_FORCING_S
382 C | o Contains problem specific forcing for merid velocity.
383 C *==========================================================*
384 C | Adds terms to gS for forcing by external sources
385 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
386 C *==========================================================*
387 C \ev
388
389 C !USES:
390 IMPLICIT NONE
391 C == Global data ==
392 #include "SIZE.h"
393 #include "EEPARAMS.h"
394 #include "PARAMS.h"
395 #include "GRID.h"
396 #include "DYNVARS.h"
397 #include "FFIELDS.h"
398 #include "SURFACE.h"
399 #ifdef ALLOW_SALT_PLUME
400 #ifdef ALLOW_SEAICE
401 #include "SEAICE.h"
402 #endif /* ALLOW_SEAICE */
403 #endif /* ALLOW_SALT_PLUME */
404
405 C !INPUT/OUTPUT PARAMETERS:
406 C == Routine arguments ==
407 C iMin,iMax :: Working range of x-index for applying forcing.
408 C jMin,jMax :: Working range of y-index for applying forcing.
409 C bi,bj :: Current tile indices
410 C kLev :: Current vertical level index
411 C myTime :: Current time in simulation
412 C myThid :: Thread Id number
413 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
414 _RL myTime
415 INTEGER myThid
416
417 C !LOCAL VARIABLES:
418 C == Local variables ==
419 C i,j :: Loop counters
420 C kSurface :: index of surface layer
421 INTEGER i, j
422 INTEGER kSurface
423 CEOP
424 #ifdef ALLOW_SALT_PLUME
425 _RL saltPlume
426 #endif /* ALLOW_SALT_PLUME */
427
428 IF ( fluidIsAir ) THEN
429 kSurface = 0
430 ELSEIF ( usingPCoords ) THEN
431 kSurface = Nr
432 ELSE
433 kSurface = 1
434 ENDIF
435
436 C-- Forcing term
437 #ifdef ALLOW_AIM
438 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
439 & iMin,iMax, jMin,jMax, bi,bj, kLev,
440 & myTime, myThid )
441 #endif /* ALLOW_AIM */
442
443 #ifdef ALLOW_FIZHI
444 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
445 & iMin,iMax, jMin,jMax, bi,bj, kLev,
446 & myTime, myThid )
447 #endif /* ALLOW_FIZHI */
448
449 #ifdef ALLOW_MYPACKAGE
450 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
451 & iMin,iMax, jMin,jMax, bi,bj, kLev,
452 & myTime, myThid )
453 #endif /* ALLOW_MYPACKAGE */
454
455 C Add fresh-water in top-layer
456 IF ( kLev .EQ. kSurface ) THEN
457 DO j=1,sNy
458 DO i=1,sNx
459 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
460 & +surfaceForcingS(i,j,bi,bj)
461 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
462 ENDDO
463 ENDDO
464 ENDIF
465
466 #ifndef ALLOW_AUTODIFF_TAMC
467 IF (linFSConserveTr) THEN
468 DO j=1,sNy
469 DO i=1,sNx
470 IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
471 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
472 & +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
473 ENDIF
474 ENDDO
475 ENDDO
476 ENDIF
477 #endif /* ndfef ALLOW_AUTODIFF_TAMC */
478
479 #ifdef ALLOW_SHELFICE
480 IF ( useShelfIce )
481 & CALL SHELFICE_FORCING_S(
482 I iMin,iMax, jMin,jMax, bi,bj, kLev,
483 I myTime, myThid )
484 #endif /* ALLOW_SHELFICE */
485
486 #ifdef ALLOW_SALT_PLUME
487 C saltPlume is the amount of salt rejected by ice while freezing;
488 C it is here redistributed to multiple vertical levels as per
489 C Duffy et al. (GRL 1999)
490 DO j=1,sNy
491 DO i=1,sNx
492 saltPlume = 0.
493 #ifdef ALLOW_SEAICE
494 IF ( saltFlux(i,j,bi,bj) .GT. 0. .AND.
495 & salt(i,j,ks,bi,bj) .GT. SEAICE_salinity ) THEN
496 saltPlume = (salt(i,j,ks,bi,bj)-SEAICE_salinity) *
497 & saltFlux(i,j,bi,bj) / salt(i,j,ks,bi,bj)
498 ENDIF
499 #endif /* ALLOW_SEAICE */
500 IF ( SaltPlumeDepth(i,j,bi,bj) .GT. -rF(kLev) ) THEN
501 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
502 & +saltPlume*horiVertRatio*recip_rhoConst
503 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
504 & *max(drF(kLev),SaltPlumeDepth(i,j,bi,bj)+rF(kLev))
505 & /SaltPlumeDepth(i,j,bi,bj)
506 ENDDO
507 ENDDO
508 #endif /* ALLOW_SALT_PLUME */
509
510 #ifdef ALLOW_RBCS
511 if (useRBCS) then
512 call RBCS_ADD_TENDENCY(bi,bj,klev, 2,
513 & myTime, myThid )
514 endif
515 #endif /* ALLOW_RBCS */
516
517 #ifdef ALLOW_OBCS
518 IF (useOBCS) THEN
519 CALL OBCS_SPONGE_S(
520 I iMin,iMax, jMin,jMax, bi,bj, kLev,
521 I myTime, myThid )
522 ENDIF
523 #endif /* ALLOW_OBCS */
524
525 RETURN
526 END

  ViewVC Help
Powered by ViewVC 1.1.22