!!  Copyright (C) Stichting Deltares, 2005-2023.
!!
!!  This file is part of iMOD.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.
!!
MODULE MOD_MSPNETRCH

USE WINTERACTER
USE RESOURCE
USE MOD_MSPINSPECTOR_PAR
USE MOD_IDF, ONLY : IDFDEALLOCATEX,IDFGETEDGE,IDFGETLOC,IDFGETDXDY,IDFWRITE,IDFCOPY,IDFEQUAL
USE IMODVAR, ONLY : IDIAGERROR
USE MOD_MEAN_PAR, ONLY : MEAN_FYR,MEAN_TYR,MEAN_NYEAR,MEAN_NPERIOD,IBATCH,MEAN_NLAYER, &
                         MEAN_ILAYER,MEAN_FMEAN,MEAN_FTOTAL,MEAN_ISEL,MEAN_RESDIR,CFUNC, &
                         MEAN_OUTFILE,MEAN_IDFNAME
USE MOD_MEAN_CLC, ONLY : MEAN1COMPUTE
USE MOD_OSD, ONLY : OSD_OPEN
USE MOD_MANAGER_UTL
USE MOD_UTL, ONLY : UTL_IDATETOGDATE,ITIMETOCDATE,UTL_IDATETOJDATE,DIFFTIME,JD


CONTAINS
 
 !###======================================================================
 LOGICAL FUNCTION MSPNETRCHCOMPUTE(IDF) 
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(INOUT) :: IDF
 CHARACTER(LEN=256) :: FNAME,ROOT
 CHARACTER(LEN=52) :: CDATE
 INTEGER :: NYR,YCNT                              ! Number of Years, YearCounter,yearnumber
 INTEGER :: FYR,FMN,FDY,TYR,TMN,TDY,JD0,JD1,JD2   ! FromYear, FromMonth, FromDay, ToYear, ToMonth, ToDay, 3 x JulianDate
 INTEGER :: I,J,IY,IM,ID,NFILES,IFILES,NU,ICOL,IROW,NCREATE,FY,FM,FD,TY,TM,TD
 INTEGER :: IHR,IMT,ISC
 REAL(KIND=DP_KIND) :: TTIME
 INTEGER(KIND=DP_KIND) :: IDATEFULL,IDATE
 CHARACTER(LEN=256),DIMENSION(:),POINTER :: IDFNAMES
 INTEGER(KIND=DP_KIND),POINTER,DIMENSION(:) :: ITIME,ITIME_BU
 INTEGER(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: LDATES
 INTEGER,ALLOCATABLE,DIMENSION(:) :: IYEAR

 MSPNETRCHCOMPUTE=.FALSE.

 !## divide full date into separate parts    
 CALL UTL_IDATETOGDATE(MSPRCH_FYR,FYR,FMN,FDY); JD1=UTL_IDATETOJDATE(MSPRCH_FYR)
 CALL UTL_IDATETOGDATE(MSPRCH_TYR,TYR,TMN,TDY); JD2=UTL_IDATETOJDATE(MSPRCH_TYR)
 
 !## get all unique files
 ALLOCATE(ITIME(1000)); ITIME=INT(0,8); NFILES=0
 IF(LEN_TRIM(SOURCEDIR).EQ.0)THEN
  ROOT=TRIM(HEDDIR)
 ELSE
  ROOT=TRIM(SOURCEDIR)//'\HEAD'
 ENDIF
 FNAME='HEAD_*_L1.IDF'
 !## get them all (the final time array ITIME)
 IF(UTL_DIRINFO_POINTER(ROOT,FNAME,IDFNAMES,'F',CORDER='N'))THEN; ENDIF
 IF(SIZE(IDFNAMES).EQ.0)THEN 
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'No Head files for '//TRIM(FNAME)//' found in '//TRIM(ROOT), 'Error')
  RETURN
 ENDIF
 DO I=1,SIZE(IDFNAMES)
  IDATE=UTL_IDFGETDATE(IDFNAMES(I),IDATEFULL=IDATEFULL) 
  IF(IDATE.NE.0D0)THEN
   NFILES=NFILES+1
   IF(NFILES.GT.SIZE(ITIME))THEN
    ALLOCATE(ITIME_BU(SIZE(ITIME)+1000))
    DO J=1,SIZE(ITIME); ITIME_BU(J)=ITIME(J); ENDDO
    DEALLOCATE(ITIME); ITIME=>ITIME_BU
   ENDIF
   ITIME(NFILES)=IDATEFULL 
  ENDIF        
 ENDDO
 DEALLOCATE(IDFNAMES)
     
 ALLOCATE(LDATES(NFILES)); DO I=1,NFILES; LDATES(I)=ITIME(I); ENDDO; DEALLOCATE(ITIME)
 
 !## get number unique dates
 CALL UTL_GETUNIQUE_DINT(LDATES,NFILES,NU,0); NFILES=NU

 !## make sure all dates are double precision
 DO I=1,NU
  LDATES(I)=UTL_COMPLETEDATE(LDATES(I))
 ENDDO
 
 !## check: HEAD files MUST exist on daily base
 TTIME=DIFFTIME(LDATES(1),LDATES(NU))+1
 IF(TTIME/NFILES.NE.1)THEN 
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,TRIM(VTOS(NU))//' Unique HEAD Files do not have all a single day difference. One file for each day expected!', 'Error')
  RETURN
 ENDIF
 
 !## check: Minimum of 2 timesteps must be given
 IF(NFILES.EQ.1)THEN 
  WRITE(*,'(/A)') 'Function is stopped. Only 1 timestep available: '//TRIM(VTOS(LDATES(1)/INT(1000000,8)))
  WRITE(*,'(A)') 'For NETRCH calculation a minimum of 2 timesteps is needed. '  
  RETURN
 ENDIF
 
 !## create NETRCH folder
 CALL UTL_CREATEDIR(TRIM(RESDIR))

 CALL IDFNULLIFY(HEAD1) ; CALL IDFNULLIFY(HEAD2) ; CALL IDFNULLIFY(SC1) ; CALL IDFNULLIFY(QMODF)
 CALL IDFNULLIFY(NETRCH)

 !## start at time position 2 while "delta lvgwmodf" is defined as "lvgmodf (t) - lvgmodf (t-1)".
 NCREATE=0; DO IFILES=1,NFILES   

 CALL UTL_IDATETOGDATE(LDATES(IFILES),IY,IM,ID,IHR,IMT,ISC)
 WRITE(CDATE,'(I4.4,2I2.2)') IY,IM,ID
 WRITE(*,*) 'Timestep '//TRIM(CDATE)//' ...'  
 JD0=JD(IY,IM,ID)
 !## check if time is within period
 IF(JD0.GE.JD1.AND.JD0.LE.JD2)THEN  

  NCREATE=NCREATE+1
  IF(NCREATE.EQ.1)THEN
   !## if not given, calc mean Storage Coefficient 
   IF(LEN_TRIM(STOAVG).EQ.0)THEN 

    !## calc average storage coefficient
    IBATCH=1
    MEAN_FYR=INT(LDATES(1)/1000000)
    MEAN_TYR=INT(LDATES(NFILES)/1000000)
    CALL UTL_IDATETOGDATE(MEAN_FYR,FY,FM,FD) 
    CALL UTL_IDATETOGDATE(MEAN_TYR,TY,TM,TD)

    NYR=TY-FY+1
    ALLOCATE(IYEAR(NYR))
    YCNT=0; DO I=FY,TY; YCNT=YCNT+1; IYEAR(YCNT)=I; ENDDO

    MEAN_NLAYER=1
    ALLOCATE(MEAN_ILAYER(MEAN_NLAYER)) ; MEAN_ILAYER(1)=1
    ALLOCATE(MEAN_FMEAN(MAX(1,MEAN_NLAYER)),MEAN_FTOTAL(MAX(1,MEAN_NLAYER)))
    IF(LEN_TRIM(SOURCEDIR).EQ.0)THEN
     MEAN_RESDIR=TRIM(MSPDIR)//'\MSW_SC1\MSW_SC1'
    ELSE
     MEAN_RESDIR=TRIM(SOURCEDIR)//'\METASWAP\MSW_SC1\MSW_SC1'
    ENDIF
    CFUNC='MEAN'
    STOAVG =TRIM(MEAN_RESDIR)//'_'//TRIM(CFUNC)//'_'// &                                 
            TRIM(VTOS(FY))//'-'//TRIM(VTOS(FM))//'-'//TRIM(VTOS(FD))//'_to_'// &
            TRIM(VTOS(TY))//'-'//TRIM(VTOS(TM))//'-'//TRIM(VTOS(TD))//'_L'//TRIM(VTOS(MEAN_NLAYER))//'.IDF'
    WRITE(*,*) 'No Mean Storage Coefficient (MSC) file given. MSC is calculated.'  
    WRITE(*,'(A/)') '.... with files: '//TRIM(MEAN_RESDIR)//'_YYYYMMDD_L1.IDF'
    MEAN_ISEL=1; MEAN_OUTFILE=''; IF(.NOT.MEAN1COMPUTE())THEN ; ENDIF
   
   ENDIF
    
   IF(.NOT.ASSOCIATED(IDF%X))THEN
    IF(.NOT.IDFREAD(IDF,TRIM(STOAVG),0))RETURN; IF(.NOT.IDFALLOCATEX(IDF))RETURN
   ENDIF

   !## initialization
   CALL IDFCOPY(IDF,SC1)
   CALL IDFCOPY(IDF,NETRCH)
   CALL IDFCOPY(IDF,HEAD1)
   CALL IDFCOPY(IDF,HEAD2)
   CALL IDFCOPY(IDF,QMODF)

   !## read storage coefficient
   IF(.NOT.IDFREADSCALE(TRIM(STOAVG),SC1,2,1,0.0D0,0))RETURN ; WRITE(*,'(A)') 'Reading '//TRIM(STOAVG)//' ...' 
 
   !## timestep=t-1 (HEAD1)
   CALL UTL_IDATETOGDATE(LDATES(IFILES),IY,IM,ID,IHR,IMT,ISC)
   WRITE(CDATE,'(I4.4,2I2.2)') IY,IM,ID

   IF(LEN_TRIM(SOURCEDIR).EQ.0)THEN
    FNAME=TRIM(HEDDIR)//'\HEAD_'//TRIM(CDATE)//'_L1.IDF'
   ELSE
    FNAME=TRIM(SOURCEDIR)//'\HEAD\HEAD_'//TRIM(CDATE)//'_L1.IDF'
   ENDIF
   IF(.NOT.IDFREADSCALE(TRIM(FNAME),HEAD1,2,1,0.0D0,0))RETURN ; WRITE(*,'(A)') 'Reading '//TRIM(FNAME)//' ...' 
  
   !## goto next timestep
   CYCLE

  ENDIF
  
  !## timestep=t  (HEAD2)
  CALL UTL_IDATETOGDATE(LDATES(IFILES),IY,IM,ID,IHR,IMT,ISC)
  WRITE(CDATE,'(I4.4,2I2.2)') IY,IM,ID
  
  IF(LEN_TRIM(SOURCEDIR).EQ.0)THEN
   FNAME=TRIM(MSPDIR)//'\BDGQMODF\BDGQMODF_'//TRIM(CDATE)//'_L1.IDF'
  ELSE
   FNAME=TRIM(SOURCEDIR)//'\METASWAP\BDGQMODF\BDGQMODF_'//TRIM(CDATE)//'_L1.IDF'
  ENDIF
  IF(.NOT.IDFREADSCALE(TRIM(FNAME),QMODF,2,1,0.0D0,0))RETURN ; WRITE(*,'(A)') 'Reading '//TRIM(FNAME)//'...' 

  IF(LEN_TRIM(SOURCEDIR).EQ.0)THEN
   FNAME=TRIM(HEDDIR)//'\HEAD_'//TRIM(CDATE)//'_L1.IDF'
  ELSE
   FNAME=TRIM(SOURCEDIR)//'\HEAD\HEAD_'//TRIM(CDATE)//'_L1.IDF'
  ENDIF
  IF(.NOT.IDFREADSCALE(TRIM(FNAME),HEAD2,2,1,0.0D0,0))RETURN ; WRITE(*,'(A)') 'Reading '//TRIM(FNAME)//'...' 
  
  !## calculating NETRCH
  DO IROW=1,HEAD1%NROW
   DO ICOL=1,HEAD1%NCOL
    IF(HEAD1%X(ICOL,IROW).NE.HEAD1%NODATA.AND.HEAD2%X(ICOL,IROW).NE.HEAD2%NODATA.AND. &
       SC1%X(  ICOL,IROW).NE.SC1%NODATA  .AND.QMODF%X(ICOL,IROW).NE.QMODF%NODATA)THEN
     !## formula is: gwa = delta lvgwmodf * sc1 - qmodf
     NETRCH%X(ICOL,IROW)= ( HEAD2%X(ICOL,IROW) - HEAD1%X(ICOL,IROW) ) * SC1%X(ICOL,IROW) - QMODF%X(ICOL,IROW)
     IF(CONVMM.EQ.1)NETRCH%X(ICOL,IROW)=1000.0D0*NETRCH%X(ICOL,IROW)
    ELSE
     NETRCH%X(ICOL,IROW)=NETRCH%NODATA
    ENDIF  
   ENDDO
  ENDDO
  
  FNAME=TRIM(RESDIR)//'\NETRCH_'//TRIM(CDATE)//'_L1.IDF' ; WRITE(*,'(A/)') 'Writing '//TRIM(FNAME)//'...' 
  IF(.NOT.IDFWRITE(NETRCH,FNAME,1))THEN ; ENDIF

  !## copy head2 -> head1
  HEAD1%X=HEAD2%X
  
 ENDIF
ENDDO    

 WRITE(*,'(A)') 'Number of NETRCH files writen: '//TRIM(VTOS(NCREATE))

 MSPNETRCHCOMPUTE=.TRUE.

END FUNCTION MSPNETRCHCOMPUTE
 
END MODULE MOD_MSPNETRCH