! "Dummy meteo- hydro- model" for Pochva scheme application

PROGRAM MODEL

USE POCHVA_SCHEME, ONLY : POCHVA

IMPLICIT NONE

INTEGER, PARAMETER :: NLON=1, NLAT=1, NLEV=1, NLEV_SOIL=22, NLEV_SNOW=10
INTEGER :: NSTEP
REAL :: DTSTEP=60., HRUN=351.0, VALMISS=-9999.
LOGICAL :: FLAG_RESTART=.false.

INTEGER, DIMENSION(NLON,NLAT) :: MASK_SOIL, IND_LEV_TBOTTOM, IND_LEV_QBOTTOM

REAL, DIMENSION(NLON,NLAT) :: SNOW, FSNOW, SNOW_ALBEDO, TSURF, QSURF, QSSURF, PSURF, FL_RAIN, FL_SNOW, &
 FRTOT, FRIRR, FRVIS, FLH_SPEC, FLH_LATENT, QFLUX, KTURBH, KTURBQ, LAIVEG, FRACVEG, ROOTVEG, &
 FLUX_HEAT_BOTTOM, FLUX_WATER_BOTTOM, RUN_OFF, RUN_OFF_TOT

REAL, DIMENSION(NLON,NLAT,NLEV_SOIL) :: TG, QG, QSOIL_MAX, QSOIL_MIN,&
 CSOIL, RHOSOIL, PSISOIL, WCONDSOIL, PARBSOIL, PARCSOIL, QSOIL_REL_WILT, QSOIL_REL_REF

REAL, DIMENSION(NLON,NLAT,0:NLEV_SOIL) :: FICE_SOIL

REAL, DIMENSION(NLON,NLAT,0:NLEV_SNOW) :: SNOWLEV, TSNOW, FICE_SNOW, SNOW_AGE, SNOW_MELT_AGE, SNOW_DENS

REAL, DIMENSION(NLON,NLAT,NLEV) :: T, Q, P, H_LEV

REAL, DIMENSION(NLEV_SOIL) :: SOILLEV

REAL :: TZER=273.15, LAIVEG_MAX=15., FRIR, SNOW_DLEV=10., LAIVEG_OLD, FRACVEG_OLD

INTEGER :: FLAG_FLUX_DEF=0, FLAG_INITIAL, FLAG_PAR_CHANGE=0, JSTEP, JSTEP_RD, &
 IND_PROC=0, I, J, K, N, NLEV_SOIL_RD

REAL :: ALBEDO1, ALBEDO2, SNOW_AGE1,SNOW_AGE_MELT1,SNOW_FICE1, FRDOWN, FRIIR, FRVIS1

 NSTEP=NINT(HRUN*3600./DTSTEP)

! Prognostic variables

 TG(:,:,:)=0.+TZER
 TSURF(:,:)=0.+TZER
 FICE_SOIL(:,:,:)=0.
! QG(:,:,:)=0.2
 QG(:,:,5)=0.35
 SNOW(:,:)=77.
 FSNOW(:,:)=1.
 TSNOW(:,:,:)=0.+TZER
 T(:,:,1)=10.+TZER

 SNOW_AGE(:,:,:)=1.
 MASK_SOIL(:,:)=1
 H_LEV(:,:,1)=15.
 PSURF(:,:)=101300.
 P(:,:,1)=100000.
 Q(:,:,1)= 1.E-4
! FL_RAIN(:,:)=0.0001 ! kg/m/m/s
! FL_RAIN(:,:)=0.0010 ! kg/m/m/s
 FL_RAIN(:,:)=0.0000 ! kg/m/m/s
 FL_SNOW(:,:)=0.0000 ! kg/m/m/s
 FRIRR(:,:)=-50.
 FRVIS(:,:)=350.
 FRTOT(:,:)=FRIRR(:,:)+FRVIS(:,:) 
! KTURBH(:,:)=0.005
! KTURBQ(:,:)=0.005
 KTURBH(:,:)=0.000
 KTURBQ(:,:)=0.000

! Soil levels depth (m): 0.005, 0.0175, 0.04375, 0.109375, 0.2734375, 0.68359375, 1.708984375, 4.272460937, 10.681152344

! SOILLEV(1)=0.5e-2
! DO K=2,NLEV_SOIL
!  SOILLEV(K)=(2.5**(K-1)+2.5**(K-2))*0.5*1.e-2
! ENDDO

! Soil levels depth (m): 0.05, 0.13, 0.32, 0.69, 1.30, 2.21, 3.48, 5.17, 7.34, 10.05, 13.36, 17.33, 22.02, 27.49, 33.8

 DO K=1,NLEV_SOIL
  SOILLEV(K)=(5.+FLOAT(K)**3)*1.e-2
 ENDDO

! Snow levels (kg/m/m)

 DO J=1,NLAT
 DO I=1,NLON
   SNOWLEV(I,J,:)=VALMISS
   IF (SNOW(I,J) > 0.) THEN
     DO WHILE (.TRUE.)
       IF (SNOW_DLEV*FLOAT(NLEV_SNOW) >= SNOW(I,J)) EXIT
       SNOW_DLEV=SNOW_DLEV+10.
     ENDDO
     SNOWLEV(I,J,0)=0.
     IF (SNOW(I,J) >= SNOW_DLEV+5.) THEN
       N=INT(SNOW(I,J)/SNOW_DLEV)
       IF (SNOW(I,J)-SNOW_DLEV*FLOAT(n) >= 5.) N=N+1
       DO K=1,N
         SNOWLEV(I,J,K)=SNOW(I,J)-FLOAT(N-K)*SNOW_DLEV
       ENDDO
       N=N-1
     ELSE
       N=0
       SNOWLEV(I,J,1)=SNOW(I,J)
     ENDIF
     TSNOW(I,J,N+1:NLEV_SNOW)=VALMISS
     FICE_SNOW(I,J,N+1:NLEV_SNOW)=VALMISS
     SNOW_AGE(I,J,N+1:NLEV_SNOW)=VALMISS
     SNOW_MELT_AGE(I,J,N+1:NLEV_SNOW)=VALMISS
   ELSE
     TSNOW(I,J,:)=VALMISS
     FICE_SNOW(I,J,:)=VALMISS
     SNOW_AGE(I,J,:)=VALMISS
     SNOW_MELT_AGE(I,J,:)=VALMISS
   ENDIF
 ENDDO
 ENDDO

! Physical parameters of soil and vegetation

! Sand over glay
 QSOIL_MAX(:,:,1:3)=0.395
 QSOIL_MIN(:,:,1:3)=0.015
 CSOIL(:,:,1:3)=948.
 RHOSOIL(:,:,1:3)=1550.
 PSISOIL(:,:,1:3)=-0.121
 WCONDSOIL(:,:,1:3)=6.0e-5
 PARBSOIL(:,:,1:3)=4.05
 PARCSOIL(:,:,1:3)=3.
 QSOIL_REL_WILT(:,:,1:3)=(0.0677-QSOIL_MIN(:,:,1:3))/(QSOIL_MAX(:,:,1:3)-QSOIL_MIN(:,:,1:3))
 QSOIL_REL_REF(:,:,1:3)=(0.2700-QSOIL_MIN(:,:,1:3))/(QSOIL_MAX(:,:,1:3)-QSOIL_MIN(:,:,1:3))
!
 QSOIL_MAX(:,:,4:NLEV_SOIL)=0.482
 QSOIL_MIN(:,:,4:NLEV_SOIL)=0.140
 CSOIL(:,:,4:NLEV_SOIL)=681.
 RHOSOIL(:,:,4:NLEV_SOIL)=1600.
 PSISOIL(:,:,4:NLEV_SOIL)=-0.405
 WCONDSOIL(:,:,4:NLEV_SOIL)=4.8e-5
 PARBSOIL(:,:,4:NLEV_SOIL)=11.4
 PARCSOIL(:,:,4:NLEV_SOIL)=3.
 QSOIL_REL_WILT(:,:,4:NLEV_SOIL)=(0.2864-QSOIL_MIN(:,:,4:NLEV_SOIL))/(QSOIL_MAX(:,:,4:NLEV_SOIL)-QSOIL_MIN(:,:,4:NLEV_SOIL))
 QSOIL_REL_REF(:,:,4:NLEV_SOIL)=(0.3700-QSOIL_MIN(:,:,4:NLEV_SOIL))/(QSOIL_MAX(:,:,4:NLEV_SOIL)-QSOIL_MIN(:,:,4:NLEV_SOIL))

!! Glacier
! QSOIL_MAX(:,:,1:NLEV_SOIL)=1.0
! QSOIL_MIN(:,:,1:NLEV_SOIL)=0.0
! CSOIL(:,:,1:NLEV_SOIL)=2093.4
! RHOSOIL(:,:,1:NLEV_SOIL)=900.
! PSISOIL(:,:,1:NLEV_SOIL)=2.
! WCONDSOIL(:,:,1:NLEV_SOIL)=0.
! PARBSOIL(:,:,1:NLEV_SOIL)=0.
! PARCSOIL(:,:,1:NLEV_SOIL)=0.
! QSOIL_REL_WILT(:,:,1:NLEV_SOIL)=1.
! QSOIL_REL_REF(:,:,1:NLEV_SOIL)=1.

! Bare soil
 LAIVEG(:,:)=0.
 FRACVEG(:,:)=0.0
 ROOTVEG(:,:)=0.0

 IND_LEV_TBOTTOM(:,:)=11
 IND_LEV_QBOTTOM(:,:)=8

 IF (.NOT.FLAG_RESTART) THEN
   FLAG_INITIAL=1
   J=1
   I=1
   OPEN (12,FILE="observ_initial.dat2",STATUS="OLD")
   READ (12,1001) NLEV_SOIL_RD, IND_LEV_TBOTTOM(I,J), IND_LEV_QBOTTOM(I,J)
   DO K=1,NLEV_SOIL_RD
     READ (12,1002) SOILLEV(K), TG(I,J,K), QG(I,J,K), &
 QSOIL_MAX(I,J,K), QSOIL_MIN(I,J,K), &
 CSOIL(I,J,K), RHOSOIL(I,J,K), PSISOIL(I,J,K), &
 WCONDSOIL(I,J,K), PARBSOIL(I,J,K), PARCSOIL(I,J,K),& 
 QSOIL_REL_WILT(I,J,K), QSOIL_REL_REF(I,J,K)
   ENDDO
   READ (12,1002) TSURF(I,J), QSSURF(I,J), SNOW(I,J), &
 LAIVEG(I,J), LAIVEG_MAX, FRACVEG(I,J), ROOTVEG(I,J) 
   DO K=0,NLEV_SNOW
     READ (12,1002) SNOWLEV(I,J,K), TSNOW(I,J,K), SNOW_AGE(I,J,K), SNOW_MELT_AGE(I,J,K)
   ENDDO
   READ (12,1002) VALMISS

! poisk
!CSOIL(I,J,1:2)=CSOIL(I,J,3)
!RHOSOIL(I,J,1:2)=RHOSOIL(I,J,3)
!PSISOIL(I,J,1:2)=PSISOIL(I,J,3)
!WCONDSOIL(I,J,1:2)=WCONDSOIL(I,J,3)
!PARBSOIL(I,J,1:2)=PARBSOIL(I,J,3)

   CLOSE (12)

 PRINT *,'        Initial'
 PRINT *,'------------------'
 PRINT *,'NLEV_SOIL_RD=',NLEV_SOIL_RD
 PRINT *
 PRINT *,'TSURF=',TSURF(I,J)-TZER,' QSSURF=',QSSURF(I,J),' SNOW=',SNOW(I,J)
 PRINT *,'              TSOIL             QSOIL'
 DO K=1,NLEV_SOIL_RD
   PRINT *,K,SOILLEV(K),TG(I,J,K)-TZER,QG(I,J,K)
 ENDDO
 PRINT *,'              QSOIL_MAX         QSOIL_MIN       QSOIL_REL_WILT    QSOIL_REL_REF'
 DO K=1,NLEV_SOIL_RD
   PRINT *,K,QSOIL_MAX(I,J,K),QSOIL_MIN(I,J,K),QSOIL_REL_WILT(I,J,K),QSOIL_REL_REF(I,J,K)
 ENDDO
 PRINT *,'              CSOIL             RHOSOIL'
 DO K=1,NLEV_SOIL_RD
   PRINT *,K,CSOIL(I,J,K),RHOSOIL(I,J,K)
 ENDDO
 PRINT *,'              PSISOIL          WCONDSOIL          PARBSOIL'
 DO K=1,NLEV_SOIL_RD
   PRINT *,K,PSISOIL(I,J,K),WCONDSOIL(I,J,K),PARBSOIL(I,J,K)
 ENDDO
 PRINT *,'LAIVEG=',LAIVEG(I,J),' LAIVEG_MAX=',LAIVEG_MAX,' FRACVEG=',FRACVEG(I,J),' ROOTVEG=',ROOTVEG(I,J)
 PRINT *,'IND_LEV_TBOTTOM=',IND_LEV_TBOTTOM(I,J),' IND_LEV_QBOTTOM=',IND_LEV_QBOTTOM(I,J)
 PRINT *,'              SNOWLEV          TSNOW            SNOW_AGE       SNOW_MELT_AGE'
 DO K=1,NLEV_SNOW
   PRINT *,K,SNOWLEV(I,J,K),TSNOW(I,J,K),SNOW_AGE(I,J,K),SNOW_MELT_AGE(I,J,K)
 ENDDO
 PRINT *,'------------------'

! STOP
! READ*

 ELSE

   FLAG_INITIAL=0
   CALL RDRF_POCHVA(NLON, NLAT)

 ENDIF

 OPEN (12,FILE="observ_atm_forcing.dat2",STATUS="OLD")

 DO JSTEP=1,NSTEP

   J=1
   I=1
   LAIVEG_OLD=LAIVEG(I,J)
   FRACVEG_OLD=FRACVEG(I,J)
   READ (12,1003) JSTEP_RD, H_LEV(I,J,1), PSURF(I,J), T(I,J,1), Q(I,J,1),&
 FL_RAIN(I,J), FL_SNOW(I,J), FRTOT(I,J), FRVIS(I,J), &
 FLH_SPEC(I,J), FLH_LATENT(I,J), QFLUX(I,J)

   FLAG_FLUX_DEF = 1 !!!!

   P(I,J,1)=PSURF(I,J)*0.995
   KTURBH(I,J)=1.0 
   KTURBQ(I,J)=1.0 

   FLAG_PAR_CHANGE=0
   IF (ABS(LAIVEG_OLD-LAIVEG(I,J)) > 1.E-3.OR.ABS(FRACVEG_OLD-FRACVEG(I,J)) > 1.E-4) FLAG_PAR_CHANGE=1

   CALL POCHVA &
 (SOILLEV, TG, TSURF, QG, QSURF, &
 SNOW, FSNOW, SNOW_ALBEDO, &
 FLUX_HEAT_BOTTOM, FLUX_WATER_BOTTOM, RUN_OFF, RUN_OFF_TOT, &
 MASK_SOIL, H_LEV(:,:,1), PSURF, P(:,:,1), T(:,:,1), Q(:,:,1), &
 FL_RAIN, FL_SNOW, FRTOT, FRVIS, FLH_SPEC, FLH_LATENT, QFLUX, KTURBH, KTURBQ, &
 QSOIL_MAX, QSOIL_MIN, CSOIL, RHOSOIL, PSISOIL, WCONDSOIL, PARBSOIL, PARCSOIL, &
 QSOIL_REL_WILT, QSOIL_REL_REF, LAIVEG, LAIVEG_MAX, FRACVEG, ROOTVEG, IND_LEV_TBOTTOM, IND_LEV_QBOTTOM, &
 VALMISS, DTSTEP, JSTEP, FLAG_FLUX_DEF, FLAG_INITIAL, FLAG_PAR_CHANGE, IND_PROC, &
 QS_SURF_EXT=QSSURF, LEV_SNOW_EXT=SNOWLEV, TSNOW_EXT=TSNOW, &
 SNOW_AGE_EXT=SNOW_AGE, SNOW_MELT_AGE_EXT=SNOW_MELT_AGE, RHO_SNOW_EXT=SNOW_DENS)

 FLAG_INITIAL=0

 ENDDO
 
1001 FORMAT (100I3) 
1002 FORMAT (100e16.8) 
1003 FORMAT (I8,100e16.8) 

STOP
END PROGRAM MODEL
!-------------------------------------------------------------------------------
SUBROUTINE RDRF_POCHVA(NPOINT_X_EXT, NPOINT_Y_EXT)
! Read restart file with pochva variables

USE MODULE_POCHVA, ONLY : DTSTEP, NPOINT_X, NPOINT_Y, NLEV_SOIL_MAX, NLEV_SNOW, &
 NLEV_SOIL_HEAT, NLEV_SOIL_WATER, LEV_SOIL, LEV_SOIL_H, D_LEV_SOIL, D_LEV_SOIL_H, SUM_D_LEV_SOIL_H, &
 TS, QS, QS_REL, PSIS, KAPPA_H, FRAC_SICE, DFRAC_SICE_DT, QSURF, TSURF, TFOREST, QFOREST, QS_REL_VEG, &
 D_LEV_SNOW, TSNOW, LEV_SNOW, FICE_SNOW, RHOSNOW, SNOW_AGE, SNOW_MELT_AGE, DLEV_SNOW_BOTTOM, FLAG_SNOW_THICK, &
 QS_OLD, QSI, DQS, RUN_OFF, FLUX_ENTROPY_ATM, FLUX_ENTROPY_ATM_SPEC, FLUX_ENTROPY_ATM_LAT, &
 FLUX_ENTROPY_ATM_RAD, FLUX_ENTROPY_ATM_SOIL, FLUX_ENTROPY_ATM_SNOW, FLUX_ENTROPY_ATM_FOREST, &
 FLUX_ENTROPY_SNOW_SOIL, FLUX_ENTROPY_FOREST_SOIL, &
 FLUX_WV, FLUX_WV_SOIL, FLUX_WV_SNOW, FLUX_WV_VEGLOW_DRY, FLUX_WV_VEGLOW_WET, &
 FLUX_WV_VEGLOW_WET, FLUX_WV_FOREST_DRY, FLUX_WV_FOREST_WET, &
 FLUX_PREC_LIQ, FLUX_PREC_ICE, FLUX_W_LIQ_ATM_SOIL, FLUX_W_LIQ_ATM_SNOW, FLUX_W_ICE_ATM_SNOW, &
 FLUX_W_LIQ_ATM_VEGLOW, FLUX_W_LIQ_ATM_FOREST, FLUX_WATER_SNOW_SOIL, &
 FLUX_SOIL_ENTROPY, FLUX_SOIL_WATER, FLUX_SOIL_WATER_VEG, FLUX_SOIL_WATER_BOTTOM_DIAG, &
 LAMBDAG, FLUX_SNOW_ENTROPY, FLUX_SNOW_WATER, &
 QSMAX, QSMIN, CG, RHOG, PSIS_SAT, KAPPA_SAT, PAR_B, PAR_C, PAR_FB, QS_REL_VEGWILT, QS_REL_VEGREF, &
 LAI_MAX, LAI_VEGLOW, LAI_FOREST, FVEGLOW, FFOREST, FVEGLOWWET, FFORESTWET, ROOT_DEPTH_VEGLOW, &
 QW_VEGLOW, QW_FOREST, ROOT_DEPTH_FOREST, QW_VEGLOW_MAX, QW_FOREST_MAX, DZSUM_VEGLOW, DZSUM_FOREST, &
 FSNCOVER, DZ_DZSUM_VEGLOW, DZ_DZSUM_FOREST, MASK_SOIL, MASK_SOIL_AXED,&
 D_LEV_SNOW_BASE, RHOSNOW_FRESH, RHOSNOW_OLD, VAL_MISSING, &
 PAR_FA, T0, T_FICE_CRIT, &
 PREC_ACCUM, FLUX_Q_ACCUM, RUN_OFF_ACCUM, FLUX_WATER_BOTTOM_ACC

IMPLICIT NONE

INTEGER :: NPOINT_X_EXT, NPOINT_Y_EXT, I, J, K, IUNIT=11
CHARACTER(LEN=30) :: FILE_IN="restart_pochva.dat2"

 OPEN (IUNIT, FILE=FILE_IN, STATUS="OLD")

 I=1
 J=1

 READ (IUNIT,*) DTSTEP, NPOINT_X, NPOINT_Y, NLEV_SOIL_MAX

 NPOINT_X=NPOINT_X_EXT
 NPOINT_Y=NPOINT_Y_EXT

 IF (ALLOCATED(LEV_SOIL)) DEALLOCATE(LEV_SOIL)
 ALLOCATE(LEV_SOIL(0:NLEV_SOIL_MAX))
 LEV_SOIL(:)=0.
 IF (ALLOCATED(LEV_SOIL_H)) DEALLOCATE(LEV_SOIL_H)
 ALLOCATE(LEV_SOIL_H(NLEV_SOIL_MAX))
 LEV_SOIL_H(:)=0.
 IF (ALLOCATED(D_LEV_SOIL)) DEALLOCATE(D_LEV_SOIL)
 ALLOCATE(D_LEV_SOIL(NLEV_SOIL_MAX))
 D_LEV_SOIL(:)=0.
 IF (ALLOCATED(D_LEV_SOIL_H)) DEALLOCATE(D_LEV_SOIL_H)
 ALLOCATE(D_LEV_SOIL_H(0:NLEV_SOIL_MAX-1))
 D_LEV_SOIL_H(:)=0.

 READ (IUNIT,*) LEV_SOIL(0:NLEV_SOIL_MAX), LEV_SOIL_H(1:NLEV_SOIL_MAX), D_LEV_SOIL(1:NLEV_SOIL_MAX), &
  D_LEV_SOIL_H(0:NLEV_SOIL_MAX-1), SUM_D_LEV_SOIL_H, LAI_MAX

! Allocations --->

 IF (ALLOCATED(NLEV_SOIL_HEAT)) DEALLOCATE(NLEV_SOIL_HEAT)
 ALLOCATE(NLEV_SOIL_HEAT(NPOINT_X,NPOINT_Y))
 NLEV_SOIL_HEAT(:,:)=0
 IF (ALLOCATED(NLEV_SOIL_WATER)) DEALLOCATE(NLEV_SOIL_WATER)
 ALLOCATE(NLEV_SOIL_WATER(NPOINT_X,NPOINT_Y))
 NLEV_SOIL_WATER(:,:)=0

 IF (ALLOCATED(TS)) DEALLOCATE(TS)
 ALLOCATE(TS(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 TS(:,:,:)=0.
 IF (ALLOCATED(QS)) DEALLOCATE(QS)
 ALLOCATE(QS(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS(:,:,:)=0.
 IF (ALLOCATED(QS_REL)) DEALLOCATE(QS_REL)
 ALLOCATE(QS_REL(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS_REL(:,:,:)=0.
 IF (ALLOCATED(PSIS)) DEALLOCATE(PSIS)
 ALLOCATE(PSIS(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 PSIS(:,:,:)=0.
 IF (ALLOCATED(KAPPA_H)) DEALLOCATE(KAPPA_H)
 ALLOCATE(KAPPA_H(NPOINT_X,NPOINT_Y,NLEV_SOIL_MAX))
 KAPPA_H(:,:,:)=0.
 IF (ALLOCATED(FRAC_SICE)) DEALLOCATE(FRAC_SICE)
 ALLOCATE(FRAC_SICE(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 FRAC_SICE(:,:,:)=0.
 IF (ALLOCATED(DFRAC_SICE_DT)) DEALLOCATE(DFRAC_SICE_DT)
 ALLOCATE(DFRAC_SICE_DT(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 DFRAC_SICE_DT(:,:,:)=0.
 IF (ALLOCATED(QSURF)) DEALLOCATE(QSURF)
 ALLOCATE(QSURF(NPOINT_X,NPOINT_Y))
 QSURF(:,:)=0.
 IF (ALLOCATED(TSURF)) DEALLOCATE(TSURF)
 ALLOCATE(TSURF(NPOINT_X,NPOINT_Y))
 TSURF(:,:)=0.
 IF (ALLOCATED(TFOREST)) DEALLOCATE(TFOREST)
 ALLOCATE(TFOREST(NPOINT_X,NPOINT_Y))
 TFOREST(:,:)=0.
 IF (ALLOCATED(QFOREST)) DEALLOCATE(QFOREST)
 ALLOCATE(QFOREST(NPOINT_X,NPOINT_Y))
 QFOREST(:,:)=0.
 IF (ALLOCATED(QS_REL_VEG)) DEALLOCATE(QS_REL_VEG)
 ALLOCATE(QS_REL_VEG(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS_REL_VEG(:,:,:)=0.

 IF (ALLOCATED(D_LEV_SNOW)) DEALLOCATE(D_LEV_SNOW)
 ALLOCATE(D_LEV_SNOW(NPOINT_X,NPOINT_Y))
 D_LEV_SNOW(:,:)=D_LEV_SNOW_BASE ! kg/m^2
 IF (ALLOCATED(TSNOW)) DEALLOCATE(TSNOW)
 ALLOCATE(TSNOW(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 TSNOW(:,:,:)=VAL_MISSING
 IF (ALLOCATED(LEV_SNOW)) DEALLOCATE(LEV_SNOW)
 ALLOCATE(LEV_SNOW(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 LEV_SNOW(:,:,:)=VAL_MISSING
 IF (ALLOCATED(FICE_SNOW)) DEALLOCATE(FICE_SNOW)
 ALLOCATE(FICE_SNOW(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 FICE_SNOW(:,:,:)=VAL_MISSING
 IF (ALLOCATED(RHOSNOW)) DEALLOCATE(RHOSNOW)
 ALLOCATE(RHOSNOW(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 RHOSNOW(:,:,:)=VAL_MISSING
 IF (ALLOCATED(SNOW_AGE)) DEALLOCATE(SNOW_AGE)
 ALLOCATE(SNOW_AGE(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 SNOW_AGE(:,:,:)=0.
 IF (ALLOCATED(SNOW_MELT_AGE)) DEALLOCATE(SNOW_MELT_AGE)
 ALLOCATE(SNOW_MELT_AGE(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 SNOW_MELT_AGE(:,:,:)=0.
 IF (ALLOCATED(DLEV_SNOW_BOTTOM)) DEALLOCATE(DLEV_SNOW_BOTTOM)
 ALLOCATE(DLEV_SNOW_BOTTOM(NPOINT_X,NPOINT_Y))
 DLEV_SNOW_BOTTOM(:,:)=0.
 IF (ALLOCATED(FLAG_SNOW_THICK)) DEALLOCATE(FLAG_SNOW_THICK)
 ALLOCATE(FLAG_SNOW_THICK(NPOINT_X,NPOINT_Y))
 FLAG_SNOW_THICK(:,:)=0

 IF (ALLOCATED(PREC_ACCUM)) DEALLOCATE(PREC_ACCUM)
 ALLOCATE(PREC_ACCUM(NPOINT_X,NPOINT_Y))
 PREC_ACCUM(:,:)=0.
 IF (ALLOCATED(FLUX_Q_ACCUM)) DEALLOCATE(FLUX_Q_ACCUM)
 ALLOCATE(FLUX_Q_ACCUM(NPOINT_X,NPOINT_Y))
 FLUX_Q_ACCUM(:,:)=0.
 IF (ALLOCATED(RUN_OFF_ACCUM)) DEALLOCATE(RUN_OFF_ACCUM)
 ALLOCATE(RUN_OFF_ACCUM(NPOINT_X,NPOINT_Y))
 RUN_OFF_ACCUM(:,:)=0.
 IF (ALLOCATED(FLUX_WATER_BOTTOM_ACC)) DEALLOCATE(FLUX_WATER_BOTTOM_ACC)
 ALLOCATE(FLUX_WATER_BOTTOM_ACC(NPOINT_X,NPOINT_Y))
 FLUX_WATER_BOTTOM_ACC(:,:)=0.

 IF (ALLOCATED(QS_OLD)) DEALLOCATE(QS_OLD)
 ALLOCATE(QS_OLD(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS_OLD(:,:,:)=0.
 IF (ALLOCATED(QSI)) DEALLOCATE(QSI)
 ALLOCATE(QSI(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QSI(:,:,:)=0.
 IF (ALLOCATED(DQS)) DEALLOCATE(DQS)
 ALLOCATE(DQS(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 DQS(:,:,:)=0.
 IF (ALLOCATED(RUN_OFF)) DEALLOCATE(RUN_OFF)
 ALLOCATE(RUN_OFF(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 RUN_OFF(:,:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM)) DEALLOCATE(FLUX_ENTROPY_ATM)
 ALLOCATE(FLUX_ENTROPY_ATM(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_SPEC)) DEALLOCATE(FLUX_ENTROPY_ATM_SPEC)
 ALLOCATE(FLUX_ENTROPY_ATM_SPEC(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_SPEC(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_LAT)) DEALLOCATE(FLUX_ENTROPY_ATM_LAT)
 ALLOCATE(FLUX_ENTROPY_ATM_LAT(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_LAT(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_RAD)) DEALLOCATE(FLUX_ENTROPY_ATM_RAD)
 ALLOCATE(FLUX_ENTROPY_ATM_RAD(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_RAD(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_SOIL)) DEALLOCATE(FLUX_ENTROPY_ATM_SOIL)
 ALLOCATE(FLUX_ENTROPY_ATM_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_SNOW)) DEALLOCATE(FLUX_ENTROPY_ATM_SNOW)
 ALLOCATE(FLUX_ENTROPY_ATM_SNOW(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_SNOW(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_ATM_FOREST)) DEALLOCATE(FLUX_ENTROPY_ATM_FOREST)
 ALLOCATE(FLUX_ENTROPY_ATM_FOREST(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_ATM_FOREST(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_SNOW_SOIL)) DEALLOCATE(FLUX_ENTROPY_SNOW_SOIL)
 ALLOCATE(FLUX_ENTROPY_SNOW_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_SNOW_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_ENTROPY_FOREST_SOIL)) DEALLOCATE(FLUX_ENTROPY_FOREST_SOIL)
 ALLOCATE(FLUX_ENTROPY_FOREST_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_ENTROPY_FOREST_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_WV)) DEALLOCATE(FLUX_WV)
 ALLOCATE(FLUX_WV(NPOINT_X,NPOINT_Y))
 FLUX_WV(:,:)=0.
 IF (ALLOCATED(FLUX_WV_SOIL)) DEALLOCATE(FLUX_WV_SOIL)
 ALLOCATE(FLUX_WV_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_WV_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_WV_SNOW)) DEALLOCATE(FLUX_WV_SNOW)
 ALLOCATE(FLUX_WV_SNOW(NPOINT_X,NPOINT_Y))
 FLUX_WV_SNOW(:,:)=0.
 IF (ALLOCATED(FLUX_WV_VEGLOW_DRY)) DEALLOCATE(FLUX_WV_VEGLOW_DRY)
 ALLOCATE(FLUX_WV_VEGLOW_DRY(NPOINT_X,NPOINT_Y))
 FLUX_WV_VEGLOW_DRY(:,:)=0.
 IF (ALLOCATED(FLUX_WV_VEGLOW_WET)) DEALLOCATE(FLUX_WV_VEGLOW_WET)
 ALLOCATE(FLUX_WV_VEGLOW_WET(NPOINT_X,NPOINT_Y))
 FLUX_WV_VEGLOW_WET(:,:)=0.
 IF (ALLOCATED(FLUX_WV_FOREST_DRY)) DEALLOCATE(FLUX_WV_FOREST_DRY)
 ALLOCATE(FLUX_WV_FOREST_DRY(NPOINT_X,NPOINT_Y))
 FLUX_WV_FOREST_DRY(:,:)=0.
 IF (ALLOCATED(FLUX_WV_FOREST_WET)) DEALLOCATE(FLUX_WV_FOREST_WET)
 ALLOCATE(FLUX_WV_FOREST_WET(NPOINT_X,NPOINT_Y))
 FLUX_WV_FOREST_WET(:,:)=0.
 IF (ALLOCATED(FLUX_PREC_LIQ)) DEALLOCATE(FLUX_PREC_LIQ)
 ALLOCATE(FLUX_PREC_LIQ(NPOINT_X,NPOINT_Y))
 FLUX_PREC_LIQ(:,:)=0.
 IF (ALLOCATED(FLUX_PREC_ICE)) DEALLOCATE(FLUX_PREC_ICE)
 ALLOCATE(FLUX_PREC_ICE(NPOINT_X,NPOINT_Y))
 FLUX_PREC_ICE(:,:)=0.
 IF (ALLOCATED(FLUX_W_LIQ_ATM_SOIL)) DEALLOCATE(FLUX_W_LIQ_ATM_SOIL)
 ALLOCATE(FLUX_W_LIQ_ATM_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_W_LIQ_ATM_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_W_LIQ_ATM_SNOW)) DEALLOCATE(FLUX_W_LIQ_ATM_SNOW)
 ALLOCATE(FLUX_W_LIQ_ATM_SNOW(NPOINT_X,NPOINT_Y))
 FLUX_W_LIQ_ATM_SNOW(:,:)=0.
 IF (ALLOCATED(FLUX_W_ICE_ATM_SNOW)) DEALLOCATE(FLUX_W_ICE_ATM_SNOW)
 ALLOCATE(FLUX_W_ICE_ATM_SNOW(NPOINT_X,NPOINT_Y))
 FLUX_W_ICE_ATM_SNOW(:,:)=0.
 IF (ALLOCATED(FLUX_W_LIQ_ATM_VEGLOW)) DEALLOCATE(FLUX_W_LIQ_ATM_VEGLOW)
 ALLOCATE(FLUX_W_LIQ_ATM_VEGLOW(NPOINT_X,NPOINT_Y))
 FLUX_W_LIQ_ATM_VEGLOW(:,:)=0.
 IF (ALLOCATED(FLUX_W_LIQ_ATM_FOREST)) DEALLOCATE(FLUX_W_LIQ_ATM_FOREST)
 ALLOCATE(FLUX_W_LIQ_ATM_FOREST(NPOINT_X,NPOINT_Y))
 FLUX_W_LIQ_ATM_FOREST(:,:)=0.
 IF (ALLOCATED(FLUX_WATER_SNOW_SOIL)) DEALLOCATE(FLUX_WATER_SNOW_SOIL)
 ALLOCATE(FLUX_WATER_SNOW_SOIL(NPOINT_X,NPOINT_Y))
 FLUX_WATER_SNOW_SOIL(:,:)=0.
 IF (ALLOCATED(FLUX_SOIL_ENTROPY)) DEALLOCATE(FLUX_SOIL_ENTROPY)
 ALLOCATE(FLUX_SOIL_ENTROPY(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 FLUX_SOIL_ENTROPY(:,:,:)=0.
 IF (ALLOCATED(FLUX_SOIL_WATER)) DEALLOCATE(FLUX_SOIL_WATER)
 ALLOCATE(FLUX_SOIL_WATER(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 FLUX_SOIL_WATER(:,:,:)=0.
 IF (ALLOCATED(FLUX_SOIL_WATER_VEG)) DEALLOCATE(FLUX_SOIL_WATER_VEG)
 ALLOCATE(FLUX_SOIL_WATER_VEG(NPOINT_X,NPOINT_Y,NLEV_SOIL_MAX))
 FLUX_SOIL_WATER_VEG(:,:,:)=0.
 IF (ALLOCATED(FLUX_SOIL_WATER_BOTTOM_DIAG)) DEALLOCATE(FLUX_SOIL_WATER_BOTTOM_DIAG)
 ALLOCATE(FLUX_SOIL_WATER_BOTTOM_DIAG(NPOINT_X,NPOINT_Y))
 FLUX_SOIL_WATER_BOTTOM_DIAG(:,:)=0.
 IF (ALLOCATED(LAMBDAG)) DEALLOCATE(LAMBDAG)
 ALLOCATE(LAMBDAG(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 LAMBDAG(:,:,:)=0.
 IF (ALLOCATED(FLUX_SNOW_ENTROPY)) DEALLOCATE(FLUX_SNOW_ENTROPY)
 ALLOCATE(FLUX_SNOW_ENTROPY(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 FLUX_SNOW_ENTROPY(:,:,:)=0.
 IF (ALLOCATED(FLUX_SNOW_WATER)) DEALLOCATE(FLUX_SNOW_WATER)
 ALLOCATE(FLUX_SNOW_WATER(NPOINT_X,NPOINT_Y,0:NLEV_SNOW))
 FLUX_SNOW_WATER(:,:,:)=0.

 IF (ALLOCATED(QSMAX)) DEALLOCATE(QSMAX)
 ALLOCATE(QSMAX(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QSMAX(:,:,:)=0.
 IF (ALLOCATED(QSMIN)) DEALLOCATE(QSMIN)
 ALLOCATE(QSMIN(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QSMIN(:,:,:)=0.
 IF (ALLOCATED(CG)) DEALLOCATE(CG)
 ALLOCATE(CG(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 CG(:,:,:)=0.
 IF (ALLOCATED(RHOG)) DEALLOCATE(RHOG)
 ALLOCATE(RHOG(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 RHOG(:,:,:)=0.
 IF (ALLOCATED(PSIS_SAT)) DEALLOCATE(PSIS_SAT)
 ALLOCATE(PSIS_SAT(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 PSIS_SAT(:,:,:)=0.
 IF (ALLOCATED(KAPPA_SAT)) DEALLOCATE(KAPPA_SAT)
 ALLOCATE(KAPPA_SAT(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 KAPPA_SAT(:,:,:)=0.
 IF (ALLOCATED(PAR_B)) DEALLOCATE(PAR_B)
 ALLOCATE(PAR_B(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 PAR_B(:,:,:)=0.
 IF (ALLOCATED(PAR_C)) DEALLOCATE(PAR_C)
 ALLOCATE(PAR_C(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 PAR_C(:,:,:)=0.
 IF (ALLOCATED(PAR_FB)) DEALLOCATE(PAR_FB)
 ALLOCATE(PAR_FB(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 PAR_FB(:,:,:)=0.
 IF (ALLOCATED(QS_REL_VEGWILT)) DEALLOCATE(QS_REL_VEGWILT)
 ALLOCATE(QS_REL_VEGWILT(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS_REL_VEGWILT(:,:,:)=0.
 IF (ALLOCATED(QS_REL_VEGREF)) DEALLOCATE(QS_REL_VEGREF)
 ALLOCATE(QS_REL_VEGREF(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX))
 QS_REL_VEGREF(:,:,:)=0.

 IF (ALLOCATED(LAI_VEGLOW)) DEALLOCATE(LAI_VEGLOW)
 ALLOCATE(LAI_VEGLOW(NPOINT_X,NPOINT_Y))
 LAI_VEGLOW(:,:)=0.
 IF (ALLOCATED(LAI_FOREST)) DEALLOCATE(LAI_FOREST)
 ALLOCATE(LAI_FOREST(NPOINT_X,NPOINT_Y))
 LAI_FOREST(:,:)=0.
 IF (ALLOCATED(FVEGLOW)) DEALLOCATE(FVEGLOW)
 ALLOCATE(FVEGLOW(NPOINT_X,NPOINT_Y))
 FVEGLOW(:,:)=0.
 IF (ALLOCATED(FFOREST)) DEALLOCATE(FFOREST)
 ALLOCATE(FFOREST(NPOINT_X,NPOINT_Y))
 FFOREST(:,:)=0.
 IF (ALLOCATED(FVEGLOWWET)) DEALLOCATE(FVEGLOWWET)
 ALLOCATE(FVEGLOWWET(NPOINT_X,NPOINT_Y))
 FVEGLOWWET(:,:)=0.
 IF (ALLOCATED(FFORESTWET)) DEALLOCATE(FFORESTWET)
 ALLOCATE(FFORESTWET(NPOINT_X,NPOINT_Y))
 FFORESTWET(:,:)=0.
 IF (ALLOCATED(ROOT_DEPTH_VEGLOW)) DEALLOCATE(ROOT_DEPTH_VEGLOW)
 ALLOCATE(ROOT_DEPTH_VEGLOW(NPOINT_X,NPOINT_Y))
 ROOT_DEPTH_VEGLOW(:,:)=0.
 IF (ALLOCATED(QW_VEGLOW)) DEALLOCATE(QW_VEGLOW)
 ALLOCATE(QW_VEGLOW(NPOINT_X,NPOINT_Y))
 QW_VEGLOW(:,:)=0.
 IF (ALLOCATED(QW_FOREST)) DEALLOCATE(QW_FOREST)
 ALLOCATE(QW_FOREST(NPOINT_X,NPOINT_Y))
 QW_FOREST(:,:)=0.
 IF (ALLOCATED(ROOT_DEPTH_FOREST)) DEALLOCATE(ROOT_DEPTH_FOREST)
 ALLOCATE(ROOT_DEPTH_FOREST(NPOINT_X,NPOINT_Y))
 ROOT_DEPTH_FOREST(:,:)=0.
 IF (ALLOCATED(QW_VEGLOW_MAX)) DEALLOCATE(QW_VEGLOW_MAX)
 ALLOCATE(QW_VEGLOW_MAX(NPOINT_X,NPOINT_Y))
 QW_VEGLOW_MAX(:,:)=0.
 IF (ALLOCATED(QW_FOREST_MAX)) DEALLOCATE(QW_FOREST_MAX)
 ALLOCATE(QW_FOREST_MAX(NPOINT_X,NPOINT_Y))
 QW_FOREST_MAX(:,:)=0.
 IF (ALLOCATED(DZSUM_VEGLOW)) DEALLOCATE(DZSUM_VEGLOW)
 ALLOCATE(DZSUM_VEGLOW(NPOINT_X,NPOINT_Y))
 DZSUM_VEGLOW(:,:)=0.
 IF (ALLOCATED(DZSUM_FOREST)) DEALLOCATE(DZSUM_FOREST)
 ALLOCATE(DZSUM_FOREST(NPOINT_X,NPOINT_Y))
 DZSUM_FOREST(:,:)=0.

 IF (ALLOCATED(FSNCOVER)) DEALLOCATE(FSNCOVER)
 ALLOCATE(FSNCOVER(NPOINT_X,NPOINT_Y))
 FSNCOVER(:,:)=0.

 IF (ALLOCATED(DZ_DZSUM_VEGLOW)) DEALLOCATE(DZ_DZSUM_VEGLOW)
 ALLOCATE(DZ_DZSUM_VEGLOW(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX-1))
 DZ_DZSUM_VEGLOW(:,:,:)=0.
 IF (ALLOCATED(DZ_DZSUM_FOREST)) DEALLOCATE(DZ_DZSUM_FOREST)
 ALLOCATE(DZ_DZSUM_FOREST(NPOINT_X,NPOINT_Y,0:NLEV_SOIL_MAX-1))
 DZ_DZSUM_FOREST(:,:,:)=0.

 IF (ALLOCATED(MASK_SOIL)) DEALLOCATE(MASK_SOIL)
 ALLOCATE(MASK_SOIL(NPOINT_X,NPOINT_Y))
 MASK_SOIL(:,:)=0
 IF (ALLOCATED(MASK_SOIL_AXED)) DEALLOCATE(MASK_SOIL_AXED)
 ALLOCATE(MASK_SOIL_AXED(NPOINT_X,NPOINT_Y))
 MASK_SOIL_AXED(:,:)=0

! <---

 READ (IUNIT,*) NLEV_SOIL_HEAT(I,J), NLEV_SOIL_WATER(I,J)
 READ (IUNIT,*) QSMAX(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QSMIN(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) CG(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) RHOG(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) PSIS_SAT(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) KAPPA_SAT(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) PAR_B(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) PAR_C(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) PAR_FB(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QS_REL_VEGWILT(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QS_REL_VEGREF(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) LAI_VEGLOW(I,J)
 READ (IUNIT,*) LAI_FOREST(I,J)
 READ (IUNIT,*) FVEGLOW(I,J)
 READ (IUNIT,*) FFOREST(I,J)
 READ (IUNIT,*) FVEGLOWWET(I,J)
 READ (IUNIT,*) FFORESTWET(I,J)
 READ (IUNIT,*) ROOT_DEPTH_VEGLOW(I,J)
 READ (IUNIT,*) ROOT_DEPTH_FOREST(I,J)
 READ (IUNIT,*) QW_VEGLOW(I,J)
 READ (IUNIT,*) QW_FOREST(I,J)
 READ (IUNIT,*) QW_VEGLOW_MAX(I,J)
 READ (IUNIT,*) QW_FOREST_MAX(I,J)
 READ (IUNIT,*) DZSUM_VEGLOW(I,J)
 READ (IUNIT,*) DZSUM_FOREST(I,J)
 READ (IUNIT,*) DZ_DZSUM_VEGLOW(I,J,0:NLEV_SOIL_MAX-1)
 READ (IUNIT,*) DZSUM_FOREST(I,J)
 READ (IUNIT,*) DZ_DZSUM_FOREST(I,J,0:NLEV_SOIL_MAX-1)
 READ (IUNIT,*) MASK_SOIL(I,J)
 READ (IUNIT,*) MASK_SOIL_AXED(I,J)
 READ (IUNIT,*) TSURF(I,J)
 READ (IUNIT,*) TS(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QS(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QS_REL(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) QSI(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) FRAC_SICE(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) DFRAC_SICE_DT(I,J,0:NLEV_SOIL_MAX)
 READ (IUNIT,*) FSNCOVER(I,J)
 READ (IUNIT,*) LEV_SNOW(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) TSNOW(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) FICE_SNOW(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) RHOSNOW(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) SNOW_AGE(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) SNOW_MELT_AGE(I,J,0:NLEV_SNOW)
 READ (IUNIT,*) D_LEV_SNOW(I,J)
 READ (IUNIT,*) FLAG_SNOW_THICK(I,J)
 READ (IUNIT,*) DLEV_SNOW_BOTTOM(I,J)
 READ (IUNIT,*) TFOREST(I,J)

 CLOSE (11) 

 DO K=0,NLEV_SNOW
   DO J=1,NPOINT_Y
     DO I=1,NPOINT_X
       IF (ABS(LEV_SNOW(I,J,K)-VAL_MISSING) < 10.) LEV_SNOW(I,J,K)=VAL_MISSING
       IF (ABS(TSNOW(I,J,K)-VAL_MISSING) < 10.) TSNOW(I,J,K)=VAL_MISSING
       IF (ABS(FICE_SNOW(I,J,K)-VAL_MISSING) < 10.) FICE_SNOW(I,J,K)=VAL_MISSING
       IF (ABS(RHOSNOW(I,J,K)-VAL_MISSING) < 10.) RHOSNOW(I,J,K)=VAL_MISSING
     ENDDO
   ENDDO
 ENDDO

 PAR_FA=-4./ALOG((T0+T_FICE_CRIT)/T0)

 PRINT *
 PRINT *," Restart, read from ",TRIM(FILE_IN)
 PRINT *

RETURN
END SUBROUTINE RDRF_POCHVA
!-------------------------------------------------------------------------------
