/[MITgcm]/MITgcm/verification/hs94.cs-32x32x5/code/apply_forcing.F
ViewVC logotype

Contents of /MITgcm/verification/hs94.cs-32x32x5/code/apply_forcing.F

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


Revision 1.1 - (show annotations) (download)
Wed Aug 20 20:23:10 2014 UTC (5 years, 9 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
convert Held & Suarez external_forcing.F (from here, with Sigma-coords
 code) to apply_forcing.F and undef USE_OLD_EXTERNAL_FORCING

1 C $Header: /u/gcmpack/MITgcm/verification/hs94.1x64x5/code/apply_forcing.F,v 1.1 2014/07/11 18:57:32 jmc Exp $
2 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 #include "SURFACE.h"
39 #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 _RL recip_P0g, termP, rFullDepth
63 _RL kV, kF, sigma_b
64
65 C-- Forcing term
66 kF = 1. _d 0/86400. _d 0
67 sigma_b = 0.7 _d 0
68 rFullDepth = rF(1)-rF(Nr+1)
69 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 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 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 #include "SURFACE.h"
132 #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 _RL recip_P0g, termP, rFullDepth
156 _RL kV, kF, sigma_b
157
158 C-- Forcing term
159 kF = 1. _d 0/86400. _d 0
160 sigma_b = 0.7 _d 0
161 rFullDepth = rF(1)-rF(Nr+1)
162 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 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 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 _RL thetaLim, kT, ka, ks, sigma_b, term1, term2, thetaEq
249 _RL termP, rFullDepth
250
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 rFullDepth = rF(1)-rF(Nr+1)
256 DO j=0,sNy+1
257 DO i=0,sNx+1
258 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 thetaEq = 315. _d 0 - term1 - term2
264 thetaEq = MAX(thetaLim,thetaEq)
265 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 kT = ka+(ks-ka)
287 & *MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
288 & *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