!! Simulate glaciers accumulation and ablation
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.2 - 29th October 2024  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 27/May/2012 | Original code implemented in the Snow module|
! | 1.1      | 19/Feb/2023 | code moved and rewritten from Snow module |
! | 1.2      | 29/Oct/2024 | fix initial condition setting |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to simulate glaciers accumulation and ablation
!
MODULE Glacier

! Modules used:

USE DataTypeSizes, ONLY : &
  ! Imported Type Definitions:
  short, float

USE GridLib, ONLY : &
  !Imported definitions:
  grid_real, grid_integer, &
  !Imported routines:
  NewGrid, GridDestroy, &
  !Imported parameters:
  NET_CDF

USE IniLib, ONLY : &
  !Imported definitions:
  IniList, &
  !Imported routines:
  IniOpen, IniClose, &
  SectionIsPresent, KeyIsPresent, &
  IniReadReal, IniReadInt

USE LogLib, ONLY: &
  !Imported routines:
  Catch

USE GridOperations, ONLY : &
  !Imported routines:
  GridByIni, CellArea,  &
  !Imported operators:
  ASSIGNMENT (=)

USE Chronos, ONLY : &
  !Imported definitions:
  DateTime, &
 !Imported operations:
  OPERATOR (+), &
  OPERATOR (==), &
!Imported routines:
 DayOfYear

USE AirTemperature, ONLY: &
    !imported variables:
    temperatureAir

USE ObservationalNetworks, ONLY : &
  !Imported routines:
  ReadMetadata, CopyNetwork, &
  WriteMetadata, DestroyNetwork, &
  AssignDataFromGrid, WriteData, &
  !Imported defined variable:
  ObservationalNetwork

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

USE Snow, ONLY : &
  !imported routines:
  DegreeDay, &
  !imported variables:
  swe, waterInSnow, dtSnow

USE Morphology, ONLY: &
  !Imported routines:
  DownstreamCell, CheckOutlet

USE MorphologicalProperties, ONLY: &
  !imported variables:
  flowDirection, dem

USE StringManipulation, ONLY : &
    !Imported routines:
    ToString


!Global variables:
INTEGER (KIND = short) :: dtIce  !!computation time step
TYPE (grid_real) :: icewe !! ice water equivalent (m)
TYPE (grid_real) :: meltCoefficientIce !! melt coefficient (mm/day/°C) 
TYPE (grid_real) :: waterInIce  !! free water in snow pack (m)
TYPE (grid_real) :: iceMelt   !! snow melt amount within the time step (m)

!Global routines:
PUBLIC :: GlacierInit
PUBLIC :: GlacierUpdate
PUBLIC :: GlacierPointInit

!local variables:
TYPE (grid_integer)        ,PRIVATE :: meltModel  !!melt model map
TYPE (grid_real)           ,PRIVATE :: meltTemperature !! threshold temperature for ablation starts (°C) 
TYPE (grid_real)           ,PRIVATE :: iceConductivity !! ice hydraulic conductivity (m/s)
TYPE (grid_real)           ,PRIVATE :: QinIce, QoutIce
TYPE (ObservationalNetwork),PRIVATE :: sites
TYPE (ObservationalNetwork),PRIVATE :: sitesICEWE
INTEGER (KIND = short)     ,PRIVATE :: fileUnitPointICEWE
TYPE (DateTime)            ,PRIVATE :: timePointExport
INTEGER (KIND = short)     ,PRIVATE :: doySnowTransform !!day of year when old snow 
                                               !!is transformed to ice

!local routines:
PRIVATE :: GlacierPointExport
PRIVATE :: GlacierParameterUpdate


!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  
!==============================================================================
!| Description: 
!  Initialize glacier model
SUBROUTINE GlacierInit   & 
  !
  (inifile, mask, time)       

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information
TYPE (grid_integer), INTENT(IN) :: mask
TYPE (DateTime),     INTENT(IN) :: time

!local declarations:
TYPE (IniList)         :: iniDB
REAL (KIND = float)    :: scalar
INTEGER (KIND = short) :: iscalar

!---------------------end of declarations--------------------------------------

!open and read configuration file
CALL IniOpen (inifile, iniDB)

!day of year for snow to ice tranformation
IF ( KeyIsPresent ('doy-snow-ice-transformation', iniDB ) ) THEN
    doySnowTransform = IniReadInt ( 'doy-snow-ice-transformation', iniDB  ) 
ELSE
     doySnowTransform = 0
END IF


!load melt model
IF (SectionIsPresent('melt-model', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-model') ) THEN
        iscalar = IniReadInt ('scalar', iniDB, 'melt-model')
        CALL NewGrid (meltModel, mask, iscalar)
    ELSE
         CALL GridByIni (iniDB, meltModel, section = 'melt-model')
    END IF
   
ELSE
    CALL Catch ('error', 'Glacier',   &
			      'melt-model not found in configuration file' )
END IF

!set melt coefficient
IF (SectionIsPresent('melt-coefficient', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-coefficient') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'melt-coefficient')
        CALL NewGrid ( meltCoefficientIce, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, meltCoefficientIce, &
                         section = 'melt-coefficient', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Glacier',   &
			      'melt-coefficient not found in configuration file' )
END IF

!set threshold temperature for ablation starts (°C)
IF (SectionIsPresent('melt-threshold-temperature', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-threshold-temperature') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'melt-threshold-temperature')
        CALL NewGrid ( meltTemperature, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, meltTemperature, &
                         section = 'melt-threshold-temperature', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Glacier',   &
        'melt-threshold-temperature not found in configuration file' )
END IF



!set ice hydraulic conductivity (m/s)
IF (SectionIsPresent('hydraulic-conductivity', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'hydraulic-conductivity') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'hydraulic-conductivity')
        CALL NewGrid ( iceConductivity, mask, scalar )
    ELSE
         CALL GridByIni ( iniDB, iceConductivity, &
                         section = 'hydraulic-conductivity', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Glacier',   &
			      'hydraulic-conductivity not found in configuration file' )
END IF


!set initial optional variables

! ice water equivalent
IF (SectionIsPresent('ice-water-equivalent', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'ice-water-equivalent') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'ice-water-equivalent')
        CALL NewGrid ( icewe, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, icewe, section = 'ice-water-equivalent')
    END IF
   
ELSE !set to default = 0
   CALL NewGrid ( icewe, mask, 0. )
END IF


!water in ice
IF (SectionIsPresent('water-in-ice', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'water-in-ice') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'water-in-ice')
        CALL NewGrid ( waterInIce, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, waterInIce, section = 'water-in-ice')
    END IF
   
ELSE !set to default = 0
   CALL NewGrid ( waterInIce, mask, 0. )
END IF


!Configuration terminated. Deallocate ini database
CALL IniClose (iniDB) 

!allocate variables
CALL NewGrid ( waterInIce, mask, 0. )
CALL NewGrid ( iceMelt,    mask, 0. )
CALL NewGrid ( QinIce,     mask, 0. )
CALL NewGrid ( QoutIce,    mask, 0. )


RETURN
END SUBROUTINE GlacierInit
  
  
!==============================================================================
!| Description:
!   Initialize export of point site data
SUBROUTINE GlacierPointInit &
!
( pointfile, path_out, time )

IMPLICIT NONE

!Arguments with intent (in):
CHARACTER (LEN = *), INTENT(IN) :: pointfile  !!file containing coordinate of points
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder
TYPE (DateTime),     INTENT(IN) :: time  !!start time

!local declarations
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: fileunit

!-------------------------end of declarations----------------------------------

timePointExport = time

!open point file

fileunit = GetUnit ()
OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), &
      status='OLD', iostat = err_io )

IF ( err_io /= 0 ) THEN
	!file does not exist
    CALL Catch ('error', 'Glacier', 'out point file does not exist')
END IF 

!Read metadata
CALL ReadMetadata (fileunit, sites)

!check dt
IF (.NOT. ( MOD ( sites % timeIncrement, dtIce ) == 0 ) ) THEN
  CALL Catch ('error', 'Glacier', 'dt out sites must be multiple of dtSnow ')
END IF

CLOSE ( fileunit )

!create virtual network and initialize file for output
fileUnitPointICEWE = GetUnit ()
OPEN ( unit = fileUnitPointICEWE, &
    file = TRIM(path_out) // 'point_glacier.fts' )
    
CALL CopyNetwork ( sites, sitesICEWE )

sitesICEWE % description = 'ice water equivalent data exported from FEST'

sitesICEWE % unit = 'mm'

sitesICEWE % offsetZ = 0.

CALL WriteMetadata ( network = sitesICEWE, &
                fileunit = fileUnitPointICEWE )

CALL WriteData (sitesICEWE, fileUnitPointICEWE, .TRUE.)    

! destroy sites
CALL DestroyNetwork ( sites )
RETURN
END SUBROUTINE GlacierPointInit


!==============================================================================
!| Description:
!   Export of point site data
SUBROUTINE GlacierPointExport &
!
( time )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
INTEGER (KIND = short) :: i

!-------------------------end of declarations----------------------------------



!set current time
sitesICEWE % time = time 
    
!populate data
CALL AssignDataFromGrid ( icewe, sitesICEWE, scaleFactor = 1000. )
    
!write data
CALL WriteData ( sitesICEWE, fileUnitPointICEWE )
    
timePointExport = timePointExport + sitesICEWE % timeIncrement


RETURN
END SUBROUTINE GlacierPointExport


!==============================================================================
!| Description: 
!  compute accumulation and ablation and update  water equivalent
SUBROUTINE GlacierUpdate   & 
  !
  ( time, mask, rainfall )       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time
TYPE (grid_integer), INTENT (IN) :: mask

!Arguments with intent(inout):
TYPE (grid_real), INTENT(INOUT) :: rainfall  !rainfall (liquid precipitation) (m/s) 

!local declarations:
INTEGER (KIND = short) :: i, j, is, js
REAL (KIND = float) :: meltRate
REAL (KIND = float) :: tAir
REAL (KIND = float) :: cmelt !!snow melt coefficient (m/s/°C)
REAL (KIND = float) :: tmelt
REAL (KIND = float) :: icewe_tminus1 !!icewe at previous time step
REAL (KIND = float) :: ddx	!!actual cell length (m)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope  !!local slope (m/m)
REAL (KIND = float) :: qin, qout  !!input and output water discharge in snowpack
REAL (KIND = float) :: area !!cell area (m2)
INTEGER (KIND = short) :: currentDOY !!current day of year
!------------------------------end of declarations-----------------------------  

!check if parameter update is required
CALL GlacierParameterUpdate (time)

!snow ice transformation
currentDOY = DayOfYear (time, 'noleap')
IF ( dtSnow > 0 .AND. currentDOY == doySnowTransform ) THEN
    DO i = 1, mask % idim
        DO j = 1, mask % jdim
            IF (mask % mat(i,j) /= mask % nodata ) THEN
                icewe % mat (i,j) = icewe % mat (i,j) + swe % mat (i,j)
                swe % mat (i,j) = 0. 
                waterInIce % mat (i,j) = waterInIce % mat (i,j) + waterInSnow % mat (i,j)
                waterInSnow % mat (i,j) = 0.
            END IF
        END DO
    END DO
    
END IF

!initialize melt grid
iceMelt = 0.


!compute inter-cell lateral water flux
QinIce = 0.
QoutIce = 0.
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
        !set downstream cell  (is,js)
        CALL DownstreamCell ( i, j, flowDirection % mat(i,j), is, js, ddx, &
                            flowDirection ) 
        
        IF ( CheckOutlet (i, j, is, js, flowDirection) ) CYCLE
      
        !thickness of the saturated depth
        ds = waterInIce % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (mask, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
        
        QoutIce % mat (i,j) = cellWidth * ds * iceConductivity % mat (i,j) *  ijSlope
      
        !output becomes input of downstream cell
        QinIce % mat (is,js) = QinIce % mat (is,js) + QoutIce % mat (i,j)
        
    END DO
END DO

!loop across cells
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        
        !skip nodata
        IF ( mask % mat (i,j) == mask % nodata ) THEN
            CYCLE
        END IF
        
        !compute potential melt rate
        tAir  = temperatureAir % mat (i,j)
        cmelt = meltCoefficientIce % mat (i,j) / 1000. / 86400.
        tmelt = meltTemperature % mat (i,j)
        SELECT CASE ( meltModel % mat (i,j) )
        CASE (1) !degree-day temperature based
            meltRate = DegreeDay ( tAir,  tmelt,  cmelt )
        CASE DEFAULT
            CALL Catch ('error', 'Glacier', 'melt model not implemented: ', &
                         argument = ToString (meltModel % mat (i,j) ) )
        END SELECT
        
        ! actual melt rate, limited by available ice 
        ! in the previous time step
        meltRate = MIN ( meltRate, icewe % mat (i,j) / dtIce )
    
        !update icewe
        IF (dtSnow > 0 ) THEN
            IF ( swe % mat (i,j) <= 0.) THEN !no protection by snow
                icewe_tminus1 = icewe % mat (i,j)
                icewe % mat (i,j) = icewe % mat (i,j) - meltRate *  dtIce
                IF ( icewe % mat (i,j) < 0. ) THEN !meltrate too high
                    icewe % mat (i,j) = 0.
                    iceMelt % mat (i,j) = icewe_tminus1
                ELSE
                    iceMelt % mat (i,j) = meltRate *  dtIce
                END IF
            ELSE  ! snow protects ice from ablation
                iceMelt % mat (i,j) = 0.
            END IF
        ELSE  !no protection by snow
            icewe_tminus1 = icewe % mat (i,j)
            icewe % mat (i,j) = icewe % mat (i,j) - meltRate *  dtIce
            IF ( icewe % mat (i,j) < 0. ) THEN !meltrate too high
                icewe % mat (i,j) = 0.
                iceMelt % mat (i,j) = icewe_tminus1
            ELSE
                iceMelt % mat (i,j) = meltRate *  dtIce
            END IF
        END IF
        
        !update water in snow
        area = CellArea (mask, i, j) 
        qin = QinIce % mat (i,j) / area
        qout = QoutIce % mat (i,j) / area
        
        !lateral water flux occurs even when ice is covered by snow
        IF ( icewe % mat (i,j) > 0. ) THEN !rainfall contributes to water in ice
            waterInIce % mat (i,j) = waterInIce % mat (i,j) + &
                                        iceMelt % mat (i,j) +  &
                                         rainfall % mat (i,j) * dtIce + &
                                        ( qin - qout ) * dtIce
             rainfall % mat (i,j) = 0.
        ELSE 
            waterInIce % mat (i,j) = waterInIce % mat (i,j) + &
                                        iceMelt % mat (i,j) +  &
                                        ( qin - qout ) * dtIce
        END IF

        
        IF ( waterInIce % mat (i,j) < 0. ) THEN
            waterInIce % mat (i,j) = 0.
        END IF
        
    END DO
END DO

!write point icewe
IF ( time == timePointExport ) THEN
   CALL GlacierPointExport ( time )
END IF

RETURN
END SUBROUTINE GlacierUpdate

!==============================================================================
!| Description: 
!  update parameter map that change in time
SUBROUTINE GlacierParameterUpdate   & 
  !
  (time)       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname

!------------------------------end of declarations-----------------------------
  
  
  
!melt coefficient
IF (  time == meltCoefficientIce % next_time ) THEN
   !destroy current grid
   filename = meltCoefficientIce % file_name
   varname = meltCoefficientIce % var_name
   CALL GridDestroy (meltCoefficientIce)
   !read new grid in netcdf file
   CALL NewGrid (meltCoefficientIce, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!melt temperature
IF (  time == meltTemperature % next_time ) THEN
   !destroy current grid
   filename = meltTemperature % file_name
   varname = meltTemperature % var_name
   CALL GridDestroy (meltTemperature)
   !read new grid in netcdf file
   CALL NewGrid (meltTemperature, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!hydraulic conductivity
IF (  time == iceConductivity % next_time ) THEN
   !destroy current grid
   filename = iceConductivity % file_name
   varname = iceConductivity % var_name
   CALL GridDestroy (iceConductivity)
   !read new grid in netcdf file
   CALL NewGrid (iceConductivity, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF
  
RETURN
END SUBROUTINE GlacierParameterUpdate
    
    
    
END MODULE Glacier