/[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.1 - (hide annotations) (download)
Tue Aug 5 21:22:43 2003 UTC (21 years, 11 months ago) by cnh
Branch: MAIN
Adding set of files for 1/12 Atlantic configuration

1 cnh 1.1 C $Header: /u/u0/gcmpack/MITgcm/model/src/external_forcing.F,v 1.19 2003/06/19 15:00:45 heimbach Exp $
2     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     & +foFacMom*surfaceTendencyU(i,j,bi,bj)*0.0100D0
73     & *_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     DO j=jMin,jMax
90     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     & +foFacMom*surfaceTendencyV(i,j,bi,bj)*0.0100D0
183     & *_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     DO j=jMin,jMax
208     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     #ifdef SHORTWAVE_HEATING
303     C Penetrating SW radiation
304     swfracb(1)=abs(rF(klev))
305     swfracb(2)=abs(rF(klev+1))
306     call SWFRAC(
307     I two,minusone,
308     I myCurrentTime,myThid,
309     U swfracb)
310     DO j=jMin,jMax
311     DO i=iMin,iMax
312     gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
313     & -maskC(i,j,klev,bi,bj)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))
314     & *recip_Cp*recip_rhoConst*recip_drF(klev)
315     ENDDO
316     ENDDO
317     #endif
318    
319     #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
320     IF (useOBCS) THEN
321     CALL OBCS_SPONGE_T(
322     I iMin, iMax, jMin, jMax,bi,bj,kLev,
323     I myCurrentTime,myThid)
324     ENDIF
325     #endif
326    
327     C-- Create a sponge layer where flow is linearly damped over entire water column
328     C Damping time scale decreases away from boundary so that
329     C tau = 1 day, 5days, 10days, 20days, 60days
330     spWidth = 5
331     recip_tauSp(1) = 1./(86400.*1. )
332     recip_tauSp(2) = 1./(86400.*5. )
333     recip_tauSp(3) = 1./(86400.*10.)
334     recip_tauSp(4) = 1./(86400.*20.)
335     recip_tauSp(5) = 1./(86400.*60.)
336     jSouthBndy = 5
337     jNorthBndy = ny-5+1
338     DO j=jMin,jMax
339     DO i=iMin,iMax
340     jG = myYGlobalLo+(bj-1)*sNy+j-1
341     jFromNBndy = jNorthBndy-jG
342     jFromSBndy = jSouthBndy-jG
343     curRecipTau=0.
344     IF ( jFromNBndy .LE. 0 ) THEN
345     curRecipTau = recip_tauSp(jFromNBndy+5)
346     ENDIF
347     IF ( jFromSBndy .GE. 0 ) THEN
348     curRecipTau = recip_tauSp(-(jFromSBndy-5))
349     ENDIF
350     gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
351     & -curRecipTau*(theta(i,j,Klev,bi,bj)-thetaRef(i,j,kLev,bi,bj))
352     C & *0.0000D0
353     ENDDO
354     ENDDO
355    
356    
357     RETURN
358     END
359     CBOP
360     C !ROUTINE: EXTERNAL_FORCING_S
361     C !INTERFACE:
362     SUBROUTINE EXTERNAL_FORCING_S(
363     I iMin, iMax, jMin, jMax,bi,bj,kLev,
364     I myCurrentTime,myThid)
365    
366     C !DESCRIPTION: \bv
367     C *==========================================================*
368     C | S/R EXTERNAL_FORCING_S
369     C | o Contains problem specific forcing for merid velocity.
370     C *==========================================================*
371     C | Adds terms to gS for forcing by external sources
372     C | e.g. fresh-water flux, climatalogical relaxation.......
373     C *==========================================================*
374     C \ev
375    
376     C !USES:
377     IMPLICIT NONE
378     C == Global data ==
379     #include "SIZE.h"
380     #include "EEPARAMS.h"
381     #include "PARAMS.h"
382     #include "GRID.h"
383     #include "DYNVARS.h"
384     #include "FFIELDS.h"
385    
386     C !INPUT/OUTPUT PARAMETERS:
387     C == Routine arguments ==
388     C iMin - Working range of tile for applying forcing.
389     C iMax
390     C jMin
391     C jMax
392     C kLev
393     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
394     _RL myCurrentTime
395     INTEGER myThid
396    
397     C !LOCAL VARIABLES:
398     C == Local variables ==
399     C Loop counters
400     INTEGER I, J
401     C number of surface interface layer
402     INTEGER kSurface
403     C Cheap sponge layer
404     _RS recip_tauSp(5)
405     INTEGER spWidth
406     _RS curRecipTau
407     INTEGER jFromNBndy, jFromSBndy,
408     & jNorthBndy, jSouthBndy, jG
409    
410     CEOP
411    
412     if ( buoyancyRelation .eq. 'OCEANICP' ) then
413     kSurface = Nr
414     else
415     kSurface = 1
416     endif
417    
418    
419     C-- Forcing term
420     C Add fresh-water in top-layer
421     IF ( kLev .EQ. kSurface ) THEN
422     DO j=jMin,jMax
423     DO i=iMin,iMax
424     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
425     & +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj)
426     ENDDO
427     ENDDO
428     ENDIF
429    
430     #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
431     IF (useOBCS) THEN
432     CALL OBCS_SPONGE_S(
433     I iMin, iMax, jMin, jMax,bi,bj,kLev,
434     I myCurrentTime,myThid)
435     ENDIF
436     #endif
437    
438     C-- Create a sponge layer where flow is linearly damped over entire water column
439     C Damping time scale decreases away from boundary so that
440     C tau = 1 day, 5days, 10days, 20days, 60days
441     spWidth = 5
442     recip_tauSp(1) = 1./(86400.*1. )
443     recip_tauSp(2) = 1./(86400.*5. )
444     recip_tauSp(3) = 1./(86400.*10.)
445     recip_tauSp(4) = 1./(86400.*20.)
446     recip_tauSp(5) = 1./(86400.*60.)
447     jSouthBndy = 5
448     jNorthBndy = ny-5+1
449     DO j=jMin,jMax
450     DO i=iMin,iMax
451     jG = myYGlobalLo+(bj-1)*sNy+j-1
452     jFromNBndy = jNorthBndy-jG
453     jFromSBndy = jSouthBndy-jG
454     curRecipTau=0.
455     IF ( jFromNBndy .LE. 0 ) THEN
456     curRecipTau = recip_tauSp(jFromNBndy+5)
457     ENDIF
458     IF ( jFromSBndy .GE. 0 ) THEN
459     curRecipTau = recip_tauSp(-(jFromSBndy-5))
460     ENDIF
461     gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
462     & -curRecipTau*(salt(i,j,Klev,bi,bj)-saltRef(i,j,kLev,bi,bj))
463     C & *0.0000D0
464     ENDDO
465     ENDDO
466    
467     RETURN
468     END

  ViewVC Help
Powered by ViewVC 1.1.22