/[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.58 - (show annotations) (download)
Sat May 14 20:01:33 2011 UTC (13 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint62y, checkpoint63, checkpoint63a
Changes since 1.57: +15 -1 lines
RBCS: add capability to apply relaxation to horizontal velocity uVel & vVel

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

  ViewVC Help
Powered by ViewVC 1.1.22