SUBROUTINE mkirpc(fnxy,fndat,nl,numsat,timage,rstride,speed,head,clat,clon,firstsat)
  !
  ! Ver 2.1.0 for WCOSS
  !
  ! This subroutine creates an analysis of IR brightness temperatures on a
  ! cylindrical grid rotated with respect to TC motion and then calculates 
  ! 2-D principle components for a single image.  These are output
  ! to a file (YYMMDDHHBBYYYY_IRPC0.dat) as scaled integers, with the first
  ! value being time (hr*100+min) before/after synoptic time followed by 9
  ! PCs.  PCs are scaled by 100.  Dividing by these scaling factors will
  ! result in the real values.  
  !
  ! INPUT: one file
  !
  ! 1) irdump1.dat (number will vary from 1 to the number of images
  !                 being processed for a given time - not exceeding 9)
  !
  ! OUTPUT:  two files
  ! 
  ! 1) varan.log (will tell you if there are any problems with the variational
  !               analysis, typically caused by not enough data to analyze)
  ! 2) YYMMDDHHBBYYYY_IRPC0.dat (file containing IRPC lines for lsdiag files)
  !
  ! Calling variables:
  !(fnxy,fndat,nl,numsat,timage,rstride,speed,head,clat,clon,firstsat)
  ! 1)fnxy     : filename for input file "irdump#.dat"
  ! 2)fndat    : filename for output file "YYMMDDHHBBYYYY_IRPC0.dat"
  !              -"YYMMDDHHBBYYYY_IRRP#.dat" is used as template
  !              -"PC0" is set in the subroutine
  ! 3)nl       : number of the satellite images currently being processed
  ! 4)numsat   : total number of images to be processed per time (<=9)
  ! 5)timage   : time in fractional hours between time and image time
  ! 6)rstride  : time in fractional hours between image times to search for
  ! 7)speed    : storm speed
  ! 8)head     : storm heading
  ! 9)clat     : storm latitude
  ! 10)clon    : storm longitude
  ! 11)firstsat : flag to determine if first image processed for time
  !
  ! Written by J. Knaff NOAA/NESDIS
  !
  !      Modified: 03 May 2016 (M. DeMaria) - to add ioper parameter to use
  !      environmental variable to define location of EVEC_cov_ra.OUT and meansd.inp files
  !
  !      Modified: 23 March 2016 (K. Musgrave CIRA) -converted to subroutine
  !      Modified: 20 November 2015 (K. Musgrave CIRA)
  !      Modified: 29 April 2015
  !      Modified: 18 June 2019 (S. Stevenson) - added logic to set output to
  !      missing if more than 10% of IR pixels in file are outside logical range
  !      to mitigate G17 LHP issues
  !      Modified: 07 April 2020 (S. Stevenson) - removed ioper flag, fixed
  !      output to write 7-day lsdiag lines
  !
  !*************************************************************************

  IMPLICIT NONE

  !Calling variables
  CHARACTER (LEN=11), INTENT(IN) :: fnxy
  CHARACTER (LEN=24), INTENT(IN) :: fndat
  INTEGER, INTENT(IN) :: nl, numsat
!  INTEGER, INTENT(IN) :: iblank
  REAL, INTENT(IN) :: timage, rstride, speed, head, clat, clon
  LOGICAL, INTENT(IN):: firstsat

  !Local variables
  INTEGER:: i,j,k
  INTEGER:: imiss=9999
  REAL:: rmiss=-99.00
  CHARACTER (LEN=256) coef_location,FN_means_W,FN_EVEC_W

  ! numpc = number of PCs produced

  INTEGER, PARAMETER:: numpc=9

  ! ix=100        ! radial points (4 km)
  ! jy=72         ! azimuthal points (5 degree)  

  INTEGER, PARAMETER :: ix=100, jy=72

  ! xmax = radius for analysis (i.e. every thing within or equal...)
  ! x0  initial radius value

  REAL, PARAMETER :: xmax=402.0, x0=2.0

  ! ymax = max azimuth for analysis (i.e. every thing within or equal...)
  ! y0  initial azimuth value

  REAL, PARAMETER :: ymax=360.0, y0=0

  REAL :: yj(0:jy), xi(0:ix)
  REAL :: axy(0:ix,0:jy),bxy(0:ix,0:jy),fxy(0:ix,0:jy)
  REAL :: dx, dy

  ! Boundary conditions for the filter parameters
  INTEGER, PARAMETER :: ibc=0, jbc=1
  !  INTEGER, PARAMETER :: iax=-80, iby=-90
  INTEGER, PARAMETER :: iax=-40, iby=-40

  INTEGER, PARAMETER ::ixx=170000  ! maximum x-dimension of the analysis grid
  INTEGER:: kk
  INTEGER::nrec, ierr

  ! Varan variables
  REAL :: gk(ixx),fk(ixx),wk(ixx),xk(ixx),yk(ixx)
  REAL :: rlat(ixx),rlon(ixx)
  REAL :: tlat, tlon, xt, yt, radt, rt, at
!  REAL :: clat,clon,head,speed,us,vs
  REAL :: us,vs

  ! File names
  CHARACTER (LEN=80) :: FN_IR00
  CHARACTER (LEN=80) :: FN_IRPC
  INTEGER:: luir00=12,luirpc=13,lulog=11

  ! IR files
  LOGICAL :: ir00=.false.
  INTEGER:: ir00yr,ir00mon,ir00day,ir00hour,ir00min,ir00sec,i00,j00
!  INTEGER:: isjd,irjd
!  REAL:: dt
  REAL, ALLOCATABLE, DIMENSION(:,:):: ir00lat,ir00lon,ir00bt

  ! Rotation with respect to motion  (set)
  LOGICAL, PARAMETER :: rotate=.true.
  LOGICAL :: ixxfailure=.false.
  REAL :: rot,xdir,ydir
  REAL, ALLOCATABLE :: frot(:,:)
  INTEGER:: IRM,irot

  ! PCs
  INTEGER:: lumean=32,lueig=33
  CHARACTER (LEN=80) :: FN_means='meansd.inp', FN_EVEC='EVEC_cov_ra.OUT'
  REAL :: xbar(0:ix,0:jy),stdev(0:ix,0:jy)
  REAL :: EVEC1(0:ix,0:jy),EVEC2(0:ix,0:jy),EVEC3(0:ix,0:jy),&
       EVEC4(0:ix,0:jy),EVEC5(0:ix,0:jy),EVEC6(0:ix,0:jy), &
       EVEC7(0:ix,0:jy),EVEC8(0:ix,0:jy),EVEC9(0:ix,0:jy)
  REAL, PARAMETER :: eval1=2325860.250, eval2=740103.438,eval3=636854.062, &
       eval4=485390.000, eval5=203416.438, eval6=120893.828, &
       eval7=119322.953, eval8=105398.062, eval9=98946.766

  INTEGER:: nummiss, counter, numhrs
  REAL, ALLOCATABLE, DIMENSION(:,:):: PC00
!  REAL :: PC(numpc),SC(numpc),slope(numpc)
!  INTEGER, DIMENSION(numpc):: countPC=0
!  REAL :: x(numsat),y(numsat)
!  REAL, ALLOCATABLE, DIMENSION(:)::xf,yf
!  REAL:: sig,a,b,sa,sb,chi2,q
!  INTEGER:: iaipc(21),isipc(21),idipc(21)
  INTEGER, DIMENSION(numsat,29):: ipc
  CHARACTER (LEN=4) :: ctech
!  INTEGER, DIMENSION(21):: iblankpc
  INTEGER:: itimage, ithr, itmin
  INTEGER:: btcount, btsize

  !***************************************************************************
  ! EXECUTABLES FOLLOW
  !***************************************************************************

!  WRITE(6,*) 'Input to mkirpc:'
!  WRITE(6,*) 'fnxy,fndat,nl,numsat,timage,rstride'
!  WRITE(6,*) fnxy,fndat,nl,numsat,timage,rstride
!  WRITE(6,*) 'speed,head,clat,clon,firstsat'
!  WRITE(*,*) speed,head,clat,clon,firstsat

  !Set up analysis options

  dx = (xmax-x0)/float(ix)
  dy = (ymax-y0)/float(jy)


  DO i=0,ix
     xi(i)=x0+dx*float(i)
  END DO

  DO j=0,jy
     yj(j)=y0+dy*float(j)
  END DO

  ! open logfile
  OPEN (UNIT=lulog,FILE='varan.log',STATUS='UNKNOWN')
  
  !read in the means and standard deviations
  call getenv( "SHIPS_COEF", coef_location )
  FN_means_W = trim( coef_location ) // FN_means
  OPEN(lumean,FILE=FN_means_W)

  DO i=0,ix

     DO j=1,jy

        READ(lumean,*) xbar(i,j),stdev(i,j)

     END DO

     xbar(i,0)  = xbar(i,jy)
     stdev(i,0) = stdev(i,jy)

  END DO

  CLOSE(lumean)

  ! read in the first 9 eigenvectors
  call getenv( "SHIPS_COEF", coef_location )
  FN_EVEC_W = trim( coef_location ) // FN_EVEC
  OPEN(lueig,FILE=FN_EVEC_W, FORM='unformatted',ACCESS='direct',&
       RECL= (ix+1)*(jy+1)   )  ! gfortran/pgf90 4*(ix+1)*(jy+1) this number
                                  ! ifort (ix+1)*(jy+1)

  READ(lueig,REC=1)EVEC1
  READ(lueig,REC=2)EVEC2
  READ(lueig,REC=3)EVEC3
  READ(lueig,REC=4)EVEC4
  READ(lueig,REC=5)EVEC5
  READ(lueig,REC=6)EVEC6
  READ(lueig,REC=7)EVEC7
  READ(lueig,REC=8)EVEC8
  READ(lueig,REC=9)EVEC9

  CLOSE (lueig)

  ! heading
  CALL hs2uv(head,speed,us,vs)
  CALL xytora(us,vs,rt,at)
  IF (rotate) THEN
     rot=at-90.0
     IF (rot .LT. 0) rot=360.0 + rot
  END IF

  ! Open the IRXY file
  FN_IR00=TRIM(fnxy)
  FN_IRPC=TRIM(fndat)
  FN_IRPC(18:19) = 'PC'
  WRITE(FN_IRPC(20:20),'(i1)') 0
  ALLOCATE (PC00(numsat,9))
  PC00=0.0

  INQUIRE(file=FN_IR00,EXIST=ir00)
  IF (ir00) THEN
     OPEN(UNIT=luir00, file=FN_IR00)

     ! Read specified IR file
     READ(luir00,'(i4,i2,i2,1x,i2,i2,i2,2i7)')ir00yr,ir00mon, &
          ir00day,ir00hour,ir00min,ir00sec,i00,j00

     IF (ALLOCATED(ir00lat)) DEALLOCATE(ir00lat)
     IF (ALLOCATED(ir00lon)) DEALLOCATE(ir00lon)
     IF (ALLOCATED(ir00bt )) DEALLOCATE(ir00bt )
     
     ALLOCATE (ir00lat(i00,j00),ir00lon(i00,j00),ir00bt(i00,j00))

     READ(luir00,'(20(f9.3,1x))')((ir00bt(i,j),j=1,j00),i=1,i00)
     READ(luir00,'(20(f9.3,1x))')((ir00lon(i,j),j=1,j00),i=1,i00)
     READ(luir00,'(20(f9.3,1x))')((ir00lat(i,j),j=1,j00),i=1,i00)

     CLOSE (luir00)

     ! Analyze the BT field for IR file
     ixxfailure=.false.
     nrec=0
     data00: DO i=1,i00
        DO j=1,j00
           tlat=ir00lat(i,j)
           tlon=ir00lon(i,j)
           CALL lltoxy(tlat,tlon,clat,clon,xt,yt)
           radt=SQRT(xt**2+yt**2)
           IF (radt .le. (xmax+20.) .AND. ir00bt(i,j).GT. 0.0) THEN
              nrec = nrec + 1
              IF (nrec .GT. ixx) THEN
                 WRITE(lulog,&
                      '("IR_cc_anal: nrec > parameter ixx... increase ixx")')
                 WRITE(lulog,'("nrec,ixx:  ")') nrec, ixx
                 
                 ixxfailure=.true.
                 ir00=.false.
                 EXIT data00
              ENDIF
              gk(nrec) = ir00bt (i,j)
              wk(nrec) = 1.0
              CALL xytora(xt,yt,rt,at)
              xk(nrec)=rt
              yk(nrec)=at
           ENDIF
        ENDDO
     ENDDO data00

     ierr=0
     fk=250.0
     fxy=250.0
     IF (.NOT. ixxfailure) THEN
        CALL varan(gk,wk,xk,yk,fk,nrec,axy,bxy,fxy,xi,yj, &
             ix,ix,jy,ibc,jbc,iax,iby,lulog,ierr)

        IF (ierr /= 0) THEN
           WRITE(lulog,'("variational anlysis failed for ",a)')TRIM(FN_IR00)
           WRITE(lulog,'("ierr = ",i3)')ierr
           WRITE(lulog,'(i10," data points used in the analysis")')nrec
           ir00=.false.
        ENDIF
        
        IF (rotate) THEN
           IF (ALLOCATED(frot)) DEALLOCATE(frot)
           ALLOCATE(frot(0:ix,0:jy))
          
           ! find the rotation in degrees
           IRM = NINT(rot/dy)
          
           DO j=0, jy
              ! find the index corresponding to the rotation
              irot=j-IRM
              IF(irot .LE. 0) THEN
                 irot = jy +irot
              ENDIF
              IF(irot .GT. jy) THEN
                 irot= irot-jy
              ENDIF

              DO i=0,ix
                 frot(i,irot)=fxy(i,j)
              ENDDO
           ENDDO
           
           DO i=0,ix
              frot(i,0)=fxy(i,jy)
           ENDDO
           
           fxy = frot
        ENDIF

        !Calculate PCs with eigenvectors and mean/sd read in earlier.
        DO i=0,ix
           DO j=1,jy
              PC00(nl,1)=PC00(nl,1)+EVEC1(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,2)=PC00(nl,2)+EVEC2(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,3)=PC00(nl,3)+EVEC3(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,4)=PC00(nl,4)+EVEC4(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,5)=PC00(nl,5)+EVEC5(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,6)=PC00(nl,6)+EVEC6(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,7)=PC00(nl,7)+EVEC7(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,8)=PC00(nl,8)+EVEC8(i,j)*(fxy(i,j)-xbar(i,j))
              PC00(nl,9)=PC00(nl,9)+EVEC9(i,j)*(fxy(i,j)-xbar(i,j))
           ENDDO
        ENDDO

        ! Normalized with unit variance...
        PC00(nl,1)=PC00(nl,1)/sqrt(eval1)
        PC00(nl,2)=PC00(nl,2)/sqrt(eval2)
        PC00(nl,3)=PC00(nl,3)/sqrt(eval3)
        PC00(nl,4)=PC00(nl,4)/sqrt(eval4)
        PC00(nl,5)=PC00(nl,5)/sqrt(eval5)
        PC00(nl,6)=PC00(nl,6)/sqrt(eval6)
        PC00(nl,7)=PC00(nl,7)/sqrt(eval7)
        PC00(nl,8)=PC00(nl,8)/sqrt(eval8)
        PC00(nl,9)=PC00(nl,9)/sqrt(eval9)
     ENDIF
  ENDIF

!  write(6,*) PC00(nl,:)

  !Prepare lsdiag lines
  ipc=imiss

  DO i=1,9
     j=i+1
     IF (ir00) THEN
        IF (i .eq. 1) THEN
           ithr=int(abs(timage))
           itmin=nint((abs(timage)-float(ithr))*60.)
           itimage=ithr*100 + itmin
!           ipc(nl,i)=int(rstride*float(nl-1)*60.)
           IF (timage .ge. 0) THEN
              ipc(nl,i)=itimage
           ELSE
              ipc(nl,i)=-1*itimage
           ENDIF
        ENDIF
        IF (PC00(nl,i) .ne. rmiss) THEN
           ipc(nl,j)=nint(PC00(nl,i)*100.)
        ELSE
           ipc(nl,j)=imiss
        ENDIF
!       Additional QC if more than 10% bad bt values
        btcount=count(ir00bt .ge. 183 .and. ir00bt .le. 313)
        btsize=size(ir00bt)
        IF ((real(btcount)/real(btsize))*100 .lt. 90) THEN
           ipc(nl,j)=imiss
        ENDIF
     ENDIF
  ENDDO

  IF (firstsat) THEN
     OPEN(unit=luirpc,file=FN_IRPC,status='replace',form='formatted')
  ELSE
     OPEN(unit=luirpc,file=FN_IRPC,status='old', &
          position='append',form='formatted')
  ENDIF
  ctech='PC00'
  IF (nl .gt. 1) THEN
     numhrs=floor(rstride*float(nl-1))
     ctech(3:3)='M'
     WRITE(ctech(4:4),'(i1)') numhrs
  ENDIF

  IF (ir00) THEN
!     WRITE(6,'(10x,21i5,1x,a4)')(ipc(nl,k),k=1,21),ctech
     WRITE(luirpc,'(10x,29i5,1x,a4)')(ipc(nl,k),k=1,29),ctech
!  ELSE
!     IF (iblank .eq. 1) THEN
!        iblankpc=imiss
!        WRITE(luirpc,'(10x,21i5,1x,a4)')(iblankpc(k),k=1,21),ctech
!     ENDIF
  ENDIF
     
  CLOSE(luirpc)
  CLOSE(lulog)
  
END SUBROUTINE mkirpc
subroutine hs2uv (h,s,u,v)
  !     this subroutine computes the u % v components
  !     from the heading and and speed (h % s)
   
  REAL:: h,s,u,v,pi
  pi= 2.0* acos(0.)
   
  u = s*sin(h*pi/180.)
  v = s*cos(h*pi/180.)
  
  return
end subroutine hs2uv