subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,          &
     qsurf,dqsurf,dqs_hyd,pcapcal,                &
     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)

  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
  USE surfdat_h
  use comdiurn_h
  USE comgeomfi_h
  USE tracer_h

  implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Calculate the surface hydrology and albedo changes.
!     
!     Authors
!     ------- 
!     Adapted from LMDTERRE by B. Charnay (2010). Further 
!     modifications by R. Wordsworth (2010).
!     
!     Called by
!     ---------
!     physiq.F
!     
!     Calls
!     -----
!     none
!
!     Notes
!     -----
!     rnat is terrain type: 0-ocean; 1-continent
!     
!==================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"

      integer ngrid,nq

!     Inputs
!     ------
      real albedoocean
      parameter (albedoocean=0.07)
      real albedoice
      save albedoice

      real snowlayer
      parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm 
      real oceantime
      parameter (oceantime=10*24*3600)

      logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
      logical activerunoff ! enable simple runoff scheme?
      logical oceanalbvary ! ocean albedo varies with the diurnal cycle?

!     Arguments
!     ---------
      integer rnat(ngrid) ! I changed this to integer (RW)
      real,dimension(:),allocatable,save :: runoff
      real totalrunoff, tsea, oceanarea
      save oceanarea

      real ptimestep
      real mu0(ngrid)
      real qsurf(ngrid,nq), tsurf(ngrid)
      real dqsurf(ngrid,nq), pdtsurf(ngrid)
      real hice(ngrid)
      real albedo0(ngrid), albedo(ngrid)

      real oceanarea2

!     Output
!     ------
      real dqs_hyd(ngrid,nq)
      real pdtsurf_hyd(ngrid)

!     Local
!     -----
      real a,b,E
      integer ig,iq, icap ! wld like to remove icap
      real fsnoi, subli, fauxo
      real twater(ngrid)
      real pcapcal(ngrid)
      real hicebis(ngrid)
      real zqsurf(ngrid,nq)
      real ztsurf(ngrid)

      integer ivap, iliq, iice
      save ivap, iliq, iice

      logical firstcall
      save firstcall

      data firstcall /.true./

      oceanbulkavg=.false.
      activerunoff=.false.
      oceanalbvary=.false.

      if(firstcall)then

         ALLOCATE(runoff(ngrid))

         ivap=igcm_h2o_vap
         iliq=igcm_h2o_vap
         iice=igcm_h2o_ice
        
         write(*,*) "hydrol: ivap=",ivap
         write(*,*) "        iliq=",iliq
         write(*,*) "        iice=",iice

!     Here's the deal: iice is used in place of igcm_h2o_ice both on the 
!                      surface and in the atmosphere. ivap is used in
!                      place of igcm_h2o_vap ONLY in the atmosphere, while
!                      iliq is used in place of igcm_h2o_vap ONLY on the 
!                      surface.

!                      Soon to be extended to the entire water cycle...

!     Ice albedo = snow albedo for now
         albedoice=albedosnow

!     Total ocean surface area
         oceanarea=0.
         do ig=1,ngrid
            if(rnat(ig).eq.0)then
               oceanarea=oceanarea+area(ig)
            endif
         enddo

         if(oceanbulkavg.and.(oceanarea.le.0.))then
            print*,'How are we supposed to average the ocean'
            print*,'temperature, when there are no oceans?'
            call abort
         endif

         if(activerunoff.and.(oceanarea.le.0.))then
            print*,'You have enabled runoff, but you have no oceans.'
            print*,'Where did you think the water was going to go?'
            call abort
         endif
         
         firstcall = .false.
      endif

!     add physical tendencies already calculated
!     ------------------------------------------

      do ig=1,ngrid
         ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
         pdtsurf_hyd(ig)=0.0
         do iq=1,nq
            zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
         enddo    
      enddo
 
      do ig=1,ngrid
         do iq=1,nq
            dqs_hyd(ig,iq) = 0.0
         enddo
      enddo

      do ig = 1, ngrid

!     Ocean
!     -----
         if(rnat(ig).eq.0)then

!     re-calculate oceanic albedo
            if(diurnal.and.oceanalbvary)then
               fauxo      = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
               albedo(ig) = 1.1*( .03 + .630/( 1. + fauxo*fauxo))
               albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
            else
               albedo(ig) = albedoocean ! modif Benjamin
            end if

!     calculate oceanic ice height including the latent heat of ice formation
!     hice is the height of oceanic ice with a maximum of maxicethick.
            hice(ig)    = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall

!            twater(ig)  = tsurf(ig) + ptimestep*zdtsurf(ig) &
            twater(ig)  = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig)
            ! this is temperature water would have if we melted entire ocean ice layer
            hicebis(ig) = hice(ig)
            hice(ig)    = 0.

            if(twater(ig) .lt. T_h2O_ice_liq)then
               E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
               hice(ig)        = E/(RLFTT*rhowater)
               hice(ig)        = max(hice(ig),0.0)
               hice(ig)        = min(hice(ig),maxicethick)
               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep              
               albedo(ig)      = albedoice

!               if (zqsurf(ig,iice).ge.snowlayer) then
!                  albedo(ig) = albedoice
!               else
!                  albedo(ig) = albedoocean &
!                       + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
!               endif

            else

               pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep                
               albedo(ig)      = albedoocean

            endif

            zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
            zqsurf(ig,iice) = hice(ig)*rhowater


!     Continent
!     ---------
         elseif (rnat(ig).eq.1) then

!     melt the snow
            if(ztsurf(ig).gt.T_h2O_ice_liq)then
               if(zqsurf(ig,iice).gt.1.0e-8)then

                  a     = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
                  b     = zqsurf(ig,iice)
                  fsnoi = min(a,b)

                  zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
                  zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi

!                 thermal effects
                  pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep  

               endif
            else

!     freeze the water
               if(zqsurf(ig,iliq).gt.1.0e-8)then

                  a     = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
                  b     = zqsurf(ig,iliq)
                  
                  fsnoi = min(a,b)

                  zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
                  zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi

!                 thermal effects
                  pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep  

               endif
            endif
            
!     deal with runoff
            if(activerunoff)then

               runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
               if(ngrid.gt.1)then ! runoff only exists in 3D
                  if(runoff(ig).ne.0.0)then
                     zqsurf(ig,iliq) = mx_eau_sol
!                    runoff is added to ocean at end
                  endif
               end if

            endif

!     re-calculate continental albedo
            albedo(ig) = albedo0(ig) ! albedo0 = base values
            if (zqsurf(ig,iice).ge.snowlayer) then
               albedo(ig) = albedosnow
            else
               albedo(ig) = albedo0(ig) &
                    + (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
            endif

         else

            print*,'Surface type not recognised in hydrol.F!'
            print*,'Exiting...'
            call abort

         endif

      end do ! ig=1,ngrid

!     perform crude bulk averaging of temperature in ocean
!     ----------------------------------------------------
      if(oceanbulkavg)then

         oceanarea2=0.       
         DO ig=1,ngrid
            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
               oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
            end if
         END DO
       
         tsea=0.
         DO ig=1,ngrid
            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then       
               tsea=tsea+ztsurf(ig)*area(ig)*pcapcal(ig)/oceanarea2
            end if
         END DO

         DO ig=1,ngrid
            if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
               pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
            end if
         END DO

         print*,'Mean ocean temperature               = ',tsea,' K'

      endif

!     shove all the runoff water into the ocean
!     -----------------------------------------
      if(activerunoff)then

         totalrunoff=0.
         do ig=1,ngrid
            if (rnat(ig).eq.1) then
               totalrunoff = totalrunoff + area(ig)*runoff(ig)
            endif
         enddo
         
         do ig=1,ngrid
            if (rnat(ig).eq.0) then
               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
                    totalrunoff/oceanarea
            endif
         enddo

      endif

!     Re-add the albedo effects of CO2 ice if necessary
!     -------------------------------------------------
      if(co2cond)then

         icap=1
         do ig=1,ngrid
            if (qsurf(ig,igcm_co2_ice).gt.0) then
               albedo(ig) = albedice(icap)
            endif
         enddo
         
      endif


      do ig=1,ngrid
         dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
         dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
      enddo

      call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)

      return
    end subroutine hydrol