/[MITgcm]/MITgcm/pkg/fizhi/step_fizhi_corr.F
ViewVC logotype

Annotation of /MITgcm/pkg/fizhi/step_fizhi_corr.F

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


Revision 1.19 - (hide annotations) (download)
Fri Jun 24 01:20:49 2011 UTC (13 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint63, 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, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.18: +19 -18 lines
avoid line longer than 72c

1 jmc 1.19 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/step_fizhi_corr.F,v 1.18 2006/08/04 17:03:58 molod Exp $
2 edhill 1.5 C $Name: $
3    
4 molod 1.8 #include "FIZHI_OPTIONS.h"
5 molod 1.14 subroutine step_fizhi_corr (myTime, myIter, myThid, dt)
6 molod 1.1 c----------------------------------------------------------------------
7     c Subroutine step_fizhi_corr - 'Wrapper' routine to advance
8 jmc 1.19 c the physics state and make the new value.
9 molod 1.1 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 jmc 1.19 c Call:phys2dyn (4) (interpolate physics state to dynamics grid
15 molod 1.1 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 molod 1.6 #include "fizhi_land_SIZE.h"
26 molod 1.1 #include "DYNVARS.h"
27     #include "fizhi_coms.h"
28     #include "gridalt_mapping.h"
29     #include "EEPARAMS.h"
30     #include "SURFACE.h"
31 molod 1.14 #ifdef ALLOW_DIAGNOSTICS
32     #include "fizhi_SHP.h"
33     #endif
34 molod 1.1
35 molod 1.13 integer myIter, myThid
36     _RL myTime
37 molod 1.14 #ifdef ALLOW_DIAGNOSTICS
38     logical diagnostics_is_on
39     external diagnostics_is_on
40     #endif
41 molod 1.1
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 molod 1.12 _RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
50 molod 1.1 _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 molod 1.7 _RL tempphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
55 molod 1.1
56     integer i, j, L, Lbotij, bi, bj
57     integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
58 molod 1.2 _RL dt
59 molod 1.1
60 molod 1.9 _RL tempij(sNx,sNy)
61 molod 1.14 _RL dtfake
62 molod 1.15 _RL dtinv
63 molod 1.9
64 molod 1.1 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 molod 1.15 dtfake = 1. _d 0
73     dtinv = 1. _d 0 / dt
74 molod 1.1
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 jmc 1.19 if(Lbotij.lt.nr+1)
90 molod 1.1 . 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 molod 1.18 . drF(L-1)* rStarExpC(i,j,bi,bj)*hfacC(i,j,L-1,bi,bj)
99 molod 1.1 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 jmc 1.19 c Create a wind magnitude field on the physics grid -
119 molod 1.15 c (Load the wind speed bottom up for use by dyn2phys)
120 molod 1.1 do L = 1,Nrphys
121     do j = 1,sNy
122     do i = 1,sNx
123 jmc 1.19 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 molod 1.15 . + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj))
126 molod 1.1 enddo
127     enddo
128     enddo
129 molod 1.3 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 molod 1.1
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 molod 1.7 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 molod 1.1 . 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,udyntemp)
148 molod 1.7 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 molod 1.1 . 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,vdyntemp)
157 molod 1.7 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 molod 1.1 . 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,thdyntemp)
166 molod 1.7 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 molod 1.1 . 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,sdyntemp)
175    
176     enddo
177     enddo
178 molod 1.3 CALL TIMER_STOP('PHYS2DYN [STEP_FIZHI_CORR]',mythid)
179 molod 1.1
180     c Second: Convert physics state on dynamics grid to C-Grid
181    
182 molod 1.3 CALL TIMER_START('ATOC [STEP_FIZHI_CORR]',mythid)
183 molod 1.1 call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr,
184     . Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
185 molod 1.3 CALL TIMER_STOP('ATOC [STEP_FIZHI_CORR]',mythid)
186 molod 1.1
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 molod 1.3 CALL TIMER_START('CTOA [STEP_FIZHI_CORR]',mythid)
207 molod 1.1 call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
208     . Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
209 molod 1.3 CALL TIMER_STOP('CTOA [STEP_FIZHI_CORR]',mythid)
210 molod 1.1
211     c Fifth: Interpolate correction terms to physics grid
212 molod 1.3 CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid)
213 molod 1.1 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 molod 1.7 . 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 molod 1.1 call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
227 molod 1.7 . 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 molod 1.1 call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
236 molod 1.7 . 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 molod 1.1 call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
245 molod 1.7 . 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 molod 1.3 enddo
254     enddo
255     CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid)
256 molod 1.2
257 molod 1.1 c Last: Increment physics state by the correction term
258 molod 1.3 do bj = myByLo(myThid), myByHi(myThid)
259     do bi = myBxLo(myThid), myBxHi(myThid)
260 molod 1.14 call step_physics(uphy,vphy,thphy,sphy,dtfake,im1,im2,jm1,jm2,
261 molod 1.2 . Nrphys,Nsx,Nsy,1,sNx,1,sNy,bi,bj,
262     . uphytemp,vphytemp,thphytemp,sphytemp)
263 molod 1.9
264 molod 1.11 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 molod 1.15 #ifdef ALLOW_DIAGNOSTICS
271 molod 1.16 call diagnostics_fill(uphy,'UAVE ',0,Nrphys,0,1,1,myThid)
272 molod 1.17 call diagnostics_fill(vphy,'VAVE ',0,Nrphys,0,1,1,myThid)
273 molod 1.16 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 molod 1.14 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 molod 1.15 . * 86400. _d 0 * dtinv
289 molod 1.14 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 molod 1.15 . * 86400. _d 0 * dtinv
301 molod 1.14 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 molod 1.15 . * 86400. _d 0 * dtinv
313 molod 1.14 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 molod 1.15 . * 86400. _d 0 * dtinv
325 molod 1.14 enddo
326     enddo
327     call diagnostics_fill(tempij,'TENDQFIZ',L,1,3,bi,bj,myThid)
328     endif
329 molod 1.11
330 molod 1.14 enddo
331     enddo
332     enddo
333 molod 1.1
334 molod 1.15 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 jmc 1.19 udyntemp(i,j,L,bi,bj) =
347     . (uvel(i,j,L,bi,bj)-udynbef(i,j,L,bi,bj))*dtinv -
348 molod 1.15 . guphy(i,j,L,bi,bj)
349 jmc 1.19 vdyntemp(i,j,L,bi,bj) =
350     . (vvel(i,j,L,bi,bj)-vdynbef(i,j,L,bi,bj))*dtinv -
351 molod 1.15 . 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 jmc 1.19 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 molod 1.15 . + 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 jmc 1.19 thdyntemp(i,j,L,bi,bj) =
429     . (theta(i,j,L,bi,bj)-thdynbef(i,j,L,bi,bj))*dtinv -
430 molod 1.15 . 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 jmc 1.19 . - dthphy(i,j,L,bi,bj) - thphytemp(i,j,L,bi,bj)
454     . ) * 86400. _d 0
455 molod 1.15 enddo
456     enddo
457     enddo
458     enddo
459     enddo
460    
461     call diagnostics_fill(thphytemp,'CORRDT ',0,Nrphys,0,1,1,myThid)
462     endif
463    
464     c Gridalt Correction Term Tendency for Q (kg/kg/day)
465     c --------------------------------------------------------
466     if(diagnostics_is_on('CORRDQ ',myThid)) then
467    
468     C gridalt correction term - first step is to compute adv+filters tendency
469     C on dynamics grid (total - physics tend)
470     do bj = myByLo(myThid), myByHi(myThid)
471     do bi = myBxLo(myThid), myBxHi(myThid)
472     do L=1,Nr
473     do j=jm1,jm2
474     do i=im1,im2
475 jmc 1.19 sdyntemp(i,j,L,bi,bj) =
476     . (salt(i,j,L,bi,bj)-sdynbef(i,j,L,bi,bj))*dtinv -
477 molod 1.15 . gsphy(i,j,L,bi,bj)
478     enddo
479     enddo
480     enddo
481     C Next step - interpolate to fizhi grid
482     CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid)
483     call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
484     . 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
485     CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid)
486     C Note: adv+filters term is now bottom up - needed in top down arrays
487     do L = 1,Nrphys
488     do j = 1,sNy
489     do i = 1,sNx
490     sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
491     enddo
492     enddo
493     enddo
494     C Last Step - subtract adv+filters and physics tend from total tend on physics grid
495     do L = 1,Nrphys
496     do j = 1,sNy
497     do i = 1,sNx
498     sphytemp(i,j,L,bi,bj)=
499     . ( (sphy(i,j,L,bi,bj)-sbef(i,j,L,bi,bj))*dtinv
500     . - dsphy(i,j,L,bi,bj) - sphytemp(i,j,L,bi,bj) ) * 86400. _d 0
501     enddo
502     enddo
503     enddo
504     enddo
505     enddo
506    
507     call diagnostics_fill(sphytemp,'CORRDQ ',0,Nrphys,0,1,1,myThid)
508     endif
509     #endif
510    
511 molod 1.14 return
512     end

  ViewVC Help
Powered by ViewVC 1.1.22