      program Noah_driver

      use module_Noahlsm_utility
      use module_sf_noahlsm_param_init
      use module_date_utilities
      use module_sf_Noahlsm_gridded_input
      use module_sf_noahlsm
      use module_sf_Noah_NC_output
      use module_rootfraction

      implicit none

      character(len=6)                     :: mdt         !model start month, day, and hour
      character(len=2)                     :: minute      !model start minute
      character(len=4)                     :: year        !model start year

      integer                              :: imonth      !only used to compute cosz
      integer                              :: iday        !only used to compute cosz
      integer                              :: itime       !only used to compute cosz

      integer                              :: nst         !number of time step during a day
      integer                              :: istep       !counting records in forcing data
      integer                              :: imstep      !for monthly output
      integer                              :: ND          !for monthly output
      integer                              :: ntimed      !for daily output
      integer                              :: idstep      !for 3-hourly output
      integer                              :: startstep   !starting time step to read forcings
      integer                              :: iyloop      !year loop index
      integer                              :: imloop      !month loop index
      integer                              :: idloop      !day loop index
      integer                              :: ihloop      !hour loop index
      integer                              :: startday    !starting day
      integer                              :: startmonth  !strating month
      integer                              :: ierr        !error message index
      integer, dimension(0:11)             :: jday        !julian day
      integer, dimension(0:11)             :: gday        !gorgian day
      integer, dimension(1:12)             :: nday        !days per month

      data (nday(imonth),imonth=1,12) /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
      data (jday(imonth),imonth=0,11) / 0, 31, 59, 90,120,151,181,212,242,273,304,334/
      data (gday(imonth),imonth=0,11) / 0, 31, 60, 91,121,152,182,213,243,274,305,335/

      integer                              :: ix          !do loop index in x-direction
      integer                              :: iy          !do loop index in y-direction
      integer                              :: iz          !do loop index in z-direction
      integer                              :: it          !do loop index for time
      integer                              :: nsoil       !number of soil layers
      integer                              :: nx          !total x-direction points
      integer                              :: ny          !total y-direction points
      integer, parameter                   :: nsnow = 3   !maximum no. of snow layers
      real, dimension(1:6)                 :: zsoil       !depth of the bottom of soil layers(m)

      integer                              :: ispin       !index for sinp-up
      integer                              :: nspin       !total spin-up times

!grided variables

      integer, allocatable, dimension(:,:) :: vegtyp    !vegetation type
      integer, allocatable, dimension(:,:) :: soiltyp   !soil type
      integer, allocatable, dimension(:,:) :: sloptyp   !sloptyp (only used for Noah runoff scheme)
      integer, allocatable, dimension(:,:) :: landsea   !land-sea mask (0-ocean; 1-land)
      integer, allocatable, dimension(:,:) :: ice       !land-sea mask (0-ocean; 1-land)
      integer, allocatable, dimension(:,:) :: isc       !soil color index
      integer, allocatable, dimension(:,:) :: ist       !surface type: 1->soil; 2-> lake
      real, allocatable, dimension(:,:)    :: lat       !latitude
      real, allocatable, dimension(:,:)    :: cosz      !solar zenith angle
      real, allocatable, dimension(:,:)    :: lon       !longitude
      real, allocatable, dimension(:,:)    :: tbot      !soil temperature at zbot
      real, allocatable, dimension(:,:,:)  :: shdfac    !greenness vegetation (shaded) fraction

      real, allocatable, dimension(:,:)    :: sfctmp    !surface air temperature
      real, allocatable, dimension(:,:)    :: q2        !surface air specific humidity
      real, allocatable, dimension(:,:)    :: lwdn      !downward longwave radiation
      real, allocatable, dimension(:,:)    :: uu        !wind speed at x-direction
      real, allocatable, dimension(:,:)    :: vv        !wind speed at y-direction
      real, allocatable, dimension(:,:)    :: sfcprs    !surfae air pressure
      real, allocatable, dimension(:,:)    :: soldn     !downward shortwave radiation
      real, allocatable, dimension(:,:)    :: prcp      !precipitation
      real, allocatable, dimension(:,:)    :: o2air     !
      real, allocatable, dimension(:,:)    :: co2air    !

      integer, allocatable, dimension(:,:) :: isnow     !actual no. of snow layers
      real, allocatable, dimension(:,:,:)  :: stc       !snow/soil tmperatures
      real, allocatable, dimension(:,:,:)  :: smc       !vol. soil moisture (m3/m3)
      real, allocatable, dimension(:,:,:)  :: sh2o      !vol. soil liquid water (m3/m3)
      real, allocatable, dimension(:,:,:)  :: sice      !vol. soil ice (m3/m3)
      real, allocatable, dimension(:,:)    :: tv        !vegetation canopy temperature
      real, allocatable, dimension(:,:)    :: tg        !ground surface temperature
      real, allocatable, dimension(:,:)    :: canice    !canopy-intercepted ice (mm)
      real, allocatable, dimension(:,:)    :: canliq    !canopy-intercepted liquid water (mm)
      real, allocatable, dimension(:,:)    :: snowh     !snow depth (m)
      real, allocatable, dimension(:,:)    :: sneqv     !snow water equivalent (mm)
      real, allocatable, dimension(:,:)    :: eah       !canopy air vapor pressure (pa)
      real, allocatable, dimension(:,:)    :: tah       !canopy air temperature (k)
      real, allocatable, dimension(:,:)    :: cm        !momentum drag coefficient
      real, allocatable, dimension(:,:)    :: ch        !sensible heat exchange coefficient

      integer, allocatable, dimension(:,:) :: nroot     !number of root layers
      real, allocatable, dimension(:,:,:)  :: qroot     !root water uptake of each layer (m/s)
      REAL, ALLOCATABLE, DIMENSION(:,:,:)  :: sadr      !root surface area density [m2/m3]
      REAL, ALLOCATABLE, DIMENSION(:,:)    :: mq        !liquid water in living plant tissues [kg/m2]

      real, allocatable, dimension(:,:)    :: fwet      !wetted or snowed fraction of the canopy (-)
      real, allocatable, dimension(:,:)    :: sneqvo    !snow mass at last time step(mm h2o)
      real, allocatable, dimension(:,:)    :: albold    !snow albedo at last time step (-)
      real, allocatable, dimension(:,:)    :: qsnow     !snowfall on the ground [mm/s]
      real, allocatable, dimension(:,:)    :: wslake    !lake water storage [mm]
      real, allocatable, dimension(:,:,:)  :: ficeold   !snow layer ice fraction at previous step [mm]
  
      real, allocatable, dimension(:,:)    :: zwt       !water table depth [m]
      real, allocatable, dimension(:,:)    :: wa        !water in the "aquifer" [mm]
      real, allocatable, dimension(:,:)    :: wt        !groundwater storage [mm]
      real, allocatable, dimension(:,:,:)  :: tsno      !snow temperature [K]
      real, allocatable, dimension(:,:,:)  :: zsnso     !snow layer depth [m]
      real, allocatable, dimension(:,:,:)  :: snice     !snow layer ice [mm]
      real, allocatable, dimension(:,:,:)  :: snliq     !snow layer liquid water [mm]
      REAL, allocatable, dimension(:,:)    :: drymass   !total dry living carbon  mass [g/m2]
      REAL, allocatable, dimension(:,:)    :: lfmass    !leaf mass [g/m2]
      REAL, allocatable, dimension(:,:,:)  :: rootms    !carbon mass of fine roots [g/m2]
      REAL, allocatable, dimension(:,:)    :: stmass    !stem mass [g/m2]
      REAL, allocatable, dimension(:,:)    :: wood      !mass of wood (incl. woody roots) [g/m2]
      REAL, allocatable, dimension(:,:)    :: soc       !soil organic carbon          [g C/m2]
      REAL, allocatable, dimension(:,:)    :: wdoc      !wet dissolved organic carbon [g C/m2]
      REAL, allocatable, dimension(:,:)    :: ddoc      !dry dissolved organic carbon [g C/m2]
      REAL, allocatable, dimension(:,:)    :: mic       !soil microbial biomass       [g C/m2]
      REAL, allocatable, dimension(:,:)    :: wenz      !wet soil enzyme              [g C/m2]
      REAL, allocatable, dimension(:,:)    :: denz      !dry soil enzyme              [g C/m2]
      REAL, allocatable, dimension(:,:)    :: soilc     !total soil carbon            [g C/m2]
      REAL, allocatable, dimension(:,:)    :: so2       ! o2 in soil pore in layers   [g O2/m2]
      REAL, allocatable, dimension(:,:)    :: kr        !
      REAL, allocatable, dimension(:,:,:)  :: froot     !

      REAL, allocatable, dimension(:,:)    :: xlai      !leaf area index
      REAL, allocatable, dimension(:,:)    :: xsai      !stem area index
      REAL, allocatable, dimension(:,:)    :: trad      !surface radiative temperature (k)
      REAL, allocatable, dimension(:,:)    :: t2m       !2 m air temperature (k)
      REAL, allocatable, dimension(:,:)    :: nee       !net ecosys exchange (g/m2/s C)
      REAL, allocatable, dimension(:,:)    :: gpp       !gross primary assimilation [g/m2/s C]
      REAL, allocatable, dimension(:,:)    :: npp       !net primary productivity [g/m2/s C]
      REAL, allocatable, dimension(:,:)    :: qco2      !soil surface co2 efflux (mg CO2/m2/s)
      REAL, allocatable, dimension(:,:)    :: fveg      !greenness vegetation fraction [-]
      REAL, allocatable, dimension(:,:)    :: ndvi      !ndvi
      real, allocatable, dimension(:,:)    :: vmax      !max SOC decomposition rate per unit
                                                        !microbial biomass [g SOC/m2/[g MIC/m2]/s]
      real, allocatable, dimension(:,:)    :: vmaxup    !max DOC uptake rate [g DOC/m2/[g CMIC/m2]/s] 
      real, allocatable, dimension(:,:)    :: km        !Michaelis-Menten constant [g/m2] SOC decomp
      real, allocatable, dimension(:,:)    :: kmup      !Michaelis-Menten constant [g/m2] DOC uptake
      real, allocatable, dimension(:,:)    :: epslon    !carbon use efficiency

      real, allocatable, dimension(:,:)    :: qin       !groundwater recharge [mm/s]
      real, allocatable, dimension(:,:)    :: runsf     !surface runoff [mm/s]
      real, allocatable, dimension(:,:)    :: runsb     !subsurface runoff [mm/s]
      real, allocatable, dimension(:,:)    :: ecan      !evaporation of intercepted water (mm/s)
      real, allocatable, dimension(:,:)    :: edir      !soil surface evaporation rate (mm/s]
      real, allocatable, dimension(:,:)    :: etran     !transpiration rate (mm/s)

      real, allocatable, dimension(:,:)    :: fsa       !total absorbed solar radiation (w/m2)
      real, allocatable, dimension(:,:)    :: fsr       !total reflected solar radiation (w/m2)
      real, allocatable, dimension(:,:)    :: fira      !total net longwave rad (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fsh       !total sensible heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fgev      !ground latent heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fctr      !transpiration latent heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fcev      !interception latent heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: flh       !total latent heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fgh       !ground heat flux (w/m2) [+ to soil]
      real, allocatable, dimension(:,:)    :: apar      !photosyn active energy by canopy (w/m2)
      real, allocatable, dimension(:,:)    :: psn       !total photosynthesis (umol co2/m2/s) [+]
      real, allocatable, dimension(:,:)    :: sav       !solar rad absorbed by veg. (w/m2)
      real, allocatable, dimension(:,:)    :: sag       !solar rad absorbed by ground (w/m2)
      real, allocatable, dimension(:,:)    :: fsnog     !snow cover fraction on the ground (-)
      real, allocatable, dimension(:,:)    :: fsnoc     !snow cover fraction of the canopy (-)
      real, allocatable, dimension(:,:)    :: scf       !snow cover fraction of the grid cell (-)
      real, allocatable, dimension(:,:,:)  :: psi       !suction head in the soil matrix [m]
      real, allocatable, dimension(:,:,:)  :: wcnd      !unsaturated hydraulic conductivity [m/s]
      real, allocatable, dimension(:,:)    :: swdown    !downward solar radiation (- value removed)
      real, allocatable, dimension(:,:)    :: qsubcan   !sublimation (+ and -) from the canopy snow [mm/s]
      real, allocatable, dimension(:,:)    :: qsubgrd   !sublimation (+ and -) from the ground snow [mm/s]

      real, allocatable, dimension(:,:)    :: bexp
      real, allocatable, dimension(:,:,:)    :: dksat
      real, allocatable, dimension(:,:)    :: dwsat
      real, allocatable, dimension(:,:)    :: f1
      real, allocatable, dimension(:,:)    :: psisat
      real, allocatable, dimension(:,:)    :: quartz 
      real, allocatable, dimension(:,:)    :: smcdry
      real, allocatable, dimension(:,:)    :: smcmax
      real, allocatable, dimension(:,:)    :: smcref
      real, allocatable, dimension(:,:)    :: smcwlt
      real, allocatable, dimension(:,:)    :: hs   
      real, allocatable, dimension(:,:)    :: rgl  
      real, allocatable, dimension(:,:)    :: rsmin
      real, allocatable, dimension(:,:)    :: slope
      real, allocatable, dimension(:,:)    :: kdt
      real, allocatable, dimension(:,:)    :: frzx

! daily mean variables

      real, allocatable, dimension(:,:)    :: snowhdm      !snow depth (m)
      real, allocatable, dimension(:,:)    :: sneqvdm      !snow water equivelant (mm or kg/m2)
      real, allocatable, dimension(:,:)    :: scfdm        !snow cover fraction (-)
      real, allocatable, dimension(:,:)    :: fsnogdm      !snow cover fraction on the ground (-)
      real, allocatable, dimension(:,:)    :: runsfdm      !surface runoff (mm/s)
      real, allocatable, dimension(:,:)    :: runsbdm      !subsurface runoff (mm/s)
      real, allocatable, dimension(:,:)    :: qsubcandm    !sublimation (+ and -) from the canopy snow [mm/s]
      real, allocatable, dimension(:,:)    :: qsubgrddm    !sublimation (+ and -) from the ground snow [mm/s]

! monthly mean variables

      real, allocatable, dimension(:,:)    :: snowhm      !snow depth (m)
      real, allocatable, dimension(:,:)    :: sneqvm      !snow water equivelant (mm or kg/m2)
      real, allocatable, dimension(:,:)    :: scfm        !snow cover fraction (-)
      real, allocatable, dimension(:,:)    :: fsnogm      !snow cover fraction on the ground (-)
      real, allocatable, dimension(:,:)    :: tgm         !ground surface temperature (K)
      real, allocatable, dimension(:,:)    :: prcpm       !precipitation (mm/s)
      real, allocatable, dimension(:,:)    :: qinm        !groundwater recharge (mm/s)
      real, allocatable, dimension(:,:)    :: runsfm      !surface runoff (mm/s)
      real, allocatable, dimension(:,:)    :: runsbm      !subsurface runoff (mm/s)
      real, allocatable, dimension(:,:)    :: ecanm       !evaporation of intercepted water (mm/s)
      real, allocatable, dimension(:,:)    :: edirm       !soil surface evaporation rate (mm/s]
      real, allocatable, dimension(:,:)    :: etranm      !transpiration rate (mm/s)
      real, allocatable, dimension(:,:)    :: zwtm        !the depth to water table (m)
      real, allocatable, dimension(:,:)    :: fsam        !total absorbed solar radiation (w/m2)
      real, allocatable, dimension(:,:)    :: firam       !total net longwave rad (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fshm        !total sensible heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: flhm        !total latent heat (w/m2) [+ to atm]
      real, allocatable, dimension(:,:)    :: fghm        !ground heat flux (w/m2) [+ to soil]
      real, allocatable, dimension(:,:,:)  :: stcm        !soil temperature (K)
      real, allocatable, dimension(:,:,:)  :: smcm        !volumetric soil moisture (m3/m3)
      real, allocatable, dimension(:,:,:)  :: sh2om       !volumetric liquid water content in soil (m3/m3)
      real, allocatable, dimension(:,:)    :: aparm       !photosyn active energy by canopy (w/m2)
      real, allocatable, dimension(:,:)    :: psnm        !total photosynthesis (umol co2/m2/s) [+]
      real, allocatable, dimension(:,:)    :: savm        !solar rad absorbed by veg. (w/m2)
      real, allocatable, dimension(:,:)    :: sagm        !solar rad absorbed by ground (w/m2)
      real, allocatable, dimension(:,:)    :: xlaim       !leaf area index (-)
      real, allocatable, dimension(:,:)    :: xsaim       !stem area index (-)
      real, allocatable, dimension(:,:)    :: tradm       !surface radiative temperature (K)
      real, allocatable, dimension(:,:)    :: t2mm        !2 m air temperature (K)
      real, allocatable, dimension(:,:)    :: neem        !net ecosystem exchange (g/m2/s C)
      real, allocatable, dimension(:,:)    :: gppm        !gross primary productivity [g/m2/s C]
      real, allocatable, dimension(:,:)    :: nppm        !net primary productivity [g/m2/s C]
      real, allocatable, dimension(:,:)    :: fvegm       !greenness vegetation fraction [-]
      real, allocatable, dimension(:,:)    :: cmm         !surface exchange coeffcient for momentum [-]
      real, allocatable, dimension(:,:)    :: chm         !surface exchange coeffcient for heat [-]
      real, allocatable, dimension(:,:)    :: qco2m       !soil surface co2 efflux (mg CO2/m2/s)
      REAL, allocatable, dimension(:,:)    :: socm        !soil organic carbon          [g C/m3]
      REAL, allocatable, dimension(:,:)    :: wdocm       !wet dissolved organic carbon [g C/m3]
      REAL, allocatable, dimension(:,:)    :: ddocm       !dry dissolved organic carbon [g C/m3]
      REAL, allocatable, dimension(:,:)    :: micm        !soil microbial biomass       [g C/m3]
      REAL, allocatable, dimension(:,:)    :: wenzm       !wet soil enzyme              [g C/m3]
      REAL, allocatable, dimension(:,:)    :: denzm       !dry soil enzyme              [g C/m3]
      real, allocatable, dimension(:,:)    :: vmaxm   !max SOC decomposition rate per unit
                                                      !microbial biomass [g SOC/m2/[g MIC/m2]/s]
      real, allocatable, dimension(:,:)    :: vmaxupm !max DOC uptake rate [g DOC/m2/[g CMIC/m2]/s] 
      real, allocatable, dimension(:,:)    :: kmm     !Michaelis-Menten constant [g/m2] SOC decomp
      real, allocatable, dimension(:,:)    :: kmupm   !Michaelis-Menten constant [g/m2] DOC uptake
      real, allocatable, dimension(:,:)    :: epslonm !carbon use efficiency
      real, allocatable, dimension(:,:)    :: wam     !groundwater storage (mm)
      real, allocatable, dimension(:,:)    :: canwatm !water on the canopy (mm)
      real, allocatable, dimension(:,:)    :: ndvim   !NDVI (-)
      real, allocatable, dimension(:,:)    :: mqm !
      real, allocatable, dimension(:,:)    :: krm
      real, allocatable, dimension(:,:,:)  :: frootm
      real, allocatable, dimension(:,:)    :: rtmassm
      real, allocatable, dimension(:,:,:)  :: sadrm
      real, allocatable, dimension(:,:,:)  :: psim
      real, allocatable, dimension(:,:,:)  :: qrootm   !root water uptake [m/s]
      real, allocatable, dimension(:,:,:)  :: qrootnm  !negative part of qroot:
                                                       !hydraulic lift or descent [m/s]
      real, allocatable, dimension(:,:)    :: sfctmpm  !2-m air temperature
      real, allocatable, dimension(:,:)    :: q2m      !2-m specific humidity
      real, allocatable, dimension(:,:)    :: windm    !surface wind
      real, allocatable, dimension(:,:)    :: sfcprsm  !surface air pressure

      integer, allocatable, dimension(:,:)    :: nndvi   !

      real :: soilcgm,socgm,wdocgm,ddocgm,micgm,wenzgm,denzgm,ndvigm
      integer :: nn,nndvig

!---------------------------------------------------------------------
! 1-d variables 
!---------------------------------------------------------------------
      real                            :: dt      !timestep (second) 
      real                            :: zlvl    !height (m) above ground of atmos. forcing
      real                            :: z0      !roughness length (m)
      real                            :: foln    !nitrogen saturation

      real, allocatable, dimension(:) :: sldpth  !the thickness (m) of each soil layer 
      real                            :: calday  !calendar day for computing solar zenith angle
!---------------------------------------------------------------------
!  declare/initialize constants
!---------------------------------------------------------------------
      character(len=14)  :: olddate
      character(len=14)  :: newdate
      character(len=256) :: DIR      !directory where input and output data files are
      character(len=256) :: fini     !directory where initial data files are
      character(len=8)   :: exp      !directory where an experiment is

      integer :: start_year,start_month,start_day,start_hour,start_min
      integer :: end_year
      namelist /noahlsm_offline/DIR,fini     ,exp    ,nsoil   ,nx      ,ny     ,& 
      zsoil  ,start_year,start_month,start_day,start_hour,start_min ,& 
      end_year  ,dt     ,zlvl

      open(30,file="noah_offline.namelist",status='old',form="formatted")
      read(30,noahlsm_offline,iostat=ierr)
      if (ierr /= 0) then
         stop
      end if
      close(30)
  
      write(mdt,'(i2.2,i2.2,i2.2)') start_month, start_day, start_hour
      write(minute,'(i2.2)') start_min
      print*, 'mdt, minute = ', mdt, minute
      write(year, '(i4.4)') start_year

      write(olddate,'(I4.4,I2.2,I2.2,I2.2,I2.2,I2.2)') & 
      start_year, start_month, start_day, start_hour, start_min, 0

      nst = int(86400./dt)
!----------------------------------------------------------------------
      call lsm_parm_init
!----------------------------------------------------------------------
! Allocate arrays for our gridded domain, now that we know the size
!----------------------------------------------------------------------
! land surface parameters

      allocate( lat    (1:nx,1:ny) )
      allocate( lon    (1:nx,1:ny) )
      allocate( cosz   (1:nx,1:ny) )
      allocate( vegtyp (1:nx,1:ny) )
      allocate( soiltyp(1:nx,1:ny) )
      allocate( sloptyp(1:nx,1:ny) )
      allocate( tbot   (1:nx,1:ny) )
      allocate( landsea(1:nx,1:ny) )
      allocate( ice    (1:nx,1:ny) )
      allocate( isc    (1:nx,1:ny) )
      allocate( ist    (1:nx,1:ny) )
      allocate( shdfac (1:nx,1:ny,1:12) )

      lat    (1:nx,1:ny)      = -999.0
      lon    (1:nx,1:ny)      = -999.0
      cosz   (1:nx,1:ny)      = -999.0
      vegtyp (1:nx,1:ny)      = -999
      soiltyp(1:nx,1:ny)      = -999
      sloptyp(1:nx,1:ny)      = -999
      landsea(1:nx,1:ny)      = -999
      ice    (1:nx,1:ny)      = -999
      isc    (1:nx,1:ny)      = -999
      ist    (1:nx,1:ny)      = -999
      tbot   (1:nx,1:ny)      = -999.0
      shdfac (1:nx,1:ny,1:12) = -999.0

! forcing data

      allocate( sfctmp (1:nx,1:ny) )
      allocate( q2     (1:nx,1:ny) )
      allocate( lwdn   (1:nx,1:ny) )
      allocate( uu     (1:nx,1:ny) )
      allocate( vv     (1:nx,1:ny) )
      allocate( sfcprs (1:nx,1:ny) )
      allocate( soldn  (1:nx,1:ny) )
      allocate( prcp   (1:nx,1:ny) )
      allocate( co2air (1:nx,1:ny) )
      allocate(  o2air (1:nx,1:ny) )

      sfctmp (1:nx,1:ny) =-999.9
      q2     (1:nx,1:ny) =-999.9
      lwdn   (1:nx,1:ny) =-999.9
      uu     (1:nx,1:ny) =-999.9
      vv     (1:nx,1:ny) =-999.9
      sfcprs (1:nx,1:ny) =-999.9
      soldn  (1:nx,1:ny) =-999.9
      prcp   (1:nx,1:ny) =-999.9
      co2air (1:nx,1:ny) =-999.9
      o2air  (1:nx,1:ny) =-999.9

! prognostic varaibles
! 1-d
      allocate( sldpth (       1:nsoil) )
      sldpth (       1:nsoil)  =-999.9

! 2d and 3d prognostic variables
      allocate( smc    (1:nx,1:ny,       1:nsoil) ) ! 1
      allocate( sh2o   (1:nx,1:ny,       1:nsoil) ) ! 3 
      allocate( sice   (1:nx,1:ny,       1:nsoil) ) ! 3
      allocate( stc    (1:nx,1:ny,-nsnow+1:nsoil) ) ! 2
      allocate( zsnso  (1:nx,1:ny,-nsnow+1:nsoil) ) ! 7
      allocate( tsno   (1:nx,1:ny,-nsnow+1:    0) ) ! 4
      allocate( snice  (1:nx,1:ny,-nsnow+1:    0) ) ! 5
      allocate( snliq  (1:nx,1:ny,-nsnow+1:    0) ) ! 6
      allocate( tv     (1:nx,1:ny) )                ! 8
      allocate( tg     (1:nx,1:ny) )                ! 9
      allocate( canliq (1:nx,1:ny) )                !10              
      allocate( canice (1:nx,1:ny) )                !11
      allocate( snowh  (1:nx,1:ny) )                !12
      allocate( sneqv  (1:nx,1:ny) )                !13
      allocate( zwt    (1:nx,1:ny) )                !14
      allocate( wa     (1:nx,1:ny) )                !15
      allocate( wt     (1:nx,1:ny) )                !16
      allocate( isnow  (1:nx,1:ny) )                !17
      allocate( lfmass (1:nx,1:ny) )                !18 
      allocate( rootms (1:nx,1:ny,1:nsoil) )        !19
      allocate( stmass (1:nx,1:ny) )                !20
      allocate( wood   (1:nx,1:ny) )                !21
      allocate( drymass(1:nx,1:ny))         !
      allocate( soc    (1:nx,1:ny))         !
      allocate( wdoc   (1:nx,1:ny))         !
      allocate( ddoc   (1:nx,1:ny))         !
      allocate( mic    (1:nx,1:ny))         !
      allocate( wenz   (1:nx,1:ny))         !  
      allocate( denz   (1:nx,1:ny))         !
      allocate( soilc  (1:nx,1:ny))                 !23
      allocate( so2    (1:nx,1:ny))                 !26
      allocate( xlai   (1:nx,1:ny) )                !24
      allocate( xsai   (1:nx,1:ny) )                !25
      allocate( eah    (1:nx,1:ny) )                !26
      allocate( tah    (1:nx,1:ny) )                !27

      allocate( fwet   (1:nx,1:ny) )                !28
      allocate( sneqvo (1:nx,1:ny) )                !29
      allocate( albold (1:nx,1:ny) )                !30
      allocate( qsnow  (1:nx,1:ny) )                !31
      allocate( wslake (1:nx,1:ny) )                !32
      allocate( ficeold(1:nx,1:ny,-nsnow+1:    0) ) ! 5

      allocate( cm     (1:nx,1:ny) )                !
      allocate( ch     (1:nx,1:ny) )                !

      allocate( kr     (1:nx,1:ny) )                !
      allocate( froot  (1:nx,1:ny,1:nsoil) )                !
      allocate( nroot  (1:nx,1:ny) )
      allocate( mq     (1:nx,1:ny) )                !
      allocate( qroot  (1:nx,1:ny,1:nsoil) )        !
      allocate( sadr   (1:nx,1:ny,1:nsoil) )        !

      smc    (1:nx,1:ny,       1:nsoil)  = -999.9 ! 1
      sh2o   (1:nx,1:ny,       1:nsoil)  = -999.9 ! 3 
      sice   (1:nx,1:ny,       1:nsoil)  = -999.9 ! 3 
      stc    (1:nx,1:ny,-nsnow+1:nsoil)  = -999.9 ! 2
      zsnso  (1:nx,1:ny,-nsnow+1:nsoil)  = -999.9 ! 7
      tsno   (1:nx,1:ny,-nsnow+1:    0)  = -999.9 ! 4
      snice  (1:nx,1:ny,-nsnow+1:    0)  = -999.9 ! 5
      snliq  (1:nx,1:ny,-nsnow+1:    0)  = -999.9 ! 6
      tv     (1:nx,1:ny)                 = -999.9 ! 8
      tg     (1:nx,1:ny)                 = -999.9 ! 9
      canliq (1:nx,1:ny)                 = -999.9 !10              
      canice (1:nx,1:ny)                 = -999.9 !11
      snowh  (1:nx,1:ny)                 = -999.9 !12
      sneqv  (1:nx,1:ny)                 = -999.9 !13
      zwt    (1:nx,1:ny)                 = -999.9 !14
      wa     (1:nx,1:ny)                 = -999.9 !15
      wt     (1:nx,1:ny)                 = -999.9 !16
      isnow  (1:nx,1:ny)                 = -999   !17
      lfmass (1:nx,1:ny)                 = -999.9 !18 
      rootms (1:nx,1:ny,1:nsoil)         = -999.9 !19
      stmass (1:nx,1:ny)                 = -999.9 !20
      wood   (1:nx,1:ny)                 = -999.9 !21
      drymass (1:nx,1:ny)                = -999.9
      soc    (1:nx,1:ny)                 = -999.9 !
      wdoc   (1:nx,1:ny)                 = -999.9 !
      ddoc   (1:nx,1:ny)                 = -999.9 !
      mic    (1:nx,1:ny)                 = -999.9 ! 
      wenz   (1:nx,1:ny)                 = -999.9 !
      denz   (1:nx,1:ny)                 = -999.9 !
      soilc  (1:nx,1:ny)                 = -999.9 !23
      so2    (1:nx,1:ny)                 = -999.9 !26
      xlai   (1:nx,1:ny)                 = -999.9 !24
      xsai   (1:nx,1:ny)                 = -999.9 !25
      eah    (1:nx,1:ny)                 = -999.9 !26
      tah    (1:nx,1:ny)                 = -999.9 !27

      fwet   (1:nx,1:ny)                 = -999.9 !28
      sneqvo (1:nx,1:ny)                 = -999.9 !29
      albold (1:nx,1:ny)                 = -999.9 !30
      qsnow  (1:nx,1:ny)                 = -999.9 !31
      wslake (1:nx,1:ny)                 = -999.9 !32
      ficeold(1:nx,1:ny,-nsnow+1:    0)  = -999.9 ! 5

      cm     (1:nx,1:ny)                 = -999.9 !
      ch     (1:nx,1:ny)                 = -999.9 !

      kr     (1:nx,1:ny)                 = -999.9 !
      froot  (1:nx,1:ny,1:nsoil)         = -999.9 !
      nroot  (1:nx,1:ny)                 = -999
      mq     (1:nx,1:ny)                 = -999.9 !
      qroot  (1:nx,1:ny,1:nsoil)         = -999.9 !
      sadr   (1:nx,1:ny,1:nsoil)         = -999.9 !

      allocate( bexp   (1:nx,1:ny) )
      allocate( dksat  (1:nx,1:ny,1:nsoil) )
      allocate( dwsat  (1:nx,1:ny) )
      allocate( f1     (1:nx,1:ny) )
      allocate( psisat (1:nx,1:ny) )
      allocate( quartz (1:nx,1:ny) )
      allocate( smcdry (1:nx,1:ny) )
      allocate( smcmax (1:nx,1:ny) )
      allocate( smcref (1:nx,1:ny) )
      allocate( smcwlt (1:nx,1:ny) )
      allocate( hs     (1:nx,1:ny) )
      allocate( rgl    (1:nx,1:ny) )
      allocate( rsmin  (1:nx,1:ny) )
      allocate( slope  (1:nx,1:ny) )
      allocate( kdt    (1:nx,1:ny) )
      allocate( frzx   (1:nx,1:ny) )

      bexp   (1:nx,1:ny) = -999.9
      dksat  (1:nx,1:ny,1:nsoil) = -999.9
      dwsat  (1:nx,1:ny) = -999.9
      f1     (1:nx,1:ny) = -999.9
      psisat (1:nx,1:ny) = -999.9
      quartz (1:nx,1:ny) = -999.9
      smcdry (1:nx,1:ny) = -999.9
      smcmax (1:nx,1:ny) = -999.9
      smcref (1:nx,1:ny) = -999.9
      smcwlt (1:nx,1:ny) = -999.9
      hs     (1:nx,1:ny) = -999.9
      rgl    (1:nx,1:ny) = -999.9
      rsmin  (1:nx,1:ny) = -999.9
      slope  (1:nx,1:ny) = -999.9
      kdt    (1:nx,1:ny) = -999.9
      frzx   (1:nx,1:ny) = -999.9

! for output

      allocate( qin    (1:nx,1:ny) )
      allocate( runsf  (1:nx,1:ny) )
      allocate( runsb  (1:nx,1:ny) )
      allocate( ecan   (1:nx,1:ny) )
      allocate( edir   (1:nx,1:ny) )  
      allocate( etran  (1:nx,1:ny) ) 
      allocate( fsa    (1:nx,1:ny) ) 
      allocate( fsr    (1:nx,1:ny) ) 
      allocate( fira   (1:nx,1:ny) )
      allocate( fsh    (1:nx,1:ny) )
      allocate( fgev   (1:nx,1:ny) )
      allocate( fctr   (1:nx,1:ny) )
      allocate( fcev   (1:nx,1:ny) )
      allocate( flh    (1:nx,1:ny) )
      allocate( fgh    (1:nx,1:ny) )
      allocate( apar   (1:nx,1:ny) )
      allocate( psn    (1:nx,1:ny) )
      allocate( sav    (1:nx,1:ny) )
      allocate( sag    (1:nx,1:ny) )
      allocate( fsnog  (1:nx,1:ny) )
      allocate( fsnoc  (1:nx,1:ny) )
      allocate( scf    (1:nx,1:ny) )
      allocate( swdown (1:nx,1:ny) )
      allocate( trad   (1:nx,1:ny) )
      allocate( t2m    (1:nx,1:ny) )
      allocate( nee    (1:nx,1:ny) )
      allocate( gpp    (1:nx,1:ny) )
      allocate( npp    (1:nx,1:ny) )
      allocate( qco2   (1:nx,1:ny) )
      allocate( vmax   (1:nx,1:ny) )
      allocate( km     (1:nx,1:ny) )
      allocate( vmaxup (1:nx,1:ny) )
      allocate( kmup   (1:nx,1:ny) )
      allocate( epslon (1:nx,1:ny) )
      allocate( fveg   (1:nx,1:ny) )
      allocate( ndvi   (1:nx,1:ny) )
      allocate( psi    (1:nx,1:ny, 1:nsoil) )
      allocate( wcnd   (1:nx,1:ny, 1:nsoil) )
      allocate( qsubcan(1:nx,1:ny) )
      allocate( qsubgrd(1:nx,1:ny) )

      qin    (1:nx,1:ny) = -999.9
      runsf  (1:nx,1:ny) = -999.9
      runsb  (1:nx,1:ny) = -999.9
      ecan   (1:nx,1:ny) = -999.9
      edir   (1:nx,1:ny) = -999.9
      etran  (1:nx,1:ny) = -999.9
      fsa    (1:nx,1:ny) = -999.9
      fsr    (1:nx,1:ny) = -999.9
      fira   (1:nx,1:ny) = -999.9
      fsh    (1:nx,1:ny) = -999.9
      fgev   (1:nx,1:ny) = -999.9
      fctr   (1:nx,1:ny) = -999.9
      fcev   (1:nx,1:ny) = -999.9
      flh    (1:nx,1:ny) = -999.9
      fgh    (1:nx,1:ny) = -999.9
      apar   (1:nx,1:ny) = -999.9
      psn    (1:nx,1:ny) = -999.9
      sav    (1:nx,1:ny) = -999.9
      sag    (1:nx,1:ny) = -999.9
      fsnog  (1:nx,1:ny) = -999.9
      fsnoc  (1:nx,1:ny) = -999.9
      scf    (1:nx,1:ny) = -999.9
      swdown (1:nx,1:ny) = -999.9
      trad   (1:nx,1:ny) = -999.9
      t2m    (1:nx,1:ny) = -999.9
      nee    (1:nx,1:ny) = -999.9
      gpp    (1:nx,1:ny) = -999.9
      npp    (1:nx,1:ny) = -999.9
      qco2   (1:nx,1:ny) = -999.9
      vmax   (1:nx,1:ny) = -999.9
      km     (1:nx,1:ny) = -999.9
      vmaxup (1:nx,1:ny) = -999.9
      kmup   (1:nx,1:ny) = -999.9
      epslon (1:nx,1:ny) = -999.9
      fveg   (1:nx,1:ny) = -999.9
      ndvi   (1:nx,1:ny) = -999.9
      psi    (1:nx,1:ny,1:nsoil) = -999.9
      wcnd   (1:nx,1:ny,1:nsoil) = -999.9
      qsubcan(1:nx,1:ny) = -999.9
      qsubgrd(1:nx,1:ny) = -999.9

! monthly mean variables for output

      allocate(  snowhdm  (1:nx,1:ny) )
      allocate(  sneqvdm  (1:nx,1:ny) )
      allocate(  scfdm    (1:nx,1:ny) )
      allocate(  fsnogdm  (1:nx,1:ny) )
      allocate(  runsfdm  (1:nx,1:ny) )
      allocate(  runsbdm  (1:nx,1:ny) )
      allocate(  qsubcandm(1:nx,1:ny) )
      allocate(  qsubgrddm(1:nx,1:ny) )

      allocate(  snowhm  (1:nx,1:ny) )
      allocate(  sneqvm  (1:nx,1:ny) )
      allocate(  scfm    (1:nx,1:ny) )
      allocate(  fsnogm  (1:nx,1:ny) )
      allocate(  tgm     (1:nx,1:ny) )
      allocate(  prcpm   (1:nx,1:ny) )
      allocate(  qinm    (1:nx,1:ny) )
      allocate(  runsfm  (1:nx,1:ny) )
      allocate(  runsbm  (1:nx,1:ny) )
      allocate(  ecanm   (1:nx,1:ny) )
      allocate(  edirm   (1:nx,1:ny) )
      allocate(  etranm  (1:nx,1:ny) )
      allocate(  zwtm    (1:nx,1:ny) )
      allocate(  fsam    (1:nx,1:ny) )
      allocate(  firam   (1:nx,1:ny) )
      allocate(  fshm    (1:nx,1:ny) )
      allocate(  flhm    (1:nx,1:ny) )
      allocate(  fghm    (1:nx,1:ny) )
      allocate(  stcm    (1:nx,1:ny,1:nsoil) )
      allocate(  smcm    (1:nx,1:ny,1:nsoil) )
      allocate(  sh2om   (1:nx,1:ny,1:nsoil) )
      allocate(  aparm   (1:nx,1:ny) )
      allocate(  psnm    (1:nx,1:ny) )
      allocate(  savm    (1:nx,1:ny) )
      allocate(  sagm    (1:nx,1:ny) )
      allocate(  xlaim   (1:nx,1:ny) )
      allocate(  xsaim   (1:nx,1:ny) )
      allocate(  tradm   (1:nx,1:ny) )
      allocate(  t2mm    (1:nx,1:ny) )
      allocate(  neem    (1:nx,1:ny) )
      allocate(  gppm    (1:nx,1:ny) )
      allocate(  nppm    (1:nx,1:ny) )
      allocate(  fvegm   (1:nx,1:ny) )
      allocate(  cmm     (1:nx,1:ny) )
      allocate(  chm     (1:nx,1:ny) )
      allocate(  qco2m   (1:nx,1:ny) )
      allocate(  socm    (1:nx,1:ny) )
      allocate(  wdocm   (1:nx,1:ny) )
      allocate(  ddocm   (1:nx,1:ny) )
      allocate(  micm    (1:nx,1:ny) )
      allocate(  wenzm   (1:nx,1:ny) )
      allocate(  denzm   (1:nx,1:ny) )
      allocate(  vmaxm   (1:nx,1:ny) )
      allocate(  kmm     (1:nx,1:ny) )
      allocate(  vmaxupm (1:nx,1:ny) )
      allocate(  kmupm   (1:nx,1:ny) )
      allocate(  epslonm (1:nx,1:ny) )
      allocate(  wam     (1:nx,1:ny) )
      allocate(  canwatm (1:nx,1:ny) )
      allocate(  ndvim   (1:nx,1:ny) )
      allocate(  nndvi   (1:nx,1:ny) )
      allocate(  sfctmpm (1:nx,1:ny) )
      allocate(  q2m     (1:nx,1:ny) )
      allocate(  windm   (1:nx,1:ny) )
      allocate(  sfcprsm (1:nx,1:ny) )
      allocate(  mqm     (1:nx,1:ny) )
      allocate(  krm     (1:nx,1:ny) )
      allocate(  frootm  (1:nx,1:ny,nsoil) )
      allocate(  rtmassm (1:nx,1:ny) )
      allocate(  sadrm   (1:nx,1:ny,nsoil) )
      allocate(  psim    (1:nx,1:ny,nsoil) )
      allocate(  qrootm  (1:nx,1:ny,nsoil) )
      allocate(  qrootnm (1:nx,1:ny,nsoil) )

      snowhdm  (1:nx,1:ny) = -999.9
      sneqvdm  (1:nx,1:ny) = -999.9
      scfdm    (1:nx,1:ny) = -999.9  !daily mean
      fsnogdm  (1:nx,1:ny) = -999.9  !daily mean
      runsfdm  (1:nx,1:ny) = -999.9  !daily mean
      runsbdm  (1:nx,1:ny) = -999.9  !daily mean
      qsubcandm(1:nx,1:ny) = -999.9  !daily mean
      qsubgrddm(1:nx,1:ny) = -999.9  !daily mean

      snowhm  (1:nx,1:ny) = -999.9
      sneqvm  (1:nx,1:ny) = -999.9
      scfm    (1:nx,1:ny) = -999.9
      fsnogm  (1:nx,1:ny) = -999.9
      tgm     (1:nx,1:ny) = -999.9
      prcpm   (1:nx,1:ny) = -999.9
      qinm    (1:nx,1:ny) = -999.9
      runsfm  (1:nx,1:ny) = -999.9
      runsbm  (1:nx,1:ny) = -999.9
      ecanm   (1:nx,1:ny) = -999.9
      edirm   (1:nx,1:ny) = -999.9
      etranm  (1:nx,1:ny) = -999.9
      zwtm    (1:nx,1:ny) = -999.9
      fsam    (1:nx,1:ny) = -999.9
      firam   (1:nx,1:ny) = -999.9
      fshm    (1:nx,1:ny) = -999.9
      flhm    (1:nx,1:ny) = -999.9
      fghm    (1:nx,1:ny) = -999.9
      stcm    (1:nx,1:ny,1:nsoil) = -999.9
      smcm    (1:nx,1:ny,1:nsoil) = -999.9
      sh2om   (1:nx,1:ny,1:nsoil) = -999.9
      aparm   (1:nx,1:ny) = -999.9
      psnm    (1:nx,1:ny) = -999.9
      savm    (1:nx,1:ny) = -999.9
      sagm    (1:nx,1:ny) = -999.9
      xlaim   (1:nx,1:ny) = -999.9
      xsaim   (1:nx,1:ny) = -999.9
      tradm   (1:nx,1:ny) = -999.9
      t2mm    (1:nx,1:ny) = -999.9
      neem    (1:nx,1:ny) = -999.9
      gppm    (1:nx,1:ny) = -999.9
      nppm    (1:nx,1:ny) = -999.9
      fvegm   (1:nx,1:ny) = -999.9
      cmm     (1:nx,1:ny) = -999.9
      chm     (1:nx,1:ny) = -999.9
      qco2m   (1:nx,1:ny) = -999.9
      socm    (1:nx,1:ny) = -999.9
      wdocm   (1:nx,1:ny) = -999.9
      ddocm   (1:nx,1:ny) = -999.9
      micm    (1:nx,1:ny) = -999.9
      wenzm   (1:nx,1:ny) = -999.9
      denzm   (1:nx,1:ny) = -999.9
      vmaxm   (1:nx,1:ny) = -999.9
      kmm     (1:nx,1:ny) = -999.9
      vmaxupm (1:nx,1:ny) = -999.9
      kmupm   (1:nx,1:ny) = -999.9
      epslonm (1:nx,1:ny) = -999.9
      wam     (1:nx,1:ny) = -999.9
      canwatm (1:nx,1:ny) = -999.9
      ndvim   (1:nx,1:ny) = -999.9
      nndvi   (1:nx,1:ny) = -999
      sfctmpm (1:nx,1:ny) = -999.9
      q2m     (1:nx,1:ny) = -999.9
      windm   (1:nx,1:ny) = -999.9
      sfcprsm (1:nx,1:ny) = -999.9
      mqm     (1:nx,1:ny) = -999.9
      krm     (1:nx,1:ny) = -999.9
      frootm  (1:nx,1:ny,1:nsoil) =-999.9
      rtmassm (1:nx,1:ny) = -999.9
      sadrm   (1:nx,1:ny,1:nsoil) =-999.9
      psim    (1:nx,1:ny,1:nsoil) =-999.9
      qrootm  (1:nx,1:ny,1:nsoil) =-999.9
      qrootnm (1:nx,1:ny,1:nsoil) =-999.9
!----------------------------------------------------------------------
! Initialize gridded domain
!----------------------------------------------------------------------
! sldpth is the thickness of each layer

      sldpth(1) = -zsoil(1)                               ! can be sent to REDPRM
      do iz = 2, nsoil
         sldpth(iz) = zsoil(iz-1)-zsoil(iz)
      enddo
 
!---------------------------------------------------------------------
! read landuse type and soil texture and other information
!----------------------------------------------------------------------

      call readland(dir,nx,ny,vegtyp,soiltyp,sloptyp, &
                   lat,lon,landsea,tbot,isc)

       foln   = 1.5        !saturated value
       z0     = 0.01       !ground roughness length

       ist(:,:) = 1

       do ix = 1,nx
       do iy = 1,ny
        if(vegtyp(ix,iy)  == 16) ist(ix,iy) = 2   ! lake points
        if(vegtyp(ix,iy) > 0) &
        call redprm(vegtyp(ix,iy),soiltyp(ix,iy),sloptyp(ix,iy),sldpth,zsoil,nsoil, & ! in
                    bexp(ix,iy)  ,dwsat(ix,iy)  ,dksat(ix,iy,:)  ,f1(ix,iy)    , & ! out
                    psisat(ix,iy),quartz(ix,iy) ,smcdry(ix,iy) ,smcmax(ix,iy), & ! out
                    smcref(ix,iy),smcwlt(ix,iy) ,hs(ix,iy)     ,rgl(ix,iy)   , & ! out
                    rsmin(ix,iy) ,nroot(ix,iy)  ,slope(ix,iy)  ,kdt(ix,iy)   , & ! out
                    frzx(ix,iy)  )
       end do
       end do

! green vegetation fraction

      call readveg(dir,nx,ny,landsea,shdfac)

!----------------------------------------------------------------------
! initialize model state
!----------------------------------------------------------------------

     call readinit(dir    ,nx      ,ny      ,nsoil   ,nsnow   ,fini  ,&
                   soiltyp,nsltype ,maxsmc  ,zsoil   ,albold  ,nroot ,&
                   fwet   ,sneqvo  ,qsnow   ,wslake  ,eah     ,tah   ,&
                   smc    ,stc     ,sh2o    ,tsno    ,snice   ,snliq ,&
                   zsnso  ,isnow   ,snowh   ,sneqv   ,canliq  ,canice,&
                   tg     ,tv      ,wa      ,wt      ,zwt     ,lfmass,&
                   rootms ,stmass  ,wood    ,soc     ,wdoc    ,ddoc  ,&
                   mic    ,wenz    ,denz    ,so2     ,xlai    ,sldpth,&
                   xsai   ,cm      ,ch      ,vegtyp  ,mq      ,kr    ,&
                   froot  )

!----------------------------------------------------------------------

      nspin = 1 !niu

! ================================================================================
! BEGIN TIME LOOP
! ================================================================================
      do ispin = 1, nspin

      it = 0
      yloop : do iyloop  = start_year, end_year
       startmonth = 1
       if(iyloop == start_year) startmonth = start_month

       if(mod(iyloop,4).eq.0) then
          nday(2) = 29
          istep   = gday(startmonth-1)*8
       else
          nday(2) = 28
          istep   = jday(startmonth-1)*8
       end if
       startstep  = istep+1

      mloop : do imloop  = startmonth,  12

       startday = 1
       if(iyloop == start_year .and. imloop == start_month) &
       startday = start_day

       imstep = 0
      dloop : do idloop  = startday, nday(imloop)
       idstep = 0
      hloop : do ihloop  = 1, nst
           
       it     = it     + 1
       istep  = istep  + 1
       imstep = imstep + 1
       idstep = idstep + 1

       read(mdt(1:2),'(i2)') imonth
       read(mdt(3:4),'(i2)') iday
       read(mdt(5:6),'(i2)') itime

       write(6,'(7i8)') it,iyloop,imloop,idloop,ihloop,itime,ispin

       !read forcing data.

       call readforc(nx   ,ny      ,dt      ,iyloop  ,imloop  ,idloop    ,&
                     itime,dir     ,sfctmp,q2    ,lwdn  ,uu     ,&
                     vv  ,sfcprs,soldn ,prcp  )

       !calday to compute cosz
       call calendr (it, dt, imloop, idloop, itime, calday)

       !compute cosz
       !$OMP PARALLEL DO
       do iy=1,ny
       do ix=1,nx
          call lsmzen (calday,lon(ix,iy),lat(ix,iy),cosz(ix,iy))
       end do
       end do
       !$OMP END PARALLEL DO

       !$OMP PARALLEL DO
       do 100 iy=1,ny
       do 200 ix=1,nx

       if(vegtyp(ix,iy) > 0) then
       !if(ix == nx-79 .and. iy == ny-40) then

         ice(ix,iy)     = abs(landsea(ix,iy)-1)
        !co2air(ix,iy)  = (372.55+it*0.000203*3.)*1.e-06 * sfcprs(ix,iy) !372.55 ppm (2000/01) -> 401 ppm (2015/12)
        !co2air(ix,iy)  = (337.+it*0.000203*3.)*1.e-06 * sfcprs(ix,iy) !337 ppm (1980/01) -> 401 ppm (2015/12)
        !co2air(ix,iy)  = (280.+it*0.000203*3.)*1.e-06 * sfcprs(ix,iy) !280 ppm (1948/01) -> 401 ppm (2015/12)
         co2air(ix,iy)  = 360.0 * 1.0e-06 * sfcprs(ix,iy) !280 ppm (1948/01) -> 401 ppm (2015/12)

         o2air(ix,iy)   = 0.209 * sfcprs(ix,iy)    !partial pressure  o2 (pa)

!         if(ix == int(nx/2) .and. iy == int(ny/2) ) &
!         write(88,'(4i5,9f10.2)')iyloop,imloop,idloop,itime,cosz(ix,iy)&
!                 ,soldn(ix,iy),lwdn(ix,iy),sfctmp(ix,iy)&
!                 ,uu(ix,iy),prcp(ix,iy)*dt,q2(ix,iy)*1000,sfcprs(ix,iy)&
!                 ,co2air(ix,iy)/sfcprs(ix,iy)*1.e6

         do iz=isnow(ix,iy)+1,0
           ficeold(ix,iy,iz) = snice(ix,iy,iz)/(snice(ix,iy,iz)+snliq(ix,iy,iz))
         end do

         call SFLX (ice(ix,iy)          ,ist(ix,iy)    ,vegtyp(ix,iy) ,&  !in
                    isc(ix,iy)          ,nsnow         ,nsoil         ,&  !in
                    zsoil(:)            ,dt            ,q2(ix,iy)     ,&  !in
                    sfctmp(ix,iy)       ,uu(ix,iy)     ,vv(ix,iy)     ,&  !in
                    soldn(ix,iy)        ,lwdn(ix,iy)   ,prcp(ix,iy)   ,&  !in
                    zlvl                ,co2air(ix,iy) ,o2air(ix,iy)  ,&  !in
                    cosz(ix,iy)         ,tbot(ix,iy)   ,foln          ,&  !in
                    sfcprs(ix,iy)       ,imloop        ,iday          ,&  !in
                    shdfac(ix,iy,imloop),lat(ix,iy)    ,z0            ,&  !in
                    ix                  ,iy            ,bexp(ix,iy)   ,&  !in
                    dwsat(ix,iy)        ,dksat(ix,iy,:)  ,f1(ix,iy)     ,&  !in
                    psisat(ix,iy)       ,quartz(ix,iy) ,smcdry(ix,iy) ,&  !in
                    smcmax(ix,iy)       ,smcref(ix,iy) ,smcwlt(ix,iy) ,&  !in
                    hs(ix,iy)           ,rgl(ix,iy)    ,rsmin(ix,iy)  ,&  !in
                    nroot(ix,iy)        ,slope(ix,iy)  ,kdt(ix,iy)    ,&  !in
                    frzx(ix,iy)         ,&  !in
                    eah(ix,iy)          ,tah(ix,iy)    ,fwet(ix,iy)   ,&  !inout
                    ficeold(ix,iy,:)    ,qsnow(ix,iy)  ,sneqvo(ix,iy) ,&  !inout
                    isnow(ix,iy)        ,zsnso(ix,iy,:),canliq(ix,iy) ,&  !inout
                    canice(ix,iy)       ,snowh(ix,iy)  ,sneqv(ix,iy)  ,&  !inout
                    snice(ix,iy,:)      ,snliq(ix,iy,:),tv(ix,iy)     ,&  !inout
                    tg(ix,iy)           ,stc(ix,iy,:)  ,sh2o(ix,iy,:) ,&  !inout
                    smc(ix,iy,:)        ,zwt(ix,iy)    ,wa(ix,iy)     ,&  !inout
                    wt(ix,iy)           ,wslake(ix,iy) ,lfmass(ix,iy) ,&  !inout
                    rootms(ix,iy,:)     ,stmass(ix,iy) ,wood(ix,iy)   ,&  !input
                    soc(ix,iy)          ,wdoc(ix,iy)   ,ddoc(ix,iy)   ,&  !inout
                    mic(ix,iy)          ,wenz(ix,iy)   ,denz(ix,iy)   ,&  !inout
                    so2(ix,iy)          ,xlai(ix,iy)   ,xsai(ix,iy)   ,&  !inout
                    albold(ix,iy)       ,cm(ix,iy)     ,ch(ix,iy)     ,&  !inout
                    sice(ix,iy,:)       ,mq(ix,iy)     ,kr(ix,iy)     ,&  !inout
                    froot(ix,iy,:)      ,                              &  !inout
                    fsa(ix,iy)          ,fsr(ix,iy)    ,fira(ix,iy)   ,&  !out
                    fsh(ix,iy)          ,fgh(ix,iy)    ,fcev(ix,iy)   ,&  !out
                    fgev(ix,iy)         ,fctr(ix,iy)   ,trad(ix,iy)   ,&  !out
                    ecan(ix,iy)         ,etran(ix,iy)  ,edir(ix,iy)   ,&  !out
                    runsf(ix,iy)        ,runsb(ix,iy)  ,apar(ix,iy)   ,&  !out
                    psn(ix,iy)          ,sav(ix,iy)    ,sag(ix,iy)    ,&  !out
                    fsnog(ix,iy)        ,nee(ix,iy)    ,gpp(ix,iy)    ,&  !out
                    npp(ix,iy)          ,t2m(ix,iy)    ,fveg(ix,iy)   ,&  !out
                    qco2(ix,iy)         ,vmax(ix,iy)   ,km(ix,iy)     ,&  !out
                    vmaxup(ix,iy)       ,kmup(ix,iy)   ,epslon(ix,iy) ,&  !out
                    qin(ix,iy)          ,ndvi(ix,iy)   ,swdown(ix,iy) ,&  !out
                    qroot(ix,iy,:)      ,sadr(ix,iy,:) ,psi(ix,iy,:)  ,&  !out
                    qsubcan(ix,iy)      ,qsubgrd(ix,iy))  !out

        endif         ! endif of land-points

  200   continue
  100   continue
        !$OMP END PARALLEL DO

        nn      = 0
        nndvig  = 0
        soilcgm = 0.
        socgm   = 0.
        wdocgm  = 0.
        ddocgm  = 0.
        micgm   = 0.
        wenzgm  = 0.
        denzgm  = 0.
        ndvigm  = 0.

        !!$OMP PARALLEL DO
        do iy=1,ny
        do ix=1,nx
           !for output
           if(vegtyp(ix,iy) > 0) then
             tsno(ix,iy,isnow(ix,iy)+1:0) = stc(ix,iy,isnow(ix,iy)+1:0)
             flh(ix,iy)   = fcev(ix,iy) + fgev(ix,iy) + fctr(ix,iy)

             fsnoc(ix,iy)    = 0.
             if(tv(ix,iy).le.273.16) fsnoc(ix,iy)     = fwet(ix,iy)
             scf(ix,iy) = (1.-fveg(ix,iy)) * fsnog(ix,iy) +   &
                              fveg(ix,iy)  * fsnoc(ix,iy)

             soilc(ix,iy) = soc(ix,iy) + wdoc(ix,iy) + ddoc(ix,iy) &
                          + mic(ix,iy) + wenz(ix,iy) + denz(ix,iy)
             soilcgm = soilcgm + soilc(ix,iy)
             socgm   = socgm   + soc  (ix,iy)
             wdocgm  = wdocgm  + wdoc (ix,iy)
             ddocgm  = ddocgm  + ddoc (ix,iy)
             micgm   = micgm   + mic  (ix,iy)
             wenzgm  = wenzgm  + wenz (ix,iy)
             denzgm  = denzgm  + denz (ix,iy)

             if(ndvi(ix,iy) > 0.0) then
                nndvig  = nndvig + 1
                ndvigm  = ndvigm  + ndvi (ix,iy)
             end if

             nn      = nn + 1
           end if
        end do
        end do
        !!$OMP END PARALLEL DO

        soilcgm = soilcgm/nn
        socgm   = socgm  /nn
        wdocgm  = wdocgm /nn
        ddocgm  = ddocgm /nn
        micgm   = micgm  /nn
        wenzgm  = wenzgm /nn
        denzgm  = denzgm /nn
        if (nndvig > 0) then
           ndvigm  = ndvigm /nndvig
        end if

        !write(*,*) 'Global SOILC,SOC,WDOC,DDOC,MIC,WENZ,DENZ (g/m2)'
        !write(6,'(7i8)') it,iyloop,imloop,idloop,ihloop,istep
        !write(400,'(2i6,2i3,2f10.1,5f9.1,2f9.3)')ispin,iyloop,imloop,idloop,soilcgm,&
        !       socgm,wdocgm,ddocgm,micgm,wenzgm,denzgm,ndvigm
!-------------------------------------------------------------------
! end of 1-d Noah processing
!-------------------------------------------------------------------

       !if(ispin == nspin) then
          !monthly output, one file per month (1 time layer)
          call nc_out(nsoil   ,nsnow   ,nroot   ,nx      ,ny      ,it      ,dt    , &
                      iyloop  ,imloop  ,idloop  ,DIR    ,EXP     ,lon   , &
                      lat     ,vegtyp  ,imstep  ,nday    ,ND      , &
                      snowh   ,sneqv   ,tg      ,stc     ,smc     , &
                      sh2o    ,prcp    ,runsf   ,runsb   ,ecan    , &
                      edir    ,etran   ,zwt     ,fsa     ,fira    , &
                      fsh     ,flh     ,fgh     ,apar    ,psn     , &
                      sav     ,sag     ,scf     ,xlai    ,xsai    , &
                      trad    ,nee     ,gpp     ,npp     ,t2m     , &
                      fveg    ,cm      ,ch      ,qin     ,qco2    , &
                      soc     ,wdoc    ,ddoc    ,mic     ,wenz    , &
                      denz    ,vmax    ,km      ,vmaxup  ,kmup    , &
                      epslon  ,wa      ,canliq  ,canice  ,ndvi    , &
                      nndvi   ,mq      ,kr      ,rootms  ,sadr    , &
                      qroot   ,psi     ,fsnog   , &
                      sfctmp  ,q2      ,uu      ,vv      ,sfcprs  , &
                      snowhm  ,sneqvm  ,tgm     ,stcm    ,smcm    , &
                      sh2om   ,prcpm   ,runsfm  ,runsbm  ,ecanm   , &
                      edirm   ,etranm  ,zwtm    ,fsam    ,firam   , &
                      fshm    ,flhm    ,fghm    ,aparm   ,psnm    , &
                      savm    ,sagm    ,scfm    ,xlaim   ,xsaim   , &
                      tradm   ,neem    ,gppm    ,nppm    ,t2mm    , &
                      fvegm   ,cmm     ,chm     ,qinm    ,qco2m   , &
                      socm    ,wdocm   ,ddocm   ,micm    ,wenzm   , &
                      denzm   ,vmaxm   ,kmm     ,vmaxupm ,kmupm   , &
                      epslonm ,wam     ,canwatm ,ndvim   ,mqm     , &
                      krm     ,rtmassm ,sadrm   ,qrootm  ,qrootnm , &
                      sfctmpm ,q2m     ,windm   ,sfcprsm ,psim    , &
                      fsnogm  )

         !daily output, one file per day

          call nc_out_daily(nx    ,ny     ,it    ,dt    , &
                           iyloop ,imloop ,idloop,DIR   , &
                           EXP    ,lon    ,lat   ,vegtyp, &
                           idstep ,ntimed ,  &
                           snowh    ,sneqv    ,scf       ,fsnog    , &
                           runsf    ,runsb    ,qsubcan   ,qsubgrd  , &
                           snowhdm  ,sneqvdm  ,scfdm     ,fsnogdm  , &
                           runsfdm  ,runsbdm  ,qsubcandm ,qsubgrddm)

         !hrly output, one file per day (24 time layers)

         if(iyloop == 2100) then 
          call nc_out_hr(nsoil      ,nx    ,ny    ,it    ,dt  , &
                      iyloop,imonth ,iday  ,DIR   ,EXP   ,lon , &
                      lat   ,vegtyp ,idstep,nday  ,        &
                      snowh ,sneqv  ,tg    ,stc   ,smc   , &
                      sh2o  ,prcp   ,runsf ,runsb ,ecan  , &
                      edir  ,etran  ,zwt   ,fsa   ,fira  , &
                      fsh   ,flh    ,fgh   ,apar  ,psn   , &
                      sav   ,sag    ,fsnog ,xlai  ,xsai  , &
                      trad  ,nee    ,gpp   ,npp   ,t2m   , &
                      fveg  ,qco2   ,soc   , &
                      wdoc  ,ddoc   ,mic   ,wenz  ,denz  )
         end if
      ! endif

!update the time 

          call geth_newdate(newdate, olddate, nint(dt))
          year    = newdate(1:4)
          mdt     = newdate(5:10)
          minute  = newdate(11:12)
          olddate = newdate

       enddo hloop
       enddo dloop

       call write_ini(DIR,EXP    ,nx     ,ny      ,nsoil   ,nsnow   , nroot, &
           iyloop   ,imloop  ,idloop-1   ,itime   ,lat   ,lon  ,vegtyp, &
           smc    ,stc   ,sh2o  ,tsno    ,snice ,snliq , &
           zsnso  ,isnow ,snowh ,sneqv   ,canliq,canice, &
           tg     ,tv    ,wa    ,wt      ,zwt   ,lfmass, &
           rootms ,stmass,wood  ,soc     ,wdoc  ,ddoc  , &
           mic    ,wenz  ,denz  ,wslake  ,mq    ,kr    , &
           froot  )

       enddo mloop   
       enddo yloop   

       enddo ! spinup
   
      end program Noah_driver
