! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE CALCSHR ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE calcshr(nrad,lm,zp,sigma, & p_pa,t3d,u3d,v3d,cape, & shr37,ustrm,vstrm,srlfl,srmfl,helicity,brn,brnu,brnu2km,blcon, & tem2,tem3,nlev,estacao,hora,dia,mes,ano,strm_opt) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate various wind shear parameters useful for gauging ! the potential for severe storms. ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Keith Brewster ! February, 1994 Based on OLAPS, hence LAPS, version of same. ! ! MODIFICATION HISTORY: ! 3/11/1996 Keith Brewster ! Added storm-relative flows, general clean-up to ! meet ARPS coding standards. ! ! 3/22/1996 (Keith Brewster) ! Fixed some bugs, added smoothing at the end. ! ! 06/20/2000 (Eric Kemp and Keith Brewster) ! Changed BRN Shear to be the denominator of BRN, instead of wind ! speed, now has units of speed squared. ! ! 04/19/2004 (Ernani L. Nascimento) ! Modified the code to fit in índices_severos3.f90, and adapted ! the computation of expected storm-motion for the Southern Hemisphere. ! !----------------------------------------------------------------------- ! ! Calculates some of the shear related variables from the ! ARPS 3D wind field. ! ! shr37 Magnitude of wind shear between 3 and 7 km AGL ! ustrm,vstrm Estimated storm motion (modified from Bob Johns) ! srlfl Low-level storm-relative wind ! srmfl Mid-level storm-relative wind ! helicity Helicity, storm relative ! brn Bulk Richardson Number (Weisman and Klemp) ! brnu Shear parameter of BRN, "U" ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nrad,lm,strm_opt ! !----------------------------------------------------------------------- ! ! "1-D" input variables (Ernani L. Nascimento) ! !----------------------------------------------------------------------- ! INTEGER :: nlev(nrad),estacao(nrad),hora(nrad),dia(nrad),mes(nrad), & ano(nrad) REAL :: sigma(nrad) REAL :: cape(nrad) ! !----------------------------------------------------------------------- ! ! "2-D" input variables (Ernani L. Nascimento) ! !----------------------------------------------------------------------- ! INTEGER :: zp(lm,nrad) REAL :: p_pa(lm,nrad) ! Pressure in Pascals REAL :: t3d(lm,nrad) ! Temperature in Kelvin REAL :: u3d(lm,nrad) REAL :: v3d(lm,nrad) ! !----------------------------------------------------------------------- ! ! Output variables ! !----------------------------------------------------------------------- ! REAL :: shr37(nrad) ! 7km - 3km wind shear REAL :: ustrm(nrad) REAL :: vstrm(nrad) ! Estimated storm motion (Bob Johns) REAL :: srlfl(nrad) REAL :: srmfl(nrad) REAL :: helicity(nrad) ! Helicity, storm relative REAL :: brn(nrad) ! Bulk Richardson Number (Weisman and Klemp) REAL :: brnu(nrad) ! Shear parameter of BRN, "U" REAL :: brnu2km(nrad) ! 2km shear parameter of BRN REAL :: blcon(nrad) ! !----------------------------------------------------------------------- ! ! Temporary variables ! !----------------------------------------------------------------------- ! REAL :: tem2(nrad) REAL :: tem3(nrad) ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! ! INTEGER :: imid,jmid REAL :: elev(nrad),hgt3d(lm,nrad) INTEGER :: i,j,k,ksfc,k2km,k3km REAL :: h3km,u3km,v3km,h7km,u7km,v7km,u2,v2 REAL :: p2km,t2km,h2km,u2km,v2km,p9km,t9km,h9km,u9km,v9km REAL :: p500m,t500m,h500m,u500m,v500m,p6km,t6km,h6km,u6km,v6km REAL :: sumu,sumv,sump,sumh,wlow,whigh,dx,dy,dz,dp,dx2,dy2 REAL :: rhohi,rholo,rhoinv,arg,new_ustrm,new_vstrm REAL :: dirmean,spmean,ushr,vshr,ddir,perc,obs_ddstrm,obs_ffstrm ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! Loop over all soundings (Ernani L. Nascimento) !----------------------------------------------------------------------- ! ! write(*,*) 'Entrada da subrotina calcshr' DO k=1,nrad if (nlev(k)==0) then shr37(k)=-999.9 ustrm(k)=-999.9 vstrm(k)=-999.9 srlfl(k)=-999.9 srmfl(k)=-999.9 helicity(k)=-999.9 brn(k)=-999.9 brnu(k)=-999.9 brnu2km(k)=-999.9 blcon(k)=-999.9 else ! print*, estacao(k),hora(k),dia(k),mes(k),ano(k) j=nlev(k) if ((p_pa(j,k)*0.01)>300.0) then cycle endif arg=0.0 elev(k)=real(zp(1,k)) ! write(*,*) 'elev(k)= ',elev(k),' para k= ',k DO j=1,nlev(k) hgt3d(j,k)=real(zp(j,k)) ENDDO ! write(*,*) 'PASSEI AQUI 3, com k= ',k ! !----------------------------------------------------------------------- ! ! Find mass weighted mean wind in first 500m ! !----------------------------------------------------------------------- ! sumu=0. sumv=0. sump=0. h500m=elev(k)+500. DO j=2,nlev(k) IF( hgt3d(j,k) < h500m ) THEN dp=p_pa(j-1,k)-p_pa(j,k) rhohi=p_pa(j,k)/t3d(j,k) rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k)) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k)) sump=sump+dp ELSE dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h500m)/dz u500m=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v500m=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow p500m=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow t500m=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow dp=p_pa(j-1,k)-p500m rhohi=p500m/t500m rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u500m) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v500m) sump=sump+dp ! print *, ' sumu,sumv,sump = ',sumu,sumv,sump,'para k= ',k EXIT END IF END DO ! 121 CONTINUE u500m=sumu/sump v500m=sumv/sump ! write(*,*) 'PASSEI AQUI 4, com k= ',k ! write(*,*) '(u500m,v500m)= ',u500m,v500m,' para k= ',k ! !----------------------------------------------------------------------- ! ! Find mass weighted mean wind sfc-2km AGL ! !----------------------------------------------------------------------- ! sumu=0. sumv=0. sump=0. h2km=elev(k)+2000. DO j=2,nlev(k) IF( hgt3d(j,k) < h2km ) THEN dp=p_pa(j-1,k)-p_pa(j,k) rhohi=p_pa(j,k)/t3d(j,k) rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k)) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k)) sump=sump+dp ELSE dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h2km)/dz u2km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v2km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow p2km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow t2km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow dp=p_pa(j-1,k)-p2km rhohi=p2km/t2km rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u2km) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v2km) sump=sump+dp EXIT END IF END DO ! 141 CONTINUE ! write(*,*) 'PASSEI AQUI 4.5, com k= ',k u2km=sumu/sump v2km=sumv/sump ! write(*,*) 'u2km,v2km: ',u2km,v2km tem2(k)=u2km tem3(k)=v2km ! write(*,*) 'tem2(k),tem3(k): ',tem2(k),tem3(k) ! write(*,*) 'PASSEI AQUI 5, com k= ',k ! !----------------------------------------------------------------------- ! ! Find mass weighted mean wind sfc-6km AGL ! !----------------------------------------------------------------------- ! sumu=0. sumv=0. sump=0. h6km=elev(k)+6000. DO j=2,nlev(k) IF( hgt3d(j,k) < h6km ) THEN dp=p_pa(j-1,k)-p_pa(j,k) rhohi=p_pa(j,k)/t3d(j,k) rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k)) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k)) sump=sump+dp ELSE dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h6km)/dz u6km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v6km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow p6km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow t6km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow dp=p_pa(j-1,k)-p6km rhohi=p6km/t6km rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u6km) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v6km) sump=sump+dp EXIT END IF END DO u6km=sumu/sump v6km=sumv/sump ! write(*,*) 'PASSEI AQUI 6, com k= ',k,' e (u6km,v6km)= ',u6km,v6km ! !----------------------------------------------------------------------- ! ! Find mass weighted mean wind 2km-9km AGL ! !----------------------------------------------------------------------- ! sumu=0. sumv=0. sump=0. h9km=elev(k)+9000. DO j=2,nlev(k) IF( hgt3d(j,k) > h2km ) EXIT END DO ! 181 CONTINUE k2km=j dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h2km)/dz u2=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v2=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow p2km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow t2km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow dp=p2km-p_pa(j,k) rholo=p2km/t2km rhohi=p_pa(j,k)/t3d(j,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u2+rhohi*u3d(j,k)) sumv=sumv+dp*rhoinv*(rholo*v2+rhohi*v3d(j,k)) sump=sump+dp DO j=k2km+1,nlev(k) IF( hgt3d(j,k) < h9km ) THEN dp=p_pa(j-1,k)-p_pa(j,k) rhohi=p_pa(j,k)/t3d(j,k) rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k)) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k)) sump=sump+dp ELSE dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h9km)/dz u9km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v9km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow p9km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow t9km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow dp=p_pa(j-1,k)-p9km rhohi=p9km/t9km rholo=p_pa(j-1,k)/t3d(j-1,k) rhoinv=1./(rhohi+rholo) sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u9km) sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v9km) sump=sump+dp EXIT END IF END DO ! 191 CONTINUE u9km=sumu/sump v9km=sumv/sump ! write(*,*) 'PASSEI AQUI 7, com k= ',k ! !----------------------------------------------------------------------- ! ! Storm motion estimation ! From Davies and Johns, 1993 ! "Some wind and instability parameters associated With ! strong and violent tornadoes." ! AGU Monograph 79, The Tornado...(Page 575) ! ! Becuase of the discontinuity produced by that method ! at the 15.5 m/s cutoff, their rules have been modified ! to provide a gradual transition, and accomodate all the ! data they mention in the article. ! ! (04/19/2004) Modified by Ernani L. Nascimento. Adapting for the ! Southern Hemisphere. ! !----------------------------------------------------------------------- ! CALL uv2ddff(u6km,v6km,dirmean,spmean) ! write(*,*) 'PASSAGEM 5.1, com k= ',k ! write(*,*) 'PASSAGEM 5.1, com k= ',k,' e spmean= ',spmean ! write(*,*) 'PASSAGEM 5.1, com k= ',k,' e dirmean= ',dirmean IF(spmean >= 20.0) THEN ! write(*,*) 'Hei, passei aqui!' dirmean=dirmean-18. IF(dirmean <= 0.) dirmean=360.+dirmean spmean=spmean*0.89 ELSE IF (spmean > 8.0) THEN ! write(*,*) 'Não! Passei foi aqui!' whigh=(spmean - 8.0)/12. wlow =1.-whigh ddir=wlow*32.0 + whigh*18.0 perc=wlow*0.75 + whigh*0.89 dirmean=dirmean-ddir IF(dirmean <= 0.) dirmean=360.+dirmean spmean=spmean*perc ELSE ! write(*,*) 'Na verdade, passei foi aqui!' dirmean=dirmean-32. IF(dirmean <= 0.) dirmean=360.+dirmean spmean=spmean*0.75 END IF ! write(*,*) 'PASSAGEM 5.2, com k= ',k arg = (dirmean * (3.141592654/180.)) ustrm(k) = -spmean * sin(arg) vstrm(k) = -spmean * cos(arg) ! Utilizando movimento estimado da célula da esquerda em 9 de outubro de 2003 IF (strm_opt == 1) then IF ((estacao(k)==83827).and.(hora(k)==00).and.(dia(k)==09).and. & (mes(k)==10).and.(ano(k)==2003)) THEN obs_ddstrm=270.0 obs_ffstrm=10.7 CALL ddff2uv(obs_ddstrm,obs_ffstrm,new_ustrm,new_vstrm) ustrm(k)=new_ustrm vstrm(k)=new_vstrm ENDIF ! Utilizando movimento estimado da célula da célula de 24 de maio de 2005 ! IF ((estacao(k)==99999).and.(hora(k)==18).and.(dia(k)==24).and. & ! (mes(k)==05).and.(ano(k)==2005)) THEN ! obs_ddstrm=297.0 ! obs_ffstrm=17.8 ! CALL ddff2uv(obs_ddstrm,obs_ffstrm,new_ustrm,new_vstrm) ! ustrm(k)=new_ustrm ! vstrm(k)=new_vstrm ! ENDIF ENDIF ! CALL ddff2uv(dirmean,spmean,ustrm(k),vstrm(k),5) ! write(*,*) 'PASSAGEM 6, com k= ',k ! !----------------------------------------------------------------------- ! ! Storm-relative low-level flow ! !----------------------------------------------------------------------- ! srlfl(k)=SQRT((ustrm(k)-u2km)*(ustrm(k)-u2km) + & (vstrm(k)-v2km)*(vstrm(k)-v2km)) ! !----------------------------------------------------------------------- ! ! Storm relative mid-level flow ! !----------------------------------------------------------------------- ! srmfl(k)=SQRT((ustrm(k)-u9km)*(ustrm(k)-u9km) + & (vstrm(k)-v9km)*(vstrm(k)-v9km)) ! !----------------------------------------------------------------------- ! ! Shear parameter for Bulk Richardson number ! !----------------------------------------------------------------------- ! ! print *, ' density-weight mean 0-500 m ',u500m,v500m ! print *, ' density-weight mean 0-6 km ',u6km,v6km ! brnu(k)=0.5*( (u6km-u500m)*(u6km-u500m) + & (v6km-v500m)*(v6km-v500m) ) ! Added the 2km BRNSHR (Ernani L. Nascimento, 03/feb/2005 brnu2km(k)=0.5*( (u2km-u500m)*(u2km-u500m) + & (v2km-v500m)*(v2km-v500m) ) ! write(*,*) 'PASSAGEM 7, com k= ',k ! !----------------------------------------------------------------------- ! ! Bulk Richardson number ! A limit of 200 is imposed, since this could ! go to inifinity. ! !----------------------------------------------------------------------- ! IF(brnu(k) > 0.) THEN brn(k)=cape(k)/brnu(k) brn(k)=AMIN1(brn(k),200.) ELSE brn(k)=200. END IF ! write(*,*) 'PASSEI AQUI 8, com k= ',k ! ! !----------------------------------------------------------------------- ! ! Calculate Helicity and 3km to 7km shear ! since both involve the 3km wind. ! ! For more efficient computation the Helicity is ! computed for zero storm motion and the storm ! motion is accounted for by adding a term at the end. ! This is mathematically equivalent to accounting ! for the storm motion at each level. ! !----------------------------------------------------------------------- ! h3km=elev(k)+3000. h7km=elev(k)+7000. ! !----------------------------------------------------------------------- ! ! Find level just above 3 km AGL ! Note, it is assumed here that there is at least ! one level between the sfc and 3 km. ! !----------------------------------------------------------------------- ! sumh=0. DO j=2,nlev(k) IF(hgt3d(j,k) > h3km) EXIT sumh=sumh + & ( u3d(j,k)*v3d(j-1,k) ) - & ( v3d(j,k)*u3d(j-1,k) ) END DO ! 240 CONTINUE k3km=j dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h3km)/dz u3km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v3km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow sumh=sumh + & ( u3km*v3d(j-1,k) ) - & ( v3km*u3d(j-1,k) ) ushr=u3km-u3d(1,k) vshr=v3km-v3d(1,k) helicity(k)=sumh + vshr*ustrm(k) - ushr*vstrm(k) ! write(*,*) 'PASSAGEM 9, com k= ',k ! !----------------------------------------------------------------------- ! ! Now Find 7km wind for 3-to-7km shear ! !----------------------------------------------------------------------- ! DO j=k3km,nlev(k) IF(hgt3d(j,k) > h7km) EXIT END DO ! 260 CONTINUE dz=hgt3d(j,k)-hgt3d(j-1,k) wlow=(hgt3d(j,k)-h7km)/dz u7km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow v7km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow shr37(k)=(SQRT( (u7km-u3km)*(u7km-u3km) + & (v7km-v3km)*(v7km-v3km) ))/4000. endif END DO ! write(*,*) 'Saída da subrotina calcshr.' !----------------------------------------------------------------------- ! ! Calculate the low-level convergence ! !----------------------------------------------------------------------- ! ! dx=(x(2)-x(1)) ! dy=(y(2)-y(1)) ! dx2=2.*dx ! dy2=2.*dy ! DO j=2,ny-2 ! DO i=2,nx-2 ! blcon(i,j)=-1000.*( & ! (tem2(i+1,j)-tem2(i-1,j))/(sigma(i,j)*dx2)+ & ! (tem3(i,j+1)-tem3(i,j-1))/(sigma(i,j)*dy2) ) ! END DO ! END DO ! ! DO j=2,ny-2 ! blcon(1,j)=-1000.*( & ! (tem2(2,j)-tem2(1,j))/(sigma(1,j)*dx)+ & ! (tem3(1,j+1)-tem3(1,j-1))/(sigma(1,j)*dy2) ) ! blcon(nx-1,j)=-1000.*( & ! (tem2(nx-1,j)-tem2(nx-2,j))/(sigma(nx-1,j)*dx)+ & ! (tem3(nx-1,j+1)-tem3(nx-1,j-1))/(sigma(nx-1,j)*dy2) ) ! END DO ! DO i=2,nx-2 ! blcon(i,1)=-1000.*( & ! (tem2(i+1,1)-tem2(i-1,1))/(sigma(i,1)*dx2)+ & ! (tem3(i,2)-tem3(i,1))/(sigma(i,1)*dy) ) ! blcon(i,ny-1)=-1000.*( & ! (tem2(i+1,ny-1)-tem2(i-1,ny-1))/(sigma(i,ny-1)*dx2)+ & ! (tem3(i,ny-1)-tem3(i,ny-2))/(sigma(i,ny-1)*dy) ) ! END DO ! blcon(1,1)=blcon(2,2) ! blcon(nx-1,ny-1)=blcon(nx-2,ny-2) ! !----------------------------------------------------------------------- ! ! Smooth the output arrays ! !----------------------------------------------------------------------- ! ! CALL smooth9p( shr37,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( ustrm,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( vstrm,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( srlfl,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( srmfl,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p(helicity,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( brn,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( brnu,nx,ny,1,nx-1,1,ny-1,tem1) ! CALL smooth9p( blcon,nx,ny,1,nx-1,1,ny-1,tem1) ! RETURN END SUBROUTINE calcshr