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

  use ioipsl_getincom 
  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,save :: oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
      logical,save :: activerunoff ! enable simple runoff scheme?
      logical,save :: 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, save :: ivap, iliq, iice

      logical, save :: firstcall

      data firstcall /.true./


      if(firstcall)then

         oceanbulkavg=.false.
         oceanalbvary=.false.
         write(*,*)"Activate runnoff into oceans?"
         activerunoff=.false.
         call getin("activerunoff",activerunoff)
         write(*,*)" activerunoff = ",activerunoff
	 
	 
	 
         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
