UPP  11.0.0
 All Data Structures Files Functions Variables Pages
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
301 30 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
ctlblk_mod
Definition: CTLBLK.f:1
calhel2
subroutine calhel2(LLOW, LUPP, DEPTH, UST, VST, HELI, CANGLE)
Subroutine that computes storm relative helicity.
Definition: CALHEL2.f:57
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
lookup_mod
Definition: LOOKUP.f:1