/[MITgcm]/MITgcm/verification/hs94.1x64x5/code_ad/apply_forcing.F
ViewVC logotype

Annotation of /MITgcm/verification/hs94.1x64x5/code_ad/apply_forcing.F

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


Revision 1.2 - (hide annotations) (download)
Wed Aug 20 20:27:43 2014 UTC (9 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65c, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.1: +81 -20 lines
use the same version of "apply_forcing.F" (with Sigma-coords code)
 from experiment hs94.cs-32x32x5.

1 jmc 1.2 C $Header: /u/gcmpack/MITgcm/verification/hs94.1x64x5/code/apply_forcing.F,v 1.1 2014/07/11 18:57:32 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     C-- File apply_forcing.F:
7     C-- Contents
8     C-- o APPLY_FORCING_U
9     C-- o APPLY_FORCING_V
10     C-- o APPLY_FORCING_T
11     C-- o APPLY_FORCING_S
12    
13     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
14     CBOP
15     C !ROUTINE: APPLY_FORCING_U
16     C !INTERFACE:
17     SUBROUTINE APPLY_FORCING_U(
18     U gU_arr,
19     I iMin,iMax,jMin,jMax, k, bi, bj,
20     I myTime, myIter, myThid )
21     C !DESCRIPTION: \bv
22     C *==========================================================*
23     C | S/R APPLY_FORCING_U
24     C | o Contains problem specific forcing for zonal velocity.
25     C *==========================================================*
26     C | Adds terms to gU for forcing by external sources
27     C | e.g. wind stress, bottom friction etc ...
28     C *==========================================================*
29     C \ev
30    
31     C !USES:
32     IMPLICIT NONE
33     C == Global data ==
34     #include "SIZE.h"
35     #include "EEPARAMS.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38 jmc 1.2 #include "SURFACE.h"
39 jmc 1.1 #include "DYNVARS.h"
40     #include "FFIELDS.h"
41    
42     C !INPUT/OUTPUT PARAMETERS:
43     C gU_arr :: the tendency array
44     C iMin,iMax :: Working range of x-index for applying forcing.
45     C jMin,jMax :: Working range of y-index for applying forcing.
46     C k :: Current vertical level index
47     C bi,bj :: Current tile indices
48     C myTime :: Current time in simulation
49     C myIter :: Current iteration number
50     C myThid :: my Thread Id number
51     _RL gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52     INTEGER iMin, iMax, jMin, jMax
53     INTEGER k, bi, bj
54     _RL myTime
55     INTEGER myIter
56     INTEGER myThid
57    
58     C !LOCAL VARIABLES:
59     C i,j :: Loop counters
60     INTEGER i, j
61     CEOP
62 jmc 1.2 _RL recip_P0g, termP, rFullDepth
63     _RL kV, kF, sigma_b
64 jmc 1.1
65     C-- Forcing term
66     kF = 1. _d 0/86400. _d 0
67     sigma_b = 0.7 _d 0
68 jmc 1.2 rFullDepth = rF(1)-rF(Nr+1)
69 jmc 1.1 c DO j=1,sNy
70     C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
71     DO j=0,sNy+1
72     DO i=1,sNx+1
73     IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN
74 jmc 1.2 IF ( selectSigmaCoord.EQ.0 ) THEN
75     recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
76     termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
77     & +rF(k+1)*recip_P0g )
78     c termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
79     ELSE
80     C-- Pressure at U.point :
81     c midP = rLowW(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
82     c & + bHybSigmC(k)
83     c & *(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
84     C-- Sigma at U.point :
85     c termP = ( midP - rLowW(i,j,bi,bj))
86     c & /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
87     C- which simplifies to:
88     termP = aHybSigmC(k)*rFullDepth
89     #ifdef NONLIN_FRSURF
90     & /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
91     #else
92     & /(rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
93     #endif
94     & + bHybSigmC(k)
95     ENDIF
96 jmc 1.1 kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
97     gU_arr(i,j) = gU_arr(i,j)
98     & - kV*uVel(i,j,k,bi,bj)
99     ENDIF
100     ENDDO
101     ENDDO
102    
103     RETURN
104     END
105    
106     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
107     CBOP
108     C !ROUTINE: APPLY_FORCING_V
109     C !INTERFACE:
110     SUBROUTINE APPLY_FORCING_V(
111     U gV_arr,
112     I iMin,iMax,jMin,jMax, k, bi, bj,
113     I myTime, myIter, myThid )
114     C !DESCRIPTION: \bv
115     C *==========================================================*
116     C | S/R APPLY_FORCING_V
117     C | o Contains problem specific forcing for merid velocity.
118     C *==========================================================*
119     C | Adds terms to gV for forcing by external sources
120     C | e.g. wind stress, bottom friction etc ...
121     C *==========================================================*
122     C \ev
123    
124     C !USES:
125     IMPLICIT NONE
126     C == Global data ==
127     #include "SIZE.h"
128     #include "EEPARAMS.h"
129     #include "PARAMS.h"
130     #include "GRID.h"
131 jmc 1.2 #include "SURFACE.h"
132 jmc 1.1 #include "DYNVARS.h"
133     #include "FFIELDS.h"
134    
135     C !INPUT/OUTPUT PARAMETERS:
136     C gV_arr :: the tendency array
137     C iMin,iMax :: Working range of x-index for applying forcing.
138     C jMin,jMax :: Working range of y-index for applying forcing.
139     C k :: Current vertical level index
140     C bi,bj :: Current tile indices
141     C myTime :: Current time in simulation
142     C myIter :: Current iteration number
143     C myThid :: my Thread Id number
144     _RL gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145     INTEGER iMin, iMax, jMin, jMax
146     INTEGER k, bi, bj
147     _RL myTime
148     INTEGER myIter
149     INTEGER myThid
150    
151     C !LOCAL VARIABLES:
152     C i,j :: Loop counters
153     INTEGER i, j
154     CEOP
155 jmc 1.2 _RL recip_P0g, termP, rFullDepth
156     _RL kV, kF, sigma_b
157 jmc 1.1
158     C-- Forcing term
159     kF = 1. _d 0/86400. _d 0
160     sigma_b = 0.7 _d 0
161 jmc 1.2 rFullDepth = rF(1)-rF(Nr+1)
162 jmc 1.1 DO j=1,sNy+1
163     c DO i=1,sNx
164     C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
165     DO i=0,sNx+1
166     IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN
167 jmc 1.2 IF ( selectSigmaCoord.EQ.0 ) THEN
168     recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
169     termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
170     & +rF(k+1)*recip_P0g )
171     c termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
172     ELSE
173     C-- Pressure at V.point :
174     c midP = rLowS(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
175     c & + bHybSigmC(k)
176     c & *(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
177     C-- Sigma at V.point :
178     c termP = ( midP - rLowS(i,j,bi,bj))
179     c & /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
180     C- which simplifies to:
181     termP = aHybSigmC(k)*rFullDepth
182     #ifdef NONLIN_FRSURF
183     & /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
184     #else
185     & /(rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
186     #endif
187     & + bHybSigmC(k)
188     ENDIF
189 jmc 1.1 kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
190     gV_arr(i,j) = gV_arr(i,j)
191     & - kV*vVel(i,j,k,bi,bj)
192     ENDIF
193     ENDDO
194     ENDDO
195    
196     RETURN
197     END
198    
199     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
200     CBOP
201     C !ROUTINE: APPLY_FORCING_T
202     C !INTERFACE:
203     SUBROUTINE APPLY_FORCING_T(
204     U gT_arr,
205     I iMin,iMax,jMin,jMax, k, bi, bj,
206     I myTime, myIter, myThid )
207     C !DESCRIPTION: \bv
208     C *==========================================================*
209     C | S/R APPLY_FORCING_T
210     C | o Contains problem specific forcing for temperature.
211     C *==========================================================*
212     C | Adds terms to gT for forcing by external sources
213     C | e.g. heat flux, climatalogical relaxation, etc ...
214     C *==========================================================*
215     C \ev
216    
217     C !USES:
218     IMPLICIT NONE
219     C == Global data ==
220     #include "SIZE.h"
221     #include "EEPARAMS.h"
222     #include "PARAMS.h"
223     #include "GRID.h"
224     #include "DYNVARS.h"
225     #include "FFIELDS.h"
226    
227     C !INPUT/OUTPUT PARAMETERS:
228     C gT_arr :: the tendency array
229     C iMin,iMax :: Working range of x-index for applying forcing.
230     C jMin,jMax :: Working range of y-index for applying forcing.
231     C k :: Current vertical level index
232     C bi,bj :: Current tile indices
233     C myTime :: Current time in simulation
234     C myIter :: Current iteration number
235     C myThid :: my Thread Id number
236     _RL gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
237     INTEGER iMin, iMax, jMin, jMax
238     INTEGER k, bi, bj
239     _RL myTime
240     INTEGER myIter
241     INTEGER myThid
242    
243     C !LOCAL VARIABLES:
244     C i,j :: Loop counters
245     C kSurface :: index of surface level
246     INTEGER i, j
247     CEOP
248 jmc 1.2 _RL thetaLim, kT, ka, ks, sigma_b, term1, term2, thetaEq
249     _RL termP, rFullDepth
250 jmc 1.1
251     C-- Forcing term
252     ka = 1. _d 0/(40. _d 0*86400. _d 0)
253     ks = 1. _d 0/(4. _d 0 *86400. _d 0)
254     sigma_b = 0.7 _d 0
255 jmc 1.2 rFullDepth = rF(1)-rF(Nr+1)
256     DO j=0,sNy+1
257     DO i=0,sNx+1
258 jmc 1.1 term1 = 60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
259     termP = 0.5 _d 0*( rF(k) + rF(k+1) )
260     term2 = 10. _d 0*LOG(termP/atm_po)
261     & *(COS(yC(i,j,bi,bj)*deg2rad)**2)
262     thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
263 jmc 1.2 thetaEq = 315. _d 0 - term1 - term2
264 jmc 1.1 thetaEq = MAX(thetaLim,thetaEq)
265 jmc 1.2 IF ( selectSigmaCoord.EQ.0 ) THEN
266     termP = 0.5 _d 0*( MIN(rF(k),Ro_surf(i,j,bi,bj))
267     & + rF(k+1) )
268     & *recip_Rcol(i,j,bi,bj)
269     ELSE
270     C-- Pressure at T.point :
271     c midP = R_low(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
272     c & + bHybSigmC(k)
273     c & *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
274     C-- Sigma at T.point :
275     c termP = ( midP - R_low(i,j,bi,bj))
276     c & /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
277     C- which simplifies to:
278     termP = aHybSigmC(k)*rFullDepth
279     #ifdef NONLIN_FRSURF
280     & /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
281     #else
282     & /(Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
283     #endif
284     & + bHybSigmC(k)
285     ENDIF
286 jmc 1.1 kT = ka+(ks-ka)
287 jmc 1.2 & *MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
288 jmc 1.1 & *COS((yC(i,j,bi,bj)*deg2rad))**4
289     gT_arr(i,j) = gT_arr(i,j)
290     & - kT*( theta(i,j,k,bi,bj)-thetaEq )
291     & *maskC(i,j,k,bi,bj)
292     ENDDO
293     ENDDO
294    
295     RETURN
296     END
297    
298     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
299     CBOP
300     C !ROUTINE: APPLY_FORCING_S
301     C !INTERFACE:
302     SUBROUTINE APPLY_FORCING_S(
303     U gS_arr,
304     I iMin,iMax,jMin,jMax, k, bi, bj,
305     I myTime, myIter, myThid )
306     C !DESCRIPTION: \bv
307     C *==========================================================*
308     C | S/R APPLY_FORCING_S
309     C | o Contains problem specific forcing for merid velocity.
310     C *==========================================================*
311     C | Adds terms to gS for forcing by external sources
312     C | e.g. fresh-water flux, climatalogical relaxation, etc ...
313     C *==========================================================*
314     C \ev
315    
316     C !USES:
317     IMPLICIT NONE
318     C == Global data ==
319     #include "SIZE.h"
320     #include "EEPARAMS.h"
321     #include "PARAMS.h"
322     #include "GRID.h"
323     #include "DYNVARS.h"
324     #include "FFIELDS.h"
325     #include "SURFACE.h"
326    
327     C !INPUT/OUTPUT PARAMETERS:
328     C gS_arr :: the tendency array
329     C iMin,iMax :: Working range of x-index for applying forcing.
330     C jMin,jMax :: Working range of y-index for applying forcing.
331     C k :: Current vertical level index
332     C bi,bj :: Current tile indices
333     C myTime :: Current time in simulation
334     C myIter :: Current iteration number
335     C myThid :: my Thread Id number
336     _RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
337     INTEGER iMin, iMax, jMin, jMax
338     INTEGER k, bi, bj
339     _RL myTime
340     INTEGER myIter
341     INTEGER myThid
342    
343     C !LOCAL VARIABLES:
344     C i,j :: Loop counters
345     c INTEGER i, j
346     CEOP
347    
348     C-- Forcing term
349    
350     RETURN
351     END

  ViewVC Help
Powered by ViewVC 1.1.22