       program cdft_cenergy

!  This program takes a full transient simulation, a reference
!  observation file
!  and bias correct data using 20-year periods moving by stes of 10
!  years
!  The program works only for data of the same month
!  The input file is used as output file and the input non bias
!  corrected
!  variable is replaced by the bias corrected one

       USE netcdf

       parameter(ndimsmax=4)

       real, allocatable, dimension(:,:,:,:) :: z500mod,z500mod1,z500era,txrthm,txrthx,txsthm,txsthx,txdym,txdyx,txrdym,txrdyx
       real, allocatable, dimension(:,:,:)   :: txmod,txmeanr,txmeans,txmeand,txmynd,txmaxr,txmaxs,txmaxd,txmn,txxn
       real, allocatable, dimension(:,:,:)   :: txera
       real, allocatable, dimension(:,:)     :: trrthm,trrdym,trsthm,trendm,trrthx,trrdyx,trsthx,trendx,trdynm,trdynx
       real, allocatable, dimension(:,:)     :: serthm,serdym,sesthm,seendm,serthx,serdyx,sesthx,seendx,sedynm,sedynx
       real, allocatable, dimension(:,:)     :: wi,we
       double precision, allocatable, dimension(:)       :: xx,xt    ! Model longitudes
       double precision, allocatable, dimension(:)       :: yy,yt    ! Model latitudes

       real cor(10000),cor1(10000),thermo(10000),trnd(10000),sed(10000),corref(10000)
       real gmstmod(1950:2022),gmstera(1950:2022),gmst(10000)
       real gm(10000),gme(10000),ww(10000)

       integer varid, ifileid, timid,timido,ofileid,obsfileid,xid,yid
       integer dimids(ndimsmax), lendim(ndimsmax),lendimo(ndimsmax)
       integer id(10000),ink(10000),iyi(10000),iyf(10000),idi(10000),idf(10000),ian(10000)
       integer nimax(1950:2022),ipatts(10000),ninc(1950:2022),nexc(1950:2022)
       integer, allocatable, dimension(:,:,:)       :: idmax

       character*256 z500modname,z500eraname,modname,cnum,fnam
       character*256 txmodname,txeraname,sgsatmodname,sgsateraname
       character*256 namedim(ndimsmax)
       character*256 filei,fileobs,fileout
       character*256 tunit,tcal,tnam,lonnam,latnam

       iread=1
       read *,modname,z500modname,z500eraname,txmodname,txeraname,sgsatmodname,sgsateraname

       nyears = 2022-1950+1
       fan = 0.99      ! Take analogs beyond the 1-fan fraction of best ones
       nna = 3
       nreg = nyears*nna

       print *,'++ USING TANALOGS with IREAD=',iread

       open(11,file=sgsatmodname)
       open(12,file=sgsateraname)
       do i=1950,2022
          read(11,*)iy,gmstmod(i)
          read(12,*)iy,gmstera(i)
       enddo
       close(11)
       close(12)

! Initializations

       j=0
       do i=1,nna
          do iy=1950,2022
             j=j+1
             gm(j)=gmstmod(iy)
             gme(j)=gmstera(iy)
          enddo
       enddo

!  Reading the z500 model anomaly netcdf file

       ierr = nf90_open(z500modname, NF90_WRITE,ifileid)
       if(modname.ne.'ERA5')ierr = nf90_inq_varid(ifileid,'stream', varid)
       if(modname.eq.'ERA5')ierr = nf90_inq_varid(ifileid,'stream', varid)
       ierr = nf90_Inquire_Variable(ifileid,varid,ndims=ndims,dimids=dimids)
!      print *,'++ MODEL File DIMS=',ndims
       do idim = 1,ndims
          ierr = nf90_Inquire_Dimension(ifileid,dimids(idim), &
    &     name = namedim(idim), len = lendim(idim))
!         print *,'GCM DIMENSION(',idim,')=',lendim(idim)
       enddo

       allocate(z500mod(lendim(1),lendim(2),1,lendim(ndims)))
       if(modname.eq.'ERA5')allocate(z500mod1(lendim(1),lendim(2),lendim(3),lendim(ndims)))
       allocate(xx(lendim(1)))
       allocate(yy(lendim(2)))

       nxx=lendim(1)
       nyy=lendim(2)
       nmod = lendim(ndims)

       ierr=nf90_inq_varid(ifileid,'lon',xid)
       ierr=nf90_get_var(ifileid,xid,xx)
       ierr=nf90_inq_varid(ifileid,'lat',yid)
       ierr=nf90_get_var(ifileid,yid,yy)
       if(modname.ne.'ERA5')ierr=nf90_get_var(ifileid,varid,z500mod,start=(/1,1,1,1/), &
     &        count=(/lendim(1),lendim(2),lendim(3),lendim(ndims)/))
       if(modname.eq.'ERA5')ierr=nf90_get_var(ifileid,varid,z500mod1,start=(/1,1,1,1/), &
     &        count=(/lendim(1),lendim(2),lendim(3),lendim(ndims)/))
       if(modname.eq.'ERA5') then
         do i=1,lendim(1)
            do j=1,lendim(2)
               do k=1,lendim(ndims)
                  z500mod(i,j,1,k)=z500mod1(i,j,1,k)
               enddo
            enddo
         enddo
       endif

       ierr = nf90_open(txmodname, NF90_NOWRITE,ifileid)
       if(modname.ne.'ERA5')ierr = nf90_inq_varid(ifileid,'tasmax', varid)
       if(modname.eq.'ERA5')ierr = nf90_inq_varid(ifileid,'tmax', varid)
       ierr = nf90_Inquire_Variable(ifileid,varid,ndims=ndims,dimids=dimids)
       do idim = 1,3
          ierr = nf90_Inquire_Dimension(ifileid,dimids(idim), &
    &     name = namedim(idim), len = lendim(idim))
       enddo
       nxt=lendim(1)
       nyt=lendim(2)
       ntimes=lendim(4)
       if(ntimes.ne.nmod) print *,'*** ERROR NTIMES, NMOD different',ntimes,nmod

       allocate(txmod(nxt,nyt,ntimes))
       allocate(txmeanr(nna,nxt,nyt))
       allocate(txmeans(nna,nxt,nyt))
       allocate(txmeand(1950:2022,nxt,nyt))
       allocate(txmaxr(nna,nxt,nyt))
       allocate(txmaxs(nna,nxt,nyt))
       allocate(txmaxd(1950:2022,nxt,nyt))
       allocate(txrthm(1950:2022,nna,nxt,nyt))
       allocate(txrthx(1950:2022,nna,nxt,nyt))
       allocate(txsthm(1950:2022,nna,nxt,nyt))
       allocate(txsthx(1950:2022,nna,nxt,nyt))
       allocate(txrdym(1950:2022,nna,nxt,nyt))
       allocate(txrdyx(1950:2022,nna,nxt,nyt))
       allocate(txdym(1950:2022,nna,nxt,nyt))
       allocate(txdyx(1950:2022,nna,nxt,nyt))
       allocate(txmn(1950:2022,nxt,nyt))
       allocate(txxn(1950:2022,nxt,nyt))
       allocate(idmax(1950:2022,nxt,nyt))
       allocate(trrthm(nxt,nyt))
       allocate(trrdym(nxt,nyt))
       allocate(trsthm(nxt,nyt))
       allocate(trendm(nxt,nyt))
       allocate(trrthx(nxt,nyt))
       allocate(trrdyx(nxt,nyt))
       allocate(trsthx(nxt,nyt))
       allocate(trendx(nxt,nyt))
       allocate(trdynm(nxt,nyt))
       allocate(trdynx(nxt,nyt))
       allocate(sedynm(nxt,nyt))
       allocate(sedynx(nxt,nyt))
       allocate(sesthm(nxt,nyt))
       allocate(sesthx(nxt,nyt))
       allocate(serthm(nxt,nyt))
       allocate(serdym(nxt,nyt))
       allocate(seendm(nxt,nyt))
       allocate(seendx(nxt,nyt))
       allocate(serthx(nxt,nyt))
       allocate(serdyx(nxt,nyt))
       allocate(xt(nxt))
       allocate(yt(nyt))
       allocate(wi(1950:2022,nna))
       allocate(we(1950:2022,nna))

       ierr=nf90_inq_varid(ifileid,'lon',xid)
       ierr=nf90_get_var(ifileid,xid,xt)
       ierr=nf90_inq_varid(ifileid,'lat',yid)
       ierr=nf90_get_var(ifileid,yid,yt)
       ierr=nf90_get_var(ifileid,varid,txmod,start=(/1,1,1/),count=(/nxt,nyt,ntimes/))

       PRINT*,'++ All Model data now read'

       if(nmod.eq.6716)ilseas=92
       if(nmod.eq.6570)ilseas=90
       ilsease=92

!  Reading the class file
!  and counting the dates in each year

       do iy=1950,2022
          nexc(iy)=0
          ninc(iy)=0
       enddo
       open(78,file='PATTERNS')
       do nm=1,nmod
          read(78,*),i,ipatts(i)
          iy=int((nm-1)/ilseas)+1950
          if(ipatts(i).eq.0)nexc(iy)=nexc(iy)+1
          if(ipatts(i).eq.1)ninc(iy)=ninc(iy)+1
       enddo
       do iy=1950,2022
          if(nexc(iy).le.2)print *,iy,nexc(iy),' NOT ENOUGH excluded dates'
          if(ninc(iy).le.2)print *,iy,ninc(iy),' NOT ENOUGH included dates'
          do i=1,nna
             wi(iy,i)=ninc(iy)
             we(iy,i)=nexc(iy)
          enddo
!         print *,iy,ninc(iy),nexc(iy)
       enddo

!  Reading the observation file

       ierr = nf90_open(z500eraname, NF90_NOWRITE,obsfileid)
       ierr = nf90_inq_varid(obsfileid,'stream',varid)
       ierr = nf90_Inquire_Variable(obsfileid,varid,dimids=dimids)
       DO idim = 1,4
          ierr = nf90_Inquire_Dimension(obsfileid,dimids(idim), &
    &     name = namedim(idim), len = lendimo(idim))
!         print *,'++ OBS DIMENSION(',idim,')=',lendimo(idim)
       ENDDO

       nera=lendimo(4)
!      print *,'++ NMOD, NERA = ',nmod,nera

       allocate(z500era(lendimo(1),lendimo(2),lendimo(3),lendimo(4)))

!  Reading Observation cordinates

       ierr = nf90_get_var(obsfileid, varid,z500era)
!      do i=1,lendimo(1)
!         do j=1,lendimo(2)
!            do k=1,lendimo(4)
!               print *,z500era(i,j,1,k)
!            enddo
!         enddo
!      enddo

       ierr = nf90_open(txeraname, NF90_NOWRITE,ifileid)
       ierr = nf90_inq_varid(ifileid,'tmax', varid)
       ierr = nf90_Inquire_Variable(ifileid,varid,ndims=ndims,dimids=dimids)
       do idim = 1,3
          ierr = nf90_Inquire_Dimension(ifileid,dimids(idim), &
    &     name = namedim(idim), len = lendim(idim))
       enddo
       nxt=lendim(1)
       nyt=lendim(2)
       ntimesera=lendim(3)

       if(ntimesera.ne.nera) print *,'*** ERROR NTIMESERA, NERA different',ntimesera,nera

       allocate(txera(nxt,nyt,ntimesera))
       ierr=nf90_get_var(ifileid,varid,txera,start=(/1,1,1/),count=(/nxt,nyt,ntimesera/))

       PRINT*,'++ All data now read'

! Initializations

       do iym=1950,2022
          do nx=1,nxt
             do ny=1,nyt
                txmn(iym,nx,ny)=0.
                txxn(iym,nx,ny)=0.
                txmeand(iym,nx,ny)=0.
                txmaxd(iym,nx,ny)=0.
                do i=1,nna
                   txrthm(iym,i,nx,ny)=0.
                   txrthx(iym,i,nx,ny)=-1.e33
                   txsthm(iym,i,nx,ny)=0.
                   txsthx(iym,i,nx,ny)=-1.e33
                   txrdym(iym,i,nx,ny)=0.
                   txrdyx(iym,i,nx,ny)=-1.e33
                   txdym(iym,i,nx,ny)=0.
                   txdyx(iym,i,nx,ny)=-1.e33
                enddo
             enddo
          enddo
       enddo
       open(71,file='PATTRENDS/cormins')

       fnam='PATTRENDS/txmx.' // trim(modname) // '.yearly.dat'
       open(61,file=trim(fnam))

       nm=0
       do iym=1950,2022
          do idm=1,ilseas
             nm=nm+1
             write(cnum,'(i0)')nm
             fnam='ancors.' // trim(modname) // '/ans.' // trim(cnum) //'.dat'
             open(51,file=trim(fnam))
             nana=int((1-fan)*nmod)
             nanap=nana
             call selan(iread,nmod,nanap,ilseas,nxx,nyy,z500mod(1,1,1,nm),z500mod,ian,iyf,idf,corm)

             write(71,189)nm,iym,idm,corm,gmstmod(iyf(1))
189          format(i6,i6,i6,f8.3,f8.3)
             do ny=1,nyt
                do nx=1,nxt
                   do i=1,nana
                      thermo(i)=txmod(nx,ny,ian(i))
                      gmst(i)=gmstmod(iyf(i))
                      ww(i)=1.
                   enddo
                   call reg(thermo,gmst,nana,ww,tthm,se)
                   i=0
! Some things in the loop are really stupid, comes from an old test, but it still works
                   do j=1,nna
                      i=i+1
                      ina=ian(j)
                      txdyn=txmod(nx,ny,ina)-tthm*gmstmod(iyf(j))
                      if(ipatts(nm).eq.1) then
                        txsthm(iym,i,nx,ny)=txsthm(iym,i,nx,ny)+txdyn/float(ninc(iym))
                        if(txdyn.gt.txsthx(iym,i,nx,ny))txsthx(iym,i,nx,ny)=txdyn
                      else
                        txdym(iym,i,nx,ny)=txdym(iym,i,nx,ny)+txdyn/float(nexc(iym))
                        if(txdyn.gt.txdyx(iym,i,nx,ny))txdyx(iym,i,nx,ny)=txdyn
                      endif
                      if(mod(iym,10).eq.1.and.i.eq.3.and.nx.eq.25.and.ny.eq.27.and.idm.eq.ilseas)print *,iym,i,txdym(iym,i,nx,ny),txdyx(iym,i,nx,ny),corm
                   enddo
7777               continue
                   if(idm.eq.ilseas.and.ninc(iym).gt.0) then
                     do i=1,nna
                        txmn(iym,nx,ny)=txmn(iym,nx,ny)+txdym(iym,i,nx,ny)/nna
                        txxn(iym,nx,ny)=txxn(iym,nx,ny)+txdyx(iym,i,nx,ny)/nna
                        txmeand(iym,nx,ny)=txmeand(iym,nx,ny)+txsthm(iym,i,nx,ny)/nna
                        txmaxd(iym,nx,ny)=txmaxd(iym,nx,ny)+txsthx(iym,i,nx,ny)/nna
                     enddo 
                     write(61,986),real(xt(nx)),real(yt(ny)),txmn(iym,nx,ny),txxn(iym,nx,ny),txmeand(iym,nx,ny),txmaxd(iym,nx,ny)
                   endif
                   if(idm.eq.ilseas.and.ninc(iym).eq.0) then
                     do i=1,nna
                        txmn(iym,nx,ny)=txmn(iym,nx,ny)+txdym(iym,i,nx,ny)/nna
                        txxn(iym,nx,ny)=txxn(iym,nx,ny)+txdyx(iym,i,nx,ny)/nna
                     enddo
                     txmeand(iym,nx,ny)=-999.
                     txmaxd(iym,nx,ny)=-999.
                     write(61,986),real(xt(nx)),real(yt(ny)),txmn(iym,nx,ny),txxn(iym,nx,ny),txmeand(iym,nx,ny),txmaxd(iym,nx,ny)
                   endif
                enddo
             enddo
          enddo
       enddo

       close(61)

       do nx=1,nxt
          do ny=1,nyt
             call reg(txsthm(1950,1,nx,ny),gm,nreg,wi(1950,1),trsthm(nx,ny),sesthm(nx,ny))
             call reg(txsthx(1950,1,nx,ny),gm,nreg,wi(1950,1),trsthx(nx,ny),sesthx(nx,ny))
             call reg(txdym(1950,1,nx,ny),gm,nreg,we(1950,1),trdynm(nx,ny),sedynm(nx,ny))
             call reg(txdyx(1950,1,nx,ny),gm,nreg,we(1950,1),trdynx(nx,ny),sedynx(nx,ny))
          enddo
       enddo

8766   continue
       
       fnam='PATTRENDS/txm.' // trim(modname) // '.trends.dat'
       open(61,file=trim(fnam))
       do ny=1,nyt
          do nx=1,nxt
             write(61,986),real(xt(nx)),real(yt(ny)),trendm(nx,ny),trrthm(nx,ny),trsthm(nx,ny),trdynm(nx,ny),trrdym(nx,ny) &
    &                                               ,seendm(nx,ny),serthm(nx,ny),sesthm(nx,ny),sedynm(nx,ny),serdym(nx,ny)
          enddo
       enddo
       close(61)
       fnam='PATTRENDS/txx.' // trim(modname) // '.trends.dat'
       open(61,file=trim(fnam))
       do ny=1,nyt
          do nx=1,nxt
             write(61,986),real(xt(nx)),real(yt(ny)),trendx(nx,ny),trrthx(nx,ny),trsthx(nx,ny),trdynx(nx,ny),trrdyx(nx,ny) &
    &                                               ,seendx(nx,ny),serthx(nx,ny),sesthx(nx,ny),sedynx(nx,ny),serdyx(nx,ny)
          enddo
       enddo
       close(61)
986    format(f9.3,f9.3,12e12.3)

       close(71)
       end

       function corr(nxx,nyy,x,y)

       real x(nxx,nyy),y(nxx,nyy)

       k=0
       zom=0.
       zem=0.
       zov=0.
       zev=0.
       zoze=0.
       do nx=1,nxx
          do ny=1,nyy
             zom=zom+x(nx,ny)
             zem=zem+y(nx,ny)
             zov=zov+x(nx,ny)*x(nx,ny)
             zev=zev+y(nx,ny)*y(nx,ny)
             zoze=zoze+x(nx,ny)*y(nx,ny)
             k=k+1
          enddo
       enddo
       zom=zom/k
       zem=zem/k
       zov=zov/k-zom*zom
       zev=zev/k-zem*zem
       zoze=zoze/k-zom*zem
       corr=zoze/sqrt(zov*zev)

       end

       subroutine reg(y,x,n,w,r,se)

       real y(n),x(n),w(n)

       k=0
       zom=0.
       zem=0.
       zov=0.
       zev=0.
       zoze=0.
       wm=0.
       do i=1,n
          zom=zom+x(i)*w(i)
          zem=zem+y(i)*w(i)
          zov=zov+x(i)*x(i)*w(i)
          zev=zev+y(i)*y(i)*w(i)
          zoze=zoze+x(i)*y(i)*w(i)
          wm=wm+w(i)
          k=k+1
       enddo
       zom=zom/wm
       zem=zem/wm
       zov=zov/wm-zom**2
       zev=zev/wm-zem**2
       zoze=zoze/wm-zom*zem
       r=zoze/zov
       se=sqrt((zev/zov-r*r)/(k-2.))

!      zom=zom/k
!      zem=zem/k
!      zov=zov/k-zom**2
!      zev=zev/k-zem**2
!      zoze=zoze/k-zom*zem
!      r=zoze/zov
!      se=sqrt((zev/zov-r*r)/(k-2.))

       end

       subroutine hpsort_eps_epw (n, ra, ind)
!---------------------------------------------------------------------
! sort an array ra(1:n) into ascending order using heapsort algorithm,
! and considering two elements being equal if their values differ
! for less than "eps".
! n is input, ra is replaced on output by its sorted rearrangement.
! create an index table (ind) by making an exchange in the index array
! whenever an exchange is made on the sorted data array (ra).
! in case of equal values in the data array (ra) the values in the
! index array (ind) are used to order the entries.
! if on input ind(1)  = 0 then indices are initialized in the routine,
! if on input ind(1) != 0 then indices are assumed to have been
!                initialized before entering the routine and these
!                indices are carried around during the sorting process
!
! no work space needed !
! free us from machine-dependent sorting-routines !
!
! adapted from Numerical Recipes pg. 329 (new edition)
!
!      use kinds, ONLY : DP
       implicit none
!-input/output variables
       integer   :: n
       real      :: eps
       integer   :: ind (n)
       real      :: ra (n)
!-local variables
       integer   :: i, ir, j, l, iind
       real      :: rra
!
! Take eps as constant eps=1e-20
       eps = 1.e-20
!
! initialize index array
!      IF (ind (1) .eq.0) then
          DO i = 1, n
             ind (i) = i
          ENDDO
!      ENDIF
! nothing to order
       IF (n.lt.2) return
! initialize indices for hiring and retirement-promotion phase
       l = n / 2 + 1

       ir = n

  sorting: do

! still in hiring phase
       IF ( l .gt. 1 ) then
          l    = l - 1
          rra  = ra (l)
          iind = ind (l)
! in retirement-promotion phase.
       ELSE
! clear a space at the end of the array
          rra  = ra (ir)
!
          iind = ind (ir)
! retire the top of the heap into it
          ra (ir) = ra (1)
!
          ind (ir) = ind (1)
! decrease the size of the corporation
          ir = ir - 1
! done with the last promotion
          IF ( ir .eq. 1 ) then
! the least competent worker at all !
             ra (1)  = rra
!
             ind (1) = iind
             exit sorting
          ENDIF
       ENDIF
! wheter in hiring or promotion phase, we
       i = l
! set up to place rra in its proper level
       j = l + l
!
       DO while ( j .le. ir )
          IF ( j .lt. ir ) then
! compare to better underling
             IF ( hslt( ra (j),  ra (j + 1) ) ) then
                j = j + 1
!else if ( .not. hslt( ra (j+1),  ra (j) ) ) then
! this means ra(j) == ra(j+1) within tolerance
!  if (ind (j) .lt.ind (j + 1) ) j = j + 1
             ENDIF
          ENDIF
! demote rra
          IF ( hslt( rra, ra (j) ) ) then
          ra (i) = ra (j)
          ind (i) = ind (j)
          i = j
          j = j + j
!else if ( .not. hslt ( ra(j) , rra ) ) then
!this means rra == ra(j) within tolerance
! demote rra
! if (iind.lt.ind (j) ) then
!    ra (i) = ra (j)
!    ind (i) = ind (j)
!    i = j
!    j = j + j
! else
! set j to terminate do-while loop
!    j = ir + 1
! endif
! this is the right place for rra
          ELSE
! set j to terminate do-while loop
             j = ir + 1
          ENDIF
       ENDDO
       ra (i) = rra
       ind (i) = iind

       END DO sorting
       contains

!  internal function
!  compare two real number and return the result

       logical function hslt( a, b )
         REAL     :: a, b
         IF( abs(a-b) <  eps ) then
           hslt = .false.
         ELSE
           hslt = ( a < b )
         end if
       end function hslt

!
       end subroutine hpsort_eps_epw

       subroutine selan(iread,nmod1,nana,ilseas,nxx,nyy,z500mod,z500mod1,ian,iyf,idf,cormin)

       real z500mod(nxx,nyy),z500mod1(nxx,nyy,nmod1),cor(10000)
       integer id(10000),iyi(10000),idi(10000),ian(nmod1),iyf(nmod1),idf(nmod1)

       do ne=1,nmod1
          iye=(ne-1)/ilseas+1950
          ide=mod(ne-1,ilseas)+1
          if(iread.eq.1) then
            read(51,*),ii,cor(ne)
          else
            cor(ne)=corr(nxx,nyy,z500mod(1,1),z500mod1(1,1,ne))
            write(51,*),ne,cor(ne)
          endif
          iyi(ne)=iye
          idi(ne)=ide
          id(ne)=ne
       enddo
       close(51)
       call hpsort_eps_epw(nmod1,cor,id)

! Selection of a set of best analogs separated by at least 5 days
! and removing cor=1 (self analogue)

       ike=nmod1
       if(cor(nmod1).gt.0.999) then
         iyf0=iyi(id(nmod1))
         idf0=idi(id(nmod1))
       else
         iyf0=-999
         idf0=-999
       endif
       do i=1,nana
          do j=ike,1,-1
             iyj=iyi(id(j))
             idj=idi(id(j))
             if(iyf0.eq.iyj.and.abs(idf0-idj).le.5) goto 6433
             do l=1,i-1
                if(iyf(l).eq.iyj.and.abs(idf(l)-idj).le.5) goto 6433
             enddo
             ian(i)=id(j)
             iyf(i)=iyj
             idf(i)=idj
             goto 6434
6433         continue
          enddo
6434      continue
          ike=j
          cormin=cor(ike)
       enddo

       end
