!==============================================================================
! This source code is part of the
! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model.
! This work is licensed under the CSIRO Open Source Software License
! Agreement (variation of the BSD / MIT License).
!
! You may not use this file except in compliance with this License.
! A copy of the License (CSIRO_BSD_MIT_License_v2.0_CABLE.txt) is located
! in each directory containing CABLE code.
!
! ==============================================================================
! Purpose: tracking climate variables for use in phenology and potential pft modules
!
! Called from: cable_serial
!
! History: Vanessa Haverd Jan 2015

! ==============================================================================

# ifndef UM_CBL
MODULE cable_climate_mod

  USE cable_def_types_mod, ONLY: met_type, climate_type, canopy_type, mp, &
       r_2, alloc_cbm_var, air_type, radiation_type
  USE TypeDef,              ONLY: i4b, dp
  USE cable_IO_vars_module, ONLY: patch
  USE casa_ncdf_module, ONLY: HANDLE_ERR
  USE CABLE_COMMON_MODULE, ONLY: CurYear, filename, cable_user

CONTAINS
  ! ==============================================================================


  SUBROUTINE cable_climate(ktau,kstart,kend,ktauday,idoy,LOY,met,climate, canopy, &
       air, rad, dels, np)


    IMPLICIT NONE

    INTEGER,      INTENT(IN) :: ktau ! integration step number
    INTEGER,      INTENT(IN) :: kstart ! starting value of ktau
    INTEGER,      INTENT(IN) :: kend ! total # timesteps in run

    INTEGER,      INTENT(IN)                  :: idoy ,LOY ! day of year (1-365) , Length oy
    INTEGER,      INTENT(IN)                  :: ktauday
    TYPE (met_type), INTENT(IN)       :: met  ! met input variables
    TYPE (climate_type), INTENT(INOUT)       :: climate  ! climate variables
    TYPE (canopy_type), INTENT(IN) :: canopy ! vegetation variables
    TYPE (air_type), INTENT(IN)       :: air
    TYPE (radiation_type), INTENT(IN)  :: rad        ! radiation variables
    REAL, INTENT(IN)               :: dels ! integration time setp (s)
    INTEGER,      INTENT(IN)                  :: np
    INTEGER :: d, y, k
    INTEGER, PARAMETER:: COLDEST_DAY_NHEMISPHERE = 355
    INTEGER, PARAMETER:: COLDEST_DAY_SHEMISPHERE = 172
    REAL, PARAMETER:: CoeffPT = 1.26
    REAL,      DIMENSION(mp)  :: mtemp_last, mmoist_last, phiEq, ppc, EpsA, RhoA
    INTEGER :: startyear
    INTEGER :: MonthDays(12)
    INTEGER::  DaysInMonth, nmonth, tmp
    LOGICAL :: IsLastDay ! last day of month?
    REAL, PARAMETER:: Gaero = 0.015  ! (m s-1) aerodynmaic conductance (for use in PT evap)
    REAL, PARAMETER:: Capp   = 29.09    ! isobaric spec heat air    [J/molA/K]
    REAL, PARAMETER:: SBoltz  = 5.67e-8  ! Stefan-Boltzmann constant [W/m2/K4]
    climate%doy = idoy

!! * Find irradiances, available energy, equilibrium latent heat flux
!PPc    = Gaero / ( Gaero + 4.0*SBoltz*((TempA+273.16)**3)/(RhoA*Capp) )
!                                                    ! PPc = Ga/(Ga+Gr)      [-]
!EpsA   = Epsif(TempA, Pmb)                          ! Epsi at TempA         [-]
!PhiSd  = SolarMJ * 1.0e6 / (DayltFrac*SecDay)       ! daylt down solar      [W/m2]
!PhiLd  = 335.97 * (((TempA + 273.16) / 293.0)**6)   ! daylt down thermal    [W/m2]
!                                                    !   (Swinbank formula)
!PhiAi  = (1.0-Albedo)*PhiSd + Emis *    &           ! daylt iso-avail engy  [W/m2]
!         (PhiLd - SBoltz*((TempA + 273.16)**4))     !   (veg + soil)
!PhiEq  = PhiAi * (PPc*EpsA) / (PPc*EpsA + 1.0)      ! equil ltnt heat flux  [W/m2]
!PhiEq  = max(PhiEq, 1.0)                            ! PhiEq > +1 W/m2, so non-negative
!                                    ! precipitation         [m/day]
!FWPT   = CoeffPT * PhiEq * ((DayltFrac*SecDay) / (RhoW*Rlat))

    ! accumulate annual evaporation and potential evaporation
    !ppc = 1.0
    RhoA = met%pmb * 100.0 / (8.314 * (met%tk)) ! air density [molA/m3]
    PPc    = Gaero / ( Gaero + 4.0*SBoltz*((met%tk**3)/(RhoA*Capp) ))

    EpSA = Epsif(met%tk - 273.16, met%pmb)
    phiEq = canopy%rniso * (PPc*EpsA) / (PPc*EpsA + 1.0)      ! equil ltnt heat flux  [W/m2]

    IF (idoy==1 .AND. MOD(ktau,ktauday)==1 ) THEN
       !  climate%evap_PT =  max(phiEq,1.0)*CoeffPT/air%rlam*dels  ! mm
       climate%evap_PT =  phiEq*CoeffPT/air%rlam*dels  ! mm
       !  climate%evap_PT = canopy%epot  ! mm
       !  climate%aevap  =   canopy%fe/air%rlam*dels ! mm
       climate%aevap = met%precip ! mm
    ELSE

       climate%evap_PT = climate%evap_PT + MAX(phiEq,1.0)*CoeffPT/air%rlam*dels  ! mm
       ! climate%evap_PT =climate%evap_PT + canopy%epot  ! mm
       ! climate%aevap =  climate%aevap + canopy%fe/air%rlam*dels ! mm
       climate%aevap = climate%aevap + met%precip ! mm
    ENDIF

    ! accumulate daily temperature, evap and potential evap
    IF(MOD(ktau,ktauday)==1) THEN
       climate%dtemp = met%tk - 273.15
       climate%dmoist = canopy%fwsoil
    ELSE
       climate%dtemp = climate%dtemp + met%tk - 273.15
       climate%dmoist = climate%dmoist + canopy%fwsoil
    ENDIF

    IF(MOD((ktau-kstart+1),ktauday)==0) THEN  ! end of day
       climate%dtemp = climate%dtemp/FLOAT(ktauday)
       climate%dmoist = climate%dmoist/FLOAT(ktauday)
    ENDIF



    IF(MOD((ktau-kstart+1),ktauday)==0) THEN  ! end of day
       ! get month and check if end of month
       IsLastDay = .FALSE.
       MonthDays = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       IF (LOY==366) MonthDays(2) = MonthDays(2) + 1
       nmonth = 1
       tmp = MonthDays(1)
       DO WHILE(tmp.LT.idoy)
          tmp = tmp +  MonthDays(nmonth+1)
          nmonth = nmonth + 1
       ENDDO

       IF (idoy == SUM(MonthDays(1:nmonth))) IsLastDay = .TRUE.

       ! On first day of year ...
       IF (idoy==1) THEN

          ! ... reset annual GDD5 counter
          climate%agdd5=0.0
          climate%agdd0=0.0
          climate%evap_PT = 0     ! annual PT evap [mm]
          climate%aevap  = 0      ! annual evap [mm]


       ENDIF

       !jhan: see CABLE Ticket#149 and above CABLE_LSM: comment
# ifndef ESM15
       WHERE ((patch%latitude>=0.0 .AND. idoy==COLDEST_DAY_NHEMISPHERE).OR. &
            (patch%latitude<0.0 .AND. idoy==COLDEST_DAY_SHEMISPHERE) )

          ! In midwinter, reset GDD counter for summergreen phenology
          climate%gdd5=0.0
          climate%gdd0=0.0
       END WHERE
# endif
       ! Update GDD counters and chill day count
       climate%gdd0 = climate%gdd0 + MAX(0.0,climate%dtemp-0.0)
       climate%agdd0= climate%agdd0 + MAX(0.0,climate%dtemp-0.0)

       climate%gdd5 = climate%gdd5 + MAX(0.0,climate%dtemp-5.0)
       climate%agdd5= climate%agdd5 + MAX(0.0,climate%dtemp-5.0)
       WHERE (climate%dtemp<5.0 .AND. climate%chilldays<=365)
          climate%chilldays = climate%chilldays + 1
       ENDWHERE

       ! Save yesterday's mean temperature for the last month
       mtemp_last=climate%mtemp

       ! Update daily temperatures, and mean overall temperature, for last 31 days

       climate%mtemp=climate%dtemp
       climate%qtemp=climate%dtemp
       climate%mmoist = climate%dmoist
       DO d=1,30
          climate%dtemp_31(:,d)=climate%dtemp_31(:,d+1)
          climate%mtemp= climate%mtemp + climate%dtemp_31(:,d)
          climate%dmoist_31(:,d)=climate%dmoist_31(:,d+1)
          climate%mmoist = climate%mmoist + climate%dmoist_31(:,d)
       ENDDO
       DO d=1,90
          climate%dtemp_91(:,d)=climate%dtemp_91(:,d+1)
          climate%qtemp= climate%qtemp + climate%dtemp_91(:,d)
       ENDDO
       climate%dtemp_31(:,31)=climate%dtemp
       climate%dtemp_91(:,91)=climate%dtemp
       climate%dmoist_31(:,31)=climate%dmoist
       climate%qtemp = climate%qtemp/91.0  ! average temperature over the last quarter
       climate%mtemp = climate%mtemp/31.0 ! average temperature over the last month
       climate%mmoist = climate%mmoist/31.0 ! average moisture index over the last month

       ! Reset GDD and chill day counter if mean monthly temperature falls below base
       ! temperature

       WHERE(mtemp_last>=5.0 .AND. climate%mtemp<5.0)
          climate%gdd5=0.0
          climate%chilldays=0
       ENDWHERE

       ! On last day of month ...

       IF (IsLastDay) THEN

          ! Update mean temperature for the last 12 months
          ! atemp_mean_new = atemp_mean_old * (11/12) + mtemp * (1/12)

          climate%atemp_mean=climate%atemp_mean*(11./12.)+climate%mtemp*(1./12.)

          ! Record minimum and maximum monthly temperatures

          IF (nmonth==1) THEN
             climate%mtemp_min=climate%mtemp;
             climate%mtemp_max=climate%mtemp;
             climate%qtemp_max_last_year =  climate%qtemp_max;
             climate%qtemp_max=climate%qtemp;
          ELSE
             WHERE (climate%mtemp<climate%mtemp_min) &
                  climate%mtemp_min=climate%mtemp
             WHERE (climate%mtemp>climate%mtemp_max) &
                  climate%mtemp_max=climate%mtemp
             WHERE (climate%qtemp>climate%qtemp_max) &
                  climate%qtemp_max=climate%qtemp
          ENDIF  ! first month of year

          ! On 31 December update records of minimum monthly temperatures for the last
          ! 20 years and find minimum monthly temperature for the last 20 years

          IF (nmonth==12) THEN
             climate%nyears = climate%nyears +1


             startyear=20-MIN(19,climate%nyears-1)
             climate%mtemp_min20=0.0
             climate%mtemp_max20=0.0
             climate%alpha_PT20 = 0.0

             IF (startyear<20) THEN
                DO y=startyear,19
                   climate%mtemp_min_20(:,y)=climate%mtemp_min_20(:,y+1)
                   climate%mtemp_min20=climate%mtemp_min20+ &
                        climate%mtemp_min_20(:,y)
                   climate%mtemp_max_20(:,y)=climate%mtemp_max_20(:,y+1)
                   climate%mtemp_max20 =climate%mtemp_max20 + &
                        climate%mtemp_max_20(:,y)
                   climate%alpha_PT_20(:,y)=climate%alpha_PT_20(:,y+1)
                   climate%alpha_PT20 =climate%alpha_PT20 + &
                        climate%alpha_PT_20(:,y)
                ENDDO

                climate%mtemp_min20=climate%mtemp_min20/REAL(20-startyear)
                climate%mtemp_max20=climate%mtemp_max20/REAL(20-startyear)
                climate%alpha_PT20=climate%alpha_PT20/REAL(20-startyear)
             ELSE
                ! only occurs when climate%nyears = 1
                climate%mtemp_min20 = climate%mtemp_min
                climate%mtemp_max20 = climate%mtemp_max
                climate%alpha_PT20 = climate%alpha_PT
             ENDIF

             climate%mtemp_min_20(:,20)=climate%mtemp_min
             climate%mtemp_max_20(:,20)=climate%mtemp_max


             climate%alpha_PT = MAX(climate%aevap/climate%evap_PT, 0.0)     ! ratio of annual evap to annual PT evap
             climate%alpha_PT_20(:,20)=climate%alpha_PT


             CALL biome1_pft(climate,np)

          ENDIF  ! last month of year



       ENDIF     ! last day of month


    ENDIF ! end of day

  END SUBROUTINE cable_climate
  !=============================================================================

  ELEMENTAL FUNCTION Epsif(TC,Pmb)
    !-------------------------------------------------------------------------------
    ! At temperature TC [deg C] and pressure Pmb [mb], return
    ! epsi = (RLAM/CAPP) * d(sat spec humidity)/dT [(kg/kg)/K], from Teten formula.
    ! MRR, xx/1987, 27-jan-94
    ! PRB, 09/1999:   Convert to F95 elemental function; works on scalars and arrays
    !                 just like intrinsic functions.
    ! MRR, 28-mar-05: use Rlat [J/molW], Capp [J/molA/K] from MODULE Constants,
    !                 to ensure consistency with other uses of Rlat, Capp
    ! MRR, 28-mar-05: Remove dependence of Rlat (latent heat vaporisation of water)
    !                 on temperature, use value at 20 C
    !-------------------------------------------------------------------------------
    !USE TypeDef
    !USE Constants
    IMPLICIT NONE
    REAL,INTENT(in):: TC, Pmb       ! temp [deg C], pressure [mb]
    REAL:: Epsif                    ! epsi
    REAL:: TCtmp, ES, dESdT         ! local
    REAL,PARAMETER:: A = 6.106      ! Teten coefficients
    REAL,PARAMETER:: B = 17.27      ! Teten coefficients
    REAL,PARAMETER:: C = 237.3      ! Teten coefficients
    REAL,PARAMETER:: Rlat      = 44140.0  ! lat heat evap H2O at 20C  [J/molW]
    REAL,PARAMETER:: Capp      = 29.09    ! isobaric spec heat air    [J/molA/K]
    !-------------------------------------------------------------------------------
    TCtmp = TC                          ! preserve TC
    IF (TCtmp.GT.100.0) TCtmp = 100.0   ! constrain TC to (-40.0,100.0)
    IF (TCtmp.LT.-40.0) TCtmp = -40.0
    ES    = A*EXP(B*TCtmp/(C+TCtmp))    ! sat vapour pressure
    dESdT = ES*B*C/(C+TCtmp)**2         ! d(sat VP)/dT: (mb/K)
    Epsif = (Rlat/Capp) * dESdT / Pmb   ! dimensionless (ES/Pmb = molW/molA)

  END FUNCTION Epsif

  !=============================================================================
  !=============================================================================
  !
  SUBROUTINE biome1_pft(climate, np)
    IMPLICIT NONE

    TYPE (climate_type), INTENT(INOUT)       :: climate  ! climate variables
    INTEGER, INTENT(IN) :: np
    INTEGER :: k, j, npft
    INTEGER, ALLOCATABLE :: pft_biome1(:,:)
    REAL, ALLOCATABLE:: alpha_PT_scaled(:)

    ! TABLE 1 , Prentice et al. J. Biogeog., 19, 117-134, 1992
    ! pft_biome1: Trees (1)tropical evergreen; (2) tropical raingreen; (3) warm temp evergreen ;
    ! (4) temperate summergreen; (5) cool-temperate conifer; (6) boreal evergreen conifer;
    ! (7) boreal summergreen;
    ! Non-trees: (8) sclerophyll/succulent; (9) warm grass/shrub; (10) cool grass/shrub;
    ! (11) cold grass/shrub; (12) hot desert shrub; (13) cold desert shrub.

    ALLOCATE(pft_biome1(np,4))
    ALLOCATE(alpha_PT_scaled(np))
    alpha_PT_scaled =  MIN(climate%alpha_PT20, 1.0)

    DO k=1,np

       pft_biome1(k,:) = 999

       IF (climate%mtemp_min20(k) .GE. 15.5) THEN
          IF (alpha_PT_scaled(k).GE.0.80) THEN
             pft_biome1(k,1) = 1
             IF (alpha_PT_scaled(k).LE.0.85) THEN
                pft_biome1(k,2) = 2
             ENDIF
          ELSEIF (alpha_PT_scaled(k).GE.0.4 .AND. alpha_PT_scaled(k).LT.0.80) THEN
             pft_biome1(k,1) = 2
          ENDIF
       ENDIF


       IF (climate%mtemp_min20(k).GE.5 .AND.alpha_PT_scaled(k).GE.0.4 &
            .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 3
       ENDIF


       IF (climate%mtemp_min20(k).GE.-15 .AND. climate%mtemp_min20(k).LE.15.5 .AND. &
            alpha_PT_scaled(k).GE.0.35 .AND. climate%agdd5(k).GT.1200 & !
            .AND. pft_biome1(k,1).GT.3 ) THEN
          pft_biome1(k,1) = 4
       ENDIF

       IF (climate%mtemp_min20(k).GE.-19 .AND. climate%mtemp_min20(k).LE.5 .AND. &
            alpha_PT_scaled(k).GE.0.35 .AND. climate%agdd5(k).GT.900 )  THEN
          IF (pft_biome1(k,1).GT.4 ) THEN
             pft_biome1(k,1) = 5
          ELSEIF (pft_biome1(k,1).EQ.4 ) THEN
             pft_biome1(k,2) = 5
          ENDIF
       ENDIF

       IF (climate%mtemp_min20(k).GE.-35 .AND. climate%mtemp_min20(k).LE.-2 .AND. &
            alpha_PT_scaled(k).GE.0.35 .AND. climate%agdd5(k).GT.350 )  THEN
          IF (pft_biome1(k,1).EQ.999 ) THEN
             pft_biome1(k,1) = 6
          ELSEIF (pft_biome1(k,2).EQ.999 ) THEN
             pft_biome1(k,2) = 6
          ELSE
             pft_biome1(k,3) = 6
          ENDIF
       ENDIF

       IF ( climate%mtemp_min20(k).LE. 5 .AND. &
            alpha_PT_scaled(k).GE.0.45 .AND. climate%agdd5(k).GT.350 )  THEN
          IF (pft_biome1(k,1).EQ.999 ) THEN
             pft_biome1(k,1) = 7
          ELSEIF (pft_biome1(k,2).EQ.999 ) THEN
             pft_biome1(k,2) = 7
          ELSEIF (pft_biome1(k,3).EQ.999 ) THEN
             pft_biome1(k,3) = 7
          ELSE
             pft_biome1(k,4) = 7
          ENDIF
       ENDIF

       IF (climate%mtemp_min20(k).GE.5 .AND.alpha_PT_scaled(k).GE.0.2 &
            .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 8
       ENDIF

       IF (climate%mtemp_max20(k).GE.22 .AND.alpha_PT_scaled(k).GE.0.1 &
            .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 9
       ENDIF

       IF (climate%agdd5(k).GE.500 .AND.alpha_PT_scaled(k).GE.0.33 &
            .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 10
       ENDIF

       IF (climate%agdd0(k).GE.100 .AND.alpha_PT_scaled(k).GE.0.33) THEN
          IF (pft_biome1(k,1).EQ.999 ) THEN
             pft_biome1(k,1) = 11
          ELSEIF (pft_biome1(k,1).EQ.10) THEN
             pft_biome1(k,2) = 11
          ENDIF
       ENDIF

       IF (climate%mtemp_max20(k).GE.22 .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 12
       ENDIF

       IF (climate%agdd0(k).GE.100 .AND. pft_biome1(k,1).EQ.999 ) THEN
          pft_biome1(k,1) = 13
       ENDIF

       ! end of evironmental constraints on pft
       npft = 0
       DO j=1,4
          IF (pft_biome1(k,j).NE.999) npft = npft+1
       ENDDO
       !     MAP to Biome1 biome and CABLE pft
       ! (1) Tropical Rainforest
       IF(pft_biome1(k,1)==1 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 1
          climate%iveg(k) = 2
       ENDIF

       ! (2) Tropical Seasonal forest
       IF(pft_biome1(k,1)==1 .AND.pft_biome1(k,2)==2.AND. npft .EQ.2 ) THEN
          climate%biome(k) = 2
          climate%iveg(k) = 2
       ENDIF

       ! (3) Tropical dry forest/savanna
       IF(pft_biome1(k,1)==2.AND. npft .EQ.1 ) THEN
          climate%biome(k) = 3
          climate%iveg(k) = 2  ! N.B. need to include c4 grass
       ENDIF


       ! (4) Broad-leaved evergreen/warm mixed-forest
       IF(pft_biome1(k,1)==3.AND. npft .EQ.1 ) THEN
          climate%biome(k) = 4
          climate%iveg(k) = 2
       ENDIF

       ! (5) Temperate deciduous forest
       IF(pft_biome1(k,1)==4.AND.pft_biome1(k,2)==5.AND. &
            pft_biome1(k,3)==7 .AND. npft .EQ.3 ) THEN
          climate%biome(k) = 5
          climate%iveg(k) = 4
       ENDIF

       ! (6) Cool mixed forest
       IF(pft_biome1(k,1)==4.AND.pft_biome1(k,2)==5.AND. &
            pft_biome1(k,3)==6 .AND.  pft_biome1(k,4)==7 &
            .AND. npft .EQ.4 ) THEN
          climate%biome(k) = 6
          climate%iveg(k) = 4
       ENDIF

       ! (7) Cool conifer forest

       IF(pft_biome1(k,1)==5.AND.pft_biome1(k,2)==6.AND. &
            pft_biome1(k,3)==7 .AND. npft .EQ.3 ) THEN
          climate%biome(k) = 7
          climate%iveg(k) = 1
       ENDIF
       ! (8) Taiga
       IF(pft_biome1(k,1)==6.AND.pft_biome1(k,2)==7 .AND. npft .EQ.2 ) THEN
          climate%biome(k) = 8
          climate%iveg(k) = 1
       ENDIF

       ! (9) Cold mixed forest
       IF(pft_biome1(k,1)==5.AND.pft_biome1(k,2)==7 .AND. npft .EQ.2 ) THEN
          climate%biome(k) = 9
          climate%iveg(k) = 1
       ENDIF

       ! (10) Cold deciduous forest
       IF(pft_biome1(k,1)==7 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 10
          climate%iveg(k) = 3
       ENDIF

       ! (11) Xerophytic woods/scrub
       IF(pft_biome1(k,1)==8 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 11
          climate%iveg(k) = 5
       ENDIF

       ! (12) Warm grass/shrub
       IF(pft_biome1(k,1)==9 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 12
          climate%iveg(k) = 5  ! include C4 grass tile ?
       ENDIF

       ! (13) Cool grass/shrub
       IF(pft_biome1(k,1)==10 .AND.pft_biome1(k,2)==11 .AND.  npft .EQ.2 ) THEN
          climate%biome(k) = 13
          climate%iveg(k) = 5  ! include C3 grass tile ?
       ENDIF

       ! (14) Tundra
       IF(pft_biome1(k,1)==11 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 14
          climate%iveg(k) = 8  !
       ENDIF

       ! (15) Hot desert
       IF(pft_biome1(k,1)==12 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 15
          climate%iveg(k) = 14  !
       ENDIF

       ! (16) Semidesert
       IF(pft_biome1(k,1)==13 .AND. npft .EQ.1 ) THEN
          climate%biome(k) = 16
          climate%iveg(k) = 5  !
       ENDIF

       ! (17) Ice/polar desert
       IF(climate%biome(k)==999) THEN
          climate%biome(k) = 17
          climate%iveg(k) = 17  !
       ENDIF

       !jhan: see CABLE Ticket#149 and above CABLE_LSM: comment
# ifndef ESM15
       ! check for DBL or NEL in SH: set to EBL instead
       IF ((climate%iveg(k)==1 .OR.climate%iveg(k)==3 .OR. climate%iveg(k)==4) &
            .AND. patch(k)%latitude<0) THEN
          climate%iveg(k) = 2
       ENDIF
# endif

       !"(/grass:1/shrub:2/woody:3"
       !1,3,Evergreen Needleleaf Forest
       !2,3,Evergreen Broadleaf Forest,,,,,,,,,,,,,,,,,,
       !3,3,Deciduous Needleleaf Forest,,,,,,,,,,,,,,,,,,
       !4,3,Deciduous Broadleaf Forest,,,,,,,,,,,,,,,,,,
       !5,2,shrub,,,,,,,,,,,,,,,,,,
       !6,1,C3 grass,,,,,,,,,,,,,,,,,,
       !7,1,C4 grass,,,,,,,,,,,,,,,,,,
       !8,1,tundra,,,,,,,,,
       ! 1.



    END DO


  END SUBROUTINE BIOME1_PFT

  ! ==============================================================================

  SUBROUTINE climate_init ( climate,np, ktauday )
    IMPLICIT NONE

    TYPE (climate_type), INTENT(INOUT)       :: climate  ! climate variables
    INTEGER, INTENT(IN) :: np, ktauday
    INTEGER :: d

    CALL alloc_cbm_var(climate,np)

    IF (cable_user%climate_fromzero .OR. .NOT.cable_user%call_climate) THEN

       ! Maciej
       !   DO d=1,31
       !climate%dtemp_31(:,d)= climate%dtemp
       ! climate%dmoist_31(:,d)= climate%dmoist
       !      climate%dtemp_31(:,d)= 0
       !      climate%dmoist_31(:,d)= 0
       !
       !   ENDDO
       climate%dtemp_31(:,:)= 0
       climate%dtemp_91(:,:)= 0
       climate%dmoist_31(:,:)= 0
       climate%atemp_mean=0



       climate%nyears = 0
       climate%chilldays = 0
       climate%mtemp = 0
       climate%mmoist = 0
       climate%mtemp_min = 0
       climate%mtemp_max=0
       climate%qtemp_max=0
       climate%qtemp_max_last_year = 0
       climate%mtemp_min20 =0
       climate%mtemp_max20=0
       climate%AGDD5=0
       climate%GDD5=0
       climate%AGDD0=0
       climate%GDD0=0
       climate%alpha_PT=0
       climate%evap_PT=0
       climate%aevap=0
       climate%mtemp_min_20=0
       climate%mtemp_max_20=0
       climate%alpha_PT_20=0
       climate%iveg = 999
       climate%biome = 999


    ELSE
       CALL READ_CLIMATE_RESTART_NC (climate, ktauday)

    ENDIF
    !else
    ! CALL READ_CLIMATE_RESTART_NC (climate)

    !endif

  END SUBROUTINE climate_init

  ! ==============================================================================

  SUBROUTINE WRITE_CLIMATE_RESTART_NC ( climate, ktauday )

    USE netcdf


    IMPLICIT NONE

    TYPE (climate_type), INTENT(IN)       :: climate  ! climate variables
    INTEGER, INTENT(IN) :: ktauday
    INTEGER*4 :: mp4
    INTEGER*4, PARAMETER   :: pmp4 =0
    INTEGER, PARAMETER   :: fmp4 = KIND(pmp4)
    INTEGER*4   :: STATUS
    INTEGER*4   :: FILE_ID, land_ID, nyear_ID, nday_ID, ndayq_ID, i
    CHARACTER :: CYEAR*4, FNAME*99,dum*50

    ! 0 dim arrays
    CHARACTER(len=20),DIMENSION(2) :: A0
    ! 1 dim arrays (npt )
    CHARACTER(len=20),DIMENSION(20) :: A1
    ! 1 dim arrays (integer) (npt )
    CHARACTER(len=20),DIMENSION(3) :: AI1
    ! 2 dim arrays (npt,20)
    CHARACTER(len=20),DIMENSION(3) :: A2
    ! 2 dim arrays (npt,31)
    CHARACTER(len=20),DIMENSION(2) :: A3
    ! 2 dim arrays (npt,91)
    CHARACTER(len=20),DIMENSION(1) :: A4


    INTEGER*4 ::  VID0(SIZE(A0)),VID1(SIZE(A1)),VIDI1(SIZE(AI1)), &
         VID2(SIZE(A2)), VID3(SIZE(A3)), VID4(SIZE(A4))

    mp4=INT(mp,fmp4)
    A0(1) = 'nyears'
    A0(2) = 'year'

    A1(1) = 'latitude'
    A1(2) = 'longitude'
    A1(3) = 'dtemp'
    A1(4) = 'mtemp'
    A1(5) = 'qtemp'
    A1(6) = 'mtemp_min'
    A1(7) = 'mtemp_max'
    A1(8) = 'qtemp_max'
    A1(9) = 'qtemp_max_last_year'
    A1(10) = 'mtemp_min20'
    A1(11) = 'mtemp_max20'
    A1(12) = 'atemp_mean'
    A1(13) = 'AGDD5'
    A1(14) = 'GDD5'
    A1(15) = 'AGDD0'
    A1(16) = 'GDD0'
    A1(17) = 'mmoist'
    A1(18) = 'evap_PT'
    A1(19) = 'aevap'
    A1(20)  = 'alpha_PT'

    AI1(1) = 'chilldays'
    AI1(2) = 'iveg'
    AI1(3) = 'biome'

    A2(1) = 'mtemp_min_20'
    A2(2) = 'mtemp_max_20'
    A2(3) = 'alpha_PT_20'

    A3(1) = 'dtemp_31'
    A3(2) = 'dmoist_31'

    A4(1) = 'dtemp_91'

# ifndef ESM15

    ! Get File-Name
    WRITE(CYEAR, FMT='(I4)') CurYear + 1
    fname = TRIM(filename%path)//'/'//TRIM( cable_user%RunIden )//&
         '_climate_rst.nc'
    ! Create NetCDF file:
    STATUS = NF90_create(fname, NF90_CLOBBER, FILE_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

    ! Put the file in define mode:
    STATUS = NF90_redef(FILE_ID)

    STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "Valid restart date", "01/01/"//CYEAR  )

    ! Define dimensions:
    ! Land (number of points)
    STATUS = NF90_def_dim(FILE_ID, 'land'   , mp4     , land_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    ! number of years (stored for 20 y running means0
    STATUS = NF90_def_dim(FILE_ID, 'nyear' , 20 , nyear_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    ! number of days (stored for 31 day monthly means)
    STATUS = NF90_def_dim(FILE_ID, 'nday' , 31 , nday_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    ! number of days (stored for 91 day quarterly means)
    STATUS = NF90_def_dim(FILE_ID, 'ndayq' , 91 , ndayq_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

    DO i = 1, SIZE(A0)
       STATUS = NF90_def_var(FILE_ID,TRIM(A0(i)) ,NF90_INT ,VID0(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO


    DO i = 1, SIZE(A1)
       STATUS = NF90_def_var(FILE_ID,TRIM(A1(i)) ,NF90_FLOAT,(/land_ID/),VID1(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO

    DO i = 1, SIZE(AI1)
       STATUS = NF90_def_var(FILE_ID,TRIM(AI1(i)) ,NF90_INT,(/land_ID/),VIDI1(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO

    DO i = 1, SIZE(A2)
       STATUS = NF90_def_var(FILE_ID,TRIM(A2(i)) ,NF90_FLOAT,(/land_ID,nyear_ID/),VID2(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO

    DO i = 1, SIZE(A3)
       STATUS = NF90_def_var(FILE_ID,TRIM(A3(i)) ,NF90_FLOAT,(/land_ID,nday_ID/),VID3(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO

    DO i = 1, SIZE(A4)
       STATUS = NF90_def_var(FILE_ID,TRIM(A4(i)) ,NF90_FLOAT,(/land_ID,ndayq_ID/),VID4(i))
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    END DO

    ! End define mode:
    STATUS = NF90_enddef(FILE_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

    ! PUT nyears and current year
    STATUS = NF90_PUT_VAR(FILE_ID, VID0(1), climate%nyears )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID0(2), CurYear )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    ! PUT LAT / LON
    STATUS = NF90_PUT_VAR(FILE_ID, VID1(1), patch%latitude )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(2), patch%longitude )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    ! PUT VARS
    STATUS = NF90_PUT_VAR(FILE_ID, VID1(3), climate%dtemp )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(4), climate%mtemp )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(5), climate%qtemp )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(6), climate%mtemp_min )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(7), climate%mtemp_max )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(8), climate%qtemp_max )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(9), climate%qtemp_max_last_year )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(10), climate%mtemp_min20 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(11), climate%mtemp_max20  )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(12), climate%atemp_mean  )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(13), climate%AGDD5  )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(14), climate%GDD5  )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(15), climate%AGDD0 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(16), climate%GDD0 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(17), climate%mmoist )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(18), climate%evap_PT )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(19), climate%aevap )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID1(20), climate%alpha_PT )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)


    STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(1), climate%chilldays )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(2), climate%iveg )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(3), climate%biome )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID2(1), climate%mtemp_min_20 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID2(2), climate%mtemp_max_20 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID2(3), climate%alpha_PT_20 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID3(1), climate%dtemp_31 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID3(2), climate%dmoist_31 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    STATUS = NF90_PUT_VAR(FILE_ID, VID4(1), climate%dtemp_91 )
    IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)

    ! Close NetCDF file:
    STATUS = NF90_close(FILE_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

# endif

  END SUBROUTINE WRITE_CLIMATE_RESTART_NC
  ! ==============================================================================

  SUBROUTINE READ_CLIMATE_RESTART_NC ( climate, ktauday )

    USE netcdf


    IMPLICIT NONE

    TYPE (climate_type), INTENT(INOUT)       :: climate  ! climate variables
    INTEGER, INTENT(IN) :: ktauday
    INTEGER*4 :: mp4
    INTEGER*4, PARAMETER   :: pmp4 =0
    INTEGER, PARAMETER   :: fmp4 = KIND(pmp4)
    INTEGER*4   :: STATUS
    INTEGER*4   :: FILE_ID, land_ID, nyear_ID, nday_ID, dID, i, land_dim
    CHARACTER :: CYEAR*4, FNAME*99,dum*50

    ! 0 dim arrays
    CHARACTER(len=20),DIMENSION(2) :: A0
    ! 1 dim arrays (npt )
    CHARACTER(len=20),DIMENSION(20) :: A1
    ! 1 dim arrays (integer) (npt )
    CHARACTER(len=20),DIMENSION(3) :: AI1
    ! 2 dim arrays (npt,20)
    CHARACTER(len=20),DIMENSION(3) :: A2
    ! 2 dim arrays (npt,31)
    CHARACTER(len=20),DIMENSION(2) :: A3
    ! 2 dim arrays (npt,91)
    CHARACTER(len=20),DIMENSION(1) :: A4

    REAL(r_2), DIMENSION(mp)          :: LAT, LON, TMP
    REAL(r_2)                         :: TMP2(mp,20),TMP3(mp,31),TMP4(mp,91)
    INTEGER*4 :: TMPI(mp), TMPI0
    LOGICAL            ::  EXISTFILE

    mp4=INT(mp,fmp4)
    A0(1) = 'nyears'
    A0(2) = 'year'

    A1(1) = 'latitude'
    A1(2) = 'longitude'
    A1(3) = 'dtemp'
    A1(4) = 'mtemp'
    A1(5) = 'qtemp'
    A1(6) = 'mtemp_min'
    A1(7) = 'mtemp_max'
    A1(8) = 'qtemp_max'
    A1(9) = 'qtemp_max_last_year'
    A1(10) = 'mtemp_min20'
    A1(11) = 'mtemp_max20'
    A1(12) = 'atemp_mean'
    A1(13) = 'AGDD5'
    A1(14) = 'GDD5'
    A1(15) = 'AGDD0'
    A1(16) = 'GDD0'
    A1(17) = 'mmoist'
    A1(18) = 'evap_PT'
    A1(19) = 'aevap'
    A1(20)  = 'alpha_PT'

    AI1(1) = 'chilldays'
    AI1(2) = 'iveg'
    AI1(3) = 'biome'

    A2(1) = 'mtemp_min_20'
    A2(2) = 'mtemp_max_20'
    A2(3) = 'alpha_PT_20'

    A3(1) = 'dtemp_31'
    A3(2) = 'dmoist_31'

    A4(1) = 'dtemp_91'


    ! Get File-Name
    WRITE(CYEAR, FMT='(I4)') CurYear + 1
    fname = TRIM(filename%path)//'/'//TRIM( cable_user%RunIden )//&
         '_climate_rst.nc'

    INQUIRE( FILE=TRIM( fname ), EXIST=EXISTFILE )

    IF ( .NOT.EXISTFILE) WRITE(*,*) fname, ' does not exist!'

# ifndef ESM15
    ! Open NetCDF file:
    STATUS = NF90_OPEN(fname, NF90_NOWRITE, FILE_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)


    ! dimensions:
    ! Land (number of points)
    STATUS = NF90_INQ_DIMID(FILE_ID, 'land'   , dID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    STATUS = NF90_INQUIRE_DIMENSION( FILE_ID, dID, LEN=land_dim )
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    ! number of years (stored for 20 y running means
    STATUS = NF90_INQ_DIMID(FILE_ID, 'nyear'   , dID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    ! number of days (stored for 31 d monthly means
    STATUS = NF90_INQ_DIMID(FILE_ID, 'nday'   , dID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

    IF ( land_dim .NE. SIZE(patch%latitude)) THEN
       WRITE(*,*) "Dimension misfit, ", fname
       WRITE(*,*) "land_dim", land_dim, SIZE(patch%latitude)
       !     STOP
    ENDIF

    ! LAT & LON
    STATUS = NF90_INQ_VARID( FILE_ID, A1(1), dID )
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    STATUS = NF90_GET_VAR( FILE_ID, dID, LAT )
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)


    STATUS = NF90_INQ_VARID( FILE_ID, A1(2), dID )
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
    STATUS = NF90_GET_VAR( FILE_ID, dID, LON )
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

    ! READ scalar fields
    DO i = 1, SIZE(A0)
       STATUS = NF90_INQ_VARID( FILE_ID, A0(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMPI0 )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(A0(i)))
       CASE ('nyears'      ) ; climate%nyears      = TMPI0
       END SELECT
    END DO


    ! READ 1-dimensional real fields
    DO i = 3, SIZE(A1)
       STATUS = NF90_INQ_VARID( FILE_ID, A1(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMP )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(A1(i)))
       CASE ('mtemp'      ) ; climate%mtemp       = TMP
       CASE ('qtemp'      ) ; climate%mtemp       = TMP
       CASE ('mtemp_min'   ) ; climate%mtemp_min   = TMP
       CASE ('mtemp_max'  ) ; climate%mtemp_max  = TMP
       CASE ('qtemp_max'  ) ; climate%qtemp_max  = TMP
       CASE ('qtemp_max_last_year'  ) ; climate%qtemp_max_last_year  = TMP
       CASE ('mtemp_min20' ) ; climate%mtemp_min20 = TMP
       CASE ('mtemp_max20'  ) ; climate%mtemp_max20  = TMP
       CASE ('atemp_mean'  ) ; climate%atemp_mean  = TMP
       CASE ('AGDD5'  ) ; climate%AGDD5  = TMP
       CASE ('GDD5'  ) ; climate%GDD5  = TMP
       CASE ('AGDD0'  ) ; climate%AGDD0  = TMP
       CASE ('GDD0'  ) ; climate%GDD0  = TMP
       CASE ('evap_PT'  ) ; climate%evap_PT  = TMP
       CASE ('aevap'  ) ; climate%aevap  = TMP
       CASE ('alpha_PT'  ) ; climate%alpha_PT  = TMP
       END SELECT
    END DO

    ! READ 1-dimensional integer fields
    DO i = 1, SIZE(AI1)

       WRITE(*,*)  TRIM(AI1(i))
       STATUS = NF90_INQ_VARID( FILE_ID, AI1(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMPI )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(AI1(i)))
       CASE ('chilldays'      ) ; climate%chilldays      = TMPI
       CASE ('iveg'      ) ; climate%iveg     = TMPI
       CASE ('biome'      ) ; climate%biome     = TMPI
       END SELECT
    END DO


    ! READ 2-dimensional fields (nyear)
    DO i = 1, SIZE(A2)
       WRITE(*,*)  TRIM(A2(i))
       STATUS = NF90_INQ_VARID( FILE_ID, A2(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMP2 )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(A2(i)))
       CASE ('mtemp_min_20' ) ; climate%mtemp_min_20 = TMP2
       CASE ('mtemp_max_20' ) ; climate%mtemp_max_20 = TMP2
       CASE ('alpha_PT_20' ) ; climate%alpha_PT_20 = TMP2
       END SELECT
    END DO

    ! READ 2-dimensional fields (nday)
    DO i = 1, SIZE(A3)
       STATUS = NF90_INQ_VARID( FILE_ID, A3(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMP3 )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(A3(i)))
       CASE ('dtemp_31' ) ; climate%dtemp_31 = TMP3
       CASE ('dmoist_31' ) ; climate%dmoist_31 = TMP3
       END SELECT
    END DO

    ! READ 2-dimensional fields (ndayq)
    DO i = 1, SIZE(A4)
       STATUS = NF90_INQ_VARID( FILE_ID, A4(i), dID )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
       STATUS = NF90_GET_VAR( FILE_ID, dID, TMP4 )
       IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)

       SELECT CASE ( TRIM(A4(i)))
       CASE ('dtemp_91' ) ; climate%dtemp_91 = TMP4
       END SELECT
    END DO


    ! Close NetCDF file:
    STATUS = NF90_close(FILE_ID)
    IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
# endif


  END SUBROUTINE  READ_CLIMATE_RESTART_NC


END MODULE cable_climate_mod
# endif
