1 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/step_fizhi_corr.F,v 1.10 2004/08/04 22:57:35 molod Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "FIZHI_OPTIONS.h" |
5 |
subroutine step_fizhi_corr (myTime, myIter, myThid) |
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 |
|
32 |
integer myTime, myIter, myThid |
33 |
|
34 |
c pe on dynamics and physics grid refers to bottom edge |
35 |
_RL pephy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy) |
36 |
_RL pedyn(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr+1,nSx,nSy) |
37 |
_RL windphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
38 |
_RL udyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
39 |
_RL vdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
40 |
_RL thdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) |
41 |
_RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
42 |
_RL uphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
43 |
_RL vphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
44 |
_RL thphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
45 |
_RL sphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
46 |
_RL tempphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy) |
47 |
|
48 |
integer i, j, L, Lbotij, bi, bj |
49 |
integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2 |
50 |
_RL dt |
51 |
|
52 |
_RL tempij(sNx,sNy) |
53 |
|
54 |
im1 = 1-OLx |
55 |
im2 = sNx+OLx |
56 |
jm1 = 1-OLy |
57 |
jm2 = sNy+OLy |
58 |
idim1 = 1 |
59 |
idim2 = sNx |
60 |
jdim1 = 1 |
61 |
jdim2 = sNy |
62 |
dt = 1. |
63 |
|
64 |
do bj = myByLo(myThid), myByHi(myThid) |
65 |
do bi = myBxLo(myThid), myBxHi(myThid) |
66 |
|
67 |
C Construct Pressures on physics and dynamics grids |
68 |
do j = 1,sNy |
69 |
do i = 1,sNx |
70 |
do L = 1,Nr |
71 |
pedyn(i,j,L,bi,bj) = 0. |
72 |
enddo |
73 |
enddo |
74 |
enddo |
75 |
do j = 1,sNy |
76 |
do i = 1,sNx |
77 |
Lbotij = ksurfC(i,j,bi,bj) |
78 |
if(Lbotij.ne.0.) |
79 |
. pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj) |
80 |
enddo |
81 |
enddo |
82 |
do j = 1,sNy |
83 |
do i = 1,sNx |
84 |
Lbotij = ksurfC(i,j,bi,bj) |
85 |
do L = Lbotij+1,Nr+1 |
86 |
pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) - |
87 |
. drF(L-1)*hfacC(i,j,L-1,bi,bj) |
88 |
enddo |
89 |
c Do not use a zero field as the top edge pressure for interpolation |
90 |
if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5) |
91 |
. pedyn(i,j,Nr+1,bi,bj) = 1.e-5 |
92 |
enddo |
93 |
enddo |
94 |
|
95 |
do j = 1,sNy |
96 |
do i = 1,sNx |
97 |
pephy(i,j,1,bi,bj)=Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj) |
98 |
do L = 2,Nrphys+1 |
99 |
pephy(i,j,L,bi,bj)=pephy(i,j,L-1,bi,bj)-dpphys(i,j,L-1,bi,bj) |
100 |
enddo |
101 |
c Do not use a zero field as the top edge pressure for interpolation |
102 |
if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5) |
103 |
. pephy(i,j,Nrphys+1,bi,bj) = 1.e-5 |
104 |
enddo |
105 |
enddo |
106 |
c |
107 |
c Create a wind magnitude field on the physics grid - |
108 |
do L = 1,Nrphys |
109 |
do j = 1,sNy |
110 |
do i = 1,sNx |
111 |
windphy(i,j,L,bi,bj) = sqrt(uphy(i,j,L,bi,bj)*uphy(i,j,L,bi,bj) |
112 |
. + vphy(i,j,L,bi,bj)*vphy(i,j,L,bi,bj)) |
113 |
enddo |
114 |
enddo |
115 |
enddo |
116 |
enddo |
117 |
enddo |
118 |
|
119 |
CALL TIMER_START('PHYS2DYN [STEP_FIZHI_CORR]',mythid) |
120 |
do bj = myByLo(myThid), myByHi(myThid) |
121 |
do bi = myBxLo(myThid), myBxHi(myThid) |
122 |
|
123 |
c Compute correction term (new dyn state-phys state to dyn) on physics grid: |
124 |
c First: interp physics state to dynamics grid |
125 |
C Note: physics field levels are numbered top down - need bottom up |
126 |
do L = 1,Nrphys |
127 |
do j = 1,sNy |
128 |
do i = 1,sNx |
129 |
tempphy(i,j,Nrphys+1-L,bi,bj) = uphy(i,j,L,bi,bj) |
130 |
enddo |
131 |
enddo |
132 |
enddo |
133 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
134 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,udyntemp) |
135 |
do L = 1,Nrphys |
136 |
do j = 1,sNy |
137 |
do i = 1,sNx |
138 |
tempphy(i,j,Nrphys+1-L,bi,bj) = vphy(i,j,L,bi,bj) |
139 |
enddo |
140 |
enddo |
141 |
enddo |
142 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
143 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,vdyntemp) |
144 |
do L = 1,Nrphys |
145 |
do j = 1,sNy |
146 |
do i = 1,sNx |
147 |
tempphy(i,j,Nrphys+1-L,bi,bj) = thphy(i,j,L,bi,bj) |
148 |
enddo |
149 |
enddo |
150 |
enddo |
151 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
152 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,thdyntemp) |
153 |
do L = 1,Nrphys |
154 |
do j = 1,sNy |
155 |
do i = 1,sNx |
156 |
tempphy(i,j,Nrphys+1-L,bi,bj) = sphy(i,j,L,bi,bj) |
157 |
enddo |
158 |
enddo |
159 |
enddo |
160 |
call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy, |
161 |
. 1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,sdyntemp) |
162 |
|
163 |
enddo |
164 |
enddo |
165 |
CALL TIMER_STOP('PHYS2DYN [STEP_FIZHI_CORR]',mythid) |
166 |
|
167 |
c Second: Convert physics state on dynamics grid to C-Grid |
168 |
|
169 |
CALL TIMER_START('ATOC [STEP_FIZHI_CORR]',mythid) |
170 |
call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr, |
171 |
. Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp) |
172 |
CALL TIMER_STOP('ATOC [STEP_FIZHI_CORR]',mythid) |
173 |
|
174 |
c Third: Subtract Phys state on dyn. grid from new dynamics state |
175 |
do bj = myByLo(myThid), myByHi(myThid) |
176 |
do bi = myBxLo(myThid), myBxHi(myThid) |
177 |
|
178 |
do L = 1,Nr |
179 |
do j = jdim1,jdim2 |
180 |
do i = idim1,idim2 |
181 |
udyntemp(i,j,L,bi,bj)=uvel(i,j,L,bi,bj)-udyntemp(i,j,L,bi,bj) |
182 |
vdyntemp(i,j,L,bi,bj)=vvel(i,j,L,bi,bj)-vdyntemp(i,j,L,bi,bj) |
183 |
thdyntemp(i,j,L,bi,bj)=theta(i,j,L,bi,bj)-thdyntemp(i,j,L,bi,bj) |
184 |
sdyntemp(i,j,L,bi,bj)=salt(i,j,L,bi,bj)-sdyntemp(i,j,L,bi,bj) |
185 |
enddo |
186 |
enddo |
187 |
enddo |
188 |
|
189 |
enddo |
190 |
enddo |
191 |
|
192 |
c Fourth: Convert correction terms to A-Grid |
193 |
CALL TIMER_START('CTOA [STEP_FIZHI_CORR]',mythid) |
194 |
call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2, |
195 |
. Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp) |
196 |
CALL TIMER_STOP('CTOA [STEP_FIZHI_CORR]',mythid) |
197 |
|
198 |
c Fifth: Interpolate correction terms to physics grid |
199 |
CALL TIMER_START('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
200 |
do bj = myByLo(myThid), myByHi(myThid) |
201 |
do bi = myBxLo(myThid), myBxHi(myThid) |
202 |
|
203 |
call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
204 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
205 |
C Note: correction term is now bottom up - needed in top down arrays |
206 |
do L = 1,Nrphys |
207 |
do j = 1,sNy |
208 |
do i = 1,sNx |
209 |
uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
210 |
enddo |
211 |
enddo |
212 |
enddo |
213 |
call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
214 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy) |
215 |
do L = 1,Nrphys |
216 |
do j = 1,sNy |
217 |
do i = 1,sNx |
218 |
vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
219 |
enddo |
220 |
enddo |
221 |
enddo |
222 |
call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
223 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
224 |
do L = 1,Nrphys |
225 |
do j = 1,sNy |
226 |
do i = 1,sNx |
227 |
thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
228 |
enddo |
229 |
enddo |
230 |
enddo |
231 |
call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx, |
232 |
. 1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy) |
233 |
do L = 1,Nrphys |
234 |
do j = 1,sNy |
235 |
do i = 1,sNx |
236 |
sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj) |
237 |
enddo |
238 |
enddo |
239 |
enddo |
240 |
enddo |
241 |
enddo |
242 |
CALL TIMER_STOP('DYN2PHYS [STEP_FIZHI_CORR]',mythid) |
243 |
|
244 |
c Last: Increment physics state by the correction term |
245 |
do bj = myByLo(myThid), myByHi(myThid) |
246 |
do bi = myBxLo(myThid), myBxHi(myThid) |
247 |
call step_physics(uphy,vphy,thphy,sphy,dt,im1,im2,jm1,jm2, |
248 |
. Nrphys,Nsx,Nsy,1,sNx,1,sNy,bi,bj, |
249 |
. uphytemp,vphytemp,thphytemp,sphytemp) |
250 |
|
251 |
if(2.eq.1 )then |
252 |
print *,' In step fizhi corr, new fizhi fields ',bi,' dt= ',dt |
253 |
do L = 1,Nrphys |
254 |
do j = jdim1,jdim2 |
255 |
do i = idim1,idim2 |
256 |
tempij(i,j) = uphy(i,j,L,bi,bj) |
257 |
enddo |
258 |
enddo |
259 |
c print *,' uphy at level ',l,' ',tempij |
260 |
enddo |
261 |
do L = 1,Nrphys |
262 |
do j = jdim1,jdim2 |
263 |
do i = idim1,idim2 |
264 |
tempij(i,j) = vphy(i,j,L,bi,bj) |
265 |
enddo |
266 |
enddo |
267 |
c print *,' vphy at level ',l,' ',tempij |
268 |
enddo |
269 |
do L = 1,Nrphys |
270 |
do j = jdim1,jdim2 |
271 |
do i = idim1,idim2 |
272 |
tempij(i,j) = thphy(i,j,L,bi,bj) |
273 |
enddo |
274 |
enddo |
275 |
print *,' thphy at level ',l,' ',tempij |
276 |
enddo |
277 |
do L = 1,Nrphys |
278 |
do j = jdim1,jdim2 |
279 |
do i = idim1,idim2 |
280 |
tempij(i,j) = sphy(i,j,L,bi,bj) |
281 |
enddo |
282 |
enddo |
283 |
print *,' sphy at level ',l,' ',tempij |
284 |
enddo |
285 |
endif |
286 |
|
287 |
call qcheck (im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,idim1,idim2, |
288 |
. jdim1,jdim2,bi,bj,dpphys,sphy) |
289 |
|
290 |
if(2.eq.1 )then |
291 |
print *,' In step fizhi corr after qcheck ',bi |
292 |
do L = 1,Nrphys |
293 |
do j = jdim1,jdim2 |
294 |
do i = idim1,idim2 |
295 |
tempij(i,j) = sphy(i,j,L,bi,bj) |
296 |
enddo |
297 |
enddo |
298 |
print *,' sphy after qcheck at level ',l,' ',tempij |
299 |
enddo |
300 |
endif |
301 |
|
302 |
|
303 |
enddo |
304 |
enddo |
305 |
|
306 |
return |
307 |
end |