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