UPP (develop)
Loading...
Searching...
No Matches
CALHEL.f
Go to the documentation of this file.
1
34! -----------------------------------------------------------------------------------
46 SUBROUTINE calhel(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6)
47
48!
49 use vrbls3d, only: zmid, uh, vh, u, v, zint
50 use vrbls2d, only: fis, u10, v10
51 use masks, only: lmv, gdlat
52 use params_mod, only: g
53 use lookup_mod, only: itb,jtb,itbq,jtbq
54 use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
55 lm, im, jm, me, spval, &
56 ista, iend, ista_m, iend_m, ista_2l, iend_2u
57 use gridspec_mod, only: gridtype
58!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 implicit none
60!
61 real,PARAMETER :: P150=15000.0,p300=30000.0,s15=15.0
62 real,PARAMETER :: D3000=3000.0,pi6=0.5235987756,pi9=0.34906585
63 real,PARAMETER :: D5500=5500.0,d6000=6000.0,d7000=7000.0
64 real,PARAMETER :: D500=500.0
65! CRA
66 real,PARAMETER :: D1000=1000.0
67 real,PARAMETER :: D1500=1500.0
68! CRA
69
70!
71! DECLARE VARIABLES
72!
73 real,intent(in) :: DEPTH(2)
74 REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) :: UST,VST
75 REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u,2),intent(out) :: HELI
76!
77 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: HTSFC, UST6, VST6, UST5, VST5, &
78 ust1, vst1, ushr1, vshr1, &
79 ushr6, vshr6, u1, v1, u2, v2, &
80 hgt1, hgt2, umean, vmean
81! REAL HTSFC(IM,JM)
82!
83! REAL UST6(IM,JM),VST6(IM,JM)
84! REAL UST5(IM,JM),VST5(IM,JM)
85! REAL UST1(IM,JM),VST1(IM,JM)
86! CRA
87! REAL USHR1(IM,JM),VSHR1(IM,JM),USHR6(IM,JM),VSHR6(IM,JM)
88! REAL U1(IM,JM),V1(IM,JM),U2(IM,JM),V2(IM,JM)
89! REAL HGT1(IM,JM),HGT2(IM,JM),UMEAN(IM,JM),VMEAN(IM,JM)
90! CRA
91
92 integer, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: COUNT6, COUNT5, COUNT1, L1, L2
93! INTEGER COUNT6(IM,JM),COUNT5(IM,JM),COUNT1(IM,JM)
94! CRA
95! INTEGER L1(IM,JM),L2(IM,JM)
96! CRA
97
98 INTEGER IVE(JM),IVW(JM)
99 integer I,J,IW,IE,JS,JN,JVN,JVS,L,N,lv
100 integer ISTART,ISTOP,JSTART,JSTOP
101 real Z2,DZABV,UMEAN5,VMEAN5,UMEAN1,VMEAN1,UMEAN6,VMEAN6, &
102 denom,z1,z3,dz,dz1,dz2,du1,du2,dv1,dv2
103
104!
105!****************************************************************
106! START CALHEL HERE
107!
108! INITIALIZE ARRAYS.
109!
110!$omp parallel do private(i,j)
111 DO j=jsta,jend
112 DO i=ista,iend
113 ust(i,j) = 0.0
114 vst(i,j) = 0.0
115 heli(i,j,1) = 0.0
116 heli(i,j,2) = 0.0
117 ust1(i,j) = 0.0
118 vst1(i,j) = 0.0
119 ust5(i,j) = 0.0
120 vst5(i,j) = 0.0
121 ust6(i,j) = 0.0
122 vst6(i,j) = 0.0
123 count6(i,j) = 0
124 count5(i,j) = 0
125 count1(i,j) = 0
126! CRA
127 ushr1(i,j) = 0.0
128 vshr1(i,j) = 0.0
129 ushr6(i,j) = 0.0
130 vshr6(i,j) = 0.0
131 u1(i,j) = 0.0
132 u2(i,j) = 0.0
133 v1(i,j) = 0.0
134 v2(i,j) = 0.0
135 umean(i,j) = 0.0
136 vmean(i,j) = 0.0
137 hgt1(i,j) = 0.0
138 hgt2(i,j) = 0.0
139 l1(i,j) = 0
140 l2(i,j) = 0
141! CRA
142
143 ENDDO
144 ENDDO
145 IF(gridtype == 'E')THEN
146 jvn = 1
147 jvs = -1
148 do j=jsta,jend
149 ive(j) = mod(j,2)
150 ivw(j) = ive(j)-1
151 enddo
152 istart = ista_m
153 istop = iend_m
154 jstart = jsta_m
155 jstop = jend_m
156 ELSE IF(gridtype == 'B')THEN
157 jvn = 1
158 jvs = 0
159 do j=jsta,jend
160 ive(j)=1
161 ivw(j)=0
162 enddo
163 istart = ista_m
164 istop = iend_m
165 jstart = jsta_m
166 jstop = jend_m
167 ELSE
168 jvn = 0
169 jvs = 0
170 do j=jsta,jend
171 ive(j) = 0
172 ivw(j) = 0
173 enddo
174 istart = ista
175 istop = iend
176 jstart = jsta
177 jstop = jend
178 END IF
179!
180! LOOP OVER HORIZONTAL GRID.
181!
182! CALL EXCH(RES(1,jsta_2l)
183! CALL EXCH(PD()
184
185! DO L = 1,LP1
186! CALL EXCH(ZINT(1,jsta_2l,L))
187! END DO
188!
189!!$omp parallel do private(htsfc,ie,iw)
190 IF(gridtype /= 'A') CALL exch(fis(ista_2l:iend_2u,jsta_2l:jend_2u))
191 DO l = 1,lm
192 IF(gridtype /= 'A') CALL exch(zmid(ista_2l:iend_2u,jsta_2l:jend_2u,l))
193 DO j=jstart,jstop
194 DO i=istart,istop
195 ie = i+ive(j)
196 iw = i+ivw(j)
197 jn = j+jvn
198 js = j+jvs
199!mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+
200!mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25
201!mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT
202 IF (gridtype=='B')THEN
203 htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(ie,jn))
204!
205! COMPUTE MASS WEIGHTED MEAN WIND IN THE 0-6 KM LAYER, THE
206! 0-0.5 KM LAYER, AND THE 5.5-6 KM LAYER
207!
208 z2 = 0.25*(zmid(iw,j,l)+zmid(ie,j,l)+zmid(i,jn,l)+zmid(ie,jn,l))
209 ELSE
210 htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(i,js))
211!
212! COMPUTE MASS WEIGHTED MEAN WIND IN THE 0-6 KM LAYER, THE
213! 0-0.5 KM LAYER, AND THE 5.5-6 KM LAYER
214!
215 z2 = 0.25*(zmid(iw,j,l)+zmid(ie,j,l)+zmid(i,jn,l)+zmid(i,js,l))
216 END IF
217 dzabv = z2-htsfc(i,j)
218
219 lv = nint(lmv(i,j))
220 IF (dzabv <= d6000 .AND. l <= lv) THEN
221 ust6(i,j) = ust6(i,j) + uh(i,j,l)
222 vst6(i,j) = vst6(i,j) + vh(i,j,l)
223 count6(i,j) = count6(i,j) + 1
224 ENDIF
225
226 IF (dzabv < d6000 .AND. dzabv >= d5500 .AND. l <= lv) THEN
227 ust5(i,j) = ust5(i,j) + uh(i,j,l)
228 vst5(i,j) = vst5(i,j) + vh(i,j,l)
229 count5(i,j) = count5(i,j) + 1
230 ENDIF
231
232 IF (dzabv < d500 .AND. l <= lv) THEN
233 ust1(i,j) = ust1(i,j) + uh(i,j,l)
234 vst1(i,j) = vst1(i,j) + vh(i,j,l)
235 count1(i,j) = count1(i,j) + 1
236 ENDIF
237! CRA
238 IF (dzabv >= d1000 .AND. dzabv <= d1500 .AND. l <= lv) THEN
239 u2(i,j) = u(i,j,l)
240 v2(i,j) = v(i,j,l)
241 hgt2(i,j) = dzabv
242 l2(i,j) = l
243 ENDIF
244
245 IF (dzabv >= d500 .AND. dzabv < d1000 .AND. &
246 l <= lv .AND. l1(i,j) <= l2(i,j)) THEN
247 u1(i,j) = u(i,j,l)
248 v1(i,j) = v(i,j,l)
249 hgt1(i,j) = dzabv
250 l1(i,j) = l
251 ENDIF
252! CRA
253
254 ENDDO
255 ENDDO
256 ENDDO
257!
258! CASE WHERE THERE IS NO LEVEL WITH HEIGHT BETWEEN 5500 AND 6000
259
260 DO j=jstart,jstop
261 DO i=istart,istop
262 IF (count5(i,j) == 0) THEN
263 lloop: DO l=lm,1,-1
264 ie=i+ive(j)
265 iw=i+ivw(j)
266 jn=j+jvn
267 js=j+jvs
268 IF (gridtype=='B')THEN
269 z2 = 0.25*(zmid(iw,j,l)+zmid(ie,j,l)+zmid(i,jn,l)+zmid(ie,jn,l))
270 ELSE
271 z2 = 0.25*(zmid(iw,j,l)+zmid(ie,j,l)+zmid(i,jn,l)+zmid(i,js,l))
272 END IF
273
274 dzabv=z2-htsfc(i,j)
275 IF (dzabv < d7000 .AND. dzabv >= d6000) THEN
276 ust5(i,j) = ust5(i,j) + uh(i,j,l)
277 vst5(i,j) = vst5(i,j) + vh(i,j,l)
278 count5(i,j) = 1
279 exit lloop
280 ENDIF
281 enddo lloop
282 ENDIF
283 ENDDO
284 ENDDO
285
286!
287!$omp parallel do private(i,j,umean6,vmean6,umean5,vmean5,umean1,vmean1,denom)
288
289 DO j=jstart,jstop
290 DO i=istart,istop
291 IF (count6(i,j) > 0 .AND. count1(i,j) > 0 .AND. count5(i,j) > 0) THEN
292 umean5 = ust5(i,j) / count5(i,j)
293 vmean5 = vst5(i,j) / count5(i,j)
294 umean1 = ust1(i,j) / count1(i,j)
295 vmean1 = vst1(i,j) / count1(i,j)
296 umean6 = ust6(i,j) / count6(i,j)
297 vmean6 = vst6(i,j) / count6(i,j)
298
299!
300! COMPUTE STORM MOTION VECTOR
301! IT IS DEFINED AS 7.5 M/S TO THE RIGHT OF THE 0-6 KM MEAN
302! WIND CONSTRAINED ALONG A LINE WHICH IS BOTH PERPENDICULAR
303! TO THE 0-6 KM MEAN VERTICAL WIND SHEAR VECTOR AND PASSES
304! THROUGH THE 0-6 KM MEAN WIND. THE WIND SHEAR VECTOR IS
305! SET AS THE DIFFERENCE BETWEEN THE 5.5-6 KM WIND (THE HEAD
306! OF THE SHEAR VECTOR) AND THE 0-0.5 KM WIND (THE TAIL).
307! THIS IS FOR THE RIGHT-MOVING CASE; WE IGNORE THE LEFT MOVER.
308
309! CRA
310 ushr6(i,j) = umean5 - umean1
311 vshr6(i,j) = vmean5 - vmean1
312
313 denom = ushr6(i,j)*ushr6(i,j)+vshr6(i,j)*vshr6(i,j)
314 IF (denom /= 0.0) THEN
315 IF (gdlat(i,j) >= 0) THEN
316 ust(i,j) = umean6 + (7.5*vshr6(i,j)/sqrt(denom))
317 vst(i,j) = vmean6 - (7.5*ushr6(i,j)/sqrt(denom))
318 ELSE IF (gdlat(i,j) < 0) THEN
319 ust(i,j) = umean6 - (7.5*vshr6(i,j)/sqrt(denom))
320 vst(i,j) = vmean6 + (7.5*ushr6(i,j)/sqrt(denom))
321 ENDIF
322 ELSE
323 ust(i,j) = 0
324 vst(i,j) = 0
325 ENDIF
326 ELSE
327 ust(i,j) = 0.0
328 vst(i,j) = 0.0
329 ushr6(i,j) = 0.0
330 vshr6(i,j) = 0.0
331 ENDIF
332
333 IF(l1(i,j) > 0 .AND. l2(i,j) > 0) THEN
334 umean(i,j) = u1(i,j) + (d1000 - hgt1(i,j))*(u2(i,j) - &
335 u1(i,j))/(hgt2(i,j) - hgt1(i,j))
336 vmean(i,j) = v1(i,j) + (d1000 - hgt1(i,j))*(v2(i,j) - &
337 v1(i,j))/(hgt2(i,j) - hgt1(i,j))
338 ELSE IF(l1(i,j) > 0 .AND. l2(i,j) == 0) THEN
339 umean(i,j) = u1(i,j)
340 vmean(i,j) = v1(i,j)
341 ELSE IF(l1(i,j) == 0 .AND. l2(i,j) > 0) THEN
342 umean(i,j) = u2(i,j)
343 vmean(i,j) = u2(i,j)
344 ELSE
345 umean(i,j) = 0.0
346 vmean(i,j) = 0.0
347 ENDIF
348
349 IF(l1(i,j) > 0 .OR. l2(i,j) > 0) THEN
350 ushr1(i,j) = umean(i,j) - u10(i,j)
351 vshr1(i,j) = vmean(i,j) - v10(i,j)
352 ELSE
353 ushr1(i,j) = 0.0
354 vshr1(i,j) = 0.0
355 ENDIF
356! CRA
357
358!tgs USHR = UMEAN5 - UMEAN1
359! VSHR = VMEAN5 - VMEAN1
360
361! UST(I,J) = UMEAN6 + (7.5*VSHR/SQRT(USHR*USHR+VSHR*VSHR))
362! VST(I,J) = VMEAN6 - (7.5*USHR/SQRT(USHR*USHR+VSHR*VSHR))
363! ELSE
364! UST(I,J) = 0.0
365! VST(I,J) = 0.0
366! ENDIF
367
368 ENDDO
369 ENDDO
370!
371! COMPUTE STORM-RELATIVE HELICITY
372!
373!!$omp parallel do private(i,j,n,l,du1,du2,dv1,dv2,dz,dz1,dz2,dzabv,ie,iw,jn,js,z1,z2,z3)
374 DO n=1,2 ! for dfferent helicity depth
375 DO l = 2,lm-1
376 if(gridtype /= 'A')then
377 call exch(zint(1,jsta_2l,l))
378 call exch(zint(1,jsta_2l,l+1))
379 end if
380 DO j=jstart,jstop
381 DO i=istart,istop
382 iw=i+ivw(j)
383 ie=i+ive(j)
384 jn=j+jvn
385 js=j+jvs
386 IF (gridtype=='B')THEN
387 z2=0.25*(zmid(iw,j,l)+zmid(ie,j,l)+ &
388 zmid(i,jn,l)+zmid(ie,jn,l))
389 ELSE
390 z2=0.25*(zmid(iw,j,l)+zmid(ie,j,l)+ &
391 zmid(i,jn,l)+zmid(i,js,l))
392 END IF
393 dzabv=z2-htsfc(i,j)
394!
395 IF(dzabv < depth(n) .AND. l <= nint(lmv(i,j)))THEN
396 IF (gridtype=='B')THEN
397 z1 = 0.25*(zmid(iw,j,l+1)+zmid(ie,j,l+1)+ &
398 zmid(i,jn,l+1)+zmid(ie,jn,l+1))
399 z3 = 0.25*(zmid(iw,j,l-1)+zmid(ie,j,l-1)+ &
400 zmid(i,jn,l-1)+zmid(ie,jn,l-1))
401 dz = 0.25*((zint(iw,j,l)+zint(ie,j,l)+ &
402 zint(i,jn,l)+zint(ie,jn,l))- &
403 (zint(iw,j,l+1)+zint(ie,j,l+1)+ &
404 zint(i,jn,l+1)+zint(ie,jn,l+1)))
405 ELSE
406 z1 = 0.25*(zmid(iw,j,l+1)+zmid(ie,j,l+1)+ &
407 zmid(i,jn,l+1)+zmid(i,js,l+1))
408 z3 = 0.25*(zmid(iw,j,l-1)+zmid(ie,j,l-1)+ &
409 zmid(i,jn,l-1)+zmid(i,js,l-1))
410 dz = 0.25*((zint(iw,j,l)+zint(ie,j,l)+ &
411 zint(i,js,l)+zint(i,jn,l))- &
412 (zint(iw,j,l+1)+zint(ie,j,l+1)+ &
413 zint(i,js,l+1)+zint(i,jn,l+1)))
414 END IF
415 dz1 = z1-z2
416 dz2 = z2-z3
417 du1 = uh(i,j,l+1)-uh(i,j,l)
418 du2 = uh(i,j,l)-uh(i,j,l-1)
419 dv1 = vh(i,j,l+1)-vh(i,j,l)
420 dv2 = vh(i,j,l)-vh(i,j,l-1)
421 IF( vh(i,j,l) <spval.and.uh(i,j,l) <spval.and. &
422 vh(i,j,l+1)<spval.and.uh(i,j,l+1)<spval.and. &
423 vh(i,j,l-1)<spval.and.uh(i,j,l-1)<spval.and. &
424 vst(i,j) <spval.and.ust(i,j) <spval) &
425 heli(i,j,n) = ((vh(i,j,l)-vst(i,j))* &
426 (dz2*(du1/dz1)+dz1*(du2/dz2)) &
427 - (uh(i,j,l)-ust(i,j))* &
428 (dz2*(dv1/dz1)+dz1*(dv2/dz2))) &
429 *dz/(dz1+dz2)+heli(i,j,n)
430
431! if(i==im/2.and.j==(jsta+jend)/2)print*,'Debug Helicity',depth(N),l,dz1,dz2,du1, &
432! du2,dv1,dv2,ust(i,j),vst(i,j)
433 ENDIF
434 ENDDO
435 ENDDO
436 ENDDO
437 END DO ! end of different helicity depth
438!
439! END OF ROUTINE.
440!
441 RETURN
442 END
calhel
subroutine calhel(depth, ust, vst, heli, ushr1, vshr1, ushr6, vshr6)
Computes storm relative helicity and storm motion.
Definition CALHEL.f:47
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27