1 |
molod |
1.15 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/step_fizhi_corr.F,v 1.14 2006/02/16 23:03:53 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 |
|
|
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 |
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 |
molod |
1.15 |
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 |
|
|
. drF(L-1)*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 |
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 |
molod |
1.15 |
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 |
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.14 |
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 |
molod |
1.15 |
. * 86400. _d 0 * dtinv |
282 |
molod |
1.14 |
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 |
molod |
1.15 |
. * 86400. _d 0 * dtinv |
294 |
molod |
1.14 |
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 |
molod |
1.15 |
. * 86400. _d 0 * dtinv |
306 |
molod |
1.14 |
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 |
molod |
1.15 |
. * 86400. _d 0 * dtinv |
318 |
molod |
1.14 |
enddo |
319 |
|
|
enddo |
320 |
|
|
call diagnostics_fill(tempij,'TENDQFIZ',L,1,3,bi,bj,myThid) |
321 |
|
|
endif |
322 |
molod |
1.11 |
|
323 |
molod |
1.14 |
enddo |
324 |
|
|
enddo |
325 |
|
|
enddo |
326 |
molod |
1.1 |
|
327 |
molod |
1.15 |
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 |
449 |
|
|
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 |
molod |
1.14 |
return |
504 |
|
|
end |