/[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.55 - (show annotations) (download)
Wed Jan 20 23:33:45 2010 UTC (14 years, 4 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint62b
Changes since 1.54: +13 -1 lines
Adding the shell of, and the hooks to, a new package that will be
used to model melting and freezing of vertical glacier ice fronts:
 Modified Files:
 	doc/tag-index model/inc/PARAMS.h model/src/do_oceanic_phys.F
 	model/src/external_forcing.F model/src/packages_boot.F
 	model/src/packages_check.F model/src/packages_init_fixed.F
 	model/src/packages_init_variables.F
 	model/src/packages_readparms.F
 Added Files:
 	pkg/icefront/ICEFRONT.h pkg/icefront/ICEFRONT_OPTIONS.h
 	pkg/icefront/icefront_check.F
 	pkg/icefront/icefront_description.tex
 	pkg/icefront/icefront_diagnostics_init.F
 	pkg/icefront/icefront_init_fixed.F
 	pkg/icefront/icefront_init_varia.F
 	pkg/icefront/icefront_readparms.F
 	pkg/icefront/icefront_tendency_apply.F
 	pkg/icefront/icefront_thermodynamics.F

1 C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.54 2008/08/24 21:46:19 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 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_OBCS
94 IF (useOBCS) THEN
95 CALL OBCS_SPONGE_U(
96 I iMin,iMax, jMin,jMax, bi,bj, kLev,
97 I myTime, myThid )
98 ENDIF
99 #endif
100
101 #ifdef ALLOW_MYPACKAGE
102 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_U(
103 & iMin,iMax, jMin,jMax, bi,bj, kLev,
104 & myTime, myThid )
105 #endif /* ALLOW_MYPACKAGE */
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 C Add windstress momentum impulse into the top-layer
179 IF ( kLev .EQ. kSurface ) THEN
180 DO j=1,sNy+1
181 c DO i=1,sNx
182 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
183 DO i=0,sNx+1
184 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
185 & +foFacMom*surfaceForcingV(i,j,bi,bj)
186 & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
187 ENDDO
188 ENDDO
189 ENDIF
190
191 #ifdef ALLOW_EDDYPSI
192 CALL TAUEDDY_EXTERNAL_FORCING_V(
193 I iMin,iMax, jMin,jMax, bi,bj, kLev,
194 I myTime, myThid )
195 #endif
196
197 #ifdef ALLOW_OBCS
198 IF (useOBCS) THEN
199 CALL OBCS_SPONGE_V(
200 I iMin,iMax, jMin,jMax, bi,bj, kLev,
201 I myTime, myThid )
202 ENDIF
203 #endif
204
205 #ifdef ALLOW_MYPACKAGE
206 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_V(
207 & iMin,iMax, jMin,jMax, bi,bj, kLev,
208 & myTime, myThid )
209 #endif /* ALLOW_MYPACKAGE */
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_ADDFLUID
291 IF ( selectAddFluid.NE.0 .AND. temp_EvPrRn.NE.UNSET_RL ) THEN
292 C- for now, use same fluid properties as for E-P-R
293 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
294 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
295 DO j=1,sNy
296 DO i=1,sNx
297 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
298 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
299 & *( temp_EvPrRn - theta(i,j,kLev,bi,bj) )
300 & *recip_rA(i,j,bi,bj)
301 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
302 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
303 ENDDO
304 ENDDO
305 ELSE
306 DO j=1,sNy
307 DO i=1,sNx
308 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
309 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
310 & *( temp_EvPrRn - tRef(kLev) )
311 & *recip_rA(i,j,bi,bj)
312 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
313 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
314 ENDDO
315 ENDDO
316 ENDIF
317 ENDIF
318 #endif /* ALLOW_ADDFLUID */
319
320 C Add heat in top-layer
321 IF ( kLev .EQ. kSurface ) THEN
322 DO j=1,sNy
323 DO i=1,sNx
324 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
325 & +surfaceForcingT(i,j,bi,bj)
326 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
327 ENDDO
328 ENDDO
329 ENDIF
330
331 IF (linFSConserveTr) THEN
332 DO j=1,sNy
333 DO i=1,sNx
334 IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
335 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
336 & +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
337 ENDIF
338 ENDDO
339 ENDDO
340 ENDIF
341
342 #ifdef ALLOW_SHELFICE
343 IF ( useShelfIce )
344 & CALL SHELFICE_FORCING_T(
345 I iMin,iMax, jMin,jMax, bi,bj, kLev,
346 I myTime, myThid )
347 #endif /* ALLOW_SHELFICE */
348
349 #ifdef ALLOW_ICEFRONT
350 IF ( useICEFRONT ) CALL ICEFRONT_TENDENCY_APPLY_T(
351 & iMin,iMax, jMin,jMax, bi,bj, kLev,
352 & myTime, myThid )
353 #endif /* ALLOW_ICEFRONT */
354
355 #ifdef SHORTWAVE_HEATING
356 C Penetrating SW radiation
357 c IF ( usePenetratingSW ) THEN
358 swfracb(1)=abs(rF(klev))
359 swfracb(2)=abs(rF(klev+1))
360 CALL SWFRAC(
361 I two, minusone,
362 U swfracb,
363 I myTime, 1, myThid )
364 kp1 = klev+1
365 IF (klev.EQ.Nr) THEN
366 kp1 = klev
367 swfracb(2)=0. _d 0
368 ENDIF
369 DO j=1,sNy
370 DO i=1,sNx
371 gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
372 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
373 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
374 & *recip_Cp*mass2rUnit
375 & *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
376 ENDDO
377 ENDDO
378 c ENDIF
379 #endif
380
381 #ifdef ALLOW_RBCS
382 IF (useRBCS) THEN
383 CALL RBCS_ADD_TENDENCY(bi,bj,klev, 1,
384 & myTime, myThid )
385 ENDIF
386 #endif
387
388 #ifdef ALLOW_OBCS
389 IF (useOBCS) THEN
390 CALL OBCS_SPONGE_T(
391 I iMin,iMax, jMin,jMax, bi,bj, kLev,
392 I myTime, myThid )
393 ENDIF
394 #endif
395
396 #ifdef ALLOW_MYPACKAGE
397 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
398 & iMin,iMax, jMin,jMax, bi,bj, kLev,
399 & myTime, myThid )
400 #endif /* ALLOW_MYPACKAGE */
401
402 RETURN
403 END
404
405 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
406 CBOP
407 C !ROUTINE: EXTERNAL_FORCING_S
408 C !INTERFACE:
409 SUBROUTINE EXTERNAL_FORCING_S(
410 I iMin,iMax, jMin,jMax, bi,bj, kLev,
411 I myTime, myThid )
412
413 C !DESCRIPTION: \bv
414 C *==========================================================*
415 C | S/R EXTERNAL_FORCING_S
416 C | o Contains problem specific forcing for merid velocity.
417 C *==========================================================*
418 C | Adds terms to gS for forcing by external sources
419 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
420 C *==========================================================*
421 C \ev
422
423 C !USES:
424 IMPLICIT NONE
425 C == Global data ==
426 #include "SIZE.h"
427 #include "EEPARAMS.h"
428 #include "PARAMS.h"
429 #include "GRID.h"
430 #include "DYNVARS.h"
431 #include "FFIELDS.h"
432 #include "SURFACE.h"
433
434 C !INPUT/OUTPUT PARAMETERS:
435 C == Routine arguments ==
436 C iMin,iMax :: Working range of x-index for applying forcing.
437 C jMin,jMax :: Working range of y-index for applying forcing.
438 C bi,bj :: Current tile indices
439 C kLev :: Current vertical level index
440 C myTime :: Current time in simulation
441 C myThid :: Thread Id number
442 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
443 _RL myTime
444 INTEGER myThid
445
446 C !LOCAL VARIABLES:
447 C == Local variables ==
448 C i,j :: Loop counters
449 C kSurface :: index of surface layer
450 INTEGER i, j
451 INTEGER kSurface
452 CEOP
453
454 IF ( fluidIsAir ) THEN
455 kSurface = 0
456 ELSEIF ( usingPCoords ) THEN
457 kSurface = Nr
458 ELSE
459 kSurface = 1
460 ENDIF
461
462 C-- Forcing term
463 #ifdef ALLOW_AIM
464 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
465 & iMin,iMax, jMin,jMax, bi,bj, kLev,
466 & myTime, myThid )
467 #endif /* ALLOW_AIM */
468
469 #ifdef ALLOW_FIZHI
470 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
471 & iMin,iMax, jMin,jMax, bi,bj, kLev,
472 & myTime, myThid )
473 #endif /* ALLOW_FIZHI */
474
475 #ifdef ALLOW_ADDFLUID
476 IF ( selectAddFluid.NE.0 .AND. salt_EvPrRn.NE.UNSET_RL ) THEN
477 C- for now, use same fluid properties as for E-P-R
478 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
479 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
480 DO j=1,sNy
481 DO i=1,sNx
482 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
483 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
484 & *( salt_EvPrRn - salt(i,j,kLev,bi,bj) )
485 & *recip_rA(i,j,bi,bj)
486 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
487 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
488 ENDDO
489 ENDDO
490 ELSE
491 DO j=1,sNy
492 DO i=1,sNx
493 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
494 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
495 & *( salt_EvPrRn - sRef(kLev) )
496 & *recip_rA(i,j,bi,bj)
497 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
498 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
499 ENDDO
500 ENDDO
501 ENDIF
502 ENDIF
503 #endif /* ALLOW_ADDFLUID */
504
505 C Add fresh-water in top-layer
506 IF ( kLev .EQ. kSurface ) THEN
507 DO j=1,sNy
508 DO i=1,sNx
509 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
510 & +surfaceForcingS(i,j,bi,bj)
511 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
512 ENDDO
513 ENDDO
514 ENDIF
515
516 IF (linFSConserveTr) THEN
517 DO j=1,sNy
518 DO i=1,sNx
519 IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
520 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
521 & +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
522 ENDIF
523 ENDDO
524 ENDDO
525 ENDIF
526
527 #ifdef ALLOW_SHELFICE
528 IF ( useShelfIce )
529 & CALL SHELFICE_FORCING_S(
530 I iMin,iMax, jMin,jMax, bi,bj, kLev,
531 I myTime, myThid )
532 #endif /* ALLOW_SHELFICE */
533
534 #ifdef ALLOW_ICEFRONT
535 IF ( useICEFRONT ) CALL ICEFRONT_TENDENCY_APPLY_S(
536 & iMin,iMax, jMin,jMax, bi,bj, kLev,
537 & myTime, myThid )
538 #endif /* ALLOW_ICEFRONT */
539
540 #ifdef ALLOW_SALT_PLUME
541 IF ( useSALT_PLUME )
542 & CALL SALT_PLUME_TENDENCY_APPLY_S(
543 I iMin,iMax, jMin,jMax, bi,bj, kLev,
544 I myTime, myThid )
545 #endif /* ALLOW_SALT_PLUME */
546
547 #ifdef ALLOW_RBCS
548 IF (useRBCS) THEN
549 CALL RBCS_ADD_TENDENCY(bi,bj,klev, 2,
550 & myTime, myThid )
551 ENDIF
552 #endif /* ALLOW_RBCS */
553
554 #ifdef ALLOW_OBCS
555 IF (useOBCS) THEN
556 CALL OBCS_SPONGE_S(
557 I iMin,iMax, jMin,jMax, bi,bj, kLev,
558 I myTime, myThid )
559 ENDIF
560 #endif /* ALLOW_OBCS */
561
562 #ifdef ALLOW_MYPACKAGE
563 IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
564 & iMin,iMax, jMin,jMax, bi,bj, kLev,
565 & myTime, myThid )
566 #endif /* ALLOW_MYPACKAGE */
567
568 RETURN
569 END

  ViewVC Help
Powered by ViewVC 1.1.22