c     proj_PM.f
c
c     Project Hs from Pt, trGt, PtPCs, trGtPCs
c
      use commondata, only : tmp_path,coef_path,eof_path,
     & era_pc_path,lambda_path

      integer itmp,cnt,i,j,k,kth,ith,jth,cnt1(9146),oHsFlg

      integer nvari(9146),ystrt,yend,lflg,yy1,yy2,yy3,mm1,mm2,mm3,
     &        dd1,dd2,dd3,hh1,hh2,hh3,Nx(9146),p(9146),idx(9146,68),
     &        idiff,sm,fstatus
      real lambda(9146),Hs_std(9146),Pt(9146),Gt(9146),
     &     Pt_std(9146),Gt_std(9146),PCs(60),PCs_std(30),
     &     coefs(9146,70),Hslambda(9146),GtPCs(60),rtmp,rtmp2(9146),
     &     rtmp3(62),transver,oHs(9146),oHs1(30,9146),PHs(9146),
     &     GtPCs_std(30),lag4(4,30),Hsmax(9146),coefs_lm(9146,70),
     &     oHs_lm(9146),Hs_mean(9146),tHs(9146),tHs_lm(9146),wtmp
      character ifPt*200, ifGt*200, ifPCs*200, iflambda*200, ifpara*200,
     &     ifidx*200, ifPtstd*200, ifGtstd*200, ifHsstd*200,ifcoef*200,
     &     ifPCsstd*200,cvari(70)*3,ifHslambda*200,ofname*200,
     &     ifGtPCs*200,ifGtPCsstd*200,c1*1,ifHsmax*200,ifpm*200,
     &     ifcoeflm*200,csea*3,cnp*5,carea*4,cmo*40,runs*3,ifHsmean*200

      call GETARG(1,cmo)
      call GETARG(2,runs)
      call GETARG(3,carea)
      call GETARG(4,csea)
      call GETARG(5,cnp)
      read(cnp,*) NP
      if(csea=="jfm") then
        sm=1
      else if(csea=="amj") then
        sm=4
      else if(csea=="jas") then
        sm=7
      else if(csea=="ond") then
        sm=10
      endif

      ifcoef=trim(coef_path)//
     &       trim(carea)//"_"//csea//"_PM_30yrs_para.txt"
      ifidx=trim(coef_path)//
     &       trim(carea)//"_"//csea//"_PM_30yrs_idx.txt"
      ifpm=trim(coef_path)//
     &       trim(carea)//"_"//csea//"_PM_30yrs_pmnum.txt"
      ifcoeflm=trim(coef_path)//
     &       trim(carea)//"_"//csea//"_PM_30yrs_paralm.txt"
      ifPCsstd=trim(era_pc_path)//
     &       trim(carea)//"_"//csea//"_PtPCs_std.txt"
      ifGtPCsstd=trim(era_pc_path)//
     &       trim(carea)//"_"//csea//"_trGtPCs_std.txt"

      ifPCs=trim(tmp_path)//trim(cmo)//"_r"//trim(runs)//
     &       "_Pt_PCs_"//trim(carea)//"_"//csea//".txt"
      ifGt=trim(tmp_path)//trim(cmo)//"_r"//trim(runs)//
     &       "_trGtonHs_"//trim(carea)//"_"//csea//".txt"
      ifPt=trim(tmp_path)//trim(cmo)//"_r"//trim(runs)//
     &       "_PtonHs_"//trim(carea)//"_"//csea//".txt"
      ifGtPCs=trim(tmp_path)//trim(cmo)//"_r"//trim(runs)//
     &        "_trGt_PCs_"//trim(carea)//"_"//csea//"_SQ.txt"

      ifHslambda=trim(lambda_path)//
     &       "Ht_BoxCox_para_"//trim(carea)//"_"//csea//".txt"
      ifHsstd=trim(coef_path)//
     &        "Hs_"//trim(carea)//"_"//csea//"_std.txt"
      ifHsmean=trim(coef_path)//
     &        "Hs_"//trim(carea)//"_"//csea//"_mean.txt"
      ifHsmax=trim(coef_path)//
     &        "Hs_"//trim(carea)//"_"//csea//"_max.txt"
      ofname=trim(tmp_path)//trim(cmo)//"_r"//trim(runs)//
     &       "_Hs_6HLY_Grid1x1_"//trim(carea)//"_"//csea//".txt"
      open(1,file=ifcoef)
      open(2,file=ifpm)
      open(3,file=ifcoeflm)
      do i=1,NP
        read(2,*) itmp,p(i)
        if(itmp.ne.i) then
          print*, "pm file",i,itmp,"error"
          stop
        endif
        read(1,*) ith, itmp, coefs(i,1:68)
        if(ith.ne.i) then
          print*,trim(ifcoef),i,ith,"idx1 error"
          stop
        endif
        do j=1,itmp
          if(abs(coefs(i,j)+999.999).lt.1E-5) then
            print*, trim(ifcoef),i, j, coefs(i,j), "missing"
            stop
          endif
        enddo
        do j=itmp+1, 68
          if(abs(coefs(i,j)+999.999).gt.1E-5) then
            print*,trim(ifcoef),i,j,itmp,coefs(i,j),"not missing"
            stop
          endif
        enddo
        NX(i)=itmp-p(i)-1
        read(3,*) ith, itmp, coefs_lm(i,1:68)
        if(ith.ne.i) then
          print*,trim(ifcoeflm),i,ith,"idx2 error"
          stop
        endif
        if(NX(i).ne.(itmp-1)) then
          print*,trim(ifcoeflm),i,itmp,"!=",NX(i),"+1"
          stop
        endif
        do j=1,itmp
          if(abs(coefs_lm(i,j)+999.999).lt.1E-5) then
            print*, trim(ifcoeflm),i, j, coefs(i,j), "missing"
            stop
          endif
        enddo
        do j=itmp+1, 68
          if(abs(coefs_lm(i,j)+999.999).gt.1E-5) then
            print*,trim(ifcoeflm),i,j,itmp,coefs(i,j),"not missing"
            stop
          endif
        enddo
      enddo
      close(1)
      close(2)
      close(3)

      open(2,file=ifidx)
      do i=1,NP
        read(2,'(i4,i3,68(1x,a3))') itmp,nvari(i),cvari(1:68)
        do j=nvari(i)+1, 68
          if(cvari(j).ne."000") then
            print*,trim(ifidx)
            print*,i,nvari(i),j,cvari(j),"error"
            stop
          endif
        enddo
        if(NX(i).ne.nvari(i)-p(i)) then
           print*,trim(ifidx),i,NX(i),nvari(i),p(i),
     &            "NX not equal nvari-p"
           stop
        endif
        do j=1,NX(i)
          if(cvari(j).eq.' Gt') then
            idx(i,j)=1
          else if(cvari(j).eq.' Pt') then
            idx(i,j)=2
          else
            read(cvari(j),'(1x,i2)') kth
            if(kth.lt.1.or.kth.gt.30) then
              print*,trim(ifidx),i,j,cvari(j),"error"
              stop
            endif
            read(cvari(j),'(a1)') c1
            if(c1.eq.'P') then
              idx(i,j)=2+kth
            else if(c1.eq.'G') then
              idx(i,j)=32+kth
            else
              print*,ifidx,i,"line",j,cvari(j),"error"
              stop
            endif
          endif
        enddo
        do j=NX(i)+1,NX(i)+p(i)
          read(cvari(j),'(1x,i2)') kth
          read(cvari(j),'(a1)') c1
          if(c1.ne.'L'.or.kth.ne.(j-NX(i))) then
            print*,ifidx,i,"line",j,cvari(j),"error"
          endif
        enddo
      enddo
      close(2)

      open(3,file=ifPCsstd)
      read(3,*) PCs_std(1:30)
      close(3)

      open(4,file=ifHslambda)
      do i=1,NP
        read(4,*) itmp,Hslambda(i)
C       if(itmp.ne.i) then
C         print*,trim(ifHslambda),i,"error"
C         stop
C       endif
        if(Hslambda(i).lt.-.5) Hslambda(i)=-.5
      enddo

      open(5,file=ifHsstd)
      read(5,*) Hs_std(1:NP)
      close(5)

      open(16,file=ifGtPCsstd)
      do i=1,30
        read(16,*) GtPCs_std(i)
      enddo
      close(16)

      open(7,file=ifHsmean)
      do i=1,NP
        read(7,*) Hs_mean(i)
        Hsmax(i)=Hs_mean(i)+35*Hs_std(i)
      enddo
      close(7)

      cnt=0
      oHs=0.
      oHs1=0.
      cnt1=0
      open(21,file=ofname)
      open(91,file=trim(ofname)//'.miss')
      open(92,file=trim(ofname)//'.flag')
C      open(93,file=trim(ofname)//'.lm_miss_flag')
      open(11,file=ifPt,iostat=fstatus)
      if(fstatus.ne.0) then
        print*,trim(ifPt)," not exist"
        stop
      endif
      open(12,file=ifGt,iostat=fstatus)
      if(fstatus.ne.0) then
        print*,trim(ifGt)," not exist"
        stop
      endif
      open(13,file=ifPCs,iostat=fstatus)
      if(fstatus.ne.0) then
        print*,trim(ifPCs)," not exist"
        stop
      endif
      open(14,file=ifGtPCs,iostat=fstatus)
      if(fstatus.ne.0) then
        print*,trim(ifGtPCs)," not exist"
        stop
      endif
      
C     print*,trim(ifPt)
C     print*,trim(ifGt)
C     print*,trim(ifPCs)
C     print*,trim(ifGtPCs)
80    read(11,*,end=100) yy1,mm1,dd1,hh1,itmp
      read(11,*) Pt(1:NP)
      if(itmp.ne.NP) then
        print*, trim(ifPt),yy1,mm1,dd1,hh1,itmp,"error"
        stop
      endif
      read(12,*,end=100) yy2,mm2,dd2,hh2,itmp
      read(12,*) Gt(1:NP)
      if(itmp.ne.NP) then
        print*, trim(ifGt),yy2,mm2,dd2,hh2,itmp,"error"
        stop
      endif
      read(13,*,end=100) yy3,mm3,dd3,hh3,PCs(1:30)
C     read(13,*) PCs(1:60)
      if(yy1*1000000+mm1*10000+dd1*100+hh1.ne.
     &   yy2*1000000+mm2*10000+dd2*100+hh2.or.
     &   yy1*1000000+mm1*10000+dd1*100+hh1.ne.
     &   yy3*1000000+mm3*10000+dd3*100+hh3) then
        print*,"PtPCs",yy1*1000000+mm1*10000+dd1*100+hh1,
     &         yy2*1000000+mm2*10000+dd2*100+hh2,
     &         yy3*1000000+mm3*10000+dd3*100+hh3,"year error"
        stop
      endif
      read(14,*,end=100) yy3,mm3,dd3,hh3,GtPCs(1:30)
      if(yy1*1000000+mm1*10000+dd1*100+hh1.ne.
     &   yy2*1000000+mm2*10000+dd2*100+hh2.or.
     &   yy1*1000000+mm1*10000+dd1*100+hh1.ne.
     &   yy3*1000000+mm3*10000+dd3*100+hh3) then
        print*,"trGtPCs",yy1*1000000+mm1*10000+dd1*100+hh1,
     &         yy2*1000000+mm2*10000+dd2*100+hh2,
     &         yy3*1000000+mm3*10000+dd3*100+hh3,"year error"
        stop
      endif
      do i=1,30
        PCs(i)=PCs(i)/PCs_std(i)
      enddo
      do i=1,30
        GtPCs(i)=GtPCs(i)/GtPCs_std(i)
      enddo
      rtmp3(3:32)=PCs(1:30)
      rtmp3(33:62)=GtPCs(1:30)
      idiff=(mm1-sm)*120+(dd1-1)*4+int(hh1/6)+1
      do ith=1,NP
        PHs(ith)=0.
        rtmp3(1)=Gt(ith)
        rtmp3(2)=Pt(ith)
C       idiff=(mm1-1)*120+(dd1-1)*4+int(hh1/6)+1
        oHs_lm(ith)=coefs_lm(ith,1)
        do i=1,NX(ith)
          oHs_lm(ith)=oHs_lm(ith)+coefs_lm(ith,i+1)*rtmp3(idx(ith,i))
        enddo
        if(Hslambda(ith)*oHs_lm(ith)*Hs_std(ith)+1..lt.0.) then
          tHs_lm(ith)=0.
        else
          tHs_lm(ith)=(Hslambda(ith)*oHs_lm(ith)*Hs_std(ith)+1.)**
     &             (1./Hslambda(ith))
        endif

        if(idiff.le.p(ith)) then
          oHs(ith)=oHs_lm(ith)
        else
          oHs(ith)=coefs(ith,1)
          do i=1,NX(ith)
            oHs(ith)=oHs(ith)+coefs(ith,(i+1))*rtmp3(idx(ith,i))
          enddo
          if(p(ith).gt.0) then
            oHsFlg=0
            do zth=1,p(ith)
              if(oHs1(zth,ith).lt.-999990.) oHsFlg=1
            enddo
            if(oHsFlg.eq.1) then
              oHs(ith)=oHs_lm(ith)
            else
              do i=1,p(ith)
                oHs(ith)=oHs(ith)+coefs(ith,NX(ith)+1+i)*oHs1(i,ith)
              enddo
            endif
          endif
        endif
        if(abs(Hslambda(ith)).le.1E-5) then
          tHs(ith)=exp(oHs(ith)*Hs_std(ith))
        else
          if(Hslambda(ith)*oHs(ith)*Hs_std(ith)+1..lt.0.) then
            tHs(ith)=0.
          else
            tHs(ith)=(Hslambda(ith)*oHs(ith)*Hs_std(ith)+1.)**
     &             (1./Hslambda(ith))
          endif
        endif

C         rtmp=transver(oHs(ith)*Hs_std(ith),Hslambda(ith))
C keep oHs1 as un-transfer-back values
        do i=1,29
          oHs1(i+1,ith)=oHs1(i,ith)
        enddo
        oHs1(1,ith)=oHs(ith)

        if(abs(tHs(ith)).gt.20.) then
          write(92,'(i5,3i3,i5,2f14.3,2f8.3)')
     &          yy1,mm1,dd1,hh1,ith,
     &          oHs(ith),tHs(ith),
     &          Hs_mean(ith),Hs_std(ith)
        endif
        PHs(ith)=tHs(ith)
      enddo
      write(21,'(i5,3i3)') yy1,mm1,dd1,hh1
      write(21,'(5f20.3)') PHs(1:NP)
C     oHs1=oHs

      goto 80

100   continue
c      print*,trim(carea)," ",trim(csea)," total missing#:", cnt1
      write(91,'(10i6)') cnt1(1:NP)
      end
