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

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

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

revision 1.2 by molod, Mon Feb 23 20:34:38 2004 UTC revision 1.15 by molod, Fri Feb 24 00:12:15 2006 UTC
# Line 1  Line 1 
1         subroutine step_fizhi_corr (myTime, myIter, myThid)  C $Header$
2    C $Name$
3    
4    #include "FIZHI_OPTIONS.h"
5           subroutine step_fizhi_corr (myTime, myIter, myThid, dt)
6  c----------------------------------------------------------------------  c----------------------------------------------------------------------
7  c  Subroutine step_fizhi_corr - 'Wrapper' routine to advance  c  Subroutine step_fizhi_corr - 'Wrapper' routine to advance
8  c        the physics state and make the new value.  c        the physics state and make the new value.
# Line 15  c      dyn2phys (4) (interpolate A-Grid Line 19  c      dyn2phys (4) (interpolate A-Grid
19  c      step_physics (advance physics state by correction term)  c      step_physics (advance physics state by correction term)
20  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
21         implicit none         implicit none
 #include "CPP_OPTIONS.h"  
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "GRID.h"  #include "GRID.h"
24  #include "fizhi_SIZE.h"  #include "fizhi_SIZE.h"
25  #include "land_SIZE.h"  #include "fizhi_land_SIZE.h"
26  #include "DYNVARS.h"  #include "DYNVARS.h"
27  #include "fizhi_coms.h"  #include "fizhi_coms.h"
28  #include "gridalt_mapping.h"  #include "gridalt_mapping.h"
29  #include "EEPARAMS.h"  #include "EEPARAMS.h"
30  #include "SURFACE.h"  #include "SURFACE.h"
31    #ifdef ALLOW_DIAGNOSTICS
32         integer myTime, myIter, myThid  #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  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)         _RL pephy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy)
# Line 35  c pe on dynamics and physics grid refers Line 46  c pe on dynamics and physics grid refers
46         _RL udyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)         _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)         _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)         _RL thdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
49         _RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)         _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)         _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)         _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)         _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)         _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         integer i, j, L, Lbotij, bi, bj
57         integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2         integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
58         _RL dt         _RL dt
59    
60           _RL tempij(sNx,sNy)
61           _RL dtfake
62           _RL dtinv
63    
64         im1 = 1-OLx         im1 = 1-OLx
65         im2 = sNx+OLx         im2 = sNx+OLx
66         jm1 = 1-OLy         jm1 = 1-OLy
# Line 53  c pe on dynamics and physics grid refers Line 69  c pe on dynamics and physics grid refers
69         idim2 = sNx         idim2 = sNx
70         jdim1 = 1         jdim1 = 1
71         jdim2 = sNy         jdim2 = sNy
72         dt = 1.         dtfake = 1. _d 0
73           dtinv = 1. _d 0 / dt
74    
75         do bj = myByLo(myThid), myByHi(myThid)         do bj = myByLo(myThid), myByHi(myThid)
76         do bi = myBxLo(myThid), myBxHi(myThid)         do bi = myBxLo(myThid), myBxHi(myThid)
# Line 69  C Construct Pressures on physics and dyn Line 86  C Construct Pressures on physics and dyn
86          do j = 1,sNy          do j = 1,sNy
87          do i = 1,sNx          do i = 1,sNx
88           Lbotij = ksurfC(i,j,bi,bj)           Lbotij = ksurfC(i,j,bi,bj)
89           if(Lbotij.ne.0.)           if(Lbotij.lt.nr+1)
90       .    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)       .    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
91          enddo          enddo
92          enddo          enddo
# Line 99  c Do not use a zero field as the top edg Line 116  c Do not use a zero field as the top edg
116          enddo          enddo
117  c  c
118  c Create a wind magnitude field on the physics grid -  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          do L = 1,Nrphys
121          do j = 1,sNy          do j = 1,sNy
122          do i = 1,sNx          do i = 1,sNx
123           windphy(i,j,L,bi,bj) = sqrt(uphy(i,j,L,bi,bj)*uphy(i,j,L,bi,bj)           windphy(i,j,L,bi,bj) =
124       .        + vphy(i,j,L,bi,bj)*vphy(i,j,L,bi,bj))       .     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          enddo
127          enddo          enddo
128          enddo          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:  c Compute correction term (new dyn state-phys state to dyn) on physics grid:
137  c    First: interp physics state to dynamics grid  c    First: interp physics state to dynamics grid
138          call phys2dyn(uphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,  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)       .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,udyntemp)
148          call phys2dyn(vphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,          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)       .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,vdyntemp)
157          call phys2dyn(thphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,          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)       .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,thdyntemp)
166          call phys2dyn(sphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,          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)       .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,sdyntemp)
175    
176         enddo         enddo
177         enddo         enddo
178           CALL TIMER_STOP('PHYS2DYN          [STEP_FIZHI_CORR]',mythid)
179    
180  c    Second: Convert physics state on dynamics grid to C-Grid  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,         call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr,
184       .                     Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)       .                     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  c    Third: Subtract Phys state on dyn. grid from new dynamics state
188         do bj = myByLo(myThid), myByHi(myThid)         do bj = myByLo(myThid), myByHi(myThid)
# Line 146  c    Third: Subtract Phys state on dyn. Line 203  c    Third: Subtract Phys state on dyn.
203         enddo         enddo
204    
205  c    Fourth: Convert correction terms to A-Grid  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,          call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
208       .     Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)       .     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  c    Fifth: Interpolate correction terms to physics grid
212           CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
213         do bj = myByLo(myThid), myByHi(myThid)         do bj = myByLo(myThid), myByHi(myThid)
214         do bi = myBxLo(myThid), myBxHi(myThid)         do bi = myBxLo(myThid), myBxHi(myThid)
215    
216          call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,          call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
217       .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,uphytemp)       .      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,          call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
227       .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,vphytemp)       .      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,          call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
236       .     1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,thphytemp)       .     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,          call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
245       .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,sphytemp)       .      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  c    Last: Increment physics state by the correction term
258          call step_physics(uphy,vphy,thphy,sphy,dt,im1,im2,jm1,jm2,         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,       .   Nrphys,Nsx,Nsy,1,sNx,1,sNy,bi,bj,
262       .                            uphytemp,vphytemp,thphytemp,sphytemp)       .                            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          do bj = myByLo(myThid), myByHi(myThid)
272          do bi = myBxLo(myThid), myBxHi(myThid)
273          do L=1,Nrphys
274    
275    c Total Tendency on Fizhi Grid for U (m/sec/day)
276    c -----------------------------------------------
277           if(diagnostics_is_on('TENDUFIZ',myThid) ) then
278            do j=jm1,jm2
279            do i=im1,im2
280             tempij(i,j) = (uphy (i,j,L,bi,bj)-ubef(i,j,L,bi,bj) )
281         .                     * 86400. _d 0 * dtinv
282            enddo
283            enddo
284            call diagnostics_fill(tempij,'TENDUFIZ',L,1,3,bi,bj,myThid)
285           endif
286    
287    c Total Tendency on Fizhi Grid for V (m/sec/day)
288    c -----------------------------------------------
289           if(diagnostics_is_on('TENDVFIZ',myThid) ) then
290            do j=jm1,jm2
291            do i=im1,im2
292             tempij(i,j) = (vphy (i,j,L,bi,bj)-vbef(i,j,L,bi,bj) )
293         .                     * 86400. _d 0 * dtinv
294            enddo
295            enddo
296            call diagnostics_fill(tempij,'TENDVFIZ',L,1,3,bi,bj,myThid)
297           endif
298    
299    c Total Tendency on Fizhi Grid for U (m/sec/day)
300    c -----------------------------------------------
301           if(diagnostics_is_on('TENDTFIZ',myThid) ) then
302            do j=jm1,jm2
303            do i=im1,im2
304             tempij(i,j) = (thphy (i,j,L,bi,bj)-thbef(i,j,L,bi,bj) )
305         .                     * 86400. _d 0 * dtinv
306            enddo
307            enddo
308            call diagnostics_fill(tempij,'TENDTFIZ',L,1,3,bi,bj,myThid)
309           endif
310    
311    c Total Tendency on Fizhi Grid for U (m/sec/day)
312    c -----------------------------------------------
313           if(diagnostics_is_on('TENDQFIZ',myThid) ) then
314            do j=jm1,jm2
315            do i=im1,im2
316             tempij(i,j) = (sphy (i,j,L,bi,bj)-sbef(i,j,L,bi,bj) )
317         .                     * 86400. _d 0 * dtinv
318            enddo
319            enddo
320            call diagnostics_fill(tempij,'TENDQFIZ',L,1,3,bi,bj,myThid)
321           endif
322    
323          enddo
324          enddo
325          enddo
326    
327    c Gridalt Correction Term Tendency for U and V (m/sec/day)
328    c --------------------------------------------------------
329          if(diagnostics_is_on('CORRDU  ',myThid) .or.
330         .    diagnostics_is_on('CORRDV  ',myThid)   ) then
331    
332    C gridalt correction term - first step is to compute adv+filters tendency
333    C                           on dynamics grid (total - physics tend)
334          do bj = myByLo(myThid), myByHi(myThid)
335          do bi = myBxLo(myThid), myBxHi(myThid)
336           do L=1,Nr
337           do j=jm1,jm2
338           do i=im1,im2
339            udyntemp(i,j,L,bi,bj) =
340         .    (uvel(i,j,L,bi,bj)-udynbef(i,j,L,bi,bj))*dtinv -
341         .                                        guphy(i,j,L,bi,bj)
342            vdyntemp(i,j,L,bi,bj) =
343         .    (vvel(i,j,L,bi,bj)-vdynbef(i,j,L,bi,bj))*dtinv -
344         .                                        gvphy(i,j,L,bi,bj)
345           enddo
346           enddo
347           enddo
348    C Next step - interpolate to fizhi grid
349    C  first put the u and v tendencies on an a-grid
350           CALL TIMER_START('CTOA              [STEP_FIZHI_CORR]',mythid)
351           call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
352         .     Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
353           CALL TIMER_STOP('CTOA              [STEP_FIZHI_CORR]',mythid)
354    C  then do vertical interpolation
355           do L = 1,Nrphys
356           do j = 1,sNy
357           do i = 1,sNx
358            windphy(i,j,L,bi,bj) =
359         .     sqrt(uphy(i,j,Nrphys+1-L,bi,bj)*uphy(i,j,Nrphys+1-L,bi,bj)
360         .        + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj))
361           enddo
362           enddo
363           enddo
364           CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
365            call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
366         .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
367           CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
368    C Note: adv+filters term is now bottom up - needed in top down arrays
369           do L = 1,Nrphys
370           do j = 1,sNy
371           do i = 1,sNx
372             uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
373           enddo
374           enddo
375           enddo
376           call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
377         .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
378           do L = 1,Nrphys
379           do j = 1,sNy
380           do i = 1,sNx
381             vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
382           enddo
383           enddo
384           enddo
385    C Last Step - subtract adv+filters and physics tend from total tend on physics grid
386           do L = 1,Nrphys
387           do j = 1,sNy
388           do i = 1,sNx
389            uphytemp(i,j,L,bi,bj)=
390         .    ( (uphy(i,j,L,bi,bj)-ubef(i,j,L,bi,bj))*dtinv
391         .     - duphy(i,j,L,bi,bj) - uphytemp(i,j,L,bi,bj) ) * 86400. _d 0
392            vphytemp(i,j,L,bi,bj)=
393         .    ( (vphy(i,j,L,bi,bj)-vbef(i,j,L,bi,bj))*dtinv
394         .     - dvphy(i,j,L,bi,bj) - vphytemp(i,j,L,bi,bj) ) * 86400. _d 0
395           enddo
396           enddo
397           enddo
398          enddo
399          enddo
400    
401          if(diagnostics_is_on('CORRDU  ',myThid)) then
402           call diagnostics_fill(uphytemp,'CORRDU  ',0,Nrphys,0,1,1,myThid)
403          endif
404          if(diagnostics_is_on('CORRDV  ',myThid)) then
405           call diagnostics_fill(vphytemp,'CORRDV  ',0,Nrphys,0,1,1,myThid)
406          endif
407    
408          endif
409    
410    c Gridalt Correction Term Tendency for TH (deg K/day)
411    c --------------------------------------------------------
412          if(diagnostics_is_on('CORRDT  ',myThid))  then
413    
414    C gridalt correction term - first step is to compute adv+filters tendency
415    C                           on dynamics grid (total - physics tend)
416          do bj = myByLo(myThid), myByHi(myThid)
417          do bi = myBxLo(myThid), myBxHi(myThid)
418           do L=1,Nr
419           do j=jm1,jm2
420           do i=im1,im2
421            thdyntemp(i,j,L,bi,bj) =
422         .    (theta(i,j,L,bi,bj)-thdynbef(i,j,L,bi,bj))*dtinv -
423         .                                        gthphy(i,j,L,bi,bj)
424           enddo
425           enddo
426           enddo
427    C Next step - interpolate to fizhi grid
428           CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
429           call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
430         .     1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
431           CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
432    C Note: adv+filters term is now bottom up - needed in top down arrays
433           do L = 1,Nrphys
434           do j = 1,sNy
435           do i = 1,sNx
436             thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
437           enddo
438           enddo
439           enddo
440    C Last Step - subtract adv+filters and physics tend from total tend on physics grid
441           do L = 1,Nrphys
442           do j = 1,sNy
443           do i = 1,sNx
444            thphytemp(i,j,L,bi,bj)=
445         .    ( (thphy(i,j,L,bi,bj)-thbef(i,j,L,bi,bj))*dtinv
446         .     - dthphy(i,j,L,bi,bj) - thphytemp(i,j,L,bi,bj) ) * 86400. _d 0
447           enddo
448         enddo         enddo
449         enddo         enddo
450          enddo
451          enddo
452    
453          call diagnostics_fill(thphytemp,'CORRDT  ',0,Nrphys,0,1,1,myThid)
454          endif
455    
456    c Gridalt Correction Term Tendency for Q (kg/kg/day)
457    c --------------------------------------------------------
458          if(diagnostics_is_on('CORRDQ  ',myThid))  then
459    
460    C gridalt correction term - first step is to compute adv+filters tendency
461    C                           on dynamics grid (total - physics tend)
462          do bj = myByLo(myThid), myByHi(myThid)
463          do bi = myBxLo(myThid), myBxHi(myThid)
464           do L=1,Nr
465           do j=jm1,jm2
466           do i=im1,im2
467            sdyntemp(i,j,L,bi,bj) =
468         .    (salt(i,j,L,bi,bj)-sdynbef(i,j,L,bi,bj))*dtinv -
469         .                                        gsphy(i,j,L,bi,bj)
470           enddo
471           enddo
472           enddo
473    C Next step - interpolate to fizhi grid
474           CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
475           call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
476         .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
477           CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
478    C Note: adv+filters term is now bottom up - needed in top down arrays
479           do L = 1,Nrphys
480           do j = 1,sNy
481           do i = 1,sNx
482             sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
483           enddo
484           enddo
485           enddo
486    C Last Step - subtract adv+filters and physics tend from total tend on physics grid
487           do L = 1,Nrphys
488           do j = 1,sNy
489           do i = 1,sNx
490            sphytemp(i,j,L,bi,bj)=
491         .    ( (sphy(i,j,L,bi,bj)-sbef(i,j,L,bi,bj))*dtinv
492         .     - dsphy(i,j,L,bi,bj) - sphytemp(i,j,L,bi,bj) ) * 86400. _d 0
493           enddo
494           enddo
495           enddo
496          enddo
497          enddo
498    
499          call diagnostics_fill(sphytemp,'CORRDQ  ',0,Nrphys,0,1,1,myThid)
500          endif
501    #endif
502    
503         return        return
504         end        end
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22