!*==TESTTEST.spg  processed by SPAG 6.72Dc at 19:44 on  5 May 2021
!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*
!                                                                      *
!     HYDRUS   - Numerical model of one-dimensional variably saturated *
!                water flow, heat transport, and transport of solutes  *
!                involved in sequential first-order decay reactions    *
!                version 7.0                                           *
!                                                                      *
!                Version coupled with user-friendly interface for      *
!                MS Windows environment                                *
!                                                                      *
!     Designed by J.Simunek, M. Sejna, and M. Th. van Genuchten (1996) *
!                                                                      *
!                                        Last modified: October, 2009  *
!                                                                      *
!                                                                      *
!     COPYRIGHT (c) 2005-9, Jiri Simunek, PC-Progress                  *
!                                                                      *
!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*
 
      SUBROUTINE HYDRUS(hydName,UZWT_elev_dhru)
	  !use parm, only: curyr,id1,idaf
	  !use smrt_parm, only UZWT_elev_dhru
      IMPLICIT NONE
!*--TESTTEST24
!*** Start of declarations inserted by SPAG
      REAL aeroterm , ah , ahw , ainterc , ak , aks , akw , albedo ,    &
         & albedon , ampl , aosm , aqh , area , ath , aths , athw ,     &
         & basegw , beta , bqh , c50
      REAL cap , captab , catm , cbot , ccuma , ccumt , chpar ,         &
         & cloudf_ac , cloudf_bc , cmin , cnew , con , conc , conlt ,   &
         & cono , conr , consat , consmax , const , contab
      REAL convh , convt , cosalf , courant , cprevo , croot ,          &
         & crootmax , cropheight , cropheightn , ct , ctemp , ctola ,   &
         & ctolr , ctop , cumch , cumq , cvbot , cvch0 , cvch1 , cvchim
      REAL cvchr , cvoli , cvtop , disp , dmoist , dmul , dmul2 ,       &
         & dsurf , dt , dtinit , dtmax , dtmaxa , dtmaxc , dtmaxt ,     &
         & dtmin , dtold , dtopt , dtsnow , dummy1 , dummy2
      REAL dz , eamean , eameann , entres , epsi , epsit , etcomb ,     &
         & evap , evapp , excesint , g0 , g1 , geofac , gwl0l ,         &
         & hbot , hcrita , hcrits , heatfl , hnew
      REAL hold , hroot , hsat , hseep , ht , htab , htemp , htop ,     &
         & omegac , omegas , omegaw , p0 , p2h , p2l , p3 ,             &
         & p3c , pard , parw , peclet
      REAL pecr , poptm , prec , precc , precd , q0 , q1 , r2h , r2l ,  &
         & rad , radius , radn , rado , radterm , ralt , rbot , retard ,&
         & rextinct , rgr , rgrowth
      REAL rhmax , rhmaxn , rhmaxo , rhmean , rhmeann , rhmeano ,       &
         & rhmin , rhminn , rhmino , rh_a , rinterc , rkhbot , rkhtop , &
         & rkm , rkvbot , rkvtop , rlamb , rlat , rlwrada , rlwrada1
      REAL rlwradb , rlwradb1 , rn , rnl , rns , rr , rroot , rrootd ,  &
         & rs , rsoil , rsoild , rspacing , rst , rtop , r_h ,          &
         & sensflux , shwrada , shwradb , sig$abort , sig$break
      REAL sig$fpe , sig$int , SIGNALQQ , sink , sinkim , snowlayer ,   &
         & snowmf , solin , sorb , sorb2 , sorbn , sorbn2 , spot ,      &
         & spot1 , ssink , strans , subvol , sunhours , sunhoursn ,     &
         & sunhourso
      REAL tauw , tbot , tconv , tdep , tempa , tempheight , tempn ,    &
         & tempo , temps , theq , thetab , thnew , thnewim , thold ,    &
         & tholdim , thr , thrr , ths , thvnew , thvold
      REAL tmax1 , tmaxa , tmaxan , tmaxao , tmaxn , tmaxo , tmin1 ,    &
         & tmina , tminan , tminao , tminn , tmino , tolh , tolth ,     &
         & tpar , tperiod , tpulse , transp , trharv , trmin
      REAL trperiod , ttop , vnew , vold , vroot , vtop , vvnew ,       &
         & vvold , watin , wc , wcuma , wcumt , wdep , wetper , wind ,  &
         & windheight , windn , windo , wtransf , wvoli
      REAL x , xconv , xlai , xlain , xlat , xrmax , xrmin , xroot ,    &
         & xrootn , xsurf , zbotdr , zintf , UZWT_elev_dhru !SYL
      INTEGER i , i4ret , icampbell , icount , icrop , idualpor ,       &
            & ienhanc , ierr , ihyst , iinterc , ikappa , ikod , ilai , &
            & ilengthpath , imethour , imodel , imoistdep , inonconv ,  &
            & iposdr , iradiation , INOD !SYL
      INTEGER irelhum , irootin , isunsh , itcum , itcumo , itemp ,     &
            & iter , iterc , iterw , itmax , itmin , itort , iunit ,    &
            & iver , jj , jprint , kappa , kappao , kbold , kbotch
      INTEGER kbott , kodbot , kodtop , ktold , ktopch , ktopt ,        &
            & laynum , m , matnum , maxal , maxit , maxitc ,            &
            & NARGS , ngrowth , nlay , nmat , NMATD , nobs , NOBSD
      INTEGER node , NPD , nprstep , ns , NSD , ntab , NTABD , ntabmod ,&
            & numnp , NUMNPD , NUNITD
      INTEGER curyr,id1,idaf
      CHARACTER(len=5) :: hydName
!*** End of declarations inserted by SPAG
 
!#######################################################################
!     Procedura pro osetreni nasledujicich signalu:
!         - signal CTRL+C
!         - signal CTRK+BREAK
!         - abnormal termination
!         - floating point error
!-----------------------------------------------------------------------
 

!#######################################################################
 
      PARAMETER (NUMNPD=1001,NMATD=20,NTABD=100,NOBSD=100,NSD=11,       &
               & NUNITD=7,NPD=1000)
 
      INTEGER plevel , alevel , tlevel , err
      LOGICAL sinkf , wlayer , qgwlf , topinf , shorto , lwat , lchem , &
            & ltemp , convgf , freed , seepf , botinf , atmbc , lroot , &
            & lupw , lwtdep , lstdep , lequil , llinear(NSD) , lartd ,  &
            & lmosink , lsolred , lsoladd , lscreen , ltable ,          &
            & lminstep , qdrain , lmssink , ltort , linitw , lvarbc ,   &
            & lprint , lfiltr , lmobim(NMATD) , lmeteo , lbact ,        &
            & lvapor , ldayvar , lenter , lenbal , lvaporout , lextrap ,&
            & lprintd , llai , ldensity , lsnow , lcentrif , leqinit ,  &
            & lsinprec , ldualneq , lflux , lhargr , lmetdaily , lend , &
            & lmassini , lstopconv , lactrsu , lomegaw , lfluxout
      DOUBLE PRECISION p , r , s , q , rtime1 , t , tinit , told ,      &
                     & RTIME , tmax , tatm , tprint1 , tatm1 , tatm2 ,  &
                     & tatmold , tatmn , tatm2o , tprint , tprintint
      CHARACTER cfilename*260 , cdatapath*260
      INTEGER*2 status
      INTEGER*2 i2 , iyear , imonth , iday , ihours , imins , isecs ,   &
              & i100th
 
      DIMENSION x(NUMNPD) , hnew(NUMNPD) , hold(NUMNPD) , htemp(NUMNPD) &
              & , matnum(NUMNPD) , pard(11,NMATD) , tprint(NPD) ,       &
              & sink(NUMNPD) , poptm(NMATD) , beta(NUMNPD) ,            &
              & laynum(NUMNPD) , cumq(12) , subvol(10) ,                &
              & contab(NTABD,NMATD) , captab(NTABD,NMATD) ,             &
              & thetab(NTABD,NMATD) , con(NUMNPD) , cap(NUMNPD) ,       &
              & consat(NMATD) , ths(NMATD) , p(NUMNPD) , r(NUMNPD) ,    &
              & s(NUMNPD) , hsat(NMATD) , parw(11,NMATD) , area(10) ,   &
              & chpar(NSD*16+4,NMATD) , conc(NSD,NUMNPD) , vnew(NUMNPD) &
              & , vold(NUMNPD) , disp(NUMNPD) , retard(NUMNPD) ,        &
              & g0(NUMNPD) , g1(NUMNPD) , q(NUMNPD) , wc(NUMNPD) ,      &
              & thnew(NUMNPD) , thold(NUMNPD) , watin(NUMNPD) ,         &
              & solin(NUMNPD) , node(NOBSD) , iunit(NUNITD) , cvch0(NSD)&
              & , cvch1(NSD) , ctop(NSD) , cbot(NSD) , croot(NSD) ,     &
              & cumch(10,NSD) , cvoli(NSD) , ccuma(NSD) , ccumt(NSD) ,  &
              & cvtop(NSD) , cvbot(NSD) , cnew(NUMNPD) , ctemp(NUMNPD) ,&
              & cprevo(NUMNPD) , tdep(NSD*16+4) , tpar(10,NMATD) ,      &
              & tempo(NUMNPD) , tempn(NUMNPD) , sorbn(NUMNPD) ,         &
              & sorb(NSD,NUMNPD) , q0(NUMNPD) , q1(NUMNPD) , aosm(NSD) ,&
              & thr(NMATD) , ssink(NUMNPD) , ah(NUMNPD) , ak(NUMNPD) ,  &
              & ath(NUMNPD) , ahw(NMATD) , athw(NMATD) , akw(NMATD) ,   &
              & kappa(NUMNPD) , aths(NUMNPD) , thrr(NUMNPD) , cvchr(NSD)&
              & , conr(NUMNPD) , cono(NUMNPD) , aks(NUMNPD) ,           &
              & kappao(NUMNPD) , crootmax(NSD) , cvchim(NSD) ,          &
              & htab(NTABD,NMATD) , ntab(NMATD) , rgrowth(1000,5) ,     &
              & conlt(NUMNPD) , convt(NUMNPD) , convh(NUMNPD) ,         &
              & theq(NUMNPD) , thvold(NUMNPD) , thvnew(NUMNPD) ,        &
              & vvold(NUMNPD) , vvnew(NUMNPD) , sorb2(NSD,NUMNPD) ,     &
              & sorbn2(NUMNPD) , thnewim(NUMNPD) , tholdim(NUMNPD) ,    &
              & sinkim(NUMNPD) , strans(NUMNPD) , dmoist(NMATD,NSD,13,6)&
              & , wdep(2+NMATD,NSD*9) , ct(NSD)
 
      DATA lstopconv , inonconv/.TRUE. , 0/
 
      icount = NARGS()
      IF ( icount.GT.1 ) THEN
         i2 = 1
         CALL GETARG(i2,cdatapath,status)
      ELSE
         cfilename = hydName//'LEVEL_01.DIR'
         OPEN (10,FILE=cfilename,STATUS='old',ERR=900)
         READ (10,99001,ERR=1100) cdatapath
 
99001    FORMAT (a)
         CLOSE (10)
      ENDIF
      WRITE (*,*) cfilename
      WRITE (*,*) cdatapath
      ilengthpath = LEN_TRIM(cdatapath)
      IF ( ilengthpath.GT.260-13 ) THEN
         ierr = 30
         GOTO 1200
      ELSE
         cfilename = cdatapath(1:ilengthpath)//'/OPTIONS.IN'
         OPEN (35,FILE=cfilename,STATUS='old',ERR=100)
      ENDIF
      WRITE (*,*) ilengthpath
 
!     Initialization
 100  CALL INIT(cosalf,ntab,itcum,tlevel,alevel,plevel,hroot,vroot,     &
              & iterw,iterc,dtmaxc,wcumt,wcuma,err,lvarbc,NSD,croot,    &
              & ccumt,ccuma,NUMNPD,sink,wc,cumq,lmeteo,lbact,lvapor,    &
              & lenbal,ldayvar,lenter,lfiltr,tauw,nprstep,ntabmod,      &
              & idualpor,dtmaxt,lextrap,lprintd,llai,rextinct,ldensity, &
              & excesint,lminstep,lprint,lsnow,snowmf,snowlayer,ctemp,  &
              & isunsh,irelhum,xroot,lcentrif,radius,gwl0l,tprintint,   &
              & itort,ienhanc,hseep,leqinit,lsinprec,imoistdep,omegac,  &
              & wtransf,ldualneq,lflux,icrop,irootin,lmassini,lmetdaily,&
              & sorb2,lactrsu,omegas,spot,lomegaw,omegaw,lend,lvaporout,&
              & lfluxout)
      DATA iunit/50 , 70 , 71 , 75 , 76 , 77 , 78/
 
      cfilename = cdatapath(1:ilengthpath)//'/SELECTOR.IN'
      OPEN (30,FILE=cfilename,STATUS='old',ERR=900)
      WRITE (*,*) cfilename
      cfilename = cdatapath(1:ilengthpath)//'/PROFILE.DAT'
      OPEN (32,FILE=cfilename,STATUS='old',ERR=900)
      cfilename = cdatapath(1:ilengthpath)//'/I_CHECK.OUT'
      OPEN (50,FILE=cfilename,STATUS='unknown',ERR=1000)
      IF ( lprint ) THEN
         cfilename = cdatapath(1:ilengthpath)//'/RUN_INF.OUT'
         OPEN (70,FILE=cfilename,STATUS='unknown',ERR=1000)
         cfilename = cdatapath(1:ilengthpath)//'/T_LEVEL.OUT'
         OPEN (71,FILE=cfilename,STATUS='unknown',ERR=1000)
         cfilename = cdatapath(1:ilengthpath)//'/NOD_INF.OUT'
         OPEN (75,FILE=cfilename,STATUS='unknown',ERR=1000)
         cfilename = cdatapath(1:ilengthpath)//'/BALANCE.OUT'
         OPEN (76,FILE=cfilename,STATUS='unknown',ERR=1000)
         cfilename = cdatapath(1:ilengthpath)//'/OBS_NODE.OUT'
         OPEN (77,FILE=cfilename,STATUS='unknown',ERR=1000)
         cfilename = cdatapath(1:ilengthpath)//'/PROFILE.OUT'
         OPEN (78,FILE=cfilename,STATUS='unknown',ERR=1000)
      ENDIF
 
!     Read input data --------------------------------------------------
 
      CALL BASINF(cosalf,maxit,tolth,tolh,topinf,botinf,shorto,lwat,    &
                & lchem,sinkf,wlayer,qgwlf,freed,seepf,atmbc,kodtop,    &
                & kodbot,rtop,rroot,rbot,hcrits,hcrita,gwl0l,aqh,bqh,   &
                & ktold,kbold,NUNITD,iunit,nmat,NMATD,nlay,lroot,ltemp, &
                & lwtdep,lequil,lscreen,qdrain,zbotdr,basegw,rspacing,  &
                & iposdr,rkhtop,rkhbot,rkvtop,rkvbot,entres,wetper,     &
                & zintf,geofac,linitw,lvarbc,xconv,tconv,lmeteo,lvapor, &
                & iver,lprint,lcentrif,lsnow,hseep,lflux,lactrsu,err)
      IF ( err.NE.0 ) THEN
         IF ( err.EQ.1 ) THEN
            ierr = 5
         ELSEIF ( err.EQ.2 ) THEN
            ierr = 6
         ELSEIF ( err.EQ.3 ) THEN
            ierr = 3
         ELSEIF ( err.EQ.4 ) THEN
            ierr = 16
         ELSE
            GOTO 200
         ENDIF
         GOTO 1200
      ENDIF
 
 200  IF ( lprint ) THEN
         IF ( lfluxout ) THEN
            cfilename = cdatapath(1:ilengthpath)//'/T_Level1.OUT'
            OPEN (44,FILE=cfilename,STATUS='unknown',ERR=1000)
         ENDIF
         IF ( lvaporout .AND. lvapor ) THEN
            cfilename = cdatapath(1:ilengthpath)//'/NOD_INF_V.OUT'
            OPEN (45,FILE=cfilename,STATUS='unknown',ERR=1000)
         ENDIF
      ENDIF
      IF ( topinf .OR. botinf .OR. atmbc ) THEN
         cfilename = cdatapath(1:ilengthpath)//'/ATMOSPH.IN'
         OPEN (31,FILE=cfilename,STATUS='old',ERR=900)
         IF ( lprint ) THEN
            cfilename = cdatapath(1:ilengthpath)//'/A_LEVEL.OUT'
            OPEN (72,FILE=cfilename,STATUS='unknown',ERR=1000)
         ENDIF
         IF ( lmeteo ) THEN
            cfilename = cdatapath(1:ilengthpath)//'/Meteo.in'
            OPEN (33,FILE=cfilename,STATUS='old',ERR=1000)
            IF ( lprint ) THEN
               cfilename = cdatapath(1:ilengthpath)//'/Meteo.out'
               OPEN (43,FILE=cfilename,STATUS='unknown',ERR=1000)
            ENDIF
         ENDIF
      ENDIF
 
      CALL NODINF(NUMNPD,numnp,NOBSD,nobs,htop,hbot,x,hnew,hold,matnum, &
                & htemp,laynum,beta,ah,ak,ath,conc,sorb,tempn,tempo,    &
                & node,NSD,ns,xsurf,lchem,ltemp,lequil,lscreen,lbact,   &
                & sorb2,err,lprint,ldualneq,lflux)
      IF ( err.NE.0 ) THEN
         IF ( err.EQ.1 ) THEN
            ierr = 12
         ELSEIF ( err.EQ.2 ) THEN
            ierr = 3
         ELSEIF ( err.EQ.3 ) THEN
            ierr = 14
         ELSEIF ( err.EQ.4 ) THEN
            ierr = 15
         ELSEIF ( err.EQ.5 ) THEN
            ierr = 29
         ELSE
            GOTO 300
         ENDIF
         GOTO 1200
      ENDIF
      WRITE (*,*) 'test1'
 300  CALL MATIN(nmat,pard,parw,htab(1,1),htab(ntab(1),1),lscreen,err,  &
               & numnp,ah,ihyst,ahw,athw,akw,matnum,hnew,kappa,aths,    &
               & thrr,conr,aks,kappao,imodel,xconv,ltable,ikappa,       &
               & ntabmod,idualpor)
      IF ( err.NE.0 ) THEN
         IF ( err.EQ.1 ) THEN
            ierr = 6
         ELSEIF ( err.EQ.2 ) THEN
            ierr = 3
         ELSE
            GOTO 400
         ENDIF
         GOTO 1200
      ENDIF
      WRITE (*,*) 'test2'
 400  IF ( ihyst.GT.0 .AND. ihyst.NE.3 )                                &
         & CALL HYSTERIN(numnp,nmat,hold,matnum,pard,parw,thnew,thold,  &
         & kappa,aths,thrr,cono,conr,aks,kappao,ah,ak,ihyst,imodel,     &
         & cdatapath)
      IF ( imodel.EQ.ntabmod ) THEN
         cfilename = cdatapath(1:ilengthpath)//'/Mater.in'
         OPEN (36,FILE=cfilename,STATUS='old',ERR=900)
      ENDIF
      IF ( linitw .AND. imodel.LT.ntabmod ) THEN
         CALL INITW(numnp,nmat,matnum,kappa,hnew,hold,htemp,pard,parw,  &
                  & imodel,htop,hbot,idualpor,thnewim,err)
         IF ( err.NE.0 ) THEN
            ierr = 24
            GOTO 1200
         ENDIF
      ENDIF
      CALL GENMAT(ntab,NTABD,nmat,thr,ths,hsat,pard,htab,contab,captab, &
                & consat,thetab,imodel,lscreen,ntabmod,consmax,xconv,   &
                & tconv,err)
      WRITE (*,*) 'test3'
      IF ( err.EQ.1 ) THEN
         ierr = 3
         GOTO 1200
      ELSE
         IF ( ihyst.NE.3 ) THEN
            CALL SETMAT(numnp,ntab,NTABD,nmat,htab,contab,captab,hnew,  &
                      & matnum,pard,con,cap,consat,ah,ak,ath,hsat,htemp,&
                      & thetab,thold,thr,ths,lwtdep,tempn,iterw,cono,   &
                      & kappa,aths,thrr,conr,aks,ahw,athw,akw,imodel,   &
                      & ltable,lvapor,thvold,conlt,convt,convh,xconv,   &
                      & tconv,ntabmod,hcrita,ldensity,conc,NSD,ienhanc)
         ELSE
           ! Bob Lenhard hysteresis
            CALL HYST(numnp,NMATD,pard,parw,matnum,kappa,hnew,hold,     &
                    & thold,con,cap,ikappa,1)
         ENDIF
         CALL INITDUALPOR(numnp,nmat,matnum,pard,thold,idualpor,thnewim,&
                        & tholdim,sinkim,hnew,strans,linitw)
         CALL TMIN(tinit,tmax,tatm,told,dt,dtmax,dmul,dmul2,dtmin,      &
                 & tprint,t,dtopt,topinf,botinf,lscreen,itmin,itmax,    &
                 & maxal,hcrits,NPD,atmbc,iver,lprintd,nprstep,         &
                 & tprintint,lenter,ldayvar,lsinprec,llai,rextinct,err)
         IF ( err.NE.0 ) THEN
            IF ( err.EQ.1 ) THEN
               ierr = 7
            ELSEIF ( err.EQ.2 ) THEN
               ierr = 28
            ELSE
               GOTO 450
            ENDIF
            GOTO 1200
         ENDIF
 450     tprint1 = tmax
         IF ( lprintd ) tprint1 = tinit + tprintint
         IF ( lmeteo ) THEN
            CALL METEOIN(rlat,ralt,shwrada,shwradb,rlwrada,rlwradb,     &
                       & rlwrada1,rlwradb1,windheight,tempheight,icrop, &
                       & ilai,cropheight,albedo,xlai,xroot,iinterc,     &
                       & ainterc,ngrowth,rgrowth,rextinct,iradiation,   &
                       & lenbal,lprint,isunsh,irelhum,cloudf_ac,        &
                       & cloudf_bc,lhargr,lmetdaily,xconv,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 32
                  GOTO 1200
               ENDIF
            ENDIF
            IF ( lenbal ) icrop = 0
         ENDIF
         tatm2 = tmax
         dtinit = dt
 
         IF ( lroot ) THEN
            CALL ROOTIN(trmin,trharv,xrmin,xrmax,rgr,lscreen,iver,      &
                      & irootin,ngrowth,rgrowth,trperiod,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 8
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 3
               ELSE
                  GOTO 500
               ENDIF
               GOTO 1200
            ENDIF
         ENDIF
 500     IF ( lroot .OR. icrop.GT.1 )                                   &
            & CALL SETRG(numnp,x,beta,t,trmin,trharv,xrmin,xrmax,rgr,   &
            & xroot,lroot,irootin,ngrowth,rgrowth,trperiod)
         IF ( ltemp ) THEN
            CALL TEMPIN(nmat,tpar,ampl,tperiod,ktopt,ttop,kbott,tbot,   &
                      & topinf,botinf,icampbell,iver,snowmf,lscreen,    &
                      & ierr)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 10
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 3
               ELSE
                  GOTO 520
               ENDIF
               GOTO 1200
            ENDIF
 520        IF ( lenbal ) ktopt = -1
         ENDIF
 
         IF ( lchem ) THEN
            CALL CHEMIN(lupw,lstdep,nmat,ns,NSD,maxitc,chpar,tdep,      &
                      & ktopch,ctop,kbotch,cbot,epsi,tpulse,cumch,ctola,&
                      & ctolr,llinear,lequil,lartd,pecr,lscreen,dsurf,  &
                      & catm,ltort,lmobim,lbact,lfiltr,imoistdep,wdep,  &
                      & NMATD,imodel,pard,iver,ldualneq,lmassini,       &
                      & leqinit,itort,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 11
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 3
               ELSEIF ( err.EQ.3 ) THEN
                  ierr = 31
               ELSE
                  GOTO 540
               ENDIF
               GOTO 1200
            ENDIF
 540        IF ( idualpor.GT.0 ) lequil = .FALSE.
            IF ( imoistdep.EQ.2 ) CALL MOISTDEPIN(cdatapath,cfilename,  &
               & nmat,NMATD,ns,NSD,dmoist,imoistdep)
            IF ( lmassini ) CALL MASSINIT(numnp,NSD,ns,nmat,matnum,tdep,&
               & tempn,chpar,conc,thold,tholdim,ths,llinear,lbact)
            IF ( .NOT.lequil .AND. (leqinit .OR. lmassini) )            &
               & CALL NONEQINIT(numnp,NSD,ns,nmat,matnum,tdep,tempo,    &
                              & chpar,conc,sorb,llinear,lmobim,idualpor,&
                              & lbact,sorb2,thold)
            IF ( lprint ) THEN
               CALL OPENSOLUTEFILES(ns,cdatapath,ilengthpath,cfilename, &
                                  & err)
               IF ( err.EQ.1 ) GOTO 1000
            ENDIF
         ENDIF
         IF ( topinf .OR. botinf .OR. atmbc ) THEN
            CALL SETBC(tmax,tatm1,rtop,rroot,rbot,hcrita,hbot,htop,     &
                     & gwl0l,topinf,botinf,ct,cbot,ns,ttop,tbot,ampl,   &
                     & ltemp,lchem,kodtop,lvarbc,err,lminstep,lmeteo,   &
                     & prec,rsoil,llai,rextinct,lcentrif,cosalf,xconv,  &
                     & tconv,imodel,hnew(numnp),irootin,xroot,wlayer,   &
                     & llinear,lactrsu,spot)
            IF ( err.EQ.1 ) THEN
               ierr = 13
               GOTO 1200
            ELSE
               IF ( lmeteo ) THEN
                  CALL METEO(1,lmetdaily,ldayvar,t,dt,tinit,tmax,tatm2, &
                           & tatmn,tatm2o,dtmax,rlat,ralt,shwrada,      &
                           & shwradb,rlwrada,rlwradb,rlwrada1,rlwradb1, &
                           & windheight,tempheight,icrop,ilai,rroot,    &
                           & xconv,tconv,rgrowth,ngrowth,iinterc,       &
                           & rinterc,ainterc,excesint,lenbal,rextinct,  &
                           & lprint,lhargr,iradiation,isunsh,irelhum,   &
                           & imethour,cloudf_ac,cloudf_bc,prec,precc,   &
                           & rsoil,evapp,transp,rns,rnl,radterm,        &
                           & aeroterm,rst,etcomb,rad,radn,rado,wind,    &
                           & windn,windo,albedo,albedon,xlai,xlain,     &
                           & xroot,xrootn,cropheight,cropheightn,ampl,  &
                           & ttop,tmaxan,tminan,tmax1,tmaxn,tmaxo,tmin1,&
                           & tminn,tmino,tempa,tmaxa,tmaxao,tmina,      &
                           & tminao,sunhours,sunhoursn,sunhourso,rhmean,&
                           & rhmeann,rhmeano,rhmax,rhmaxn,rhmaxo,rhmin, &
                           & rhminn,rhmino,rh_a,eamean,eameann,rtop,err)
                  IF ( err.NE.0 ) THEN
                     IF ( err.EQ.1 ) THEN
                        ierr = 32
                     ELSEIF ( err.EQ.2 ) THEN
                        ierr = 13
                     ELSEIF ( err.EQ.3 ) THEN
                        ierr = 33
                     ELSE
                        GOTO 550
                     ENDIF
                     GOTO 1200
                  ENDIF
               ENDIF
 550           IF ( lchem .AND. snowlayer.LE.0 )                        &
                  & CALL SETCHEMBC(prec,rsoil,ns,ctop,ct,wlayer,        &
                  & hnew(numnp),kodtop,ktopch)
               IF ( lsnow ) THEN
                  dtsnow = SNGL(tatm1-tinit)
                  CALL SNOW(prec,dtsnow,ttop,snowmf,snowlayer,rsoil,    &
                          & xconv,lminstep,ctop,ct,ns)
               ENDIF
               tatm = DMIN1(tatm1,tatm2)
               IF ( ldayvar ) THEN
                  rrootd = rroot
                  rsoild = rsoil
                  CALL DAILYVAR(tconv,t,rroot,rrootd)
                  CALL DAILYVAR(tconv,t,rsoil,rsoild)
               ENDIF
               IF ( lsinprec ) THEN
                  precd = prec
                  tatmold = SNGL(tinit)
                  CALL SINPREC(t,tatmold,tatm1,prec,precd)
               ENDIF
               IF ( kodtop.EQ.-4 ) rtop = ABS(rsoil) - ABS(prec)
            ENDIF
         ENDIF
         IF ( sinkf ) THEN
            CALL SINKIN(nmat,lchem,lmosink,lsolred,lsoladd,p0,poptm,p2h,&
                      & p2l,p3,r2h,r2l,aosm,c50,p3c,ns,lmssink,crootmax,&
                      & iver,omegac,lactrsu,omegas,spot1,rkm,cmin,      &
                      & lomegaw,lscreen,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 9
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 3
               ELSE
                  GOTO 560
               ENDIF
               GOTO 1200
            ENDIF
 560        IF ( .NOT.topinf .AND. .NOT.botinf .AND. .NOT.atmbc )       &
               & spot = spot1
            CALL SETSNK(numnp,nmat,matnum,x,hroot,vroot,sink,rroot,hnew,&
                      & lmosink,lsolred,lsoladd,p0,poptm,p2h,p2l,p3,r2h,&
                      & r2l,aosm,c50,p3c,beta,lchem,ns,NSD,conc,croot,  &
                      & lmssink,thold,pard,dt,omegac,imodel,con,lomegaw,&
                      & omegaw,rbot)
         ENDIF
         CLOSE (30)
         CLOSE (32)
         CLOSE (50)
 
         IF ( lprint ) THEN
            CALL PROFIL(numnp,nmat,x,matnum,xsurf,beta,ah,ak,ath,thr,   &
                      & ths,consat,hsat,lscreen,err)
            IF ( err.EQ.1 ) THEN
               ierr = 23
               GOTO 1200
            ELSE
               CALL NODOUT(numnp,nmat,hnew,thold,con,x,xsurf,cosalf,    &
                         & tinit,matnum,cap,ak,sink,consat,ns,NSD,conc, &
                         & tempo,sorb,kappa,lbact,sorb2,lvapor,lwtdep,  &
                         & conlt,convt,convh,thold(numnp),dt,idualpor,  &
                         & thnewim,sinkim,strans,ldensity,lcentrif,     &
                         & radius,lvaporout,ldualneq,err)
               IF ( err.EQ.1 ) THEN
                  ierr = 20
                  GOTO 1200
               ENDIF
            ENDIF
         ENDIF
         CALL SUBREG(numnp,nmat,nlay,hnew,thold,thold,x,matnum,laynum,  &
                   & t-dt,dt,cosalf,con,lchem,conc,chpar,0,ths,wcumt,   &
                   & wcuma,ccumt,ccuma,wvoli,cvoli,watin,solin,lwat,    &
                   & ltemp,tempn,tpar,tdep,ns,NSD,sorb,llinear,lequil,  &
                   & lmobim,err,subvol,area,lprint,lbact,sorb2,lvapor,  &
                   & thvold,thvnew,lwtdep,conlt,convh,convt,idualpor,   &
                   & thnewim,tholdim,ldensity,lcentrif,radius,ldualneq, &
                   & cprevo)
         IF ( err.EQ.1 ) THEN
            ierr = 21
            GOTO 1200
         ELSE
            IF ( lchem .OR. ltemp .OR. lflux )                          &
               & CALL VELOC(numnp,hold,con,x,cosalf,vold,thold,thold,   &
               & sink,dt,lvapor,lwtdep,conlt,convt,convh,tempo,vvold,   &
               & thvold,thvold,ldensity,conc,NSD,lcentrif,radius)
            DO i = 1 , numnp
               vnew(i) = vold(i)
               thnew(i) = thold(i)
               IF ( lvapor ) thvnew(i) = thvold(i)
               IF ( lvapor ) vvnew(i) = vvold(i)
            ENDDO
            IF ( lbact .AND. .NOT.lwat )                                &
               & CALL EXCLUSION(numnp,nmat,NSD,pard,chpar,thnew,vnew,   &
               & thold,vold)
 
            IF ( lscreen ) WRITE (*,*) 'beginning of numerical solution'
            CALL GETDAT(iyear,imonth,iday)
            CALL GETTIM(ihours,imins,isecs,i100th)
            rtime1 = RTIME(imonth,iday,ihours,imins,isecs,i100th)
 
!     Time stepping calculations ---------------------------------------
 
!     Loop between water flow and heat and vapor transport
            itemp = 0
         ENDIF
      ENDIF



 600  IF ( lenbal ) THEN
         epsit = 0.5
         temps = epsit*tempn(numnp) + (1.-epsit)*tempo(numnp)
         CALL EVAPOR(t,temps,tmaxa,tmina,rad,hnew(numnp),tempheight,    &
                   & windheight,wind,rhmean,heatfl,rtop,prec,tperiod,   &
                   & rlat,albedo,sunhours,thnew(numnp),xconv,tconv,     &
                   & iradiation,rns,rnl,rn,sensflux,evap,xlat,const,    &
                   & isunsh,r_h,lmetdaily,rst,tempa,rh_a,shwrada,       &
                   & shwradb,rlwrada,rlwradb,imethour)
      ENDIF
	  
       !!!!!!!!SYL, update bottom pressure head and the pressure head of all the nodes in the soil profile.
	  curyr = 1
	  id1 = 1
	  idaf = 1
      IF (curyr.eq.1 .and. id1.eq.idaf) THEN
        hbot = hold(1)
        kodbot = ABS(kodbot)
      ELSE
        hbot = UZWT_elev_dhru-x(1) !UZWT_elev_dhru = water table elevation , x(1) = evelation of the bottom of the soil profile
        hold(1) = UZWT_elev_dhru-x(1)
        kodbot = ABS(kodbot)
        DO INOD=1,numnp-1
          IF ((x(INOD+1)-x(1)).GT.hold(1)) EXIT
          IF (kodtop.eq.1) THEN
            IF (INOD.eq.(numnp-1)) THEN
			  EXIT
			ELSE
              hold(INOD+1)=hold(INOD)-x(INOD+1)+x(INOD)
			ENDIF
          ELSE
            hold(INOD+1)=hold(INOD)-x(INOD+1)+x(INOD)
          ENDIF
        END DO
      ENDIF

      hnew(:)=hold(:)
      htemp(:)=hold(:)
      CALL SETMAT(numnp,ntab,NTABD,nmat,htab,contab,captab,hnew,  &
                      & matnum,pard,con,cap,consat,ah,ak,ath,hsat,htemp,&
                      & thetab,thold,thr,ths,lwtdep,tempn,iterw,cono,   &
                      & kappa,aths,thrr,conr,aks,ahw,athw,akw,imodel,   &
                      & ltable,lvapor,thvold,conlt,convt,convh,xconv,   &
                      & tconv,ntabmod,hcrita,ldensity,conc,NSD,ienhanc)
      DO i = 1 , Numnp
        !WRITE(*,*) thnew(i)
        !WRITE(*,*) thold(i) 
        thold(i)=thnew(i)
      END DO
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!SYL
	  
!     Solve water movement ---------------------------------------------
      IF ( lwat ) THEN
         itcumo = itcum
         CALL WATFLOW(numnp,ntab,NTABD,nmat,htab,contab,captab,hnew,    &
                    & hold,matnum,pard,parw,con,cap,consat,ah,ak,ath,   &
                    & hsat,htemp,kodtop,kodbot,rtop,rbot,cosalf,t,dt,x, &
                    & sink,p,r,s,freed,seepf,qgwlf,aqh,bqh,gwl0l,htop,  &
                    & hbot,hcrita,hcrits,wlayer,iterw,itcum,topinf,     &
                    & ktold,kbold,tolth,tolh,maxit,dtmin,told,dtopt,    &
                    & convgf,thetab,thnew,thold,thr,ths,lwtdep,tempn,   &
                    & kappa,kappao,aths,thrr,cono,conr,aks,ahw,athw,akw,&
                    & ihyst,imodel,qdrain,zbotdr,basegw,rspacing,iposdr,&
                    & rkhtop,rkhbot,rkvtop,rkvbot,entres,wetper,zintf,  &
                    & geofac,ltable,lvapor,xconv,tconv,conlt,convt,     &
                    & convh,tauw,theq,thvnew,thvold,ntabmod,idualpor,   &
                    & thnewim,tholdim,sinkim,vtop,tempo,itemp,wtransf,  &
                    & ldensity,conc,NSD,ienhanc,lcentrif,radius,hseep)
         IF ( .NOT.convgf ) THEN
            inonconv = inonconv + 1
         ELSE
            inonconv = 0
         ENDIF
         IF ( lstopconv .AND. .NOT.convgf .AND. inonconv.GE.10 ) THEN
            CALL CLOSEOUTPUT(rtime1,ns,topinf,botinf,lchem,lscreen,     &
                           & lmeteo,lprint)
            IF ( lenter ) THEN
               WRITE (*,*) 'Press Enter to continue'
               READ (*,*)
            ENDIF
            STOP
         ENDIF
         IF ( lmeteo .AND. lmetdaily .AND. .NOT.lenbal )                &
            & CALL DAYMETEOOUT(t,etcomb,evapp,transp,rns,rnl,radterm,   &
                             & aeroterm,precc,rinterc,excesint,tempa,   &
                             & rh_a,rst,lprint)  ! Output for daily variated meteo information
         IF ( lenbal ) THEN
            IF ( itcum.GT.itcumo+maxit .AND. .NOT.lmetdaily )           &
               & CALL METEOINT(3,t,tatm2o,tatm2,rad,rado,radn,tmaxa,    &
                             & tmaxao,tmaxan,tmina,tminao,tminan,wind,  &
                             & windo,windn,rhmean,rhmeano,rhmeann,      &
                             & sunhours,sunhourso,sunhoursn,lenbal)
            temps = epsit*tempn(numnp) + (1.-epsit)*tempo(numnp)
            m = matnum(numnp)
            rlamb = AMAX1(0.,tpar(4,m)+tpar(5,m)*thnew(numnp)+tpar(6,m) &
                  & *SQRT(thnew(numnp)))
            dz = x(numnp) - x(numnp-1)
            CALL UPDATEENERGY(t,vtop,rtop,heatfl,temps,rns,rnl,rn,evap, &
                            & xlat,sensflux,xconv,tconv,const,tlevel,   &
                            & nprstep,r_h,dz,rlamb,itemp,lprint,        &
                            & lmetdaily,tempa,rh_a,rst)
         ENDIF
      ELSE
         iterw = 1
         itcum = itcum + 1
      ENDIF
 
!     To calculate the velocities --------------------------------------
      IF ( lwat .AND. (ltemp .OR. lchem .OR. lflux) )                   &
         & CALL VELOC(numnp,hnew,con,x,cosalf,vnew,thnew,thold,sink,dt, &
         & lvapor,lwtdep,conlt,convt,convh,tempn,vvnew,thvnew,thvold,   &
         & ldensity,conc,NSD,lcentrif,radius)
 
!     Root zone calculations
      IF ( lroot .OR. icrop.GT.1 )                                      &
         & CALL SETRG(numnp,x,beta,t,trmin,trharv,xrmin,xrmax,rgr,xroot,&
         & lroot,irootin,ngrowth,rgrowth,trperiod)
      IF ( sinkf ) CALL SETSNK(numnp,nmat,matnum,x,hroot,vroot,sink,    &
                             & rroot,hnew,lmosink,lsolred,lsoladd,p0,   &
                             & poptm,p2h,p2l,p3,r2h,r2l,aosm,c50,p3c,   &
                             & beta,lchem,ns,NSD,conc,croot,lmssink,    &
                             & thnew,pard,dt,omegac,imodel,con,lomegaw, &
                             & omegaw,rbot)
 
!     Calculation of heat transport ------------------------------------
      IF ( ltemp ) THEN
         CALL TEMPER(numnp,nmat,x,dt,t,matnum,tempo,tempn,tpar,ampl,p,r,&
                   & s,q,vold,vnew,thold,thnew,retard,disp,sink,tperiod,&
                   & ktopt,ttop,kbott,tbot,lvapor,thvold,thvnew,vvold,  &
                   & vvnew,g0,lenbal,heatfl,xconv,tconv,dtmaxt,         &
                   & icampbell,itemp)
         IF ( (lvapor .OR. lenbal) .AND. itemp.LT.1 ) THEN
            itemp = itemp + 1
            GOTO 600
         ENDIF
      ENDIF
 
!     Calculations of the solute transport -----------------------------
      IF ( lchem ) THEN
         ikod = 0
         IF ( lflux ) ikod = 2
         CALL SOLUTE(numnp,nmat,ns,NSD,x,dt,t,tpulse,chpar,matnum,thold,&
                   & thnew,vold,vnew,disp,epsi,ktopch,ctop,kbotch,cbot, &
                   & conc,p,r,s,q,g0,g1,retard,cvtop,cvbot,cvch0,cvch1, &
                   & lupw,wc,peclet,courant,dtmaxc,tempo,tempn,cnew,    &
                   & cprevo,ctemp,tdep,ths,ctola,ctolr,iterc,maxitc,    &
                   & htemp,sorb,sorbn,llinear,lequil,lartd,pecr,q0,q1,  &
                   & dsurf,catm,ltort,sink,crootmax,ssink,cvchr,lmobim, &
                   & cvchim,tlevel,lbact,sorb2,sorbn2,dtmin,dtopt,lwat, &
                   & lfiltr,idualpor,tholdim,thnewim,sinkim,strans,     &
                   & itort,xconv,tconv,lvapor,rbot,err,imoistdep,NMATD, &
                   & dmoist,wdep,ikod,beta,ldualneq,atmbc,sinkf,lactrsu,&
                   & omegas,omegaw,spot,rkm,cmin,ldensity)
         IF ( err.NE.0 ) THEN
            ierr = 27
            GOTO 1200
         ENDIF
      ENDIF
 
!     Output ------------------------------------------------------------
!     T-level information
      IF ( ABS(t-tmax).LE.0.5*dtmin .OR. t.GT.tmax ) lend = .TRUE.
      jprint = 0
      IF ( (ABS(tprint(plevel)-t).LT.0.001*dt .OR.                      &
         & (lprintd .AND. (ABS(tprint1-t).LT.0.001*dt)) .OR.            &
         & (.NOT.shorto .AND.                                           &
         & ABS(FLOAT((tlevel+nprstep-1)/nprstep)-(tlevel+nprstep-1)     &
         & /FLOAT(nprstep)).LT.0.0001)) ) jprint = 1
 
      CALL TLINF(numnp,con,x,cosalf,t,dt,iterw,iterc,tlevel,rtop,rroot, &
               & vroot,hnew,hroot,cumq,itcum,kodtop,kodbot,convgf,lwat, &
               & lchem,croot,ns,NSD,conc,cvtop,cvbot,cvch0,cvch1,peclet,&
               & courant,wcumt,wcuma,ccumt,ccuma,cumch,thnew,thold,sink,&
               & lscreen,err,cvchr,cvchim,lprint,lvapor,lwtdep,conlt,   &
               & convh,convt,tempn,rsoil,prec,nprstep,thvold,thvnew,    &
               & xconv,idualpor,sinkim,wtransf,ldensity,snowlayer,      &
               & lcentrif,radius,thnewim,wlayer,hcrits,lend,lfluxout,   &
               & jprint,dummy1,dummy2,lflux,nobs,node,vnew,cnew,ctop)
      IF ( err.NE.0 ) THEN
         IF ( err.EQ.1 ) THEN
            ierr = 19
         ELSEIF ( err.EQ.2 ) THEN
!917   ierr=17
!      goto 1000
            ierr = 18
         ELSEIF ( err.EQ.3 ) THEN
            ierr = 26
         ELSE
            GOTO 700
         ENDIF
         GOTO 1200
      ENDIF
 
 700  IF ( nobs.GT.0 .AND. lprint .AND. jprint.EQ.1 ) THEN
         CALL OBSNOD(t,numnp,nobs,ns,NSD,node,conc,hnew,thnew,tempn,    &
                   & lchem,thnewim,vnew,vvnew,lflux,err)
         IF ( err.EQ.1 ) THEN
            ierr = 22
            GOTO 1200
         ENDIF
      ENDIF
      IF ( lprintd .AND. DABS(tprint1-t).LT.0.001*dt )                  &
         & tprint1 = tprint1 + tprintint
 
!     P-level information ----------------------------------------------
      IF ( ABS(tprint(plevel)-t).LT.0.001*dt ) THEN
         IF ( lprint ) THEN
            CALL NODOUT(numnp,nmat,hnew,thnew,con,x,xsurf,cosalf,       &
                      & tprint(plevel),matnum,cap,ak,sink,consat,ns,NSD,&
                      & conc,tempn,sorb,kappa,lbact,sorb2,lvapor,lwtdep,&
                      & conlt,convt,convh,thold(numnp),dt,idualpor,     &
                      & thnewim,sinkim,strans,ldensity,lcentrif,radius, &
                      & lvaporout,ldualneq,err)
            IF ( err.EQ.1 ) THEN
               ierr = 20
               GOTO 1200
            ENDIF
         ENDIF
         CALL SUBREG(numnp,nmat,nlay,hnew,thnew,thold,x,matnum,laynum,t,&
                   & dt,cosalf,con,lchem,conc,chpar,plevel,ths,wcumt,   &
                   & wcuma,ccumt,ccuma,wvoli,cvoli,watin,solin,lwat,    &
                   & ltemp,tempn,tpar,tdep,ns,NSD,sorb,llinear,lequil,  &
                   & lmobim,err,subvol,area,lprint,lbact,sorb2,lvapor,  &
                   & thvold,thvnew,lwtdep,conlt,convh,convt,idualpor,   &
                   & thnewim,tholdim,ldensity,lcentrif,radius,ldualneq, &
                   & cprevo)
         IF ( err.EQ.1 ) THEN
            ierr = 21
            GOTO 1200
         ELSE
            plevel = plevel + 1
         ENDIF
      ENDIF
 
!     A-level information ----------------------------------------------
      IF ( DABS(t-tatm).LE.0.001*dt .AND.                               &
         & (topinf .OR. botinf .OR. atmbc) ) THEN
         IF ( lprint ) THEN
            CALL ALINF(t,cumq,hnew(numnp),hroot,hnew(1),alevel,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 25
                  GOTO 1200
               ENDIF
            ENDIF
         ENDIF
         IF ( DABS(t-tatm1).LE.0.001*dt ) THEN
            tatmold = tatm1
            CALL SETBC(tmax,tatm1,rtop,rr,rbot,hcrita,hbot,htop,gwl0l,  &
                     & topinf,botinf,ct,cbot,ns,ttop,tbot,ampl,ltemp,   &
                     & lchem,kodtop,lvarbc,err,lminstep,lmeteo,prec,rs, &
                     & llai,rextinct,lcentrif,cosalf,xconv,tconv,imodel,&
                     & hnew(numnp),irootin,xroot,wlayer,llinear,lactrsu,&
                     & spot)
            IF ( err.EQ.1 ) THEN
               ierr = 13
               GOTO 1200
            ELSEIF ( .NOT.lmeteo ) THEN
               rroot = rr
               rsoil = rs
            ENDIF
         ENDIF
         IF ( lmeteo .AND. DABS(t-tatm2).LE.0.001*dt ) THEN
            CALL METEO(2,lmetdaily,ldayvar,t,dt,tinit,tmax,tatm2,tatmn, &
                     & tatm2o,dtmax,rlat,ralt,shwrada,shwradb,rlwrada,  &
                     & rlwradb,rlwrada1,rlwradb1,windheight,tempheight, &
                     & icrop,ilai,rroot,xconv,tconv,rgrowth,ngrowth,    &
                     & iinterc,rinterc,ainterc,excesint,lenbal,rextinct,&
                     & lprint,lhargr,iradiation,isunsh,irelhum,imethour,&
                     & cloudf_ac,cloudf_bc,prec,precc,rsoil,evapp,      &
                     & transp,rns,rnl,radterm,aeroterm,rst,etcomb,rad,  &
                     & radn,rado,wind,windn,windo,albedo,albedon,xlai,  &
                     & xlain,xroot,xrootn,cropheight,cropheightn,ampl,  &
                     & ttop,tmaxan,tminan,tmax1,tmaxn,tmaxo,tmin1,tminn,&
                     & tmino,tempa,tmaxa,tmaxao,tmina,tminao,sunhours,  &
                     & sunhoursn,sunhourso,rhmean,rhmeann,rhmeano,rhmax,&
                     & rhmaxn,rhmaxo,rhmin,rhminn,rhmino,rh_a,eamean,   &
                     & eameann,rtop,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 32
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 13
               ELSEIF ( err.EQ.3 ) THEN
                  ierr = 33
               ELSE
                  GOTO 750
               ENDIF
               GOTO 1200
            ENDIF
         ENDIF
 750     tatm = DMIN1(tatm1,tatm2)
         IF ( lchem .AND. snowlayer.LE.0 )                              &
            & CALL SETCHEMBC(prec,rsoil,ns,ctop,ct,wlayer,hnew(numnp),  &
            & kodtop,ktopch)
         IF ( lsnow ) THEN
            dtsnow = SNGL(tatm1-tatmold)
            CALL SNOW(prec,dtsnow,ttop,snowmf,snowlayer,rsoil,xconv,    &
                    & lminstep,ctop,ct,ns)
         ENDIF
         IF ( ldayvar ) THEN
            rrootd = rroot
            rsoild = rsoil
         ENDIF
         IF ( lsinprec ) precd = prec
         IF ( .NOT.lvarbc .AND. kodtop.EQ.-4 ) rtop = ABS(rsoil)        &
            & - ABS(prec)
         IF ( .NOT.lvarbc ) rtop = ABS(rsoil) - ABS(prec)
         alevel = alevel + 1
      ENDIF
 
      IF ( wlayer .AND. hnew(numnp).GT.0. ) THEN
                                            ! mass balance in the surface layer
         ht = hnew(numnp)
         DO jj = 1 , ns
            IF ( (ht+dt*(prec-rsoil)).GT.0. ) ctop(jj)                  &
               & = (ht*ctop(jj)+dt*prec*ct(jj))/(ht+dt*(prec-rsoil))
         ENDDO
      ENDIF
 
!     Time governing ---------------------------------------------------
      IF ( ABS(t-tmax).LE.0.5*dtmin .OR. t.GT.tmax ) THEN
         CALL CLOSEOUTPUT(rtime1,ns,topinf,botinf,lchem,lscreen,lmeteo, &
                        & lprint)
         GOTO 601
         ! IF ( lenter ) THEN
            ! WRITE (*,*) 'Press Enter to continue'
            ! READ (*,*)
         ! ENDIF
         ! STOP
      ELSE
         told = t
         dtold = dt
         ktold = kodtop
         kbold = kodbot
         IF ( .NOT.lwat ) iterw = 1
         iter = MAX0(iterw,iterc)
         dtmaxa = MIN(dtmaxc,dtmaxt)
         IF ( lsinprec .AND. precd.GT.0 )                               &
            & dtmaxa = MIN(dtmaxa,SNGL(tatm1-tatmold)/20.)
         CALL TMCONT(dt,dtmax,dtopt,dmul,dmul2,dtmin,iter,              &
                   & MIN(tprint(plevel),tprint1),tatm,t,tmax,dtmaxa,    &
                   & itmin,itmax,lminstep,dtinit)
         t = t + dt
         IF ( ldayvar ) THEN
            CALL DAILYVAR(tconv,t,rroot,rrootd)
            CALL DAILYVAR(tconv,t,rsoil,rsoild)
            rtop = ABS(rsoil) - ABS(prec)
         ENDIF
         IF ( lsinprec ) THEN
            CALL SINPREC(t,tatmold,tatm1,prec,precd)
            rtop = ABS(rsoil) - ABS(prec)
         ENDIF
         IF ( lmeteo ) THEN
            CALL METEO(3,lmetdaily,ldayvar,t,dt,tinit,tmax,tatm2,tatmn, &
                     & tatm2o,dtmax,rlat,ralt,shwrada,shwradb,rlwrada,  &
                     & rlwradb,rlwrada1,rlwradb1,windheight,tempheight, &
                     & icrop,ilai,rroot,xconv,tconv,rgrowth,ngrowth,    &
                     & iinterc,rinterc,ainterc,excesint,lenbal,rextinct,&
                     & lprint,lhargr,iradiation,isunsh,irelhum,imethour,&
                     & cloudf_ac,cloudf_bc,prec,precc,rsoil,evapp,      &
                     & transp,rns,rnl,radterm,aeroterm,rst,etcomb,rad,  &
                     & radn,rado,wind,windn,windo,albedo,albedon,xlai,  &
                     & xlain,xroot,xrootn,cropheight,cropheightn,ampl,  &
                     & ttop,tmaxan,tminan,tmax1,tmaxn,tmaxo,tmin1,tminn,&
                     & tmino,tempa,tmaxa,tmaxao,tmina,tminao,sunhours,  &
                     & sunhoursn,sunhourso,rhmean,rhmeann,rhmeano,rhmax,&
                     & rhmaxn,rhmaxo,rhmin,rhminn,rhmino,rh_a,eamean,   &
                     & eameann,rtop,err)
            IF ( err.NE.0 ) THEN
               IF ( err.EQ.1 ) THEN
                  ierr = 32
               ELSEIF ( err.EQ.2 ) THEN
                  ierr = 13
               ELSEIF ( err.EQ.3 ) THEN
                  ierr = 33
               ELSE
                  GOTO 800
               ENDIF
               GOTO 1200
            ENDIF
         ENDIF
 800     tlevel = tlevel + 1
         IF ( tlevel.GT.1000000 ) tlevel = 2
      ENDIF
 
!     New updated values
      CALL UPDATE(numnp,lwat,lchem,ltemp,lvapor,idualpor,lextrap,dt,    &
                & dtold,htemp,hnew,hold,thold,thnew,vold,vnew,thvold,   &
                & thvnew,vvold,vvnew,tholdim,thnewim,tempo,tempn,rtop,  &
                & xconv,consmax,kodtop,kodbot)
 
      itemp = 0
      GOTO 600
 
! --- End of time loop -------------------------------------------------
 
!     Error messages
 900  ierr = 1
      GOTO 1200
 1000 ierr = 2
      GOTO 1200
 1100 ierr = 4
 
 1200 CALL ERROROUT(ierr,cfilename,cdatapath,ilengthpath,lscreen)
      IF ( lenter ) THEN
         WRITE (*,*) ierr
         WRITE (*,*) 'Press Enter to continue'
         READ (*,*)
      ENDIF
      STOP
601   END
!*==ERROROUT.spg  processed by SPAG 6.72Dc at 19:44 on  5 May 2021
 
!***********************************************************************
 
      SUBROUTINE ERROROUT(Ierr,Cfilename,Cdatapath,Ilengthpath,Lscreen)
      IMPLICIT NONE
!*--ERROROUT976
!*** Start of declarations inserted by SPAG
      INTEGER Ierr , Ilengthpath
!*** End of declarations inserted by SPAG
 
      CHARACTER*260 cerr(33) , Cfilename , Cdatapath , cfilenameerr
      LOGICAL Lscreen
 
      cerr(1) = 'Open file error in file :'
      cerr(2) =                                                         &
     &'File already exists or hard disk is full ! Open file error in out&
     &put file : '
      cerr(3) = 'Error when writing to an output file !'
      cerr(4) =                                                         &
     &'Error when reading from an input file Level_01.dir data pathway !&
     &'
      cerr(5) =                                                         &
     &'Error when reading from an input file Selector.in Basic Informati&
     &ons !'
      cerr(6) =                                                         &
     &'Error when reading from an input file Selector.in Water Flow Info&
     &rmations !'
      cerr(7) =                                                         &
     &'Error when reading from an input file Selector.in Time Informatio&
     &ns !'
      cerr(8) =                                                         &
     &'Error when reading from an input file Selector.in Root Growth Inf&
     &ormations !'
      cerr(9) =                                                         &
     &'Error when reading from an input file Selector.in Sink Informatio&
     &ns !'
      cerr(10) =                                                        &
     &'Error when reading from an input file Selector.in Heat Transport &
     &Informations !'
      cerr(11) =                                                        &
     &'Error when reading from an input file Selector.in Solute Transpor&
     &t Informations !'
      cerr(12) = 'Error when reading from an input file Profile.dat !'
      cerr(13) = 'Error when reading from an input file Atmosph.in !'
      cerr(14) = 'Dimension in NumNPD is exceeded !'
      cerr(15) = 'Dimension in NObsD is exceeded !'
      cerr(16) = 'Dimension in NMatD or NLay is exceeded !'
      cerr(17) = 'Error when writing into an output file I_CHECK.OUT !'
      cerr(18) = 'Error when writing into an output file RUN_INF.OUT !'
      cerr(19) = 'Error when writing into an output file T_LEVEL.OUT !'
      cerr(20) = 'Error when writing into an output file NOD_INF.OUT !'
      cerr(21) = 'Error when writing into an output file BALANCE.OUT !'
      cerr(22) = 'Error when writing into an output file OBS_NODE.OUT !'
      cerr(23) = 'Error when writing into an output file PROFILE.OUT !'
      cerr(24) = 'Initial water content condition is lower than Qr !'
      cerr(25) = 'Error when writing into an output file A_LEVEL.OUT !'
      cerr(26) = 'Error when writing into an output file SOLUTE.OUT !'
      cerr(27) = 'Does not converge in the solute transport module !'
      cerr(28) = 'Number of Print-Times is exceeded !'
      cerr(29) = 'Dimension in NSD is exceeded !'
      cerr(30) = 'The path to the project is too long !!!'
      cerr(31) = 'Bulk density can not be equal to zero !'
      cerr(32) = 'Error when reading from an input file Meteo.in !'
      cerr(33) =                                                        &
     &'Crop Height must be smaller than the height of wind and temperatu&
     &re measurements!'
 
      cfilenameerr = Cdatapath(1:Ilengthpath)//'/Error.msg'
      OPEN (99,FILE=cfilenameerr,STATUS='unknown',ERR=100)
      IF ( Ierr.LE.2 ) THEN
         IF ( Lscreen ) WRITE (*,*) cerr(Ierr) , Cfilename
         WRITE (99,*) cerr(Ierr) , Cfilename
      ELSE
         IF ( Lscreen ) WRITE (*,*) cerr(Ierr)
         WRITE (99,*) cerr(Ierr)
      ENDIF
      CLOSE (99)
      RETURN
 
 100  WRITE (*,*)                                                       &
     &'Folder with input data of the specified project does not exist or&
     & pathway is too long or corrupted'
      WRITE (*,*) Cfilename
      END
!*==CLOSEOUTPUT.spg  processed by SPAG 6.72Dc at 19:44 on  5 May 2021
 
!***********************************************************************
 
      SUBROUTINE CLOSEOUTPUT(Rtime1,Ns,Topinf,Botinf,Lchem,Lscreen,     &
                           & Lmeteo,Lprint)
      IMPLICIT NONE
!*--CLOSEOUTPUT1063
!*** Start of declarations inserted by SPAG
      INTEGER jj , MSFLIB , Ns
!*** End of declarations inserted by SPAG
 
      LOGICAL Topinf , Botinf , Lchem , Lscreen , Lmeteo , Lprint
      INTEGER*2 iyear , imonth , iday , ihours , imins , isecs , i100th
      DOUBLE PRECISION Rtime1 , rtime2 , RTIME
 
      CALL GETDAT(iyear,imonth,iday)
      CALL GETTIM(ihours,imins,isecs,i100th)
      rtime2 = RTIME(imonth,iday,ihours,imins,isecs,i100th)
      IF ( Lscreen ) WRITE (*,*) 'Real time [sec]' , rtime2 - Rtime1
      IF ( Lprint ) THEN
         WRITE (70,'(''end'')')
         WRITE (76,*)
         WRITE (76,*) 'Calculation time [sec]' , rtime2 - Rtime1
         WRITE (71,'(''end'')')
         IF ( Topinf .OR. Botinf ) THEN
            WRITE (72,'(''end'')')
            IF ( Lmeteo ) WRITE (43,'(''end'')')
         ENDIF
         WRITE (77,'(''end'')')
         IF ( Lchem ) THEN
            DO jj = 1 , Ns
               WRITE (80+jj,'(''end'')')
            ENDDO
         ENDIF
      ENDIF
      END

!*==CLOSEFILES.spg  processed by SPAG 6.72Dc at 19:44 on  5 May 2021
 
!***********************************************************************
 
      SUBROUTINE CLOSEFILES
      IMPLICIT NONE
!*--CLOSEFILES1164
      LOGICAL lopen
 
      INQUIRE (UNIT=43,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (43,'(''end'')')
         CLOSE (43)
      ENDIF
      INQUIRE (UNIT=44,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (44,'(''end'')')
         CLOSE (44)
      ENDIF
      INQUIRE (UNIT=70,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (70,'(''end'')')
         CLOSE (70)
      ENDIF
      INQUIRE (UNIT=71,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (71,'(''end'')')
         CLOSE (71)
      ENDIF
      INQUIRE (UNIT=72,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (72,'(''end'')')
         CLOSE (72)
      ENDIF
      INQUIRE (UNIT=77,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (77,'(''end'')')
         CLOSE (77)
      ENDIF
      INQUIRE (UNIT=75,OPENED=lopen)
      IF ( lopen ) CLOSE (75)
      INQUIRE (UNIT=76,OPENED=lopen)
      IF ( lopen ) CLOSE (76)
      INQUIRE (UNIT=78,OPENED=lopen)
      IF ( lopen ) CLOSE (78)
      INQUIRE (UNIT=31,OPENED=lopen)
      IF ( lopen ) CLOSE (31)
      INQUIRE (UNIT=81,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (81,'(''end'')')
         CLOSE (81)
      ENDIF
      INQUIRE (UNIT=82,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (82,'(''end'')')
         CLOSE (82)
      ENDIF
      INQUIRE (UNIT=83,OPENED=lopen)
      IF ( lopen ) THEN
         WRITE (83,'(''end'')')
         CLOSE (83)
      ENDIF
      INQUIRE (UNIT=41,OPENED=lopen)
      IF ( lopen ) CLOSE (41)
      INQUIRE (UNIT=42,OPENED=lopen)
      IF ( lopen ) CLOSE (42)
 
      END
 
!***********************************************************************
