SUBROUTINE CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: CALHEL COMPUTES STORM RELATIVE HELICITY
! PRGRMMR: BALDWIN ORG: W/NP2 DATE: 94-08-22
!
! ABSTRACT:
! THIS ROUTINE COMPUTES ESTIMATED STORM MOTION AND
! STORM-RELATIVE ENVIRONMENTAL HELICITY.
! (DAVIES-JONES ET AL 1990) THE ALGORITHM PROCEEDS AS
! FOLLOWS.
!
! THE STORM MOTION COMPUTATION NO LONGER EMPLOYS THE DAVIES AND
! JOHNS (1993) METHOD WHICH DEFINED STORM MOTION AS 30 DEGREES TO
! THE RIGHT OF THE 0-6 KM MEAN WIND AT 75% OF THE SPEED FOR MEAN
! SPEEDS LESS THAN 15 M/S AND 20 DEGREES TO THE RIGHT FOR SPEEDS
! GREATER THAN 15 M/S. INSTEAD, WE NOW USE THE DYNAMIC METHOD
! (BUNKERS ET AL. 1998) WHICH HAS BEEN FOUND TO DO BETTER IN
! CASES WITH 'NON-CLASSIC' HODOGRAPHS (SUCH AS NORTHWEST-FLOW
! EVENTS) AND DO AS WELL OR BETTER THAN THE OLD METHOD IN MORE
! CLASSIC SITUATIONS.
!
! PROGRAM HISTORY LOG:
! 94-08-22 MICHAEL BALDWIN
! 97-03-27 MICHAEL BALDWIN - SPEED UP CODE
! 98-06-15 T BLACK - CONVERSION FROM 1-D TO 2-D
! 00-01-04 JIM TUCCILLO - MPI VERSION
! 00-01-10 G MANIKIN - CHANGED TO BUNKERS METHOD
! 02-05-22 G MANIKIN - NOW ALLOW CHOICE OF COMPUTING
! HELICITY OVER TWO DIFFERENT
! (0-1 and 0-3 KM) DEPTHS
! 03-03-25 G MANIKIN - MODIFIED CODE TO COMPUTE MEAN WINDS
! USING ARITHMETIC AVERAGES INSTEAD OF
! MASS WEIGHTING; DIFFERENCES ARE MINOR
! BUT WANT TO BE CONSISTENT WITH THE
! BUNKERS METHOD
! 04-04-16 M PYLE - MINIMAL MODIFICATIONS, BUT PUT INTO
! NMM WRFPOST CODE
! 05=02-25 H CHUANG - ADD COMPUTATION FOR ARW A GRID
! 05-07-07 BINBIN ZHOU - ADD RSM FOR A GRID
!
! USAGE: CALHEL(UST,VST,HELI)
! INPUT ARGUMENT LIST:
! DPTH - DEPTH IN METERS OVER WHICH HELICITY SHOULD BE COMPUTED;
! ALLOWS ONE TO DISTINGUISH 0-3 KM AND 0-1 KM VALUES
!
! OUTPUT ARGUMENT LIST:
! UST - ESTIMATED U COMPONENT (M/S) OF STORM MOTION.
! VST - ESTIMATED V COMPONENT (M/S) OF STORM MOTION.
! HELI - STORM-RELATIVE HELICITY (M**2/S**2)
! CRA
! USHR1 - U COMPONENT (M/S) OF 0-1 KM SHEAR
! VSHR1 - V COMPONENT (M/S) OF 0-1 KM SHEAR
! USHR6 - U COMPONENT (M/S) OF 0-0.5 to 5.5-6.0 KM SHEAR
! VSHR6 - V COMPONENT (M/S) OF 0-0.5 to 5.5-6.0 KM SHEAR
! CRA
!
! OUTPUT FILES:
! NONE
!
! SUBPROGRAMS CALLED:
! UTILITIES:
!
! LIBRARY:
! COMMON - VRBLS
! LOOPS
! PHYS
! EXTRA
! MASKS
! OPTIONS
! INDX
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE : IBM SP
!$$$
!
use vrbls3d, only: zmid, uh, vh, u, v, zint
use vrbls2d, only: fis, u10, v10
use masks, only: lmv
use params_mod, only: g
use lookup_mod, only: ITB,JTB,ITBQ,JTBQ
use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
lm, im, jm, me
use gridspec_mod, only: gridtype
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
implicit none
!
real,PARAMETER :: P150=15000.0,P300=30000.0,S15=15.0
real,PARAMETER :: D3000=3000.0,PI6=0.5235987756,PI9=0.34906585
real,PARAMETER :: D5500=5500.0,D6000=6000.0,D7000=7000.0
real,PARAMETER :: D500=500.0
! CRA
real,PARAMETER :: D1000=1000.0
real,PARAMETER :: D1500=1500.0
! CRA
!
! DECLARE VARIABLES
!
real,intent(in) :: DEPTH(2)
REAL,dimension(IM,jsta_2l:jend_2u), intent(out) :: UST,VST
REAL,dimension(IM,jsta_2l:jend_2u,2),intent(out) :: HELI
!
real, dimension(im,jsta_2l:jend_2u) :: HTSFC, UST6, VST6, UST5, VST5, &
UST1, VST1, USHR1, VSHR1, &
USHR6, VSHR6, U1, V1, U2, V2, &
HGT1, HGT2, UMEAN, VMEAN
! REAL HTSFC(IM,JM)
!
! REAL UST6(IM,JM),VST6(IM,JM)
! REAL UST5(IM,JM),VST5(IM,JM)
! REAL UST1(IM,JM),VST1(IM,JM)
! CRA
! REAL USHR1(IM,JM),VSHR1(IM,JM),USHR6(IM,JM),VSHR6(IM,JM)
! REAL U1(IM,JM),V1(IM,JM),U2(IM,JM),V2(IM,JM)
! REAL HGT1(IM,JM),HGT2(IM,JM),UMEAN(IM,JM),VMEAN(IM,JM)
! CRA
integer, dimension(im,jsta_2l:jend_2u) :: COUNT6, COUNT5, COUNT1, L1, L2
! INTEGER COUNT6(IM,JM),COUNT5(IM,JM),COUNT1(IM,JM)
! CRA
! INTEGER L1(IM,JM),L2(IM,JM)
! CRA
INTEGER IVE(JM),IVW(JM)
integer I,J,IW,IE,JS,JN,JVN,JVS,L,N,lv
integer ISTART,ISTOP,JSTART,JSTOP
real Z2,DZABV,UMEAN5,VMEAN5,UMEAN1,VMEAN1,UMEAN6,VMEAN6, &
DENOM,Z1,Z3,DZ,DZ1,DZ2,DU1,DU2,DV1,DV2
!
!****************************************************************
! START CALHEL HERE
!
! INITIALIZE ARRAYS.
!
!$omp parallel do private(i,j)
DO J=JSTA,JEND
DO I=1,IM
UST(I,J) = 0.0
VST(I,J) = 0.0
HELI(I,J,1) = 0.0
HELI(I,J,2) = 0.0
UST1(I,J) = 0.0
VST1(I,J) = 0.0
UST5(I,J) = 0.0
VST5(I,J) = 0.0
UST6(I,J) = 0.0
VST6(I,J) = 0.0
COUNT6(I,J) = 0
COUNT5(I,J) = 0
COUNT1(I,J) = 0
! CRA
USHR1(I,J) = 0.0
VSHR1(I,J) = 0.0
USHR6(I,J) = 0.0
VSHR6(I,J) = 0.0
U1(I,J) = 0.0
U2(I,J) = 0.0
V1(I,J) = 0.0
V2(I,J) = 0.0
UMEAN(I,J) = 0.0
VMEAN(I,J) = 0.0
HGT1(I,J) = 0.0
HGT2(I,J) = 0.0
L1(I,J) = 0
L2(I,J) = 0
! CRA
ENDDO
ENDDO
IF(gridtype == 'E')THEN
JVN = 1
JVS = -1
do J=JSTA,JEND
IVE(J) = MOD(J,2)
IVW(J) = IVE(J)-1
enddo
ISTART = 2
ISTOP = IM-1
JSTART = JSTA_M
JSTOP = JEND_M
ELSE IF(gridtype == 'B')THEN
JVN = 1
JVS = 0
do J=JSTA,JEND
IVE(J)=1
IVW(J)=0
enddo
ISTART = 2
ISTOP = IM-1
JSTART = JSTA_M
JSTOP = JEND_M
ELSE
JVN = 0
JVS = 0
do J=JSTA,JEND
IVE(J) = 0
IVW(J) = 0
enddo
ISTART = 1
ISTOP = IM
JSTART = JSTA
JSTOP = JEND
END IF
!
! LOOP OVER HORIZONTAL GRID.
!
! CALL EXCH(RES(1,jsta_2l)
! CALL EXCH(PD()
! DO L = 1,LP1
! CALL EXCH(ZINT(1,jsta_2l,L))
! END DO
!
!!$omp parallel do private(htsfc,ie,iw)
IF(gridtype /= 'A') CALL EXCH(FIS(1:IM,JSTA_2L:JEND_2U))
DO L = 1,LM
IF(gridtype /= 'A') CALL EXCH(ZMID(1:IM,JSTA_2L:JEND_2U,L))
DO J=JSTART,JSTOP
DO I=ISTART,ISTOP
IE = I+IVE(J)
IW = I+IVW(J)
JN = J+JVN
JS = J+JVS
!mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+
!mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25
!mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT
IF (gridtype=='B')THEN
HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(IE,JN))
!
! COMPUTE MASS WEIGHTED MEAN WIND IN THE 0-6 KM LAYER, THE
! 0-0.5 KM LAYER, AND THE 5.5-6 KM LAYER
!
Z2 = 0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ZMID(I,JN,L)+ZMID(IE,JN,L))
ELSE
HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(I,JS))
!
! COMPUTE MASS WEIGHTED MEAN WIND IN THE 0-6 KM LAYER, THE
! 0-0.5 KM LAYER, AND THE 5.5-6 KM LAYER
!
Z2 = 0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ZMID(I,JN,L)+ZMID(I,JS,L))
END IF
DZABV = Z2-HTSFC(I,J)
lv = NINT(LMV(I,J))
IF (DZABV <= D6000 .AND. L <= lv) THEN
UST6(I,J) = UST6(I,J) + UH(I,J,L)
VST6(I,J) = VST6(I,J) + VH(I,J,L)
COUNT6(I,J) = COUNT6(I,J) + 1
ENDIF
IF (DZABV < D6000 .AND. DZABV >= D5500 .AND. L <= lv) THEN
UST5(I,J) = UST5(I,J) + UH(I,J,L)
VST5(I,J) = VST5(I,J) + VH(I,J,L)
COUNT5(I,J) = COUNT5(I,J) + 1
ENDIF
IF (DZABV < D500 .AND. L <= lv) THEN
UST1(I,J) = UST1(I,J) + UH(I,J,L)
VST1(I,J) = VST1(I,J) + VH(I,J,L)
COUNT1(I,J) = COUNT1(I,J) + 1
ENDIF
! CRA
IF (DZABV >= D1000 .AND. DZABV <= D1500 .AND. L <= lv) THEN
U2(I,J) = U(I,J,L)
V2(I,J) = V(I,J,L)
HGT2(I,J) = DZABV
L2(I,J) = L
ENDIF
IF (DZABV >= D500 .AND. DZABV < D1000 .AND. &
L <= lv .AND. L1(I,J) <= L2(I,J)) THEN
U1(I,J) = U(I,J,L)
V1(I,J) = V(I,J,L)
HGT1(I,J) = DZABV
L1(I,J) = L
ENDIF
! CRA
ENDDO
ENDDO
ENDDO
!
! CASE WHERE THERE IS NO LEVEL WITH HEIGHT BETWEEN 5500 AND 6000
!
DO J=JSTART,JSTOP
DO I=ISTART,ISTOP
IF (COUNT5(I,J) == 0) THEN
DO L=LM,1,-1
IE=I+IVE(J)
IW=I+IVW(J)
JN=J+JVN
JS=J+JVS
IF (gridtype=='B')THEN
Z2 = 0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ZMID(I,JN,L)+ZMID(IE,JN,L))
ELSE
Z2 = 0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ZMID(I,JN,L)+ZMID(I,JS,L))
END IF
DZABV=Z2-HTSFC(I,J)
IF (DZABV < D7000 .AND. DZABV >= D6000) THEN
UST5(I,J) = UST5(I,J) + UH(I,J,L)
VST5(I,J) = VST5(I,J) + VH(I,J,L)
COUNT5(I,J) = 1
GOTO 30
ENDIF
ENDDO
ENDIF
30 CONTINUE
ENDDO
ENDDO
!
!$omp parallel do private(i,j,umean6,vmean6,umean5,vmean5,umean1,vmean1,denom)
DO J=JSTART,JSTOP
DO I=ISTART,ISTOP
IF (COUNT6(I,J) > 0 .AND. COUNT1(I,J) > 0 .AND. COUNT5(I,J) > 0) THEN
UMEAN5 = UST5(I,J) / COUNT5(I,J)
VMEAN5 = VST5(I,J) / COUNT5(I,J)
UMEAN1 = UST1(I,J) / COUNT1(I,J)
VMEAN1 = VST1(I,J) / COUNT1(I,J)
UMEAN6 = UST6(I,J) / COUNT6(I,J)
VMEAN6 = VST6(I,J) / COUNT6(I,J)
!
! COMPUTE STORM MOTION VECTOR
! IT IS DEFINED AS 7.5 M/S TO THE RIGHT OF THE 0-6 KM MEAN
! WIND CONSTRAINED ALONG A LINE WHICH IS BOTH PERPENDICULAR
! TO THE 0-6 KM MEAN VERTICAL WIND SHEAR VECTOR AND PASSES
! THROUGH THE 0-6 KM MEAN WIND. THE WIND SHEAR VECTOR IS
! SET AS THE DIFFERENCE BETWEEN THE 5.5-6 KM WIND (THE HEAD
! OF THE SHEAR VECTOR) AND THE 0-0.5 KM WIND (THE TAIL).
! THIS IS FOR THE RIGHT-MOVING CASE; WE IGNORE THE LEFT MOVER.
! CRA
USHR6(I,J) = UMEAN5 - UMEAN1
VSHR6(I,J) = VMEAN5 - VMEAN1
DENOM = USHR6(I,J)*USHR6(I,J)+VSHR6(I,J)*VSHR6(I,J)
IF (DENOM .NE. 0.0) THEN
UST(I,J) = UMEAN6 + (7.5*VSHR6(I,J)/SQRT(DENOM))
VST(I,J) = VMEAN6 - (7.5*USHR6(I,J)/SQRT(DENOM))
ELSE
UST(I,J) = 0
VST(I,J) = 0
ENDIF
ELSE
UST(I,J) = 0.0
VST(I,J) = 0.0
USHR6(I,J) = 0.0
VSHR6(I,J) = 0.0
ENDIF
IF(L1(I,J) > 0 .AND. L2(I,J) > 0) THEN
UMEAN(I,J) = U1(I,J) + (D1000 - HGT1(I,J))*(U2(I,J) - &
U1(I,J))/(HGT2(I,J) - HGT1(I,J))
VMEAN(I,J) = V1(I,J) + (D1000 - HGT1(I,J))*(V2(I,J) - &
V1(I,J))/(HGT2(I,J) - HGT1(I,J))
ELSE IF(L1(I,J) > 0 .AND. L2(I,J) == 0) THEN
UMEAN(I,J) = U1(I,J)
VMEAN(I,J) = V1(I,J)
ELSE IF(L1(I,J) == 0 .AND. L2(I,J) > 0) THEN
UMEAN(I,J) = U2(I,J)
VMEAN(I,J) = U2(I,J)
ELSE
UMEAN(I,J) = 0.0
VMEAN(I,J) = 0.0
ENDIF
IF(L1(I,J) > 0 .OR. L2(I,J) > 0) THEN
USHR1(I,J) = UMEAN(I,J) - U10(I,J)
VSHR1(I,J) = VMEAN(I,J) - V10(I,J)
ELSE
USHR1(I,J) = 0.0
VSHR1(I,J) = 0.0
ENDIF
! CRA
!tgs USHR = UMEAN5 - UMEAN1
! VSHR = VMEAN5 - VMEAN1
! UST(I,J) = UMEAN6 + (7.5*VSHR/SQRT(USHR*USHR+VSHR*VSHR))
! VST(I,J) = VMEAN6 - (7.5*USHR/SQRT(USHR*USHR+VSHR*VSHR))
! ELSE
! UST(I,J) = 0.0
! VST(I,J) = 0.0
! ENDIF
ENDDO
ENDDO
!
! COMPUTE STORM-RELATIVE HELICITY
!
!!$omp parallel do private(i,j,n,l,du1,du2,dv1,dv2,dz,dz1,dz2,dzabv,ie,iw,jn,js,z1,z2,z3)
DO N=1,2 ! for dfferent helicity depth
DO L = 2,LM-1
if(GRIDTYPE /= 'A')then
call exch(ZINT(1,jsta_2l,L))
call exch(ZINT(1,jsta_2l,L+1))
end if
DO J=JSTART,JSTOP
DO I=ISTART,ISTOP
IW=I+IVW(J)
IE=I+IVE(J)
JN=J+JVN
JS=J+JVS
IF (gridtype=='B')THEN
Z2=0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ &
ZMID(I,JN,L)+ZMID(IE,JN,L))
ELSE
Z2=0.25*(ZMID(IW,J,L)+ZMID(IE,J,L)+ &
ZMID(I,JN,L)+ZMID(I,JS,L))
END IF
DZABV=Z2-HTSFC(I,J)
!
IF(DZABV < DEPTH(N) .AND. L <= NINT(LMV(I,J)))THEN
IF (gridtype=='B')THEN
Z1 = 0.25*(ZMID(IW,J,L+1)+ZMID(IE,J,L+1)+ &
ZMID(I,JN,L+1)+ZMID(IE,JN,L+1))
Z3 = 0.25*(ZMID(IW,J,L-1)+ZMID(IE,J,L-1)+ &
ZMID(I,JN,L-1)+ZMID(IE,JN,L-1))
DZ = 0.25*((ZINT(IW,J,L)+ZINT(IE,J,L)+ &
ZINT(I,JN,L)+ZINT(IE,JN,L))- &
(ZINT(IW,J,L+1)+ZINT(IE,J,L+1)+ &
ZINT(I,JN,L+1)+ZINT(IE,JN,L+1)))
ELSE
Z1 = 0.25*(ZMID(IW,J,L+1)+ZMID(IE,J,L+1)+ &
ZMID(I,JN,L+1)+ZMID(I,JS,L+1))
Z3 = 0.25*(ZMID(IW,J,L-1)+ZMID(IE,J,L-1)+ &
ZMID(I,JN,L-1)+ZMID(I,JS,L-1))
DZ = 0.25*((ZINT(IW,J,L)+ZINT(IE,J,L)+ &
ZINT(I,JS,L)+ZINT(I,JN,L))- &
(ZINT(IW,J,L+1)+ZINT(IE,J,L+1)+ &
ZINT(I,JS,L+1)+ZINT(I,JN,L+1)))
END IF
DZ1 = Z1-Z2
DZ2 = Z2-Z3
DU1 = UH(I,J,L+1)-UH(I,J,L)
DU2 = UH(I,J,L)-UH(I,J,L-1)
DV1 = VH(I,J,L+1)-VH(I,J,L)
DV2 = VH(I,J,L)-VH(I,J,L-1)
HELI(I,J,N) = ((VH(I,J,L)-VST(I,J))* &
(DZ2*(DU1/DZ1)+DZ1*(DU2/DZ2)) &
- (UH(I,J,L)-UST(I,J))* &
(DZ2*(DV1/DZ1)+DZ1*(DV2/DZ2))) &
*DZ/(DZ1+DZ2)+HELI(I,J,N)
! if(i==im/2.and.j==(jsta+jend)/2)print*,'Debug Helicity',depth(N),l,dz1,dz2,du1, &
! du2,dv1,dv2,ust(i,j),vst(i,j)
ENDIF
ENDDO
ENDDO
ENDDO
END DO ! end of different helicity depth
!
! END OF ROUTINE.
!
RETURN
END