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