/[MITgcm]/MITgcm_contrib/jscott/code_rafmod/external_forcing.F
ViewVC logotype

Annotation of /MITgcm_contrib/jscott/code_rafmod/external_forcing.F

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


Revision 1.1 - (hide annotations) (download)
Tue Aug 21 16:34:17 2007 UTC (17 years, 11 months ago) by jscott
Branch: MAIN
code directory for crude ML horizontal mixing scheme

1 jscott 1.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    
400     C !INPUT/OUTPUT PARAMETERS:
401     C == Routine arguments ==
402     C iMin,iMax :: Working range of x-index for applying forcing.
403     C jMin,jMax :: Working range of y-index for applying forcing.
404     C bi,bj :: Current tile indices
405     C kLev :: Current vertical level index
406     C myTime :: Current time in simulation
407     C myThid :: Thread Id number
408     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
409     _RL myTime
410     INTEGER myThid
411    
412     C !LOCAL VARIABLES:
413     C == Local variables ==
414     C i,j :: Loop counters
415     C kSurface :: index of surface layer
416     INTEGER i, j
417     INTEGER kSurface
418     cjrs .04 Sv outflow, 2nd number is surface area of points 88,31:32)
419     _RL outflow
420     PARAMETER( outflow = 4.0d4/3.197336771082545d11)
421     CEOP
422    
423     IF ( fluidIsAir ) THEN
424     kSurface = 0
425     ELSEIF ( usingPCoords ) THEN
426     kSurface = Nr
427     ELSE
428     kSurface = 1
429     ENDIF
430    
431     C-- Forcing term
432     #ifdef ALLOW_AIM
433     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
434     & iMin,iMax, jMin,jMax, bi,bj, kLev,
435     & myTime, myThid )
436     #endif /* ALLOW_AIM */
437    
438     #ifdef ALLOW_FIZHI
439     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
440     & iMin,iMax, jMin,jMax, bi,bj, kLev,
441     & myTime, myThid )
442     #endif /* ALLOW_FIZHI */
443    
444     #ifdef ALLOW_MYPACKAGE
445     IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
446     & iMin,iMax, jMin,jMax, bi,bj, kLev,
447     & myTime, myThid )
448     #endif /* ALLOW_MYPACKAGE */
449    
450     if (kLev.eq.5) then
451     gS(88,31,kLev,bi,bj)=gS(88,31,kLev,bi,bj) +
452     & outflow/190.d0*convertFW2Salt
453     gS(88,32,kLev,bi,bj)=gS(88,32,kLev,bi,bj) +
454     & outflow/190.d0*convertFW2Salt
455     endif
456     C Add fresh-water in top-layer
457     IF ( kLev .EQ. kSurface ) THEN
458     DO j=1,sNy
459     DO i=1,sNx
460     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
461     & +surfaceForcingS(i,j,bi,bj)
462     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
463     ENDDO
464     ENDDO
465     ENDIF
466    
467     #ifndef ALLOW_AUTODIFF_TAMC
468     IF (linFSConserveTr) THEN
469     DO j=1,sNy
470     DO i=1,sNx
471     IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
472     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
473     & +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
474     ENDIF
475     ENDDO
476     ENDDO
477     ENDIF
478     #endif /* ndfef ALLOW_AUTODIFF_TAMC */
479    
480     #ifdef ALLOW_SHELFICE
481     IF ( useShelfIce )
482     & CALL SHELFICE_FORCING_S(
483     I iMin,iMax, jMin,jMax, bi,bj, kLev,
484     I myTime, myThid )
485     #endif /* ALLOW_SHELFICE */
486    
487     #ifdef ALLOW_RBCS
488     if (useRBCS) then
489     call RBCS_ADD_TENDENCY(bi,bj,klev, 2,
490     & myTime, myThid )
491     endif
492     #endif
493    
494     #ifdef ALLOW_OBCS
495     IF (useOBCS) THEN
496     CALL OBCS_SPONGE_S(
497     I iMin,iMax, jMin,jMax, bi,bj, kLev,
498     I myTime, myThid )
499     ENDIF
500     #endif
501    
502     RETURN
503     END

  ViewVC Help
Powered by ViewVC 1.1.22