      subroutine zhatus(iseas,ilag,nlag,ilagf,rdir,
     +zreal,zpers,zpred,QKT2,nyear,ndnd,nled,nmod,
     +acgeoc,acgeop,accca,acper)
C
      parameter(LEDBEG=1,LEDEND=1)
C
      dimension zpers(nyear,ndnd),zreal(nyear,nled,ndnd)
      dimension zpred(nyear,nled,ndnd)
      DIMENSION ZZ(ndnd,nmod)
      DIMENSION S(nled,ndnd)
      DIMENSION LEADX(nled),Z(nyear,nled,ndnd)
      DIMENSION VRBL1(nyear),VRBL2(nyear),NYER(nled)
      DIMENSION RMSER(nled,ndnd),RMSER2(nled,ndnd)
      DIMENSION SCORM(nled,ndnd),VRBL3(nyear)
      DIMENSION SCORME(nled,ndnd),CORSKL(nled,ndnd)
      DIMENSION ERBR(nled,ndnd),ERBR2(nled,ndnd)
      DIMENSION SCORCR(nled,ndnd),CORLAG(nled,ndnd)
      DIMENSION QKT2(nled,nyear)
      DIMENSION Q(nyear),QCOR1(nled,ndnd)
      DIMENSION QCOR2(nled,ndnd),CORLGC(nled,ndnd),VBL1CI(nyear)
      DIMENSION CORJK(nyear),VBLJ1(nyear)
      DIMENSION VBLJ2(nyear),VRBL1C(nyear),VRBL2Z(nyear)
      DIMENSION RMSPD(nled,ndnd),RMSPC(nled,ndnd)
      DIMENSION RMSPCI(nled,ndnd),RMSPR(nled,ndnd)
      DIMENSION HVCOR(nled,ndnd)
      DIMENSION HEIDKE(nled,ndnd)
      DIMENSION IVBL1(nyear),IVBL2(nyear),KNTFC(3),KNTOB(3)
      DIMENSION V1HKE(nyear),V2HKE(nyear)
      DIMENSION acgeoc(ndnd,nlag),acgeop(ndnd,nlag)
      DIMENSION accca(nlag),acper(nlag)
      character*80 rdir
C
      DATA IPREPC/0/
      DATA NODAMP/0/
C
C   --WHEN IPREPC=1 (JCL), A WRITE(93 STMT WRITES PERSISTENCES TO FILE
C   --FOR RUN WHERE PREDICTAND IS CHANGE FROM FINAL P'DOR TO P'DAND.
C   --IF IW9=1, SKILL RESULTS ARE WRITTEN TO IU=9, REGARDLESS OF VALUE
C   --OF IPREPC; IF IW9=0 WRITING TO IU=9 STILL MAY OCCUR DUE TO IPRECP.
C   --IF NODAMP=1 (WITH IPREPC=1), HVCOR IS SET TO 1 (NO DAMPING OF
C   --PERSIS FCSTS IN HYBRID CCA/PERSISTENCE FORECAST SKILL TESTS).
C
      DATA IW9/1/
C
C--------Output files---------------------------------------------------
C
      open(9,file=rdir(1:nfill(rdir)) // '/ccaskl',form='formatted')
C
      open(10,file=rdir(1:nfill(rdir)) // '/geocca',form='formatted')
      open(11,file=rdir(1:nfill(rdir)) // '/geoper',form='formatted')
C
      open(40,file=rdir(1:nfill(rdir)) // '/skill.1d.ac.cca',
     +form='formatted')
      open(41,file=rdir(1:nfill(rdir)) // '/skill.1d.ac.per',
     +form='formatted')
C
      open(51,file=rdir(1:nfill(rdir)) // '/sdtmpfct',form='formatted')
C
      open(60,file=rdir(1:nfill(rdir)) // '/ccaskl2',form='formatted')
C
      open(81,file=rdir(1:nfill(rdir)) // '/zhatus.output',
     +form='formatted')
C
      do lead=1,nled
      NYER(lead)=0
      enddo
      NYER(5)=1
C
      do i=1,nyear
      VRBL3(i)=0.
      enddo
C
      nyearM1=nyear-1
      nyearM2=nyear-2
C
      do LEAD=1,nled
      NYER(LEAD)=nyear-NYER(LEAD)
      enddo
C
C-----------------------------------------------------------------------
C
C         COMPUTE STATISTICS
C
C-----------------------------------------------------------------------
C
C    SKILL OF UNINFLATED FORECASTS OUT TO 5 SEAS LEAD FOR
C    THE SFC STNS  BASED ON ROOT MEAN SQUARED ERROR (RMSE)
C    FOR LEADING PREDICTOR SEASON iseas OVER NYER(1) YRS
C    (BUT NYER(5) YRS FOR 5 SEASON LEAD) SEASON LEAD
C    RMSE BY STATION IMPROVEMENT OVER CLIMATOLOGY FCST IN PARENS
C    FOR THE ndnd PREDICTANDS
C
      DO 201 LEAD=LEDBEG,LEDEND
      nyearX=NYER(LEAD)
C
      DO 202 istn=1,ndnd
C
      do IYROUT=1,nyearX
      Q(IYROUT)=QKT2(LEAD,IYROUT)
      VRBL1(IYROUT)=zpred(IYROUT,LEAD,istn)
      VRBL2(IYROUT)=ZREAL(IYROUT,LEAD,istn)
      enddo
C
      CALL RMSERR(VRBL1,VRBL2,nyearX,RMSE,ER1)
      CALL RMSERR(VRBL3,VRBL2,nyearX,RMSE2,ER2)
      CALL QCOR(VRBL1,VRBL2,Q,nyearX,QC1,QC2)
      QCOR1(LEAD,istn)=QC1
      QCOR2(LEAD,istn)=QC2
      RMSER(LEAD,istn)=RMSE
      RMSER2(LEAD,istn)=RMSE2
      ERBR(LEAD,istn)=ER1
      ERBR2(LEAD,istn)=ER2
      SCORM(LEAD,istn)=100.-(RMSER(LEAD,istn)/RMSER2(LEAD,istn)*100.)
      SCORME(LEAD,istn)=100.-(ERBR(LEAD,istn)/ERBR2(LEAD,istn)*100.)
C
  202 CONTINUE
C
  201 continue
C
C   --GO THROUGH TMP STATN 1 THROUGH ndnd, EACH FOR 5 SEASONAL PDANDS
C   --CALCULATE Z-SCORE FOR THE TMP PREDICTIONS FOR FINAL (YEAR=nyear) YR
C
      DO 16 I=1,ndnd
      DO 17 LEAD=LEDBEG,LEDEND
      nyearX=NYER(LEAD)
C
      CALL STATS(zpred(1,LEAD,I),nyearX,-999.,NX,XBAR,STDEV,XMAX,XMIN)
C
      S(LEAD,I)=STDEV
C
C   --INFLATE AMPLITUDE OF FCSTS TO THAT OF THE OBSERVATIONS (SD OF 1).
C
      do IYROUT=1,nyearX
      Z(IYROUT,LEAD,I)=zpred(IYROUT,LEAD,I)/STDEV
      enddo
C
   17 CONTINUE
   16 CONTINUE
C
      DO 301 LEAD=LEDBEG,LEDEND
      nyearX=NYER(LEAD)
      DO 302 istn=1,ndnd
C
      do IYROUT=1,nyearX
      VRBL1(IYROUT)=Z(IYROUT,LEAD,istn)
      VRBL2(IYROUT)=ZREAL(IYROUT,LEAD,istn)
      enddo
C
      CALL RMSERR(VRBL1,VRBL2,nyearX,RMSE,ER1)
      CALL RMSERR(VRBL3,VRBL2,nyearX,RMSE2,ER2)
      CALL QCOR(VRBL1,VRBL2,Q,nyearX,QC1,QC2)
      QCOR1(LEAD,istn)=QC1
      QCOR2(LEAD,istn)=QC2
      RMSER(LEAD,istn)=RMSE
      RMSER2(LEAD,istn)=RMSE2
      ERBR(LEAD,istn)=ER1
      ERBR2(LEAD,istn)=ER2
      SCORM(LEAD,istn)=100.-(RMSER(LEAD,istn)/RMSER2(LEAD,istn)*100.)
      SCORME(LEAD,istn)=100.-(ERBR(LEAD,istn)/ERBR2(LEAD,istn)*100.)
  302 CONTINUE
  301 continue
C
      DO 116 LEAD=LEDBEG,LEDEND
      nyearX=NYER(LEAD)
C
C ...... NEXT CARD IS DEFAULT IF NOT OVERWRITTEN BY MORE EXACT CALC BELOW
C
      EXPD=FLOAT(nyearX)/3.
      AVSKL=0.
      SKL50=0.
      AVSKL0=0.
      HIT50=1./FLOAT(ndnd)
      AVSKLH=0.
C
      DO 101 istn=1,ndnd
      HITS=0.
      do ICAT=1,3
      KNTFC(ICAT)=0
      KNTOB(ICAT)=0
      enddo
C
C   --IN PRINCIPLE, THE VRBL2 SHOULD HAVE ITS ELEMENTS JACKKNIFED (C-V)
C   --BY SEPARATE RESTANDARDIZATIONS FOR EACH YR WITHHELD (VAN DEN DOOL,
C   --1987) BUT HERE IT IS FOUND TO MAKE VIRTUALLY NO DIFFERENCE.
C
      do IYROUT=1,nyearX
      VRBL1(IYROUT)=Z(IYROUT,LEAD,istn)
      VRBL2(IYROUT)=ZREAL(IYROUT,LEAD,istn)
      enddo
C
      CALL CORREL(VRBL1,VRBL2,nyearX,AVG1,AVG2,SD1,SD2,COR)
C
C   --COMPUTE HEIDKE SKILL SCORE
C
      do IYROUT=1,nyearX
C
      V1HKE(IYROUT)=(VRBL1(IYROUT)-AVG1)/SD1
      V2HKE(IYROUT)=(VRBL2(IYROUT)-AVG2)/SD2
      IVBL1(IYROUT)=2
      IVBL2(IYROUT)=2
      IF(V1HKE(IYROUT).GT.0.4307)IVBL1(IYROUT)=3
      IF(V1HKE(IYROUT).LT.-0.4307)IVBL1(IYROUT)=1
      IF(V2HKE(IYROUT).GT.0.4307)IVBL2(IYROUT)=3
      IF(V2HKE(IYROUT).LT.-0.4307)IVBL2(IYROUT)=1
      IF(IVBL1(IYROUT).EQ.IVBL2(IYROUT))HITS=HITS+1.
C
C   --COUNT TOTAL # FCSTS AND #OBS FOR CALC OF EXPECTED FRACTION OF HITS
C     KNTFC(IVBL2(IYROUT))=KNTFC(IVBL2(IYROUT))+1.
C     KNTOB(IVBL2(IYROUT))=KNTOB(IVBL2(IYROUT))+1.
C
      enddo
C
C     SUMEX=(KNTFC(1)*KNTOB(1)+KNTFC(2)*KNTOB(2)+KNTFC(3)*KNTOB(3))
C    $  /FLOAT(nyearX*nyearX)
C     EXPD=SUMEX*FLOAT(nyearX)
C
      HEIDKE(LEAD,istn)=(HITS-EXPD)/(FLOAT(nyearX)-EXPD)
C
C   --ADJUST FOR DEGENERATE (HIGHLY NEGATIVE) CCA C-V CORREL SKL SCORES;
C   --AMPLITUDE OF FCSTS WAS DIVIDED BY S (TO INFLATE) SO MUST MULT BACK
C
      IF(COR.LT.0.)COR=COR*(SD1*S(LEAD,istn)/SD2)
C
C   --DO THE SAME FOR HEIDKE SKILL, WHICH IS DEGENERATE WHEN INFLATED
C
      IF(HEIDKE(LEAD,istn).LT.0.)
     $ HEIDKE(LEAD,istn)=HEIDKE(LEAD,istn)*(SD1*S(LEAD,istn)/SD2)
      CORSKL(LEAD,istn)=COR
      AVSKL=AVSKL+COR/FLOAT(ndnd)
      AVSKLH=AVSKLH+HEIDKE(LEAD,istn)/FLOAT(ndnd)
      AVSKL0=AVSKL0+MAX(COR,0.)/FLOAT(ndnd)
      IF(COR.GE.0.5)SKL50=SKL50+HIT50
C
  101 CONTINUE
C
C-----------------------------------------------------------------------
C
C        Write out geographical distribution of CV AC skill score
C        for CCA forecasts
C
C-----------------------------------------------------------------------
C
      mon=iseas+ilag+(LEAD-1)*3+2
      monpl=mon
      if ((monpl.gt.12).and.(monpl.le.24)) mon=mon-12
      if (monpl.gt.24) mon=mon-24
      monpl=mon
      leadpl=ilag+(LEAD-1)*3
      write(10,9010) MONPL,LEADPL,(CORSKL(LEAD,istn),istn=1,ndnd)
 9010 format(2X,'SEASON = ',I4,4X,'LEAD = ',I4,
     +/,6(10(F7.3),/))
C
      do istn=1,ndnd
      acgeoc(istn,ilagf)=CORSKL(LEAD,istn)
      enddo
C
C-----------------------------------------------------------------------
C
C        Write out area averaged CV AC skill score
C        for CCA forecast
C
C-----------------------------------------------------------------------
C
      mon=iseas+ilag+(LEAD-1)*3+2
      monpl=mon
      if ((monpl.gt.12).and.(monpl.le.24)) mon=mon-12
      if (monpl.gt.24) mon=mon-24
      monpl=mon
      leadpl=ilag+(LEAD-1)*3
      write(40,9040) MONPL,LEADPL,avskl
 9040 format(i4,i4,f8.3)
C
      accca(ilagf)=avskl
      print *,' for ilag ',ilag,' ilagf ',ilagf,' accca ',accca(ilagf)
C
C-----------------------------------------------------------------------
C
  116 CONTINUE
C
      DO 118 LEAD=LEDBEG,LEDEND
      nyearX=NYER(LEAD)
      NYRXM1=nyearX-1
      NYRXM2=nyearX-2
      AVPDIR=0.
      AVPABD=0.
      AVPCV=0.
      AVPCV0=0.
C
      DO 107 istn=1,ndnd
C
      do IYROUT=1,nyearX
      VRBL1(IYROUT)=zpers(IYROUT,istn)
      VRBL2(IYROUT)=ZREAL(IYROUT,LEAD,istn)
      enddo
C
C   --CALCULATE PERSISTENCE IN A CROSS-VALIDATION FRAMEWORK TO MATCH
C   --THE CCA FRAMEWORK, IN ADDITION TO THE DIRECT PERSISTENCE CALCUL.
C
      DO 500 IYROUT=1,nyearX
      KNTT=0
C
      do JK=1,nyearX
      IF(JK.ne.IYROUT) then
      KNTT=KNTT+1
      VBLJ1(KNTT)=VRBL1(JK)
      VBLJ2(KNTT)=VRBL2(JK)
      endif
      enddo
C
      CALL CORREL(VBLJ1,VBLJ2,NYRXM1,AV1,AV2,SD1,SD2,CORR)
C
      CORJK(IYROUT)=CORR
C
C ...... NEXT CARD TO RESTANDARDIZE VERIFICATION TMP WRT GROUP W/O TARGET
C ...... YEAR INCLUDED, TO AVOID VAN DEN DOOL BIAS.
C
      VRBL2Z(IYROUT)=(VRBL2(IYROUT)-AV2)/SD2
      VRBL1C(IYROUT)=VRBL1(IYROUT)*CORJK(IYROUT)
C
  500 CONTINUE
C
      CALL CORREL(VRBL1,VRBL2,nyearX,AVG1,AVG2,SD1,SD2,COR)
      CALL CORREL(VRBL1C,VRBL2Z,nyearX,AVGG1,AVGG2,SDD1,SDD2,CORCV)
C
C ...... LAST CARD USES VERIF WRT JACKKNIFE GROUP (VRBL2Z, NOT VRBL2), AND
C ...... PREVIOUS CARD VERIFIES DIRECTLY (NOT C-V). BOTH FOR PERSISTENCE.
C
      CORLAG(LEAD,istn)=COR
C
C ...... ADJUST FOR DEGENERATE (HIGHLY NEGATIVE) C-V PERSIS CORREL SCORES
C
      IF(CORCV.LT.0.)CORCV=CORCV*(SDD1/SDD2)
      CORLGC(LEAD,istn)=CORCV
      AVPDIR=AVPDIR+COR/FLOAT(ndnd)
      AVPABD=AVPABD+ABS(COR)/FLOAT(ndnd)
      AVPCV=AVPCV+CORCV/FLOAT(ndnd)
      AVPCV0=AVPCV0+MAX(CORCV,0.)/FLOAT(ndnd)
C
C   --COMPUTE RMSE OF THE TWO VERSIONS OF THE PERSISTENCE FORECAST
C   --(DIRECT-RMSPD, C-V -RMSPC), AND FOR DIRECT PERSISTENCE REVERSED
C   --(RMSPR) WHICH FARES BETTER THAN NONREVERSED FOR ANTIPERSISTENCE.
C   --ALSO COMPUTE RMSE FOR INFLATED C-V PERSIS (VBL1CI) FOR COMPAR-
C   --ABILITY TO THE INFLATED CCA FCST FOR RMSE (STANDARDIZE VRBL1C).
C
      CALL STATS(VRBL1C,nyearX,-999.9,NX,XBAR,STDEV,XMAX,XMIN)
C
      do IYR=1,nyearX
      VBL1CI(IYR)=(VRBL1C(IYR)-XBAR)/STDEV
      enddo
C
      SUMD=0.
      SUMR=0.
      SUMC=0.
      SUMCI=0.
      do IYR=1,nyearX
      SUMD=SUMD+(VRBL1(IYR)-VRBL2(IYR))**2
      SUMR=SUMR+((-1.*VRBL1(IYR))-VRBL2(IYR))**2
      SUMC=SUMC+(VRBL1C(IYR)-VRBL2Z(IYR))**2
      SUMCI=SUMCI+(VBL1CI(IYR)-VRBL2Z(IYR))**2
      enddo
C
      RMSPD(LEAD,istn)=SQRT(SUMD/FLOAT(nyearX))
      RMSPR(LEAD,istn)=SQRT(SUMR/FLOAT(nyearX))
      RMSPC(LEAD,istn)=SQRT(SUMC/FLOAT(nyearX))
      RMSPCI(LEAD,istn)=SQRT(SUMCI/FLOAT(nyearX))
C
  107 CONTINUE
C
C-----------------------------------------------------------------------
C
C        Write out area averaged CV AC skill score
C        for persistence forecast
C
C-----------------------------------------------------------------------
C
      mon=iseas+ilag+(LEAD-1)*3+2
      monpl=mon
      if ((monpl.gt.12).and.(monpl.le.24)) mon=mon-12
      if (monpl.gt.24) mon=mon-24
      monpl=mon
      leadpl=ilag+(LEAD-1)*3
      write(41,9040) MONPL,LEADPL,avpcv
C
      acper(ilagf)=avpcv
      print *,' for ilag ',ilag,' ilagf ',ilagf,' acper ',acper(ilagf)
C
C-----------------------------------------------------------------------
C
      IF(IPREPC.EQ.1)THEN
C
C   --WRITE DIRECT PERSIS DATA OUT TO FILE FOR DOOL EXPT (HYBRID EXPT)
C
        do IA=1,ndnd
        IF(NODAMP.NE.1)HVCOR(LEAD,IA)=CORLAG(LEAD,IA)
        IF(NODAMP.EQ.1)HVCOR(LEAD,IA)=1.
        enddo
C
      END IF
C
C   --CONVERT TO MAP PRINT ORDER FOR CROSS-VALID. PERSIS MAP
C
C-----------------------------------------------------------------------
C
C        Write out geographical distribution of CV AC skill score
C        for persistence forecast
C
C-----------------------------------------------------------------------
C
      mon=iseas+ilag+(LEAD-1)*3+2
      monpl=mon
      if ((monpl.gt.12).and.(monpl.le.24)) mon=mon-12
      if (monpl.gt.24) mon=mon-24
      monpl=mon
      leadpl=ilag+(LEAD-1)*3
      write(11,9010) MONPL,LEADPL,(CORLGC(LEAD,istn),istn=1,ndnd)
c 
      do istn=1,ndnd
      acgeop(istn,ilagf)=CORLGC(LEAD,istn)
      enddo
C
C-----------------------------------------------------------------------
C
  118 CONTINUE
C
      do LEAD=LEDBEG,LEDEND
      do istn=1,ndnd
      SCORCR(LEAD,istn)=(CORSKL(LEAD,istn)-ABS(CORLGC(LEAD,istn)))/
     $ (1.-ABS(CORLGC(LEAD,istn)))*100.
      enddo
      enddo
C
C     WRITE FCST STDEV (S) TO EXTERNAL FILE, NOT IN MAP PRINT ORDER.
C
      close(40)
      close(41)
      close(42)
      close(43)
      close(45)
c
      close(60)
c
      close(82)
      close(83)
      close(84)
c
      RETURN
      END
