1 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/step_fizhi_corr.F,v 1.17 2006/06/09 16:05:28 molod Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "FIZHI_OPTIONS.h" |
5 |
subroutine step_fizhi_corr (myTime, myIter, myThid, dt) |
6 |
c---------------------------------------------------------------------- |
7 |
c Subroutine step_fizhi_corr - 'Wrapper' routine to advance |
8 |
c the physics state and make the new value. |
9 |
c At this point, increment with the "correction term" |
10 |
c which includes the dynamics tendency and the integral |
11 |
c constraint to enforce agreement with the dynamics state |
12 |
c Also: Set up "bi, bj loop" and some timers and clocks here. |
13 |
c |
14 |
c Call:phys2dyn (4) (interpolate physics state to dynamics grid |
15 |
c for use in the correction terms) |
16 |
c AtoC (convert physics state on dynamics grid to C-Grid) |
17 |
c CtoA (convert correction term on dynamics grid to A-Grid) |
18 |
c dyn2phys (4) (interpolate A-Grid correction term to physics grid) |
19 |
c step_physics (advance physics state by correction term) |
20 |
c----------------------------------------------------------------------- |
21 |
implicit none |
22 |
#include "SIZE.h" |
23 |
#include "GRID.h" |
24 |
#include "fizhi_SIZE.h" |
25 |
#include "fizhi_land_SIZE.h" |
26 |
#include "DYNVARS.h" |
27 |
#include "fizhi_coms.h" |
28 |
#include "gridalt_mapping.h" |
29 |
#include "EEPARAMS.h" |
30 |
#include "SURFACE.h" |
31 |
#ifdef ALLOW_DIAGNOSTICS |
32 |
#include "fizhi_SHP.h" |
33 |
#endif |
34 |
|
35 |
integer myIter, myThid |
36 |
_RL myTime |
37 |
#ifdef ALLOW_DIAGNOSTICS |
38 |
logical diagnostics_is_on |
39 |
external diagnostics_is_on |
40 |
#endif |
41 |
|
42 |
c pe on dynamics and physics grid refers to bottom edge |
43 |
_RL pephy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy) |
44 |
_RL pedyn(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr+1,nSx,nSy) |
45 |
_RL windphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
46 |
_RL udyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
47 |
_RL vdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
48 |
_RL thdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
49 |
_RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
50 |
_RL uphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
51 |
_RL vphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
52 |
_RL thphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
53 |
_RL sphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
54 |
_RL tempphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
55 |
|
56 |
integer i, j, L, Lbotij, bi, bj |
57 |
integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2 |
58 |
_RL dt |
59 |
|
60 |
_RL tempij(sNx,sNy) |
61 |
_RL dtfake |
62 |
_RL dtinv |
63 |
|
64 |
im1 = 1-OLx |
65 |
im2 = sNx+OLx |
66 |
jm1 = 1-OLy |
67 |
jm2 = sNy+OLy |
68 |
idim1 = 1 |
69 |
idim2 = sNx |
70 |
jdim1 = 1 |
71 |
jdim2 = sNy |
72 |
dtfake = 1. _d 0 |
73 |
dtinv = 1. _d 0 / dt |
74 |
|
75 |
do bj = myByLo(myThid), myByHi(myThid) |
76 |
do bi = myBxLo(myThid), myBxHi(myThid) |
77 |
|
78 |
C Construct Pressures on physics and dynamics grids |
79 |
do j = 1,sNy |
80 |
do i = 1,sNx |
81 |
do L = 1,Nr |
82 |
pedyn(i,j,L,bi,bj) = 0. |
83 |
enddo |
84 |
enddo |
85 |
enddo |
86 |
do j = 1,sNy |
87 |
do i = 1,sNx |
88 |
Lbotij = ksurfC(i,j,bi,bj) |
89 |
if(Lbotij.lt.nr+1) |
90 |
. pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj) |
91 |
enddo |
92 |
enddo |
93 |
do j = 1,sNy |
94 |
do i = 1,sNx |
95 |
Lbotij = ksurfC(i,j,bi,bj) |
96 |
do L = Lbotij+1,Nr+1 |
97 |
pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) - |
98 |
. drF(L-1)* rStarExpC(i,j,bi,bj)*hfacC(i,j,L-1,bi,bj) |
99 |
enddo |
100 |
c Do not use a zero field as the top edge pressure for interpolation |
101 |
if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5) |
102 |
. pedyn(i,j,Nr+1,bi,bj) = 1.e-5 |
103 |
enddo |
104 |
enddo |
105 |
|
106 |
do j = 1,sNy |
107 |
do i = 1,sNx |
108 |
pephy(i,j,1,bi,bj)=Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj) |
109 |
do L = 2,Nrphys+1 |
110 |
pephy(i,j,L,bi,bj)=pephy(i,j,L-1,bi,bj)-dpphys(i,j,L-1,bi,bj) |
111 |
enddo |
112 |
c Do not use a zero field as the top edge pressure for interpolation |
113 |
if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5) |
114 |
. pephy(i,j,Nrphys+1,bi,bj) = 1.e-5 |
115 |
enddo |
116 |
enddo |
117 |
c |
118 |
c Create a wind magnitude field on the physics grid - |
119 |
c (Load the wind speed bottom up for use by dyn2phys) |
120 |
do L = 1,Nrphys |
121 |
do j = 1,sNy |
122 |
do i = 1,sNx |
123 |
windphy(i,j,L,bi,bj) = |
124 |
. sqrt(uphy(i,j,Nrphys+1-L,bi,bj)*uphy(i,j,Nrphys+1-L,bi,bj) |
125 |
. + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj)) |
126 |
enddo |
127 |
enddo |
128 |
enddo |
129 |
enddo |
130 |
enddo |
131 |
|
132 |
CALL TIMER_START('PHYS2DYN [STEP_FIZHI_CORR]',mythid) |
133 |
do bj = myByLo(myThid), myByHi(myThid) |
134 |
do bi = myBxLo(myThid), myBxHi(myThid) |
135 |
|
136 |
c Compute correction term (new dyn state-phys state to dyn) on physics grid: |
137 |
c First: interp physics state to dynamics grid |
138 |
C Note: physics field levels are numbered top down - need bottom up |
139 |
do L = 1,Nrphys |
140 |
do j = 1,sNy |
141 |
do i = 1,sNx |
142 |
tempphy(i,j,Nrphys+1-L,bi,bj) = uphy(i,j,L,bi,bj) |
143 |
enddo |
144 |
enddo |
145 |
enddo |
146 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
147 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,udyntemp) |
148 |
do L = 1,Nrphys |
149 |
do j = 1,sNy |
150 |
do i = 1,sNx |
151 |
tempphy(i,j,Nrphys+1-L,bi,bj) = vphy(i,j,L,bi,bj) |
152 |
enddo |
153 |
enddo |
154 |
enddo |
155 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
156 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,vdyntemp) |
157 |
do L = 1,Nrphys |
158 |
do j = 1,sNy |
159 |
do i = 1,sNx |
160 |
tempphy(i,j,Nrphys+1-L,bi,bj) = thphy(i,j,L,bi,bj) |
161 |
enddo |
162 |
enddo |
163 |
enddo |
164 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
165 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,thdyntemp) |
166 |
do L = 1,Nrphys |
167 |
do j = 1,sNy |
168 |
do i = 1,sNx |
169 |
tempphy(i,j,Nrphys+1-L,bi,bj) = sphy(i,j,L,bi,bj) |
170 |
enddo |
171 |
enddo |
172 |
enddo |
173 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
174 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,sdyntemp) |
175 |
|
176 |
enddo |
177 |
enddo |
178 |
CALL TIMER_STOP('PHYS2DYN [STEP_FIZHI_CORR]',mythid) |
179 |
|
180 |
c Second: Convert physics state on dynamics grid to C-Grid |
181 |
|
182 |
CALL TIMER_START('ATOC [STEP_FIZHI_CORR]',mythid) |
183 |
call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr, |
184 |
. Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp) |
185 |
CALL TIMER_STOP('ATOC [STEP_FIZHI_CORR]',mythid) |
186 |
|
187 |
c Third: Subtract Phys state on dyn. grid from new dynamics state |
188 |
do bj = myByLo(myThid), myByHi(myThid) |
189 |
do bi = myBxLo(myThid), myBxHi(myThid) |
190 |
|
191 |
do L = 1,Nr |
192 |
do j = jdim1,jdim2 |
193 |
do i = idim1,idim2 |
194 |
udyntemp(i,j,L,bi,bj)=uvel(i,j,L,bi,bj)-udyntemp(i,j,L,bi,bj) |
195 |
vdyntemp(i,j,L,bi,bj)=vvel(i,j,L,bi,bj)-vdyntemp(i,j,L,bi,bj) |
196 |
thdyntemp(i,j,L,bi,bj)=theta(i,j,L,bi,bj)-thdyntemp(i,j,L,bi,bj) |
197 |
sdyntemp(i,j,L,bi,bj)=salt(i,j,L,bi,bj)-sdyntemp(i,j,L,bi,bj) |
198 |
enddo |
199 |
enddo |
200 |
enddo |
201 |
|
202 |
enddo |
203 |
enddo |
204 |
|
205 |
c Fourth: Convert correction terms to A-Grid |
206 |
CALL TIMER_START('CTOA [STEP_FIZHI_CORR]',mythid) |
207 |
call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2, |
208 |
. Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp) |
209 |
CALL TIMER_STOP('CTOA [STEP_FIZHI_CORR]',mythid) |
210 |
|
211 |
c Fifth: Interpolate correction terms to physics grid |
212 |
CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
213 |
do bj = myByLo(myThid), myByHi(myThid) |
214 |
do bi = myBxLo(myThid), myBxHi(myThid) |
215 |
|
216 |
call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
217 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
218 |
C Note: correction term is now bottom up - needed in top down arrays |
219 |
do L = 1,Nrphys |
220 |
do j = 1,sNy |
221 |
do i = 1,sNx |
222 |
uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
223 |
enddo |
224 |
enddo |
225 |
enddo |
226 |
call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
227 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
228 |
do L = 1,Nrphys |
229 |
do j = 1,sNy |
230 |
do i = 1,sNx |
231 |
vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
232 |
enddo |
233 |
enddo |
234 |
enddo |
235 |
call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
236 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
237 |
do L = 1,Nrphys |
238 |
do j = 1,sNy |
239 |
do i = 1,sNx |
240 |
thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
241 |
enddo |
242 |
enddo |
243 |
enddo |
244 |
call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
245 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
246 |
do L = 1,Nrphys |
247 |
do j = 1,sNy |
248 |
do i = 1,sNx |
249 |
sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
250 |
enddo |
251 |
enddo |
252 |
enddo |
253 |
enddo |
254 |
enddo |
255 |
CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
256 |
|
257 |
c Last: Increment physics state by the correction term |
258 |
do bj = myByLo(myThid), myByHi(myThid) |
259 |
do bi = myBxLo(myThid), myBxHi(myThid) |
260 |
call step_physics(uphy,vphy,thphy,sphy,dtfake,im1,im2,jm1,jm2, |
261 |
. Nrphys,Nsx,Nsy,1,sNx,1,sNy,bi,bj, |
262 |
. uphytemp,vphytemp,thphytemp,sphytemp) |
263 |
|
264 |
call qcheck (im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,idim1,idim2, |
265 |
. jdim1,jdim2,bi,bj,dpphys,sphy) |
266 |
|
267 |
enddo |
268 |
enddo |
269 |
|
270 |
#ifdef ALLOW_DIAGNOSTICS |
271 |
call diagnostics_fill(uphy,'UAVE ',0,Nrphys,0,1,1,myThid) |
272 |
call diagnostics_fill(vphy,'VAVE ',0,Nrphys,0,1,1,myThid) |
273 |
call diagnostics_fill(thphy,'TAVE ',0,Nrphys,0,1,1,myThid) |
274 |
call diagnostics_fill(sphy,'QAVE ',0,Nrphys,0,1,1,myThid) |
275 |
#endif |
276 |
|
277 |
#ifdef ALLOW_DIAGNOSTICS |
278 |
do bj = myByLo(myThid), myByHi(myThid) |
279 |
do bi = myBxLo(myThid), myBxHi(myThid) |
280 |
do L=1,Nrphys |
281 |
|
282 |
c Total Tendency on Fizhi Grid for U (m/sec/day) |
283 |
c ----------------------------------------------- |
284 |
if(diagnostics_is_on('TENDUFIZ',myThid) ) then |
285 |
do j=jm1,jm2 |
286 |
do i=im1,im2 |
287 |
tempij(i,j) = (uphy (i,j,L,bi,bj)-ubef(i,j,L,bi,bj) ) |
288 |
. * 86400. _d 0 * dtinv |
289 |
enddo |
290 |
enddo |
291 |
call diagnostics_fill(tempij,'TENDUFIZ',L,1,3,bi,bj,myThid) |
292 |
endif |
293 |
|
294 |
c Total Tendency on Fizhi Grid for V (m/sec/day) |
295 |
c ----------------------------------------------- |
296 |
if(diagnostics_is_on('TENDVFIZ',myThid) ) then |
297 |
do j=jm1,jm2 |
298 |
do i=im1,im2 |
299 |
tempij(i,j) = (vphy (i,j,L,bi,bj)-vbef(i,j,L,bi,bj) ) |
300 |
. * 86400. _d 0 * dtinv |
301 |
enddo |
302 |
enddo |
303 |
call diagnostics_fill(tempij,'TENDVFIZ',L,1,3,bi,bj,myThid) |
304 |
endif |
305 |
|
306 |
c Total Tendency on Fizhi Grid for U (m/sec/day) |
307 |
c ----------------------------------------------- |
308 |
if(diagnostics_is_on('TENDTFIZ',myThid) ) then |
309 |
do j=jm1,jm2 |
310 |
do i=im1,im2 |
311 |
tempij(i,j) = (thphy (i,j,L,bi,bj)-thbef(i,j,L,bi,bj) ) |
312 |
. * 86400. _d 0 * dtinv |
313 |
enddo |
314 |
enddo |
315 |
call diagnostics_fill(tempij,'TENDTFIZ',L,1,3,bi,bj,myThid) |
316 |
endif |
317 |
|
318 |
c Total Tendency on Fizhi Grid for U (m/sec/day) |
319 |
c ----------------------------------------------- |
320 |
if(diagnostics_is_on('TENDQFIZ',myThid) ) then |
321 |
do j=jm1,jm2 |
322 |
do i=im1,im2 |
323 |
tempij(i,j) = (sphy (i,j,L,bi,bj)-sbef(i,j,L,bi,bj) ) |
324 |
. * 86400. _d 0 * dtinv |
325 |
enddo |
326 |
enddo |
327 |
call diagnostics_fill(tempij,'TENDQFIZ',L,1,3,bi,bj,myThid) |
328 |
endif |
329 |
|
330 |
enddo |
331 |
enddo |
332 |
enddo |
333 |
|
334 |
c Gridalt Correction Term Tendency for U and V (m/sec/day) |
335 |
c -------------------------------------------------------- |
336 |
if(diagnostics_is_on('CORRDU ',myThid) .or. |
337 |
. diagnostics_is_on('CORRDV ',myThid) ) then |
338 |
|
339 |
C gridalt correction term - first step is to compute adv+filters tendency |
340 |
C on dynamics grid (total - physics tend) |
341 |
do bj = myByLo(myThid), myByHi(myThid) |
342 |
do bi = myBxLo(myThid), myBxHi(myThid) |
343 |
do L=1,Nr |
344 |
do j=jm1,jm2 |
345 |
do i=im1,im2 |
346 |
udyntemp(i,j,L,bi,bj) = |
347 |
. (uvel(i,j,L,bi,bj)-udynbef(i,j,L,bi,bj))*dtinv - |
348 |
. guphy(i,j,L,bi,bj) |
349 |
vdyntemp(i,j,L,bi,bj) = |
350 |
. (vvel(i,j,L,bi,bj)-vdynbef(i,j,L,bi,bj))*dtinv - |
351 |
. gvphy(i,j,L,bi,bj) |
352 |
enddo |
353 |
enddo |
354 |
enddo |
355 |
C Next step - interpolate to fizhi grid |
356 |
C first put the u and v tendencies on an a-grid |
357 |
CALL TIMER_START('CTOA [STEP_FIZHI_CORR]',mythid) |
358 |
call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2, |
359 |
. Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp) |
360 |
CALL TIMER_STOP('CTOA [STEP_FIZHI_CORR]',mythid) |
361 |
C then do vertical interpolation |
362 |
do L = 1,Nrphys |
363 |
do j = 1,sNy |
364 |
do i = 1,sNx |
365 |
windphy(i,j,L,bi,bj) = |
366 |
. sqrt(uphy(i,j,Nrphys+1-L,bi,bj)*uphy(i,j,Nrphys+1-L,bi,bj) |
367 |
. + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj)) |
368 |
enddo |
369 |
enddo |
370 |
enddo |
371 |
CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
372 |
call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
373 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
374 |
CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
375 |
C Note: adv+filters term is now bottom up - needed in top down arrays |
376 |
do L = 1,Nrphys |
377 |
do j = 1,sNy |
378 |
do i = 1,sNx |
379 |
uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
380 |
enddo |
381 |
enddo |
382 |
enddo |
383 |
call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
384 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
385 |
do L = 1,Nrphys |
386 |
do j = 1,sNy |
387 |
do i = 1,sNx |
388 |
vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
389 |
enddo |
390 |
enddo |
391 |
enddo |
392 |
C Last Step - subtract adv+filters and physics tend from total tend on physics grid |
393 |
do L = 1,Nrphys |
394 |
do j = 1,sNy |
395 |
do i = 1,sNx |
396 |
uphytemp(i,j,L,bi,bj)= |
397 |
. ( (uphy(i,j,L,bi,bj)-ubef(i,j,L,bi,bj))*dtinv |
398 |
. - duphy(i,j,L,bi,bj) - uphytemp(i,j,L,bi,bj) ) * 86400. _d 0 |
399 |
vphytemp(i,j,L,bi,bj)= |
400 |
. ( (vphy(i,j,L,bi,bj)-vbef(i,j,L,bi,bj))*dtinv |
401 |
. - dvphy(i,j,L,bi,bj) - vphytemp(i,j,L,bi,bj) ) * 86400. _d 0 |
402 |
enddo |
403 |
enddo |
404 |
enddo |
405 |
enddo |
406 |
enddo |
407 |
|
408 |
if(diagnostics_is_on('CORRDU ',myThid)) then |
409 |
call diagnostics_fill(uphytemp,'CORRDU ',0,Nrphys,0,1,1,myThid) |
410 |
endif |
411 |
if(diagnostics_is_on('CORRDV ',myThid)) then |
412 |
call diagnostics_fill(vphytemp,'CORRDV ',0,Nrphys,0,1,1,myThid) |
413 |
endif |
414 |
|
415 |
endif |
416 |
|
417 |
c Gridalt Correction Term Tendency for TH (deg K/day) |
418 |
c -------------------------------------------------------- |
419 |
if(diagnostics_is_on('CORRDT ',myThid)) then |
420 |
|
421 |
C gridalt correction term - first step is to compute adv+filters tendency |
422 |
C on dynamics grid (total - physics tend) |
423 |
do bj = myByLo(myThid), myByHi(myThid) |
424 |
do bi = myBxLo(myThid), myBxHi(myThid) |
425 |
do L=1,Nr |
426 |
do j=jm1,jm2 |
427 |
do i=im1,im2 |
428 |
thdyntemp(i,j,L,bi,bj) = |
429 |
. (theta(i,j,L,bi,bj)-thdynbef(i,j,L,bi,bj))*dtinv - |
430 |
. gthphy(i,j,L,bi,bj) |
431 |
enddo |
432 |
enddo |
433 |
enddo |
434 |
C Next step - interpolate to fizhi grid |
435 |
CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
436 |
call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
437 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
438 |
CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
439 |
C Note: adv+filters term is now bottom up - needed in top down arrays |
440 |
do L = 1,Nrphys |
441 |
do j = 1,sNy |
442 |
do i = 1,sNx |
443 |
thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
444 |
enddo |
445 |
enddo |
446 |
enddo |
447 |
C Last Step - subtract adv+filters and physics tend from total tend on physics grid |
448 |
do L = 1,Nrphys |
449 |
do j = 1,sNy |
450 |
do i = 1,sNx |
451 |
thphytemp(i,j,L,bi,bj)= |
452 |
. ( (thphy(i,j,L,bi,bj)-thbef(i,j,L,bi,bj))*dtinv |
453 |
. - dthphy(i,j,L,bi,bj) - thphytemp(i,j,L,bi,bj) ) * 86400. _d 0 |
454 |
enddo |
455 |
enddo |
456 |
enddo |
457 |
enddo |
458 |
enddo |
459 |
|
460 |
call diagnostics_fill(thphytemp,'CORRDT ',0,Nrphys,0,1,1,myThid) |
461 |
endif |
462 |
|
463 |
c Gridalt Correction Term Tendency for Q (kg/kg/day) |
464 |
c -------------------------------------------------------- |
465 |
if(diagnostics_is_on('CORRDQ ',myThid)) then |
466 |
|
467 |
C gridalt correction term - first step is to compute adv+filters tendency |
468 |
C on dynamics grid (total - physics tend) |
469 |
do bj = myByLo(myThid), myByHi(myThid) |
470 |
do bi = myBxLo(myThid), myBxHi(myThid) |
471 |
do L=1,Nr |
472 |
do j=jm1,jm2 |
473 |
do i=im1,im2 |
474 |
sdyntemp(i,j,L,bi,bj) = |
475 |
. (salt(i,j,L,bi,bj)-sdynbef(i,j,L,bi,bj))*dtinv - |
476 |
. gsphy(i,j,L,bi,bj) |
477 |
enddo |
478 |
enddo |
479 |
enddo |
480 |
C Next step - interpolate to fizhi grid |
481 |
CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
482 |
call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
483 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
484 |
CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
485 |
C Note: adv+filters term is now bottom up - needed in top down arrays |
486 |
do L = 1,Nrphys |
487 |
do j = 1,sNy |
488 |
do i = 1,sNx |
489 |
sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
490 |
enddo |
491 |
enddo |
492 |
enddo |
493 |
C Last Step - subtract adv+filters and physics tend from total tend on physics grid |
494 |
do L = 1,Nrphys |
495 |
do j = 1,sNy |
496 |
do i = 1,sNx |
497 |
sphytemp(i,j,L,bi,bj)= |
498 |
. ( (sphy(i,j,L,bi,bj)-sbef(i,j,L,bi,bj))*dtinv |
499 |
. - dsphy(i,j,L,bi,bj) - sphytemp(i,j,L,bi,bj) ) * 86400. _d 0 |
500 |
enddo |
501 |
enddo |
502 |
enddo |
503 |
enddo |
504 |
enddo |
505 |
|
506 |
call diagnostics_fill(sphytemp,'CORRDQ ',0,Nrphys,0,1,1,myThid) |
507 |
endif |
508 |
#endif |
509 |
|
510 |
return |
511 |
end |