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