UPP (develop)
Loading...
Searching...
No Matches
CALHEL3.f
Go to the documentation of this file.
1
3!
48!-----------------------------------------------------------------------
49 SUBROUTINE calhel3(LLOW,LUPP,UST,VST,HELI)
50
51!
52 use vrbls3d, only: zmid, uh, vh, u, v, zint
53 use vrbls2d, only: fis, u10, v10
54 use masks, only: lmv
55 use params_mod, only: g
56 use lookup_mod, only: itb,jtb,itbq,jtbq
57 use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
58 lm, im, jm, me, spval, &
59 ista, iend, ista_m, iend_m, ista_2l, iend_2u
60 use gridspec_mod, only: gridtype
61!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62 implicit none
63!
64 real,PARAMETER :: P150=15000.0,p300=30000.0,s15=15.0
65 real,PARAMETER :: D3000=3000.0,pi6=0.5235987756,pi9=0.34906585
66 real,PARAMETER :: D5500=5500.0,d6000=6000.0,d7000=7000.0
67 real,PARAMETER :: D500=500.0
68! CRA
69 real,PARAMETER :: D1000=1000.0
70 real,PARAMETER :: D1500=1500.0
71! CRA
72 REAL, PARAMETER :: pi = 3.1415927
73
74!
75! DECLARE VARIABLES
76!
77 integer,dimension(ista_2l:iend_2u,jsta_2l:jend_2u),intent(in) :: LLOW, LUPP
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),intent(out) :: HELI
80!
81 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: HTSFC, UST6, VST6, UST5, VST5, &
82 ust1, vst1, ushr1, vshr1, &
83 ushr6, vshr6, u1, v1, u2, v2, &
84 hgt1, hgt2, umean, vmean
85 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: USHR05,VSHR05,ELT,ELB
86
87! REAL HTSFC(IM,JM)
88!
89! REAL UST6(IM,JM),VST6(IM,JM)
90! REAL UST5(IM,JM),VST5(IM,JM)
91! REAL UST1(IM,JM),VST1(IM,JM)
92! CRA
93! REAL USHR1(IM,JM),VSHR1(IM,JM),USHR6(IM,JM),VSHR6(IM,JM)
94! REAL U1(IM,JM),V1(IM,JM),U2(IM,JM),V2(IM,JM)
95! REAL HGT1(IM,JM),HGT2(IM,JM),UMEAN(IM,JM),VMEAN(IM,JM)
96! CRA
97
98 integer, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: COUNT6, COUNT5, COUNT1, L1, L2
99! INTEGER COUNT6(IM,JM),COUNT5(IM,JM),COUNT1(IM,JM)
100! CRA
101! INTEGER L1(IM,JM),L2(IM,JM)
102! CRA
103
104 INTEGER IVE(JM),IVW(JM)
105 integer I,J,IW,IE,JS,JN,JVN,JVS,L,N,lv
106 integer ISTART,ISTOP,JSTART,JSTOP
107 real Z2,DZABV,UMEAN5,VMEAN5,UMEAN1,VMEAN1,UMEAN6,VMEAN6, &
108 denom,z1,z3,dz,dz1,dz2,du1,du2,dv1,dv2
109!
110!****************************************************************
111! START CALHEL HERE
112!
113! INITIALIZE ARRAYS.
114!
115!$omp parallel do private(i,j)
116 DO j=jsta,jend
117 DO i=ista,iend
118 ust(i,j) = 0.0
119 vst(i,j) = 0.0
120 heli(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
282 IF (dzabv < d7000 .AND. dzabv >= d6000) THEN
283 ust5(i,j) = ust5(i,j) + uh(i,j,l)
284 vst5(i,j) = vst5(i,j) + vh(i,j,l)
285 count5(i,j) = 1
286 GOTO 30
287 ENDIF
288 ENDDO
289 ENDIF
29030 CONTINUE
291 ENDDO
292 ENDDO
293
294!
295!$omp parallel do private(i,j,umean6,vmean6,umean5,vmean5,umean1,vmean1,denom)
296
297 DO j=jstart,jstop
298 DO i=istart,istop
299 IF (count6(i,j) > 0 .AND. count1(i,j) > 0 .AND. count5(i,j) > 0) THEN
300 umean5 = ust5(i,j) / count5(i,j)
301 vmean5 = vst5(i,j) / count5(i,j)
302 umean1 = ust1(i,j) / count1(i,j)
303 vmean1 = vst1(i,j) / count1(i,j)
304 umean6 = ust6(i,j) / count6(i,j)
305 vmean6 = vst6(i,j) / count6(i,j)
306
307!
308! COMPUTE STORM MOTION VECTOR
309! IT IS DEFINED AS 7.5 M/S TO THE RIGHT OF THE 0-6 KM MEAN
310! WIND CONSTRAINED ALONG A LINE WHICH IS BOTH PERPENDICULAR
311! TO THE 0-6 KM MEAN VERTICAL WIND SHEAR VECTOR AND PASSES
312! THROUGH THE 0-6 KM MEAN WIND. THE WIND SHEAR VECTOR IS
313! SET AS THE DIFFERENCE BETWEEN THE 5.5-6 KM WIND (THE HEAD
314! OF THE SHEAR VECTOR) AND THE 0-0.5 KM WIND (THE TAIL).
315! THIS IS FOR THE RIGHT-MOVING CASE; WE IGNORE THE LEFT MOVER.
316
317! CRA
318 ushr6(i,j) = umean5 - umean1
319 vshr6(i,j) = vmean5 - vmean1
320
321 denom = ushr6(i,j)*ushr6(i,j)+vshr6(i,j)*vshr6(i,j)
322 IF (denom /= 0.0) THEN
323 ust(i,j) = umean6 + (7.5*vshr6(i,j)/sqrt(denom))
324 vst(i,j) = vmean6 - (7.5*ushr6(i,j)/sqrt(denom))
325 ELSE
326 ust(i,j) = 0
327 vst(i,j) = 0
328 ENDIF
329 ELSE
330 ust(i,j) = 0.0
331 vst(i,j) = 0.0
332 ushr6(i,j) = 0.0
333 vshr6(i,j) = 0.0
334 ENDIF
335
336 IF(l1(i,j) > 0 .AND. l2(i,j) > 0) THEN
337 umean(i,j) = u1(i,j) + (d1000 - hgt1(i,j))*(u2(i,j) - &
338 u1(i,j))/(hgt2(i,j) - hgt1(i,j))
339 vmean(i,j) = v1(i,j) + (d1000 - hgt1(i,j))*(v2(i,j) - &
340 v1(i,j))/(hgt2(i,j) - hgt1(i,j))
341 ELSE IF(l1(i,j) > 0 .AND. l2(i,j) == 0) THEN
342 umean(i,j) = u1(i,j)
343 vmean(i,j) = v1(i,j)
344 ELSE IF(l1(i,j) == 0 .AND. l2(i,j) > 0) THEN
345 umean(i,j) = u2(i,j)
346 vmean(i,j) = u2(i,j)
347 ELSE
348 umean(i,j) = 0.0
349 vmean(i,j) = 0.0
350 ENDIF
351
352 IF(l1(i,j) > 0 .OR. l2(i,j) > 0) THEN
353 ushr05(i,j) = umean1 - u10(i,j)
354 vshr05(i,j) = vmean1 - v10(i,j)
355 ushr1(i,j) = umean(i,j) - u10(i,j)
356 vshr1(i,j) = vmean(i,j) - v10(i,j)
357 ELSE
358 ushr05(i,j) = 0.0
359 vshr05(i,j) = 0.0
360 ushr1(i,j) = 0.0
361 vshr1(i,j) = 0.0
362 ENDIF
363! CRA
364
365!tgs USHR = UMEAN5 - UMEAN1
366! VSHR = VMEAN5 - VMEAN1
367
368! UST(I,J) = UMEAN6 + (7.5*VSHR/SQRT(USHR*USHR+VSHR*VSHR))
369! VST(I,J) = VMEAN6 - (7.5*USHR/SQRT(USHR*USHR+VSHR*VSHR))
370! ELSE
371! UST(I,J) = 0.0
372! VST(I,J) = 0.0
373! ENDIF
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 l = 2,lm-1
381 if(gridtype /= 'A')then
382 call exch(zint(1,jsta_2l,l))
383 call exch(zint(1,jsta_2l,l+1))
384 end if
385 DO j=jstart,jstop
386 DO i=istart,istop
387 iw=i+ivw(j)
388 ie=i+ive(j)
389 jn=j+jvn
390 js=j+jvs
391 IF (gridtype=='B')THEN
392 z2=0.25*(zmid(iw,j,l)+zmid(ie,j,l)+ &
393 zmid(i,jn,l)+zmid(ie,jn,l))
394 ELSE
395 z2=0.25*(zmid(iw,j,l)+zmid(ie,j,l)+ &
396 zmid(i,jn,l)+zmid(i,js,l))
397 END IF
398 dzabv=z2-htsfc(i,j)
399 elt(i,j) = zint(i,j,lupp(i,j))-htsfc(i,j)
400 elb(i,j) = zint(i,j,llow(i,j))-htsfc(i,j)
401
402!
403 IF(dzabv <= elt(i,j) .AND. dzabv >= elb(i,j) .AND. l <= nint(lmv(i,j)))THEN
404 IF (gridtype=='B')THEN
405 z1 = 0.25*(zmid(iw,j,l+1)+zmid(ie,j,l+1)+ &
406 zmid(i,jn,l+1)+zmid(ie,jn,l+1))
407 z3 = 0.25*(zmid(iw,j,l-1)+zmid(ie,j,l-1)+ &
408 zmid(i,jn,l-1)+zmid(ie,jn,l-1))
409 dz = 0.25*((zint(iw,j,l)+zint(ie,j,l)+ &
410 zint(i,jn,l)+zint(ie,jn,l))- &
411 (zint(iw,j,l+1)+zint(ie,j,l+1)+ &
412 zint(i,jn,l+1)+zint(ie,jn,l+1)))
413 ELSE
414 z1 = 0.25*(zmid(iw,j,l+1)+zmid(ie,j,l+1)+ &
415 zmid(i,jn,l+1)+zmid(i,js,l+1))
416 z3 = 0.25*(zmid(iw,j,l-1)+zmid(ie,j,l-1)+ &
417 zmid(i,jn,l-1)+zmid(i,js,l-1))
418 dz = 0.25*((zint(iw,j,l)+zint(ie,j,l)+ &
419 zint(i,js,l)+zint(i,jn,l))- &
420 (zint(iw,j,l+1)+zint(ie,j,l+1)+ &
421 zint(i,js,l+1)+zint(i,jn,l+1)))
422 END IF
423 dz1 = z1-z2
424 dz2 = z2-z3
425 du1 = uh(i,j,l+1)-uh(i,j,l)
426 du2 = uh(i,j,l)-uh(i,j,l-1)
427 dv1 = vh(i,j,l+1)-vh(i,j,l)
428 dv2 = vh(i,j,l)-vh(i,j,l-1)
429 IF( l >= lupp(i,j) .AND. l <= llow(i,j) ) THEN
430 IF( vh(i,j,l) <spval.and.uh(i,j,l) <spval.and. &
431 vh(i,j,l+1)<spval.and.uh(i,j,l+1)<spval.and. &
432 vh(i,j,l-1)<spval.and.uh(i,j,l-1)<spval.and. &
433 vst(i,j) <spval.and.ust(i,j) <spval) &
434 heli(i,j) = ((vh(i,j,l)-vst(i,j))* &
435 (dz2*(du1/dz1)+dz1*(du2/dz2)) &
436 - (uh(i,j,l)-ust(i,j))* &
437 (dz2*(dv1/dz1)+dz1*(dv2/dz2))) &
438 *dz/(dz1+dz2)+heli(i,j)
439 ENDIF
440 IF(lupp(i,j) == llow(i,j)) heli(i,j) = 0.
441
442! if(i==im/2.and.j==(jsta+jend)/2)print*,'Debug Helicity',depth(N),l,dz1,dz2,du1, &
443! du2,dv1,dv2,ust(i,j),vst(i,j)
444 ENDIF
445 ENDDO
446 ENDDO
447 ENDDO
448
449! END OF ROUTINE.
450!
451 RETURN
452 END
calhel3
subroutine calhel3(llow, lupp, ust, vst, heli)
This routine computes estimated storm motion and storm-relative environmental helicity.
Definition CALHEL3.f:50
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27