!! Elaborate and export raster map 
!|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.0 - 27th March 2024    
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 27/Mar/2024 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! relevant internal variables are aggregated in time
! and space by changing raster resolution and/or 
! spatial reference system, and written to output.
! Raster maps can be used for subsequent elaboration 
! like for example for computing indexes like  
! Standardized Precipitation Index (SPI), or for 
! using them as input to difefrent models. 
! The user provides a list of variables to be processed 
! and exported in the configuration file
! like in the following example
!
!```
! # configure variables for raster export
!
! time = 0 23 * * *
!
! folder = ./results/raster_maps/
!
! [map-template]
!   file = ./data/map_template.asc
!   format = esri-ascii
!   epsg = 32633
!
! [soil-balance]
!   soil-moisture = 1
!   runoff = 1
!   infiltration = 1
!   percolation = 1
!   actual-ET = 1
!   potential-ET = 1
!
! [meteo]
!   precipitation = 1
!   temperature = 1
!   relative-humidity = 0
!   solar-radiation = 0
!   net-radiation = 0
!   wind-speed = 0
!
! [snow]
!   snow-water-equivalent = 1
!
!```
!   
! All variables marked by 1 are elaborated and exported. 
! When one variable is marked by 1 but it is not allocated 
! because not computed by the FEST model according to options
! defined in the configuration files, value is not exported. 
! For example, if user set to export wind-speed 
! but windspeed is not used in the current simulation, raster 
! maps of windspeed are not written to output folder.
! Currently maps are exported in `esri-ascii` format,
! one map for each variable and for each time step.
!
! Variables that can be exported, description, and unit are listed 
! in the following table.
!
! | variable               | Description                                                        | Unit           |
! |------------------------|--------------------------------------------------------------------|----------------|
! | precipitation          | Precipitation fallen in the time step from the last exportation    | mm             |
! | temperature            | Air temperature of the current time step fallen in 24 hours        | Celsius degree |
! | relative-humidity      | Air relative humidity of the current time step                     | % (0-100)      |
! | solar-radiation        | Solar radiation of the current time step                           | w/m²           |
! | net-radiation          | Net radiation of the current time step                             | w/m²           |
! | wind-speed             | Wind speed of the current time step                                | m/s            |
! | snow-water-equivalent  | Snow water equivalent of the current time step                     | mm             |
! | soil-moisture          | Soil moisture of the current time step                             | \-             |
! | runoff                 | Runoff of the current time step                                    | mm             |
! | infiltration           | Infiltration into soil of the current time step                    | mm             |
! | percolation            | Deep percolation out of transmission zone of the current time step | mm             |
! | actual-ET              | Actual evapotranspiration of the current time step                 | mm             |
! | potential-ET           | Potential evapotranspiration of the current time step              | mm             |
!
!
! The name of output files is the concatenation of result 
! folder name  <folder>, a suffix that reminds date and time 
! in the form `YYYY-MM-DDThh-mm`  and the name of variable, 
! as listed in the following table.
!
! | variable                    | Output file name                                         |
! |-----------------------------|----------------------------------------------------------|
! | precipitation               | `<folder>` `YYYY-MM-DDThh-mm` `_precipitation.asc`       |
! | temperature                 | `<folder>` `YYYY-MM-DDThh-mm` `_temperature.asc`         |
! | relative-humidity           | `<folder>` `YYYY-MM-DDThh-mm` `_rh.asc`                  |
! | solar-radiation             | `<folder>` `YYYY-MM-DDThh-mm` `_rad.asc`                 |
! | net-radiation               | `<folder>` `YYYY-MM-DDThh-mm` `_netrad.asc`              |
! | wind-speed                  | `<folder>` `YYYY-MM-DDThh-mm` `_windspeed.asc`           |
! | snow-water-equivalent       | `<folder>` `YYYY-MM-DDThh-mm` `_swe.asc`                 |
! | soil-moisture               | `<folder>` `YYYY-MM-DDThh-mm` `_soil-moisture.asc`       |
! | runoff                      | `<folder>` `YYYY-MM-DDThh-mm` `_runoff.asc`              |
! | infiltration                | `<folder>` `YYYY-MM-DDThh-mm` `_infiltration.asc`        |
! | percolation                 | `<folder>` `YYYY-MM-DDThh-mm` `_percolation.asc`         |
! | actual-ET                   | `<folder>` `YYYY-MM-DDThh-mm` `_et.asc`                  |
! | potential-ET                | `<folder>` `YYYY-MM-DDThh-mm` `_pet.asc`                 |
!
MODULE RasterExport

! Modules used: 

USE DataTypeSizes, ONLY : &  
    ! Imported Parameters:
    float, &
    short 

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE IniLib, ONLY: &
    !Imported derived types:
    IniList, &
    !Imported routines:
    IniOpen, &
    IniClose, &
    IniReadInt,&
    IniReadString,&
    KeyIsPresent, &
    SectionIsPresent

USE GridLib, ONLY : &
    !imported definitions:
    grid_real, &
    !Imported routines:
    NewGrid , &
    ExportGrid, &
    !Imported parameters:
    ESRI_ASCII

USE GridOperations, ONLY : &
    !Imported routines:
    GridByIni, &
    GridConvert, &
    GridResample, &
    !Imported operands:
    ASSIGNMENT( = )

USE Loglib, ONLY : &
    !Imported routines:
    Catch

USE CronSchedule, ONLY : &
    !Imported types:
    CronTab, &
    !Imported routines:
    CronParseString, &
    CronIsTime
        
USE Chronos, ONLY : &
    !Imported types:
    DateTime, &
    !Imported variables:
    timeString, &
    !Imported operands:
    ASSIGNMENT( = )
    

IMPLICIT NONE

!Public routines
PUBLIC :: InitRasterExport
PUBLIC :: ExportMaps

!private declarations

CHARACTER (LEN = 1000)   :: pathout 
LOGICAL                  :: useTemplate
INTEGER (KIND = short), PRIVATE :: countVar  !!count number of variables active for output 
INTEGER (KIND = short), PRIVATE :: countSteps !!number of steps cumulated before exporting
TYPE (CronTab), PRIVATE :: cron

!active output
LOGICAL, PRIVATE :: varOut (13)  !1 = precipitation, 
                                 !2 = air-temperature, 
                                 !3 = relative-humidity
                                 !4 = solar-radiation,
                                 !5 = net-radiation
                                 !6 = wind-speed
                                 !7 = snow-water-equivalent
                                 !8 = soil-moisture
                                 !9 = runoff
                                 !10 = infiltration
                                 !11 = percolation
                                 !12 = actual-ET
                                 !13 = potential-ET

TYPE (grid_real), PRIVATE :: rasterTemplate
TYPE (grid_real), PRIVATE :: rasterPrecipitation
TYPE (grid_real), PRIVATE :: rasterTemperature
TYPE (grid_real), PRIVATE :: rasterRH
TYPE (grid_real), PRIVATE :: rasterRad
TYPE (grid_real), PRIVATE :: rasterNetRad
TYPE (grid_real), PRIVATE :: rasterWS
TYPE (grid_real), PRIVATE :: rasterSWE
TYPE (grid_real), PRIVATE :: rasterSM
TYPE (grid_real), PRIVATE :: rasterRunoff
TYPE (grid_real), PRIVATE :: rasterInfiltration
TYPE (grid_real), PRIVATE :: rasterPercolation
TYPE (grid_real), PRIVATE :: rasterET
TYPE (grid_real), PRIVATE :: rasterPET
TYPE (grid_real), PRIVATE :: gridTemp  !!temporary grid with the same coordinate reference system of template
TYPE (grid_real), PRIVATE :: gridTemp2 !! temporary grid with the same coordinate reference system 
                          !!and spatial extent and resolution of template


!=======
CONTAINS
!=======
! Define procedures contained in this module. 
    
    
!==============================================================================
!| Description: 
!  Initialization of raster export
SUBROUTINE InitRasterExport   & 
!
 (fileini, temp, precipitation, &
  rh, radiation, netradiation, windspeed, &
  swe, sm, runoff, infiltration, percolation, et, pet)  

IMPLICIT NONE

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN)    :: fileini   
TYPE (grid_real), INTENT(IN) :: temp !!air temperarure (°C)
TYPE (grid_real), INTENT(IN) :: precipitation !!precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: rh !!air relative humidity (0-100)
TYPE (grid_real), INTENT(IN) :: radiation !!solar radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: netradiation !!net radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: windspeed !!wind speed (m/s)
TYPE (grid_real), INTENT(IN) :: swe !!snow water equivalent (m)
TYPE (grid_real), INTENT(IN) :: sm !!soil mositure (-)
TYPE (grid_real), INTENT(IN) :: runoff !!runoff (m/s)
TYPE (grid_real), INTENT(IN) :: infiltration !!infiltration (m/s)
TYPE (grid_real), INTENT(IN) :: percolation !!deep percolation (m/s)
TYPE (grid_real), INTENT(IN) :: et !!actual evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration (m/s)

 

!local declarations
TYPE (IniList)          :: iniDB
CHARACTER (LEN = 300)  :: string
!-------------------------------end of declaration-----------------------------

!initialize counter
countSteps = 0

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

! configure time to export data
IF (KeyIsPresent ('time', iniDB) ) THEN
    string =  IniReadString ('time', iniDB)
    CALL CronParseString (string, cron) 
ELSE
    CALL Catch ('error', 'RasterExport', &
            'missing time ' )
END IF

! set template for exported raster
IF (SectionIsPresent ('map-template', iniDB) ) THEN
    useTemplate = .TRUE.
    CALL GridByIni (iniDB, rasterTemplate, section = 'map-template')
    gridTemp % grid_mapping = rasterTemplate % grid_mapping 
    CALL NewGrid ( gridTemp2, rasterTemplate )
ELSE
    useTemplate = .FALSE.
    CALL NewGrid (rasterTemplate, mask) 
END IF

! set out folder
IF (KeyIsPresent ('folder', iniDB) ) THEN
    pathout =  IniReadString ('folder', iniDB)
ELSE
    CALL Catch ('error', 'RasterExport', &
            'missing folder for output ' )
END IF

! search for active variable for output
CALL Catch ('info', 'RasterExport', 'checking for active variables ')

countVar = 0

!precipitation
IF ( IniReadInt ('precipitation', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (temp % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'air-temperature not allocated, &
                                        forced to not export raster ')
       varOut (1) = .FALSE.
   ELSE
       varOut (1) = .TRUE.
       CALL NewGrid (rasterPrecipitation, rasterTemplate)
       
   END IF
ELSE
   varOut (1) = .FALSE.
END IF

!air-temperature
IF ( IniReadInt ('temperature', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (precipitation % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'precipitation not allocated, &
                                        forced to not export raster ')
       varOut (2) = .FALSE.
   ELSE
       varOut (2) = .TRUE.
       CALL NewGrid (rasterTemperature, rasterTemplate)
       
   END IF
ELSE
   varOut (2) = .FALSE.
END IF

!relative-humidity
IF ( IniReadInt ('relative-humidity', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (rh % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'relative humidity not allocated, &
                                        forced to not export raster ')
       varOut (3) = .FALSE.
   ELSE
       varOut (3) = .TRUE.
       CALL NewGrid (rasterRH, rasterTemplate)
       
   END IF
ELSE
   varOut (3) = .FALSE.
END IF


!solar-radiation
IF ( IniReadInt ('solar-radiation', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (radiation % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'solar radiation not allocated, &
                                        forced to not export raster ')
       varOut (4) = .FALSE.
   ELSE
       varOut (4) = .TRUE.
       CALL NewGrid (rasterRad, rasterTemplate)
       
   END IF
ELSE
   varOut (4) = .FALSE.
END IF


!net-radiation
IF ( IniReadInt ('net-radiation', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (netradiation % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'net radiation not allocated, &
                                        forced to not export raster ')
       varOut (5) = .FALSE.
   ELSE
       varOut (5) = .TRUE.
       CALL NewGrid (rasterNetRad, rasterTemplate)
       
   END IF
ELSE
   varOut (5) = .FALSE.
END IF


!wind-speed
IF ( IniReadInt ('wind-speed', iniDB, section = 'meteo') == 1) THEN
   IF ( .NOT. ALLOCATED (windspeed % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'wind speed not allocated, &
                                        forced to not export raster ')
       varOut (6) = .FALSE.
   ELSE
       varOut (6) = .TRUE.
       CALL NewGrid (rasterWS, rasterTemplate)
       
   END IF
ELSE
   varOut (6) = .FALSE.
END IF

!snow-water-equivalent
IF ( IniReadInt ('snow-water-equivalent', iniDB, section = 'snow') == 1) THEN
   IF ( .NOT. ALLOCATED (swe % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'snow water equivalent not allocated, &
                                        forced to not export raster ')
       varOut (7) = .FALSE.
   ELSE
       varOut (7) = .TRUE.
       CALL NewGrid (rasterSWE, rasterTemplate)
       
   END IF
ELSE
   varOut (7) = .FALSE.
END IF

!soil-moisture
IF ( IniReadInt ('soil-moisture', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (sm % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'soil moisture not allocated, &
                                        forced to not export raster ')
       varOut (8) = .FALSE.
   ELSE
       varOut (8) = .TRUE.
       CALL NewGrid (rasterSM, rasterTemplate)
       
   END IF
ELSE
   varOut (8) = .FALSE.
END IF

!runoff
IF ( IniReadInt ('runoff', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (runoff % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'runoff not allocated, &
                                        forced to not export raster ')
       varOut (9) = .FALSE.
   ELSE
       varOut (9) = .TRUE.
       CALL NewGrid (rasterRunoff, rasterTemplate)
       
   END IF
ELSE
   varOut (9) = .FALSE.
END IF

!infiltration
IF ( IniReadInt ('infiltration', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (infiltration % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'infiltration not allocated, &
                                        forced to not export raster ')
       varOut (10) = .FALSE.
   ELSE
       varOut (10) = .TRUE.
       CALL NewGrid (rasterInfiltration, rasterTemplate)
       
   END IF
ELSE
   varOut (10) = .FALSE.
END IF

!percolation
IF ( IniReadInt ('percolation', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (percolation % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'percolation not allocated, &
                                        forced to not export raster ')
       varOut (11) = .FALSE.
   ELSE
       varOut (11) = .TRUE.
       CALL NewGrid (rasterPercolation, rasterTemplate)
       
   END IF
ELSE
   varOut (11) = .FALSE.
END IF

!actual-ET
IF ( IniReadInt ('actual-ET', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (et % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'ET not allocated, &
                                        forced to not export raster ')
       varOut (12) = .FALSE.
   ELSE
       varOut (12) = .TRUE.
       CALL NewGrid (rasterET, rasterTemplate)
       
   END IF
ELSE
   varOut (12) = .FALSE.
END IF

!potential-ET
IF ( IniReadInt ('potential-ET', iniDB, section = 'soil-balance') == 1) THEN
   IF ( .NOT. ALLOCATED (pet % mat) ) THEN
       CALL Catch ('warning', 'RasterExport', 'PET not allocated, &
                                        forced to not export raster ')
       varOut (13) = .FALSE.
   ELSE
       varOut (13) = .TRUE.
       CALL NewGrid (rasterPET, rasterTemplate)
       
   END IF
ELSE
   varOut (13) = .FALSE.
END IF


CALL IniClose (iniDB) 


!Initialize times
!timeNewTemp = time

RETURN
END SUBROUTINE InitRasterExport
  
  
  
!==============================================================================
!| Description: 
!  Update and export maps
SUBROUTINE ExportMaps   & 
!
 (time, dt, temp, precipitation, rh, radiation, netradiation, windspeed, &
  swe, sm, runoff, infiltration, percolation, et, pet) 

IMPLICIT NONE

!arguments with intent in:
TYPE (DateTime),  INTENT(IN) :: time
INTEGER (KIND = short), INTENT(IN) :: dt !! time step (s)
TYPE (grid_real), INTENT(IN) :: temp !!air temperarure (°C)
TYPE (grid_real), INTENT(IN) :: precipitation !!precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: rh !!air relative humidity (0-100)
TYPE (grid_real), INTENT(IN) :: radiation !!solar radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: netradiation !!net radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: windspeed !!wind speed (m/s)
TYPE (grid_real), INTENT(IN) :: swe !!snow water equivalent (m)
TYPE (grid_real), INTENT(IN) :: sm !!soil moisture (-)
TYPE (grid_real), INTENT(IN) :: runoff !!runoff (m/s)
TYPE (grid_real), INTENT(IN) :: infiltration !!infiltration (m/s)
TYPE (grid_real), INTENT(IN) :: percolation !!deep percolation (m/s)
TYPE (grid_real), INTENT(IN) :: et !!actual evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration (m/s)

!local declarations:
INTEGER (KIND = short) :: i, j
CHARACTER (LEN = 300) :: string
CHARACTER (LEN = 16)  :: string16
!--------------------------end of declarations---------------------------------

!update precipitation
IF ( varOut (1) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (precipitation, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterPrecipitation % mat (i,j) = rasterPrecipitation % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterPrecipitation % mat (i,j) =  rasterPrecipitation % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterPrecipitation % mat (i,j) =   rasterPrecipitation % mat (i,j) + &
                          precipitation % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF


!update temperature
IF ( varOut (2) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (temp, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterTemperature % mat (i,j) =   &
                          ( countSteps * rasterTemperature % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterTemperature % mat (i,j) =  rasterTemperature % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterTemperature % mat (i,j) =   &
                          ( countSteps * rasterTemperature % mat (i,j) + &
                          temp % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF


!update relative humidity
IF ( varOut (3) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (rh, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterRH % mat (i,j) =   &
                          ( countSteps * rasterRH % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterRH % mat (i,j) =  rasterRH % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterRH % mat (i,j) =   &
                          ( countSteps * rasterRH % mat (i,j) + &
                          rh % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF


!update radiation
IF ( varOut (4) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (radiation, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterRad % mat (i,j) =   &
                          ( countSteps * rasterRad % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterRad % mat (i,j) =  rasterRad % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterRad % mat (i,j) =   &
                          ( countSteps * rasterRad % mat (i,j) + &
                          radiation % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF


!update net radiation
IF ( varOut (5) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (netradiation, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterNetRad % mat (i,j) =   &
                          ( countSteps * rasterNetRad % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterNetRad % mat (i,j) =  rasterNetRad % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterNetRad % mat (i,j) =   &
                          ( countSteps * rasterNetRad % mat (i,j) + &
                          netradiation % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF


!update wind speed
IF ( varOut (6) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (windspeed, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterWS % mat (i,j) =   &
                          ( countSteps * rasterWS % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterWS % mat (i,j) =  rasterWS % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterWS % mat (i,j) =   &
                          ( countSteps * rasterWS % mat (i,j) + &
                          windspeed % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF

!update snow water equivalent
IF ( varOut (7) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (swe, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterSWE % mat (i,j) =   &
                          ( countSteps * rasterSWE % mat (i,j) + &
                          gridTemp2 % mat (i,j) * 1000. ) / ( countSteps + 1)
                ELSE
                    rasterSWE % mat (i,j) =  rasterSWE % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterSWE % mat (i,j) =   &
                          ( countSteps * rasterSWE % mat (i,j) + &
                          swe % mat (i,j) * 1000. ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF

!update soil moisture
IF ( varOut (8) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (sm, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterSM % mat (i,j) =   &
                          ( countSteps * rasterSM % mat (i,j) + &
                          gridTemp2 % mat (i,j) ) / ( countSteps + 1)
                ELSE
                    rasterSM % mat (i,j) =  rasterSM % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterSM % mat (i,j) =   &
                          ( countSteps * rasterSM % mat (i,j) + &
                          sm % mat (i,j) ) / ( countSteps + 1)
                END IF
            END DO
        END DO
    END IF
END IF

!update runoff
IF ( varOut (9) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (runoff, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterRunoff % mat (i,j) = rasterRunoff % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterRunoff % mat (i,j) =  rasterRunoff % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterRunoff % mat (i,j) =   rasterRunoff % mat (i,j) + &
                          runoff % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF

!update infiltration
IF ( varOut (10) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (infiltration, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterInfiltration % mat (i,j) = rasterInfiltration % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterInfiltration % mat (i,j) =  rasterInfiltration % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterInfiltration % mat (i,j) =   rasterInfiltration % mat (i,j) + &
                          infiltration % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF

!update percolation
IF ( varOut (11) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (percolation, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterPercolation % mat (i,j) = rasterPercolation % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterPercolation % mat (i,j) =  rasterPercolation % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterPercolation % mat (i,j) =   rasterPercolation % mat (i,j) + &
                          percolation % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF

!update actual evapotranspiration
IF ( varOut (12) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (et, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterET % mat (i,j) = rasterET % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterET % mat (i,j) =  rasterET % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterET % mat (i,j) =   rasterET % mat (i,j) + &
                          et % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF

!update potential evapotranspiration
IF ( varOut (13) ) THEN
    IF ( useTemplate ) THEN !need to convert maps
        CALL GridConvert (pet, gridTemp)
        CALL GridResample (gridTemp, gridTemp2)
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata .AND. &
                     gridTemp2 % mat (i,j) /= gridTemp2 % nodata ) THEN
                    rasterPET % mat (i,j) = rasterPET % mat (i,j) + &
                          gridTemp2 % mat (i,j) * dt * 1000.
                ELSE
                    rasterPET % mat (i,j) =  rasterPET % nodata
                END IF
            END DO
        END DO
    ELSE
        DO  j = 1, rasterTemplate % jdim
            DO i = 1, rasterTemplate % idim
                IF ( rasterTemplate % mat (i,j) /= rasterTemplate % nodata) THEN
                    rasterPET % mat (i,j) =   rasterPET % mat (i,j) + &
                          pet % mat (i,j)  * dt * 1000.
                END IF
            END DO
        END DO
    END IF
END IF


countSteps = countSteps + 1

IF (CronIsTime (time, cron) ) THEN
    !set path
    timeString = time  !convert to string  'YYYY-MM-DDThh:mm:ssTZD'
    string16 = timeString (1:16)
    string16 (14:14) = '-'
    
    !precipitation
    IF (varOut (1) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_precipitation.asc'
        CALL ExportGrid (rasterPrecipitation, string, ESRI_ASCII)
        !reset raster
        rasterPrecipitation = 0.
    END IF
    
    !air temperature
    IF (varOut (2) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_temperature.asc'
        CALL ExportGrid (rasterTemperature, string, ESRI_ASCII)
        !reset raster
        rasterTemperature = 0.
    END IF
    
    !air relative humidity
    IF (varOut (3) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_rh.asc'
        CALL ExportGrid (rasterRH, string, ESRI_ASCII)
        !reset raster
        rasterRH = 0.
    END IF
    
    !radiation
    IF (varOut (4) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_rad.asc'
        CALL ExportGrid (rasterRad, string, ESRI_ASCII)
        !reset raster
        rasterRad = 0.
    END IF
    
    !net radiation
    IF (varOut (5) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_netrad.asc'
        CALL ExportGrid (rasterNetRad, string, ESRI_ASCII)
        !reset raster
        rasterNetRad = 0.
    END IF
    
    !wind speed
    IF (varOut (6) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_windspeed.asc'
        CALL ExportGrid (rasterWS, string, ESRI_ASCII)
        !reset raster
        rasterWS = 0.
    END IF
    
    !snow water equivalent
    IF (varOut (7) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_swe.asc'
        CALL ExportGrid (rasterSWE, string, ESRI_ASCII)
        !reset raster
        rasterSWE = 0.
    END IF
    
    !soil moisture
    IF (varOut (8) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_soil-moisture.asc'
        CALL ExportGrid (rasterSM, string, ESRI_ASCII)
        !reset raster
        rasterSM = 0.
    END IF
    
    !runoff
    IF (varOut (9) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_runoff.asc'
        CALL ExportGrid (rasterRunoff, string, ESRI_ASCII)
        !reset raster
        rasterRunoff = 0.
    END IF
    
    !infiltration
    IF (varOut (10) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_infiltration.asc'
        CALL ExportGrid (rasterInfiltration, string, ESRI_ASCII)
        !reset raster
        rasterInfiltration = 0.
    END IF
    
    !percolation
    IF (varOut (11) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_percolation.asc'
        CALL ExportGrid (rasterPercolation, string, ESRI_ASCII)
        !reset raster
        rasterPercolation = 0.
    END IF
    
    !actual evapotranspiration
    IF (varOut (12) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_et.asc'
        CALL ExportGrid (rasterET, string, ESRI_ASCII)
        !reset raster
        rasterET = 0.
    END IF
    
    !potential evapotranspiration
    IF (varOut (13) ) THEN 
        string = TRIM (pathout) // TRIM (string16) // '_pet.asc'
        CALL ExportGrid (rasterPET, string, ESRI_ASCII)
        !reset raster
        rasterPET = 0.
    END IF
    
    !reset counter
    countSteps = 0
 
END IF

RETURN
END SUBROUTINE ExportMaps

END MODULE RasterExport