1 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_fillnegs.F,v 1.6 2004/08/04 22:50:41 molod Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "FIZHI_OPTIONS.h" |
5 |
SUBROUTINE QCHECK ( idim1,idim2,jdim1,jdim2,ldim,im1,im2,jm1,jm2, |
6 |
. dp,pz,qz) |
7 |
C*********************************************************************** |
8 |
C Purpose |
9 |
C Check Specific Humidity Field for Negative values |
10 |
C |
11 |
C Argument Description |
12 |
C IDIM1 .... Start Zonal Dimension |
13 |
C IDIM2 .... End Zonal Dimension |
14 |
C JDIM1 .... Start Meridional Dimension |
15 |
C JDIM2 .... End Meridional Dimension |
16 |
C IM1 ...... Start Zonal Span |
17 |
C IM2 ...... End Zonal Span |
18 |
C JM1 ...... Start Meridional Span |
19 |
C JM2 ...... End Meridional Span |
20 |
C LDIM ..... Vertical Dimension |
21 |
C DP ....... Delta Pressure |
22 |
C PZ ........Surface Pressure |
23 |
C QZ ........Specific Humidity (g/g) |
24 |
C*********************************************************************** |
25 |
implicit none |
26 |
|
27 |
integer idim1,idim2,jdim1,jdim2,Ldim,im1,im2,jm1,jm2 |
28 |
_RL delt |
29 |
_RL qz(idim1:idim2,jdim1:jdim2,ldim) |
30 |
_RL dp(idim1:idim2,jdim1:jdim2,ldim) |
31 |
_RL pz(idim1:idim2,jdim1:jdim2) |
32 |
|
33 |
integer i,j,L,LM1 |
34 |
_RL getcon,grav,ddsig |
35 |
|
36 |
grav = getcon('GRAVITY') |
37 |
|
38 |
c Fill Negative Specific Humidities |
39 |
c --------------------------------- |
40 |
DO L=2,Ldim |
41 |
LM1 = L-1 |
42 |
do j=jm1,jm2 |
43 |
do i=im1,im2 |
44 |
ddsig = dp(i,j,LM1)/dp(i,j,L) |
45 |
if( qz(i,j,LM1).lt.0.0 ) then |
46 |
qz(i,j,L ) = qz(i,j,L) + qz(i,j,LM1)*ddsig |
47 |
qz(i,j,LM1) = 0.0 |
48 |
endif |
49 |
enddo |
50 |
enddo |
51 |
enddo |
52 |
|
53 |
do j=jm1,jm2 |
54 |
do i=im1,im2 |
55 |
if( qz(i,j,Ldim).lt.0.0 ) qz(i,j,Ldim) = 0.0 |
56 |
enddo |
57 |
enddo |
58 |
|
59 |
return |
60 |
end |
61 |
|
62 |
subroutine tracer_fill ( pq,im,jm,lm,dlam,dphi,dp) |
63 |
C*********************************************************************** |
64 |
C PURPOSE |
65 |
C Fill negative tracer values using local borrowing |
66 |
C |
67 |
C INPUT |
68 |
C pq ..... Mass-weighted (PI) Tracer |
69 |
C im ..... Zonal Dimension |
70 |
C jm ..... Meridional Dimension |
71 |
C lm ..... Vertical Dimension |
72 |
C dlam ... Zonal Grid Increment |
73 |
C dphi ... Meridional Grid Increment |
74 |
C dp ..... Vertical Grid Increment |
75 |
C |
76 |
C Note: |
77 |
C If no immediate surrounding value is large enough to fill negative |
78 |
C value, |
79 |
C the sum of immediate surrounding positive values is tried. |
80 |
C If sum is not large enough, tracer is simply set to zero. |
81 |
C |
82 |
C*********************************************************************** |
83 |
|
84 |
implicit none |
85 |
|
86 |
c Input Variables |
87 |
c --------------- |
88 |
integer im,jm,lm |
89 |
_RL pq(im,jm,lm),dlam(im),dphi(jm),dp(im,jm,lm) |
90 |
|
91 |
c Local Variables |
92 |
c --------------- |
93 |
integer i,j,l,im1,ip1,imax,m |
94 |
_RL lam(im), phi(jm) |
95 |
_RL array(6) |
96 |
_RL pi,a,getcon,undef |
97 |
_RL qmax,qval,sum,fact |
98 |
|
99 |
_RL dxu(im,jm) |
100 |
_RL dxv(im,jm) |
101 |
_RL dxp(im,jm) |
102 |
_RL dyv(im,jm) |
103 |
_RL dyp(im,jm) |
104 |
|
105 |
_RL d2p(im,jm) |
106 |
|
107 |
C ********************************************************* |
108 |
C **** Initialization **** |
109 |
C ********************************************************* |
110 |
|
111 |
pi = 4.0*atan(1.0) |
112 |
a = getcon('EARTH RADIUS') |
113 |
|
114 |
c Compute Longitudes |
115 |
c ------------------ |
116 |
lam(1) = -pi |
117 |
do i=2,im |
118 |
lam(i) = lam(i-1) + dlam(i-1) |
119 |
enddo |
120 |
|
121 |
c Compute Latitudes |
122 |
c ----------------- |
123 |
phi(1) = -pi/2. |
124 |
do j=2,jm-1 |
125 |
phi(j) = phi(j-1) + dphi(j-1) |
126 |
enddo |
127 |
phi(jm) = pi/2. |
128 |
|
129 |
c Compute DXU and DYV |
130 |
c ------------------- |
131 |
do j=2,jm-1 |
132 |
do i=1,im |
133 |
dxu(i,j) = a*cos(phi(j))*dlam(i) |
134 |
enddo |
135 |
enddo |
136 |
|
137 |
do j=2,jm-2 |
138 |
do i=1,im |
139 |
dyv(i,j) = a*dphi(j) |
140 |
enddo |
141 |
enddo |
142 |
do i=1,im |
143 |
dyv(i,1) = a*(dphi(1) +0.5*dphi(2) ) |
144 |
dyv(i,jm-1) = a*(dphi(jm-1)+0.5*dphi(jm-2)) |
145 |
enddo |
146 |
|
147 |
c Compute DXP and DXV |
148 |
c ------------------- |
149 |
do j=2,jm-1 |
150 |
im1 = im |
151 |
do i=1,im |
152 |
dxp(i,j) = ( dxu(i,j)+dxu(im1,j) )*0.5 |
153 |
im1 = i |
154 |
enddo |
155 |
enddo |
156 |
|
157 |
do j=2,jm-2 |
158 |
do i=1,im |
159 |
dxv(i,j) = ( dxp(i,j)+dxp(i,j+1) )*0.5 |
160 |
enddo |
161 |
enddo |
162 |
|
163 |
c Compute DYP |
164 |
c ----------- |
165 |
do j=3,jm-2 |
166 |
do i=1,im |
167 |
dyp(i,j) = ( dyv(i,j)+dyv(i,j-1) )*0.5 |
168 |
enddo |
169 |
enddo |
170 |
do i=1,im |
171 |
dyp(i,2) = dyv(i,1) |
172 |
dyp(i,jm-1) = dyv(i,jm-1) |
173 |
enddo |
174 |
|
175 |
c Compute Area Factor D2P |
176 |
c ----------------------- |
177 |
do j=3,jm-2 |
178 |
do i=1,im |
179 |
d2p(i,j) = 0.5*( dxv(i,j)+dxv(i,j-1) )*dyp(i,j) |
180 |
enddo |
181 |
enddo |
182 |
do i=1,im |
183 |
d2p(i,2) = dxv(i,2) *dyp(i,2) |
184 |
d2p(i,jm-1) = dxv(i,jm-2)*dyp(i,jm-1) |
185 |
enddo |
186 |
|
187 |
undef = getcon('UNDEF') |
188 |
|
189 |
C ********************************************************* |
190 |
C **** Fill Negative Values **** |
191 |
C ********************************************************* |
192 |
|
193 |
do l=1,lm |
194 |
do j=2,jm-1 |
195 |
|
196 |
im1 = im-1 |
197 |
i = im |
198 |
do ip1=1,im |
199 |
|
200 |
if( pq(i,j,L).lt.0.0 ) then |
201 |
|
202 |
qval = pq(i ,j,L)*d2p(i ,j)*dp(i,j,L) |
203 |
array(1) = pq(ip1,j,L)*d2p(ip1,j)*dp(i,j,L) |
204 |
array(2) = pq(im1,j,L)*d2p(im1,j)*dp(i,j,L) |
205 |
|
206 |
if( j.eq.jm-1 ) then |
207 |
array(3) = -undef |
208 |
else |
209 |
array(3) = pq(i,j+1,L)*d2p(i,j+1)*dp(i,j,L) |
210 |
endif |
211 |
if( j.eq.2 ) then |
212 |
array(4) = -undef |
213 |
else |
214 |
array(4) = pq(i,j-1,L)*d2p(i,j-1)*dp(i,j,L) |
215 |
endif |
216 |
if( L.eq.1 ) then |
217 |
array(5) = -undef |
218 |
else |
219 |
array(5) = pq(i,j,L-1)*d2p(i,j)*dp(i,j,L) |
220 |
endif |
221 |
if( L.eq.lm ) then |
222 |
array(6) = -undef |
223 |
else |
224 |
array(6) = pq(i,j,L+1)*d2p(i,j)*dp(i,j,L) |
225 |
endif |
226 |
|
227 |
call maxval1 (array,6,-qval,qmax,imax) |
228 |
|
229 |
if( imax.eq.0 ) then |
230 |
sum = 0.0 |
231 |
do m=1,6 |
232 |
if( array(m).gt.0.0 ) sum = sum + array(m) |
233 |
enddo |
234 |
if( sum.gt.-qval ) then |
235 |
fact = 1.0 + qval/sum |
236 |
if( array(1).gt.0 ) pq(ip1,j,L) = pq(ip1,j,L) * fact |
237 |
if( array(2).gt.0 ) pq(im1,j,L) = pq(im1,j,L) * fact |
238 |
if( array(3).gt.0 ) pq(i,j+1,L) = pq(i,j+1,L) * fact |
239 |
if( array(4).gt.0 ) pq(i,j-1,L) = pq(i,j-1,L) * fact |
240 |
if( array(5).gt.0 ) pq(i,j,L-1) = pq(i,j,L-1) * fact |
241 |
if( array(6).gt.0 ) pq(i,j,L+1) = pq(i,j,L+1) * fact |
242 |
pq(i,j,L) = 0.0 |
243 |
else |
244 |
pq(i,j,L) = 0.0 |
245 |
endif |
246 |
else |
247 |
if( imax.eq.1 ) pq(ip1,j,L) = pq(ip1,j,L) + |
248 |
. pq(i,j,L)*d2p(i,j)/d2p(ip1,j) |
249 |
if( imax.eq.2 ) pq(im1,j,L) = pq(im1,j,L) + |
250 |
. pq(i,j,L)*d2p(i,j)/d2p(im1,j) |
251 |
if( imax.eq.3 ) pq(i,j+1,L) = pq(i,j+1,L) + |
252 |
. pq(i,j,L)*d2p(i,j)/d2p(i,j+1) |
253 |
if( imax.eq.4 ) pq(i,j-1,L) = pq(i,j-1,L) + |
254 |
. pq(i,j,L)*d2p(i,j)/d2p(i,j-1) |
255 |
if( imax.eq.5 ) pq(i,j,L-1) = pq(i,j,L-1) + |
256 |
. pq(i,j,L)*dp(i,j,L) /dp(i,j,L-1) |
257 |
if( imax.eq.6 ) pq(i,j,L+1) = pq(i,j,L+1) + |
258 |
. pq(i,j,L)*dp(i,j,L) /dp(i,j,L+1) |
259 |
pq(i,j,L) = 0.0 |
260 |
endif |
261 |
|
262 |
endif ! End pq<0 Test |
263 |
|
264 |
im1 = i |
265 |
i = ip1 |
266 |
enddo |
267 |
enddo |
268 |
enddo |
269 |
|
270 |
return |
271 |
end |
272 |
|
273 |
subroutine maxval1 (q,im,qval,qmax,imax) |
274 |
C*********************************************************************** |
275 |
C PURPOSE |
276 |
C Find the location and value of the array element which is greater |
277 |
C than a prescribed value. |
278 |
C |
279 |
C INPUT |
280 |
C q ...... Array Elements |
281 |
C im ..... Dimension of Array q |
282 |
C qval ... Prescribed Value |
283 |
C |
284 |
C OUTPUT |
285 |
C qmax ... Largest Array element which is greater than qval |
286 |
C imax ... Location of Largest Array Element |
287 |
C |
288 |
C Note: |
289 |
C If no array element is larger than qval, then imax = 0 |
290 |
C |
291 |
C*********************************************************************** |
292 |
implicit none |
293 |
integer im, i, imax |
294 |
_RL q(im), qmax, qval |
295 |
qmax = qval |
296 |
imax = 0 |
297 |
do i=1,im |
298 |
if( q(i).gt.qmax ) then |
299 |
qmax = q(i) |
300 |
imax = i |
301 |
endif |
302 |
enddo |
303 |
return |
304 |
end |