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

Annotation 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 - (hide annotations) (download)
Thu Aug 7 13:03:31 2003 UTC (20 years, 10 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 cnh 1.4 C $Header: /u/u0/gcmpack/MITgcm/model/src/external_forcing.F,v 1.19 2003/06/19 15:00:45 heimbach Exp $
2 cnh 1.1 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 cnh 1.2 & +foFacMom*surfaceTendencyU(i,j,bi,bj)
73 cnh 1.1 & *_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 cnh 1.4 DO j=1,sNy
90 cnh 1.1 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 cnh 1.2 & +foFacMom*surfaceTendencyV(i,j,bi,bj)
183 cnh 1.1 & *_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 cnh 1.4 DO j=1,sNy
208 cnh 1.1 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 cnh 1.2 C-- Forcing term
303     C Add heat in top-layer ( 90 day climatalogical average relaxation )
304     IF ( kLev .EQ. kSurface ) THEN
305 cnh 1.3 curRecipTau=1./(86400.*90.)
306 cnh 1.2 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 cnh 1.1 #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 cnh 1.4 DO j=1,sNy
353 cnh 1.1 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 cnh 1.2 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 cnh 1.4 curRecipTau=1./(86400.*90.)
448 cnh 1.2 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 cnh 1.1 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 cnh 1.4 DO j=1,sNy
478 cnh 1.1 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