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