! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE HIS2VERGRID ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE HIS2VERGRID(nx,ny,nz,nzsoil,x,y,z,zp,zpsoil,hterain, & 1,23 uprt,vprt,wprt,ptprt,pprt,qvprt, & qc,qr,qi,qs,qh,tke,kmh,kmv, & ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, & nstyps,soiltyp,stypfrct,vegtyp,lai,roufns, & veg,tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate,radfrc,radsw,rnflx, & radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & hinfmt,fgrdbasfn,hisfile,nhisfile, & mRefopt,mReflist,mRefdir,mRefrunname, & refopt,refcomp,nreflvl,mxscorelvl,threshold) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Performs read in the ARPS' history dump data and ! the grided observations, and then calclulate ! ETS and bias ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Hu ! ! Original Coding: 01/29/2006 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! DATA ARRAYS READ IN: ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp z coordinate of grid points in physical space (m) ! zpsoil z coordinate of soil model ! ! uprt Perturbation x component of velocity (m/s) ! vprt Perturbation y component of velocity (m/s) ! wprt Perturbation z component of velocity (m/s) ! ! ptprt Perturbation potential temperature (K) ! pprt Perturbation pressure (Pascal) ! ! qc Cloud water mixing ratio (kg/kg) ! qr Rainwater mixing ratio (kg/kg) ! qi Cloud ice mixing ratio (kg/kg) ! qs Snow mixing ratio (kg/kg) ! qh Hail mixing ratio (kg/kg) ! ! tke Turbulent Kinetic Energy ((m/s)**2) ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! wbar Base state z velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil soil temperature (K) ! qsoil Soil moisture ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation, SWin - SWout ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2 * s )) ! qvsflx Surface moisture flux of (kg/(m**2 * s)) ! ! CALCULATED DATA ARRAYS: ! ! su Sounding x component of velocity (m/s) ! sv Sounding y component of velocity (m/s) ! stheta Sounding potential temperature (K) ! spres Sounding pressure (mb) ! stemp Sounding temperature (C) ! sdewp Sounding dew-point (C) ! sdrct Sounding wind direction (degrees) ! ssped Sounding wind speed (m/s) ! shght Sounding height (m) ! ! WORK ARRAYS: ! ! tem1 Temporary work array. ! ! Temporary arrays are defined and used differently by each ! subroutine. ! !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! IMPLICIT NONE ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'vericst.inc' INCLUDE 'mp.inc' ! !----------------------------------------------------------------------- ! ! Arrays to be read in: ! !----------------------------------------------------------------------- ! INTEGER :: nx, ny, nz, nzsoil REAL :: x (nx) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL :: y (ny) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL :: z (nz) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL :: zpsoil(nx,ny,nzsoil) ! The physical height coordinate of soil model ! REAL :: hterain (nx,ny) ! Terrain height ! REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor specific ! humidity (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar (nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific INTEGER :: nstyps INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type INTEGER :: vegtyp (nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulative precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s) INTEGER :: nhisfile ! ! !----------------------------------------------------------------------- ! ! Map variables, declared here for use in subroutines... ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: xs(:) REAL, ALLOCATABLE :: ys(:) REAL, ALLOCATABLE :: xmap(:) REAL, ALLOCATABLE :: ymap(:) REAL, ALLOCATABLE :: latgr(:,:) REAL, ALLOCATABLE :: longr(:,:) ! !----------------------------------------------------------------------- ! ! Work Arrays ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: tem1(:,:,:) REAL, ALLOCATABLE :: tem2(:,:,:) REAL, ALLOCATABLE :: tem3(:,:,:) REAL, ALLOCATABLE :: tem4(:,:) ! REAL, ALLOCATABLE :: mosaic3d(:,:,:) REAL, ALLOCATABLE :: threatscore(:,:),biasscore(:,:) ! !----------------------------------------------------------------------- ! ! options ! !----------------------------------------------------------------------- ! INTEGER :: mRefopt CHARACTER (LEN=256) :: mReflist CHARACTER (LEN=256) :: mRefdir CHARACTER (LEN=256) :: mRefrunname INTEGER :: refopt INTEGER :: refcomp INTEGER :: nreflvl INTEGER :: mxscorelvl REAL :: threshold(20) !----------------------------------------------------------------------- ! ! Misc. internal variables ! !----------------------------------------------------------------------- ! INTEGER :: len1, istatus REAL :: time INTEGER :: i,j,k INTEGER :: ireturn,lengbf,lenfil,nchin,lengbf00Z,lengbf12Z INTEGER :: hinfmt INTEGER :: nhisfile_max PARAMETER (nhisfile_max=200) CHARACTER (LEN=256) :: filename CHARACTER (LEN=256) :: fgrdbasfn CHARACTER (LEN=256) :: grdbasfn CHARACTER (LEN=256) :: hisfile(nhisfile_max) CHARACTER (LEN=256) :: ftemp CHARACTER (LEN=256) :: fmosaic(nhisfile_max) ! CHARACTER (LEN=256) :: refdmpfn LOGICAL :: needmRef INTEGER :: nf CHARACTER (LEN=8) :: the_date CHARACTER (LEN=10) :: the_time CHARACTER (LEN=5) :: the_zone INTEGER :: the_values(8) INTEGER :: sd_id,istat INTEGER :: nxm,nym,nzm INTEGER :: nchanl ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! CALL DATE_AND_TIME(the_date,the_time,the_zone,the_values) ! DO nf=1,nhisfile ! lenfil =256 ! CALL slength( hisfile(nf), lenfil) ! WRITE(6,'(1x,a,a)') 'History file is ',hisfile(nf)(1:lenfil) ! END DO OPEN(12, file=trim(mReflist)) DO nf=1,nhisfile read(12,*) fmosaic(nf) write(*,*) 'the mosaic file is:',trim(fmosaic(nf)) enddo close(12) ! nstyp = nstyps ! ALLOCATE ! ALLOCATE( mosaic3d(nx,ny,nz), STAT=istatus ) CALL check_alloc_status(istatus, "mosaic3d") mosaic3d=0.0 ALLOCATE( threatscore(mxscorelvl,nhisfile), STAT=istatus ) CALL check_alloc_status(istatus, "threatscore") threatscore=0.0 ALLOCATE( biasscore(mxscorelvl,nhisfile), STAT=istatus ) CALL check_alloc_status(istatus, "biasscore") biasscore=0.0 ALLOCATE(xs (nx),STAT=istatus) CALL check_alloc_status(istatus, "arps:xs") xs = 0.0 ALLOCATE(ys (ny),STAT=istatus) CALL check_alloc_status(istatus, "arps:ys") ys = 0.0 ALLOCATE(tem1 (nx,ny,nz),STAT=istatus) CALL check_alloc_status(istatus, "arps:tem1") tem1 = 0.0 ALLOCATE(tem2 (nx,ny,nz),STAT=istatus) CALL check_alloc_status(istatus, "arps:tem2") tem2 = 0.0 ALLOCATE(tem3 (nx,ny,nz),STAT=istatus) CALL check_alloc_status(istatus, "arps:tem3") tem3 = 0.0 ALLOCATE(tem4 (nx,ny),STAT=istatus) CALL check_alloc_status(istatus, "arps:tem4") tem4 = 0.0 ALLOCATE(xmap (nx),STAT=istatus) CALL check_alloc_status(istatus, "arps:xmap") xmap = 0.0 ALLOCATE(ymap (ny),STAT=istatus) CALL check_alloc_status(istatus, "arps:ymap") ymap = 0.0 ALLOCATE(latgr (nx,ny),STAT=istatus) CALL check_alloc_status(istatus, "arps:latgr") latgr = 0.0 ALLOCATE(longr (nx,ny),STAT=istatus) CALL check_alloc_status(istatus, "arps:longr") longr = 0.0 !----------------------------------------------------------------------- ! ! Loop through all history dumps ! !----------------------------------------------------------------------- ! grdbasfn = fgrdbasfn lengbf=256 CALL slength(grdbasfn,lengbf) grdbasfn(1:lengbf)=grdbasfn(1:lengbf) IF ( mp_opt > 0 .AND. readsplit <= 0 ) THEN WRITE(grdbasfn(lengbf+1:lengbf+5),'(a,i2.2,i2.2)') '_', & loc_x,loc_y lengbf = lengbf + 5 END IF DO nf = 1, nhisfile ! !----------------------------------------------------------------------- ! ! Read all input data arrays ! !----------------------------------------------------------------------- ! filename=hisfile(nf) lenfil =256 CALL slength( hisfile(nf), lenfil) IF ( mp_opt > 0 .AND. readsplit <= 0 ) THEN WRITE(filename(lenfil+1:lenfil+5),'(a,i2.2,i2.2)') '_', & loc_x,loc_y lenfil = lenfil + 5 END IF IF ( myproc == 0 ) & WRITE(6,'(1x,a,a)') 'History file is ',filename(1:lenfil) CALL dtaread(nx,ny,nz,nzsoil,nstyps, & hinfmt,nchin,grdbasfn(1:lengbf),lengbf, & filename(1:lenfil),lenfil,time, & x,y,z,zp,zpsoil,uprt,vprt,wprt,ptprt,pprt, & qvprt,qc,qr,qi,qs,qh,tke,kmh,kmv, & ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn, tem1,tem2,tem3) DO j=1,ny DO i=1,nx hterain(i,j)=zp(i,j,2) END DO END DO ! !----------------------------------------------------------------------- ! ! ireturn = 0 for a successful read ! !----------------------------------------------------------------------- ! IF( ireturn == 0 ) THEN ! successful read ! !----------------------------------------------------------------------- ! get predicted reflectivity ! !----------------------------------------------------------------------- IF( mRefopt == 1 ) then ! ! Calculate temperature (K) at ARPS grid points ! DO k=1,nz DO j=1,ny DO i=1,nx tem2(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k) END DO END DO END DO CALL temper (nx,ny,nz,tem2, pprt ,pbar,tem1) ! ! Calculate reflectivity according to precipitation mixing ratio ! tem2 = 3d reflectiivty ! tem3 = composite reflectivity ! tem2=0 tem3=0 IF (refopt == 2) THEN CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tem1, tem2) ELSE CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem2) ENDIF if(refcomp == 1) CALL cal_rfc(nx, ny, nz, tem2, tem3) ENDIF ! !----------------------------------------------------------------------- ! ! get reflectivity observation ! !----------------------------------------------------------------------- ! IF( mRefopt == 1 ) then ftemp=trim(mRefdir)//trim(fmosaic(nf)) WRITE(6,'(a,a)') ' Opening radar file ',trim(ftemp) if(2==1) then nchanl=38 OPEN(UNIT=nchanl,FILE=trim(ftemp),ERR=400, & FORM='unformatted',STATUS='old') ! READ(nchanl) nxm,nym,nzm if( (nxm .ne. nx) .or. (nym .ne. ny) .or. (nzm .ne. nz) ) then write(*,*) 'Mosaic has a different dimension from the predictio!' write(*,*) 'nx,ny,nz=',nx,ny,nz write(*,*) 'nxm,nym,nzm=',nxm,nym,nzm stop 123 endif read(nchanl) mosaic3d CLOSE(nchanl) else call readadcol_Mosaic(nx,ny,nz,ftemp,mosaic3d) endif ENDIF ! !----------------------------------------------------------------------- ! ! calculate ETS ! !----------------------------------------------------------------------- IF( mRefopt == 1 ) then call ETS_Ref(nx,ny,nz,mosaic3d,tem3,refcomp,nreflvl, & mxscorelvl,threshold,threatscore(:,nf),biasscore(:,nf)) ENDIF ! ! mRefdir,mRefrunname, mReflist needmRef=.false. IF (needmRef) THEN len1 =256 CALL slength( mRefrunname, len1) refdmpfn=mRefrunname(1:len1)//'.ref'// & filename(lenfil-5:lenfil) END IF ! ! !----------------------------------------------------------------------- ! ! Extract model data for MOS. ! !----------------------------------------------------------------------- ! ELSE WRITE(6,'(a)') ' Error reading data. HIS2VERGRID ends' STOP END IF ! !----------------------------------------------------------------------- ! ! Go to next history dump ! !----------------------------------------------------------------------- ! END DO ! nf !----------------------------------------------------------------------- ! ! save ETS ! !----------------------------------------------------------------------- ! refdmpfn=trim(mRefrunname)//'_ref.ETS' OPEN(12,file=trim(refdmpfn) ) write(12,*) ' Bias Score Equitable threat score ' DO nf=1,nhisfile write(12,'(8f8.4,8f8.4)') & (biasscore(k,nf),k=1,mxscorelvl), & (threatscore(k,nf), k=1,mxscorelvl) ENDDO close(12) RETURN 400 continue write(*,*) 'Error in open Mosaic file!!!!' stop 123 END SUBROUTINE his2vergrid ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE TEMPER ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE temper ( nx,ny,nz,theta, ppert, pbar, t ) 9 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Using a version of Poisson's formula, calculate temperature. ! !----------------------------------------------------------------------- ! ! AUTHOR: Joe Bradley ! 12/05/91 ! ! MODIFICATIONS: ! Modified by Ming Xue so that arrays are only defined at ! one time level. ! 6/09/92 Added full documentation and phycst include file for ! rddcp=Rd/Cp (K. Brewster) ! !----------------------------------------------------------------------- ! ! INPUT: ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! theta Potential temperature (degrees Kelvin) ! ppert Perturbation pressure (Pascals) ! pbar Base state pressure (Pascals) ! ! OUTPUT: ! ! t Temperature (degrees Kelvin) ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE ! INTEGER :: nx,ny,nz ! REAL :: theta(nx,ny,nz) ! potential temperature (degrees Kelvin) REAL :: ppert(nx,ny,nz) ! perturbation pressure (Pascals) REAL :: pbar (nx,ny,nz) ! base state pressure (Pascals) ! REAL :: t (nx,ny,nz) ! temperature (degrees Kelvin) ! !----------------------------------------------------------------------- ! ! Include file ! !----------------------------------------------------------------------- ! INCLUDE 'phycst.inc' ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! Calculate the temperature using Poisson's formula. ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 t(i,j,k) = theta(i,j,k) * & (((ppert(i,j,k) + pbar(i,j,k)) / p0) ** rddcp) END DO END DO END DO RETURN END SUBROUTINE temper ! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINES cal_rfc ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! ! PURPOSE: ! Calculate rfc value. ! !----------------------------------------------------------------------- ! ! AUTHOR: Min Zou ! 3/2/98 ! ! MODIFICATION HISTORY: ! ! Ming Xue (10/16/2001) ! Now passing in precalculated reflectivity field instead of calculating ! it inside. ! !----------------------------------------------------------------------- ! SUBROUTINE cal_rfc(nx, ny, nz, ref, refc) 3 IMPLICIT NONE INTEGER :: nx,ny,nz REAL, INTENT(IN ) :: ref (nx,ny,nz) ! Reflectivity REAL, INTENT(OUT) :: refc(nx,ny,nz) ! Composite reflectivity INTEGER :: i,j,k ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO j=1,ny DO i=1,nx refc(i,j,1)= ref(i,j,1) DO k=2,nz-1 refc(i,j,1) = MAX(refc(i,j,1),ref(i,j,k)) END DO END DO END DO DO j=1,ny DO i=1,nx DO k=2,nz-1 refc(i,j,k) = refc(i,j,1) END DO END DO END DO RETURN END SUBROUTINE cal_rfc ! !################################################################## !################################################################## !###### ###### !###### PROGRAM ARPSDAS ###### !###### ARPS Data Analysis System ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE ETS_Ref(nx,ny,nz,mosaic3d,fcst3d, & 1 refcomp,nreflvl,mxscorelvl,lvls, & threatscore,biasscore) ! ! PURPOSE: calculate ETS of predicted reflectivity against observations ! ! ! AUTHOR: ! ! Ming Hu, CAPS, May, 2006 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! IMPLICIT NONE ! INTEGER :: refcomp,nreflvl ! ! predicted and observed reflectiivty ! REAL :: mosaic3d(nx,ny,nz) REAL :: fcst3d(nx,ny,nz) ! !----------------------------------------------------------------------- ! ! Temporary work arrays ! !----------------------------------------------------------------------- ! REAL :: tem2d(nx,ny) REAL :: fcst2d(nx,ny) !----------------------------------------------------------------------- ! ! for ETS calculation !----------------------------------------------------------------------- ! ! INTEGER :: mxscorelvl REAL :: hits(mxscorelvl), misses(mxscorelvl) REAL :: falsealarms(mxscorelvl), corneg(mxscorelvl) REAL :: forecastpoints(mxscorelvl), observationpoints(mxscorelvl) REAL :: hitsrandom REAL :: lvls(20) REAL :: threatscore(mxscorelvl),biasscore(mxscorelvl) ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: nx,ny,nz INTEGER :: nxf,nyf,nzf ! integer :: num integer :: i,j,k ,ilng,nchanl real :: rmax ! !----------------------------------------------------------------------- ! if( mxscorelvl > 20 ) then write(*,*) 'Too many threshold for statistic' stop 123 endif IF(refcomp == 1 ) then tem2d=0.0 DO j=1,ny DO i=1,nx rmax=-9999.9 DO k=2,nz-1 if( mosaic3d(i,j,k) > rmax) rmax=mosaic3d(i,j,k) ENDDO tem2d(i,j) = rmax ENDDO ENDDO fcst2d=0.0 DO j=1,ny DO i=1,nx rmax=-9999.9 DO k=2,nz-1 if( fcst3d(i,j,k) > rmax) rmax=fcst3d(i,j,k) ENDDO fcst2d(i,j) = rmax ENDDO ENDDO ELSEIF( refcomp == 0 ) then DO j=1,ny DO i=1,nx tem2d(i,j) = mosaic3d(i,j,nreflvl) fcst2d(i,j) = fcst3d(i,j,nreflvl) ENDDO ENDDO ELSE write(*,*) ' Unknown choice for refcomp' write(*,*) ' refcomp should be 0 or 1' stop 123 ENDIF ! Equitable threat score and bias score num=0 hits=0 misses=0 falsealarms=0 corneg=0 forecastpoints=0 observationpoints=0 DO i=2, nx-1 DO j=2, ny-1 num=num+1 DO k=1,mxscorelvl IF(fcst2d(i,j) >= lvls(k) .and. tem2d(i,j) >= lvls(k) ) THEN hits(k)=hits(k) + 1 ELSE IF ( fcst2d(i,j) < lvls(k) .and. tem2d(i,j) >= lvls(k) ) THEN misses(k)=misses(k)+1 ELSE IF ( fcst2d(i,j) >= lvls(k) .and. tem2d(i,j) < lvls(k) ) THEN falsealarms(k)=falsealarms(k)+1 ELSE corneg(k)=corneg(k)+1 END IF IF( fcst2d(i,j) >= lvls(k) ) forecastpoints(k) = forecastpoints(k) + 1 IF( tem2d(i,j) >= lvls(k) ) observationpoints(k) = observationpoints(k) + 1 ENDDO ENDDO ENDDO IF( num > 0 ) THEN write(*,*) ' The total number parting statistic is:',num DO k=1,mxscorelvl hitsrandom= forecastpoints(k) * observationpoints(k) / float(num) threatscore(k)=(hits(k)-hitsrandom)/(hits(k)+misses(k)+falsealarms(k)-hitsrandom) biasscore(k)=(hits(k)+falsealarms(k))/(hits(k)+misses(k)) ENDDO ELSE write(*,*) ' Error: no data were found in statistic !' stop 234 ENDIF ! write(*,*) ' Bias Score Equitable threat score ' ! write(*,'(f10.2,8f8.4,8f8.4)') (biasscore(k),k=1,mxscorelvl), & ! (threatscore(k), k=1,mxscorelvl) return 400 write(*,*) ' Error in read in Radar file ! ' stop 8888 102 write(*,*) ' Error in read in radar data !' END Subroutine ETS_Ref ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE READADCOL ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE readadcol_Mosaic(nx,ny,nz,rfname,gridref) 2,8 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! read in gridded radar data to a file as columns with ! individual lat,lons. ! !----------------------------------------------------------------------- ! ! 01/28/06 MIng HU ! Modified to fit the need to write NSSL mosaic reflectiivty ! !----------------------------------------------------------------------- ! ! INPUT: ! dmpfmt file format (1:binary, 2:hdf) ! iradfmt binary format ! hdf4cmpr hdf4 compression level ! rfname radar file name (character*80) ! iyr year ! imon month ! iday day ! ihr hour ! imin min ! isec sec ! isource) source number ! 1= WSR-88D raw ! 2= WSR-88D NIDS ! 3= NSSL mosaic reflectivity ! ! OUTPUT: ! data are written to file ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'grid.inc' ! INTEGER :: nx,ny,nz ! INTEGER :: dmpfmt INTEGER :: iradfmt INTEGER :: hdf4cmpr CHARACTER (LEN=256) :: rfname INTEGER :: iyr,imon,iday,ihr,imin,isec INTEGER :: isource ! REAL :: zs(nz) REAL :: zpsc(nx,ny,nz) REAL :: gridref(nx,ny,nz) ! REAL :: readk(nz) REAL :: readhgt(nz) REAL :: readref(nz) ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! ! Radar output descriptors ! !----------------------------------------------------------------------- ! ! INTEGER :: mxradvr,nradvr ! PARAMETER(mxradvr=10,nradvr=6) ! INTEGER :: iradvr(mxradvr) ! DATA iradvr /1,2,3,4,5,6,0,0,0,0/ ! !----------------------------------------------------------------------- ! ! Radar output thresholds ! !----------------------------------------------------------------------- ! REAL :: refmin,refmax,velmin,velmax PARAMETER(refmin=-5.0, refmax=100., & velmin=-200.,velmax=200.) REAL :: misval PARAMETER(misval=-999.0) ! !----------------------------------------------------------------------- ! ! Radar output variables ! !----------------------------------------------------------------------- ! REAL :: grdlatc(nx,ny) REAL :: grdlonc(nx,ny) ! !----------------------------------------------------------------------- ! ! Misc local variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=4) :: radid CHARACTER (LEN=12) :: runname INTEGER :: iunit,myr,itime INTEGER :: i,j,k,klev,kk,kntcol,nn INTEGER :: idummy INTEGER :: istat,sd_id INTEGER :: ipt REAL :: gridlat,gridlon,elev,rdummy INTEGER(2), allocatable :: itmp(:,:,:) ! Temporary array REAL, allocatable :: hmax(:), hmin(:) ! Temporary array ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! print *, ' inside readadcol_mosaic' idummy=-999 rdummy=-999. dmpfmt=1 iradfmt=1 hdf4cmpr=1 ! IF (dmpfmt > 1 .AND. hdf4cmpr > 3) THEN write(*,*) 'HDF is not available!' stop 123 END IF IF(dmpfmt == 1)THEN CALL getunit(iunit) ! !----------------------------------------------------------------------- ! ! Open file for output ! !----------------------------------------------------------------------- ! write(*,*) rfname OPEN(iunit,FILE=TRIM(rfname),STATUS='UNKNOWN',FORM='UNFORMATTED') ! !----------------------------------------------------------------------- ! ! Write radar description variables ! !----------------------------------------------------------------------- ! read(iunit) radid read(iunit) idummy,itime,idummy,isource,idummy, & idummy,idummy,idummy,idummy,idummy ! CALL abss2ctim(itime,myr,imon,iday,ihr,imin,isec) ! write(*,*) 'The time of the data is: ',myr,imon,iday,ihr,imin,isec !----------------------------------------------------------------------- ! ! Write grid description variables ! This should provide enough info to verify that the ! proper grid has been chosen. To recreate the grid, ! icluding elevation information, ! the reading program should get a grid-base-file ! named runname.grdbasfil ! !----------------------------------------------------------------------- ! idummy=0 rdummy=0. read(iunit) runname read(iunit) iradfmt,strhopt,mapproj,idummy,idummy, & idummy,idummy,idummy,idummy,idummy read(iunit) dx,dy,dz,dzmin,ctrlat, & ctrlon,trulat1,trulat2,trulon,sclfct, & rdummy,rdummy,rdummy,rdummy,rdummy read(iunit) idummy,idummy ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! !----------------------------------------------------------------------- ! ! For each horizontal grid point form a column of remapped ! data containing the non-missing grid points ! !----------------------------------------------------------------------- ! IF(dmpfmt==1)THEN DO ipt=1,(nx*ny) read(iunit,END=51) i,j,rdummy,rdummy, & gridlat,gridlon,elev,klev read(iunit,END=51) (readk(k),k=1,klev) read(iunit,END=51) (readhgt(k),k=1,klev) read(iunit,END=51) (readref(k),k=1,klev) IF(i <= nx.AND.i >= 1 .AND. j <= ny.AND.j >= 1) THEN DO kk=1,klev k=nint(readk(kk)) IF(k <= nz.AND.k >= 1) THEN gridref(i,j,k)=readref(kk) END IF ! 1 < k < nz END DO ! kk = 1, klev END IF ! 1 < i < nx & 1 < j < ny END DO ! ipt = 1, nx*ny 51 continue ipt=ipt-1 WRITE(6,'(a,i6,a)') ' End of file reached after reading', & ipt,' columns' ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! RETURN END SUBROUTINE readadcol_Mosaic