cc This program is used to look up the equal water level paired dataset and prepare para.tpl, obs.tpl, *pst, input.dat, and periods.dat .
cc In the 1.4 version the monthly discharge values are no longer repeated in different years. Instead there are one independent 
cc DM value for each month within the time series.

c define variables 
        USE DFLIB  !调用IFPORT
        USE DFWIN  !调用IFPORT

        implicit none
        real,ALLOCATABLE:: P(:),WL(:,:),R(:,:),PTOTAL(:),PWL(:)
        real weightwl,weightp,tempwl,pp1,pp2,pp3,pp4,pp5,interval,diff
        integer i,j,k,m,n,q,site,gauge
        integer datanum,wellnum,periodsnum,gaugenum
        integer,ALLOCATABLE:: PSTART(:),PEND(:),duration(:)
        integer,ALLOCATABLE:: date1(:),month(:)
        integer,ALLOCATABLE:: active(:)
        integer GAUGEDIR,WELLDIR,MONTHCOUNT
        integer obsnum,paranum,active1
        LOGICAL fexist
        Character CASENAM*58,dummy*158,DIR*100
        Character tt,currentline*79,obsnam*10,dmnam*10,currentline2*79
        Character,ALLOCATABLE::  wellnam(:)*100
        Character,ALLOCATABLE::  RAINGAUGE(:)*100,RECESSION(:)*100
        Character,ALLOCATABLE::  GAUGENAM(:)*100,CONFIG(:)*100



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  READ PRICIPITATION AND WATER LEVEL DATA FILE AND Recession data file       C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
        write (*,*) "请输入NAM文件的文件名"
        read *, CASENAM

c        Call GETARG(1,CASENAM)
c        read(CASENAM,*)
        open (5,file=trim(CASENAM))
            read (5,*) dummy          
            read (5,*) gaugenum
            read (5,*) dummy          
            Allocate(GAUGENAM(gaugenum))
            Allocate(RAINGAUGE(gaugenum))
            Allocate(RECESSION(gaugenum))
            Allocate(CONFIG(gaugenum))
            read (5,*) (GAUGENAM(i),i=1,gaugenum)
          do i=1,gaugenum
            read (5,*) dummy          
            read (5,*) RAINGAUGE(i)
            read (5,*) RECESSION(i)
            read (5,*) CONFIG(i)
          end do
        close (5)


C LOOP FOR GAUGE I
       DO gauge=1,gaugenum
C Create directory for current raigauge
           PRINT *,GAUGENAM(gauge)
           GAUGEDIR=MAKEDIRQQ(TRIM(GAUGENAM(gauge)))

C Read in precipitation and water level data
        open (10,FILE=trim(RAINGAUGE(gauge)))
           read (10,*) datanum,wellnum,monthcount
           print *,datanum,wellnum

            Allocate(P(datanum))
            Allocate(WL(wellnum,datanum))
            Allocate(R(wellnum,datanum))
            Allocate(PTOTAL(datanum))
            Allocate(PWL(datanum))
            Allocate(PSTART(datanum))
            Allocate(PEND(datanum))
            Allocate(duration(datanum))
            Allocate(date1(datanum))
            Allocate(month(datanum))
            Allocate(active(wellnum))
            Allocate(wellnam(wellnum))
            read (10,*) tt,tt,tt,(wellnam(i),i=1,wellnum)
          do i=1,datanum
            read (10,*) date1(i),month(i),P(i),(WL(j,i),j=1,wellnum)
          end do
        close (10)
 
        open (20,FILE=trim(RECESSION(gauge)))
            read (20,*) datanum,wellnum
            read (20,*) tt,tt,(wellnam(i),i=1,wellnum)
          do i=1,datanum
            read (20,*) date1(i),month(i),(R(j,i),j=1,wellnum)
          end do
        close (20)
     
C Judge whether control file is exit.
       INQUIRE( File = trim(GAUGENAM(gauge))//'.ctr',EXIST = FEXIST)
        
      IF(.not.FEXIST) THEN  
        open (25,FILE=trim(GAUGENAM(gauge))//'.ctr')
            write (25,*) "ID   WELLNAME     RUN-CTR"
          do i=1,wellnum
            active1=1
            write (25,600) i,adjustr(trim(wellnam(i))),active1
	      active(i)=active1
          end do
        close (25)
       ELSE 
        open (25,FILE= trim(GAUGENAM(gauge))//'.ctr')
            read (25,*) dummy
          do i=1,wellnum
            read (25,*) tt,wellnam(i),active(i)
          end do
        close (25)
      ENDIF 

c LOOK UP DATA PAIRS FOR WELLnum is site
       DO site=1,wellnum
C Create directory for current well
      WELLDIR=MAKEDIRQQ(TRIM(GAUGENAM(gauge))//'\'//TRIM(wellnam(site)))

       OPEN (95,FILE='currentw.lst',status='replace')
	    write (95,*) GAUGENAM(gauge),wellnam(site)
	 close (95)

C judge whether to active simulation for wellnum(site)
        if (active(site).eq.1) then

c define the wl period interval and diff standard

       OPEN (30,FILE=trim(CONFIG(gauge)))
        read(30,500,advance='YES') currentline 
        read(30,*)interval,diff
          q=1
          do m=1,datanum-interval-1
            do n=m+interval,datanum ! the interval between two WL values
               if (abs(WL(site,n)-WL(site,m)).lt.diff) then
                  PSTART(q)=date1(m)
                  PEND(q)=date1(n)
                  duration(q)=date1(n)-date1(m)
                  PWL(q)=WL(site,m)
                  PTOTAL(q)=0
                  do i=m,n-1
                     PTOTAL(q)=PTOTAL(q)+P(i)
                  end do
                     q=q+1
                  goto 1000
               endif
            end do
1000      end do
          periodsnum=q-1
          obsnum=datanum+periodsnum
	 PRINT *,wellnam(site),obsnum


c sort paired WL data according to the lengthen of duration
        do i=1,periodsnum-1
           do j=i+1,periodsnum-1
          if (duration(i) .gt. duration(j)) then
                 pp1=duration(i)
                 pp2=PSTART(i)
                 pp3=PEND(i)
                 pp4=PTOTAL(i)

                 duration(i)=duration(j)
                 PSTART(i)=PSTART(j)
                 PEND(i)=PEND(j)
                 PTOTAL(i)=PTOTAL(j)

                 duration(j)=pp1
                 PSTART(j)=pp2
                 PEND(j)=pp3
                 PTOTAL(j)=pp4
              endif
            end do
         end do

500     format(a79)
100     format(24x,100(3x,5x))
200     format(2i8,100f8.0)
400     format(2i4)
600     format(i8,a30,2i8)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C                    WRITE PEST INPUTS FILES                                  C  
C    case.pst   input.dat  periods.dat   obs.tpl   para.tpl                   C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C WRITE PST FILE
c OPEN (30,FILE=trim(CONFIG(gauge)))
C           read(30,500,advance='YES') currentline 
C           read(30,*)paranum
		 paranum=MONTHCOUNT+1
	     read(30,500,advance='YES') currentline
      DIR = trim(GAUGENAM(gauge))//'\'//trim(wellnam(site))
      OPEN (40,FILE=trim(DIR)//'\case.pst')
           write(40,'(A)')'pcf'
           write(40,'(A)')'* control data'
           write(40,*)'restart  estimation'
           write(40,*)paranum,obsnum,paranum,'           0           2'        
           read(30,500,advance='YES') currentline 
	     write(40,*) currentline
	     read(30,500,advance='YES') currentline 
	     write(40,*) currentline
           read(30,500,advance='YES') currentline 
	     write(40,*) currentline
	     read(30,500,advance='YES') currentline 
	     write(40,*) currentline
	     read(30,500,advance='YES') currentline 
	     write(40,*) currentline
	     read(30,500,advance='YES') currentline 
	     write(40,500,advance='YES') currentline

cc write parameter group and parameter data
           write(40,'(A)')'* parameter groups'
		 currentline="    relative 0.01  0.0  switch  2.0 parabolic "
		 do i=1,MONTHCOUNT
             write(dmnam,"(i5.5)")i  
		   dmnam=trim('DM'//dmnam)
             write(40,*) trim(dmnam)//trim(currentline)
           end do
           write(40,'(A)')'beta  '//trim(currentline)
cc write DM optimization parameters
	     read(30,500,advance='YES') currentline 
	     read(30,500,advance='YES') currentline 
           write(40,'(A)')'* parameter data'
		 currentline2="   1.	0.	1 "
		 do i=1,MONTHCOUNT
             write(dmnam,"(i5.5)")i  
		   dmnam=trim('DM'//dmnam)
             write(40,*) trim(dmnam)//"  "//trim(currentline)//"  "//
     $trim(dmnam)//trim(currentline2)
           end do
		   currentline="beta    log factor	5	1.E-10	1.E+3	beta"
             write(40,*) trim(currentline)//trim(currentline2)
cc write observation data
             read(30,500,advance='YES') currentline 
		   read(30,*)weightwl,weightp 
             write(40,'(A)') '* observation groups'  
             write(40,*)'obs_wl'
		   write(40,*)'obs_p'
             write(40,'(A)')'* observation data'   
           do i=1,datanum
              write(obsnam,"(i5.5)")i           
              obsnam=adjustl('o'//obsnam)
              write(40,*) obsnam,WL(site,i),weightwl," obs_wl"
           end do

           do i=datanum+1,obsnum
              write(obsnam,"(i5.5)")i           
              obsnam=adjustl('o'//obsnam)
              write(40,*) obsnam,PTOTAL(i-datanum),weightp," obs_p"
           end do
             write(40,'(A)')'* model command line'
             write(40,*)'simusl.exe'       
             write(40,'(A)')'* model input/output' 
             write(40,*)trim(DIR)//'\para.tpl '//trim(DIR)//'\para.dat'           
             write(40,*)trim(DIR)//'\out.tpl '//trim(DIR)//'\out.dat'    
             write(40,'(A)')'* prior information' 
        close(30)
	  close(40)
c write observation data template file
       OPEN (50,FILE=trim(DIR)//'\out.tpl')
           write(50,'(A)')'pif #'
         do i=1,obsnum
	     write(obsnam,"(i5.5)")i           
           obsnam=trim('o'//obsnam)
           write(50,*) "l1	("//trim(obsnam)//")4:11"
         end do
       Close (50)

C WRTIE individual data file
      OPEN (60,FILE=trim(DIR)//'\input.dat')
           write (60,*) datanum,MONTHCOUNT
         do i=1,datanum
           write (60,*) date1(i),Month(i),P(i),R(site,i),WL(site,i)
        end do
      Close (60)

C WRTIE PERIODS FILE
      OPEN (70,FILE=trim(DIR)//'\periods.dat')
           write (70,*) periodsnum
        do i=1,periodsnum
           write (70,*) PSTART(i),PEND(i)
        end do
      Close (70)

C WRTIE PERIODS FILE
c      OPEN (66,FILE=trim(DIR)//'\periodsall.dat')
c           write (66,*) periodsnum
c        do i=1,periodsnum
c           write (66,*) PSTART(i),PEND(i),PTOTAL(i),duration(i),PWL(i)
c        end do
c      Close (66)

C WRTIE PARAMETER TEMPLATE FILE
      OPEN (90,FILE=trim(DIR)//'\para.tpl')
      write(90,'(A)')'ptf #'
	do i=1,MONTHCOUNT
         write(dmnam,"(i5.5)")i  
	   dmnam=trim('DM'//dmnam)
         write(90,'(A)') '# '//trim(dmnam)//' #'
      end do

      write(90,'(A)')'# beta #'
      Close (90)
       
cc call PEST excution to run simulations for each well       

	 call system('pest '//trim(DIR)//'\case')
       
cc end for the well numbere site
       ENDIF      

	 END DO
            DEAllocate(P)
            DEAllocate(WL)
            DEAllocate(R)
            DEAllocate(PTOTAL)
            DEAllocate(PWL)
            DEAllocate(PSTART)
            DEAllocate(PEND)
            DEAllocate(duration)
            DEAllocate(date1)
            DEAllocate(month)
            DEAllocate(active)
            DEAllocate(wellnam)
cc end for the gauge numbere i
       END DO
            DEAllocate(GAUGENAM)
            DEAllocate(RAINGAUGE)
            DEAllocate(RECESSION)
            DEAllocate(CONFIG)
	    
		call system('del currentw.lst')
       STOP
       END

