source: trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90 @ 2176

subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,  &
     qsurf,dqsurf,dqs_hyd,pcapcal,                &
     albedo,albedo_bareground,                    &
     albedo_snow_SPECTV,albedo_co2_ice_SPECTV,    &
     mu0,pdtsurf,pdtsurf_hyd,hice,                &
     pctsrf_sic,sea_ice)

  use ioipsl_getin_p_mod, only: getin_p
  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
  USE surfdat_h
  use comdiurn_h
  USE geometry_mod, only: cell_area
  USE tracer_h
  use slab_ice_h
  use callkeys_mod, only: albedosnow,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
  use radinc_h, only : L_NSPECTV

  implicit none

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

      integer ngrid,nq

!     Inputs
!     ------
      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?
!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)

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

      real ptimestep
      real mu0(ngrid)
      real qsurf(ngrid,nq), tsurf(ngrid)
      real dqsurf(ngrid,nq), pdtsurf(ngrid)
      real hice(ngrid)
      real albedo(ngrid,L_NSPECTV)
      real albedo_bareground(ngrid)
      real albedo_snow_SPECTV(L_NSPECTV)
      real albedo_co2_ice_SPECTV(L_NSPECTV)
      real pctsrf_sic(ngrid), sea_ice(ngrid)

      real oceanarea2

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

!     Local
!     -----
      real a,b,E
      integer ig,iq, nw
      real fsnoi, subli, fauxo
      real twater(ngrid)
      real pcapcal(ngrid)
      real hicebis(ngrid)
      real zqsurf(ngrid,nq)
      real ztsurf(ngrid)
      real albedo_sic, alb_ice
      real zfra

      integer, save :: ivap, iliq, iice
!$OMP THREADPRIVATE(ivap,iliq,iice)

      logical, save :: firstcall
!$OMP THREADPRIVATE(firstcall)

      data firstcall /.true./


      if(firstcall)then

         oceanbulkavg=.false.
         oceanalbvary=.false.
         write(*,*)"Activate runnoff into oceans?"
         activerunoff=.false.
         call getin_p("activerunoff",activerunoff)
         write(*,*)" activerunoff = ",activerunoff
         
         
         
         if (activerunoff) then
           ALLOCATE(runoff(ngrid))
           runoff(1:ngrid)=0
         endif

         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...

!     Total ocean surface area
         oceanarea=0.
         do ig=1,ngrid
            if(nint(rnat(ig)).eq.0)then
               oceanarea=oceanarea+cell_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(nint(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
               do nw=1,L_NSPECTV
                  albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in slab_ice_h.F90. Later we could introduce spectral dependency for alb_ocean.
               enddo
!            end if


            if(ok_slab_ocean) then
          
               zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0))     ! Snow Fraction (Critical height 45kg/m2~15cm)
               alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) &   ! Ice Albedo
               *exp(-sea_ice(ig)/h_alb_ice)
               ! Albedo final calculation :
               do nw=1,L_NSPECTV
                  albedo(ig,nw) = pctsrf_sic(ig)*                                        &
                                 (albedo_snow_SPECTV(nw)*zfra + alb_ice*(1.0-zfra))      &
                               + (1.-pctsrf_sic(ig))*alb_ocean
               enddo

               ! Oceanic ice height, just for diagnostics
               hice(ig)    = MIN(10.,sea_ice(ig)/rhowater)
            else !ok_slab_ocean


!     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              
                  do nw=1,L_NSPECTV
                     albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
                  enddo  

!                 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
                  DO nw=1,L_NSPECTV
                     albedo(ig,nw) = alb_ocean
                  ENDDO               

               endif

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

            endif!(ok_slab_ocean)


!     Continent
!     ---------
         elseif (nint(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
            DO nw=1,L_NSPECTV
               albedo(ig,nw) = albedo_bareground(ig)
            ENDDO
            if (zqsurf(ig,iice).ge.snowlayer) then
               DO nw=1,L_NSPECTV
                  albedo(ig,nw) = albedo_snow_SPECTV(nw)
               ENDDO
            else
               DO nw=1,L_NSPECTV
                  albedo(ig,nw) = albedo_bareground(ig)                            & 
                               + (albedo_snow_SPECTV(nw) - albedo_bareground(ig))  &
                                 *zqsurf(ig,iice)/snowlayer
               ENDDO
            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((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
               oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
            end if
         END DO
       
         tsea=0.
         DO ig=1,ngrid
            if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then       
               tsea=tsea+ztsurf(ig)*cell_area(ig)*pcapcal(ig)/oceanarea2
            end if
         END DO

         DO ig=1,ngrid
            if((nint(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 (nint(rnat(ig)).eq.1) then
               totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
            endif
         enddo
         
         do ig=1,ngrid
            if (nint(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

         do ig=1,ngrid
            if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
               DO nw=1,L_NSPECTV
                  albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
               ENDDO
            endif
         enddo ! ngrid
         
      endif ! co2cond


      do ig=1,ngrid ! We calculate here the tracer tendencies. Don't forget that we have to retrieve the dqsurf tendencies we added at the beginning of the routine !
         dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
         dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
      enddo

      if (activerunoff) then
        call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
      endif

      return
    end subroutine hydrol