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