module process_cat
  
  ! reichle, 13 Jun 2005
  ! qliu+reichle, 14 Aug 2008 - major overhaul, use subroutine louissurface()
  ! reichle, 29 Nov 2010 - added Helfand Monin-Obukhov surface layer turbulence scheme
  ! reichle, 20 Dec 2011 - reinstated "PARdrct" and "PARdffs" for MERRA-Land file specs
  ! reichle, 28 Dec 2011 - removed field "totalb" from "cat_diagn" structure
  ! reichle, 15 Mar 2013 - minor changes to adjust to infrastructure changes 
  !                        in MAPL_Base.F90 of GEOSadas-5_9_1_p8
  ! reichle, 30 Oct 2013 - removed field rzeq from cat_diagn structure
  ! reichle, 31 Oct 2013 - split "cat_diagn" structure into "cat_diagS" and "cat_diagF"
  ! reichle,  8 Aug 2014 - use Ganymed-4_0 version of catchment.F90 and StieglitzSnow.F90
  ! reichle,  5 Dec 2014 - use Ganymed-4_1 revisions of roughness length
  
  use MAPL_SatVaporMod,                 ONLY:     &
       MAPL_EQsat
  
  use catch_types,                      ONLY:     &
       cat_param_type,                            &
       cat_progn_type,                            &
       cat_diagF_type,                            &
       cat_diagS_type,                            &
       catprogn2wesn,                             &
       assignment (=)
       
  use driver_types,                     ONLY:     &
       met_force_type,                            &
       veg_param_type,                            &
       alb_param_type                           

  use catchment_model,                  ONLY:     &
       catchment,                                 &
       sibalb,                                    &
       catch_calc_soil_moist,                     &
       subtile2tile => catch_calc_subtile2tile
  
  use catch_constants,                  ONLY :    &
       N_snow         => CATCH_N_SNOW,            &
       N_gt           => CATCH_N_GT,              &
       RHOFS          => CATCH_SNWALB_RHOFS,      &
       SNWALB_VISMAX  => CATCH_SNWALB_VISMAX,     &
       SNWALB_NIRMAX  => CATCH_SNWALB_NIRMAX,     &
       SLOPE          => CATCH_SNWALB_SLOPE,      &
       MAPL_NumVegTypes,                          &
       Tzero            => MAPL_TICE,             &
       stefan_boltzmann => MAPL_STFBOL,           &
       rgas             => MAPL_rgas,             &
       grav             => MAPL_GRAV,             &
       cpair            => MAPL_CP,               &
       VIREPS           => MAPL_VIREPS,           &
       MAPL_USMIN,                                &
       MAPL_PI,                                   &
       nodata_generic,                            &
       nodata_tol_generic,                        &
       sfc_turb_scheme,                           &
       GEOS5_coupled_style,                       &
       incl_Louis_ww,                             &
       incl_Louis_extra_derivs

  use StieglitzSnow,                    ONLY:     &
       snowrt,                                    &
       snow_albedo,                               &
       StieglitzSnow_calc_asnow,                  &
       StieglitzSnow_calc_tpsnow
  
  use sfclayer,                         ONLY:     &
       louissurface,                              &
       helfsurface
  
  implicit none
  
  private
  
  public :: propagate_cat
  
  ! parameters for computing zthick, taken from GEOS_CatchGridComp.F90 

  real,    parameter :: MIN_VEG_HEIGHT = 0.01           

  real,    parameter :: SCALE4Z0    = 2.        ! default from Ganymed-4_1
    
  real,    parameter :: Z0_BY_ZVEG  = 0.13 
  real,    parameter :: D0_BY_ZVEG  = 0.66 
  real,    parameter :: HPBL        = 1000.                 

  real,    parameter :: DZE_MIN     = 10.
  
  real,    parameter :: TOTALB_MIN  = 0.01
  real,    parameter :: TOTALB_MAX  = 0.90

  ! parameters for computing land surface emissivity, taken from GEOS_CatchGridComp.F90 
  
  ! Emissivity values from Wilber et al (1999, NATA-TP-1999-209362)
  ! Fu-Liou bands have been combined to Chou bands (though these are broadband only)
  ! IGBP veg types have been mapped to Sib-Mosaic types 
  ! Details in ~suarez/Emiss on cerebus
  
  real,   parameter :: EMSVEG(MAPL_NumVegTypes) = (/  &
       0.99560, 0.99000, 0.99560, 0.99320,            &
       0.99280, 0.99180 /)
  real,   parameter :: EMSBARESOIL   =    0.94120
  real,   parameter :: EMSSNO        =    0.99999
  

  ! definition of land in sfclayer module  
  
  integer, parameter :: SFCTYPE     = 3     

  ! parameters for the Helfand Monin-Obuhkhov surface layer routine
  
  integer, parameter :: HELF_NITER  = 6  ! number of internal iterations
  integer, parameter :: HELF_Z0     = 1  ! CHOOSEZ0 in GEOS_CatchGridComp.F90 (ocean only)
  
contains
  
  subroutine propagate_cat( &
       N_cat, dtstep, tile_id,                                       &  
       cat_param,                                                    &
       met_force, veg_param, sunang, alb_param,                      &  ! ntp values
       cat_progn,                                                    &
       cat_diagS, cat_diagF )
    
    implicit none
    
    integer, intent(in) :: N_cat

    integer, dimension(N_cat), intent(in) :: tile_id

    real, intent(in) :: dtstep             
    
    type(cat_param_type), dimension(N_cat), intent(in) :: cat_param
    
    type(met_force_type), dimension(N_cat), intent(in) :: met_force

    type(veg_param_type), dimension(N_cat), intent(in) :: veg_param

    real,                 dimension(N_cat), intent(in) :: sunang

    type(alb_param_type), dimension(N_cat), intent(in) :: alb_param    
    
    type(cat_progn_type), dimension(N_cat), intent(inout) :: cat_progn
    
    type(cat_diagS_type), dimension(N_cat), intent(out)   :: cat_diagS

    type(cat_diagF_type), dimension(N_cat), intent(out)   :: cat_diagF
    
    ! --------------------------------------------------------------

    real :: SFRAC

    real, dimension(N_cat) :: &
         SQSCAT, Z2, DZM, RSOIL1, RSOIL2, SATCAP, RDC, U2FAC,            &
         TRAINL, PARDIR, PARDIF,                                         &
         SWNETF, SWNETS,                                                 &
         RA1, ETURB1, DEDTC1, HSTURB1, DHSDTC1,                          &
         RA2, ETURB2, DEDTC2, HSTURB2, DHSDTC2,                          &
         RA4, ETURB4, DEDTC4, HSTURB4, DHSDTC4,                          &
         RAS, ETURBS, DEDTCS, HSTURBS, DHSDTCS,                          &
         QSAT1, DQS1, ALW1, BLW1,                                        &
         QSAT2, DQS2, ALW2, BLW2,                                        &
         QSAT4, DQS4, ALW4, BLW4,                                        &
         QSATS, DQSS, ALWS, BLWS,                                        &
         DEDQA1, DHSDQA1,                                                &
         DEDQA2, DHSDQA2,                                                &
         DEDQA4, DHSDQA4,                                                &
         DEDQAS, DHSDQAS,                                                &
         ALBVR, ALBNR, ALBVF, ALBNF,                                     & 
         SNOVR, SNONR, SNOVF, SNONF,                                     &
         SWLAND, RZEQ,                                                   &
         GHFLXSNO, GHFLXTSKIN,                                           &
         ENTOT, WTOT, WCHANGE, ECHANGE,                                  &
         SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1,                 &
         LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW,                       &
         TCSORIG1, TPSN1IN1, TPSN1OUT1,                                  &
         TC1_0, TC2_0, TC4_0, QA1_0, QA2_0, QA4_0, EACC_0   
                         

    real, dimension(N_cat) :: zol, zol_sc, d0, dze, tmprealvec, uuu, zvg, ar4

    real, dimension(:),   allocatable :: ch_tile, cq_tile, qa_tile
    real, dimension(:,:), allocatable :: qa
    
    real, dimension(N_gt,  N_cat) :: GHTCNT
    real, dimension(N_snow,N_cat) :: WESN, HTSN, SNDZ
    
    real :: albedo, albsn, asnow, totalb
    real :: qsattc, tmpt, tmpswnet, emis
    
    integer :: k, n

    logical :: snow_dominates_alb
    
    character(len=*), parameter :: Iam = 'propagate_cat'
    character(len=400) :: err_msg

    ! ----------------------------------------------------------------
    
    ! initialize cat_diagS and cat_diagF fields that are needed here 
    !  and in sfc_turb 
    
    do n=1,N_cat
       
       cat_diagS(n) = nodata_generic
       cat_diagF(n) = nodata_generic
       
    end do
    
    ! diagnose top layer snow temperature in K
    
    call StieglitzSnow_calc_tpsnow( N_cat, cat_progn%htsn(1), cat_progn%wesn(1),    &
         cat_diagS%tpsn(1), tmprealvec )
    
    cat_diagS%tpsn(1) = cat_diagS%tpsn(1) + Tzero   ! convert to Kelvin
    
    ! diagnose sub-tile area fractions

    ! updated to new interfaces - reichle, 3 Apr 2012
    
    call catch_calc_soil_moist(                                     &
         N_cat,                                                     &
         cat_param%vegcls, cat_param%dzsf,  cat_param%vgwmax,       &
         cat_param%cdcr1,  cat_param%cdcr2, cat_param%psis,         &
         cat_param%bee,    cat_param%poros, cat_param%wpwet,        &
         cat_param%ars1,   cat_param%ars2,  cat_param%ars3,         &
         cat_param%ara1,   cat_param%ara2,  cat_param%ara3,         &
         cat_param%ara4,   cat_param%arw1,  cat_param%arw2,         &
         cat_param%arw3,   cat_param%arw4,                          &
         cat_progn%srfexc, cat_progn%rzexc, cat_progn%catdef,       &
         cat_diagS%ar1,    cat_diagS%ar2,   ar4                )
        
    call StieglitzSnow_calc_asnow( N_snow, N_cat, catprogn2wesn(N_cat,cat_progn),   &
         cat_diagS%asnow )
    
    ! ---------------------------------------------------------------

    ! assorted Catchment model inputs
    
    call pmonth(                                                 &
         N_cat, cat_param%vegcls, veg_param%grn, veg_param%lai,  &
         zol, d0, SQSCAT, Z2, zvg, RSOIL1, RSOIL2, RDC )
    
    ! ---------------------------------------------------------------
    
    ! precipitation
    
    TRAINL = met_force%Rainf - met_force%Rainf_C
    
    SFRAC = 1.
    
    ! ---------------------------------------------------------------
    
    ! turbulence

    do n=1,N_cat
       
       ! Ganymed-4_1 revision of roughness length.  Intentionally applies only to 
       ! roughness that is input into surface turbulence calculations, and not to 
       ! wind-related input into subroutine catchment(), see below.
       !
       ! NOTE: Changes in this subroutine (propagate_cat()) to accommodate the
       !       clean-up of subroutine pmonth() result in roundoff differences
       !       w.r.t. previous LDASsa version even if SCALE4Z0=1.
       !
       ! - reichle, 8 Dec 2014
       
       zol_sc(n) = zol(n)*SCALE4Z0                

       dze(n) = max(met_force(n)%RefH-d0(n), DZE_MIN)   
              
    end do
    
    ! compute surface exchange coefficients etc BEFORE possibly resetting
    ! profile of Qair-QAx-Qsat(surf) -- for consistency with two-stage 
    ! run-method in GEOS_CatchGridComp.F90
    ! reichle+qliu,  9 Oct 2008

    select case (sfc_turb_scheme)
       
    case (0)     ! Louis

       allocate(ch_tile(N_cat))
       allocate(cq_tile(N_cat))
       
       call sfc_turb(N_cat, met_force, cat_progn, cat_diagS, zol_sc, dze, veg_param%lai, &
            RA1, ETURB1, DEDQA1, DEDTC1, HSTURB1, DHSDQA1, DHSDTC1,           &
            RA2, ETURB2, DEDQA2, DEDTC2, HSTURB2, DHSDQA2, DHSDTC2,           &
            RA4, ETURB4, DEDQA4, DEDTC4, HSTURB4, DHSDQA4, DHSDTC4,           &
            RAS, ETURBS, DEDQAS, DEDTCS, HSTURBS, DHSDQAS, DHSDTCS,           & 
            ch_tile=ch_tile, cq_tile=cq_tile )

    case (1)     ! Helfand Monin-Obukhov
       
       call sfc_turb(N_cat, met_force, cat_progn, cat_diagS, zol_sc, dze, veg_param%lai, &
            RA1, ETURB1, DEDQA1, DEDTC1, HSTURB1, DHSDQA1, DHSDTC1,           &
            RA2, ETURB2, DEDQA2, DEDTC2, HSTURB2, DHSDQA2, DHSDTC2,           &
            RA4, ETURB4, DEDQA4, DEDTC4, HSTURB4, DHSDQA4, DHSDTC4,           &
            RAS, ETURBS, DEDQAS, DEDTCS, HSTURBS, DHSDQAS, DHSDTCS,           &
            t2m_tile=cat_diagF%t2m, q2m_tile=cat_diagF%q2m  )
       
    case default
       
       err_msg = 'Unknown selection for sfc_turb_scheme.'
       stop
       
    end select

       

    do n=1,N_cat

       ! make sure that Qair-QAx-Qsat(surf) profile is ok
       ! (in a nutshell, need Qair<=QAx<=Qsat(surf) or Qair>=QAx>=Qsat(surf))
       ! - reichle, 5 Mar 2008 
       ! - qliu+reichle, 12 Aug 2008
       
       !qsattc = qsat(cat_progn(n)%tc1,PSUR(n),ALHE)
       qsattc = MAPL_EQsat(cat_progn(n)%tc1,met_force(n)%Psurf)
       if (cat_progn(n)%qa1>qsattc .and. cat_progn(n)%qa1>met_force(n)%Qair ) then
          cat_progn(n)%qa1 = max( met_force(n)%Qair, qsattc )
       endif
       if (cat_progn(n)%qa1<qsattc .and. cat_progn(n)%qa1<met_force(n)%Qair ) then
          cat_progn(n)%qa1 = min( met_force(n)%Qair, qsattc )
       endif

       !qsattc = qsat(cat_progn(n)%tc2,PSUR(n),ALHE)
       qsattc = MAPL_EQsat(cat_progn(n)%tc2,met_force(n)%Psurf)
       if (cat_progn(n)%qa2>qsattc .and. cat_progn(n)%qa2>met_force(n)%Qair ) then
          cat_progn(n)%qa2 = max( met_force(n)%Qair, qsattc )
       endif
       if (cat_progn(n)%qa2<qsattc .and. cat_progn(n)%qa2<met_force(n)%Qair ) then
          cat_progn(n)%qa2 = min( met_force(n)%Qair, qsattc )
       endif

       !qsattc = qsat(cat_progn(n)%tc4,PSUR(n),ALHE)
       qsattc = MAPL_EQsat(cat_progn(n)%tc4,met_force(n)%Psurf)
       if (cat_progn(n)%qa4>qsattc .and. cat_progn(n)%qa4>met_force(n)%Qair ) then
          cat_progn(n)%qa4 = max( met_force(n)%Qair, qsattc )
       endif
       if (cat_progn(n)%qa4<qsattc .and. cat_progn(n)%qa4<met_force(n)%Qair ) then
          cat_progn(n)%qa4 = min( met_force(n)%Qair, qsattc )
       endif
       
    end do
    
    ! convert ght and snow structures into arrays needed by sibalb() and catchment()
    
    do k=1,N_gt
       
       GHTCNT(k,:) = cat_progn%ght(k)
       
    end do
    
    do k=1,N_snow
       
       WESN(k,:) = cat_progn%wesn(k)
       HTSN(k,:) = cat_progn%htsn(k)
       SNDZ(k,:) = cat_progn%sndz(k)
       
    end do

    ! radiation
    
    call SIBALB(N_cat, cat_param%vegcls,                       &
         veg_param%lai, veg_param%grn, sunang,                 & 
         alb_param%sc_albvf, alb_param%sc_albvf,               & ! use only diffuse
         alb_param%sc_albnf, alb_param%sc_albnf,               & ! use only diffuse
         ALBVR, ALBNR, ALBVF, ALBNF                )             ! inst snow-free albedos
    
    call SNOW_ALBEDO(N_cat, N_snow, 0, cat_param%vegcls,       &
         veg_param%lai, sunang,                                &
         RHOFS,                                                &   
         SNWALB_VISMAX, SNWALB_NIRMAX, SLOPE,                  & 
         WESN, HTSN, SNDZ,                                     & 
         ALBVR, ALBNR, ALBVF, ALBNF,                           & ! inst snow-free albedos
         SNOVR, SNONR, SNOVF, SNONF )                            ! inst snow albedos 
        
    do n=1,N_cat
       
       ! tentative albedo for snow-free and snow-covered portion
       ! (naively assumes that each SW radiation component contributes 0.25 to 
       !  total SWdown -- these tentative albedos will NOT be used as such
       !  if SWnet is provided in forcing data)
       ! reichle+qliu,  9 Oct 2008

       albedo = min(1., 0.25*(ALBVR(n) + ALBNR(n) + ALBVF(n) + ALBNF(n)))
       albsn =  min(1., 0.25*(SNOVR(n) + SNONR(n) + SNOVF(n) + SNONF(n)))
       
       !SWNETF(n) = (1.-albedo) * met_force(n)%SWdown
       !SWNETS(n) = (1.-albsn ) * met_force(n)%SWdown
       
       ! diagnose total albedo and upward shortwave
       ! (removed duplicate call to calc_asnow - reichle, 3 Apr 2012)
       
       asnow = cat_diagS(n)%asnow    ! will also be needed for emissivity calcs below
       
       ! total albedo calculation requires a threshold "SWDN_THRESHOLD" for minimum 
       ! allowable positive SWdown as implemented in subroutine repair_forcing()
       ! qliu+reichle, 21 Aug 2008
       
       ! improved treatment of albedo and SWnet calculations that take advantage
       ! of SWnet forcing to obviate the need for individual SW radiation components
       ! reichle+qliu,  9 Oct 2008
       
       if (met_force(n)%SWdown>0. .and. sunang(n)>0.) then
          
          ! sun is up and SWdown>0

          ! determine whether snow or snow-free albedo dominates

          snow_dominates_alb = ( (1-asnow)*(1-albedo)<=asnow*(1-albsn) )
          
          ! compute SWNET for snow-free and snow conditions
          
          if ( (abs(met_force(n)%SWnet-nodata_generic)<nodata_tol_generic) &
               .or.(snow_dominates_alb) )  then
             
             ! SWnet is not available or wanted, use simplistic albedos (see above)
             
             SWNETF(n) = (1.-albedo) * met_force(n)%SWdown
             SWNETS(n) = (1.-albsn ) * met_force(n)%SWdown
             
             tmpswnet = (1-asnow)*SWNETF(n) + asnow* SWNETS(n)
             
          else
             
             ! use knowlegdge of SWnet to back out better estimate of albedo
             ! (works well if *all* tiles in forcing grid cell are either
             ! completely snow-covered or completely snow-free)
             
             if (.not. snow_dominates_alb) then
                
                ! snow-free area dominates, use simplistic snow albedo
                ! and back out better snow-free albedo
                
                SWNETS(n) = (1.-albsn ) * met_force(n)%SWdown
                SWNETF(n) = (met_force(n)%SWnet-asnow*SWNETS(n))/(1.-asnow)
                
             else
                
                ! snow-covered area dominates, use simplistic snow-free albedo
                ! and back out better snow-covered albedo                
                
                SWNETF(n) = (1.-albedo ) * met_force(n)%SWdown
                SWNETS(n) = (met_force(n)%SWnet-(1.-asnow)*SWNETF(n))/asnow
                
             end if
             
             tmpswnet = met_force(n)%SWnet
             
          end if
          
          totalb = min(max(1. - tmpswnet/met_force(n)%SWdown, TOTALB_MIN), TOTALB_MAX)
          
          cat_diagF(n)%swup   = totalb*met_force(n)%SWdown
          
       else
          
          SWNETS(n)           = 0.
          SWNETF(n)           = 0.
          
          cat_diagF(n)%swup   = 0.

       end if
       
                     
       ! photo-synthetically active radiation
              
       if(abs(met_force(n)%PARdrct-nodata_generic)<nodata_tol_generic) then
          
          ! if PAR components are not available:
          !
          !   assume half of SWdown is photosynthetically active
          !   assume half of PAR is direct, half diffuse
       
          PARDIR(n) = 0.5*0.5*met_force(n)%SWdown          
          PARDIF(n) = PARDIR(n)
       
       else

          ! assume that PARdffs is available whenever PARdrct is

          PARDIR(n) = met_force(n)%PARdrct
          PARDIF(n) = met_force(n)%PARdffs
          
       end if
       
       ! compute QSATx, DQSx
       
       QSAT1(n)   = MAPL_EQsat(cat_progn(n)%tc1, met_force(n)%Psurf, DQS1(n) )       
       QSAT2(n)   = MAPL_EQsat(cat_progn(n)%tc2, met_force(n)%Psurf, DQS2(n) )
       QSAT4(n)   = MAPL_EQsat(cat_progn(n)%tc4, met_force(n)%Psurf, DQS4(n) )
       
       QSATS(n)   = MAPL_EQsat(cat_diagS(n)%tpsn(1), met_force(n)%Psurf, DQSS(n), &
            OverIce=.true. )
       

       ! compute ALWx, BLWx
       
       ! tile-average surface emissivity 
       
       emis = EMSVEG(cat_param(n)%vegcls) + &
            (EMSBARESOIL - EMSVEG(cat_param(n)%vegcls))*exp(-veg_param(n)%lai)
       
       emis = emis*(1.-asnow) + EMSSNO*asnow
       
       ! ALWx and BLWx are computed separately for each sub-tile 
       ! (unlike in GEOScatch_GridComp.F90)
       
       tmpt    = cat_progn(n)%tc1
       BLW1(n) = emis*stefan_boltzmann*tmpt*tmpt*tmpt
       ALW1(n) = -3.*tmpt*BLW1(n)
       BLW1(n) =  4.*BLW1(n)
       
       tmpt    = cat_progn(n)%tc2
       BLW2(n) = emis*stefan_boltzmann*tmpt*tmpt*tmpt
       ALW2(n) = -3.*tmpt*BLW2(n)
       BLW2(n) =  4.*BLW2(n)
       
       tmpt    = cat_progn(n)%tc4
       BLW4(n) = emis*stefan_boltzmann*tmpt*tmpt*tmpt
       ALW4(n) = -3.*tmpt*BLW4(n)
       BLW4(n) =  4.*BLW4(n)
       
       tmpt    = cat_diagS(n)%tpsn(1)
       BLWS(n) = emis*stefan_boltzmann*tmpt*tmpt*tmpt
       ALWS(n) = -3.*tmpt*BLWS(n)
       BLWS(n) =  4.*BLWS(n)
       
    end do
    
    
    ! get wind-related variable for input to catchment()

    ! Ganymed-4_1 intentionally uses *unscaled* zol here (per Andrea Molod).
    ! - reichle, 5 Dec 2014

    uuu = max(met_force%Wind, MAPL_USMIN) *(log((max(zvg-d0,0.)+zol)/zol) &
         / log((max(met_force%RefH-d0,10.)+zol)/zol))
    
    
    ! call to Catchment model 
    
    call catchment( N_cat, dtstep, sfrac,                                &
         tile_id,cat_param%vegcls,                                       &
         cat_param%dzsf, met_force%Rainf_C, TRAINL, met_force%Snowf, uuu,&
         ETURB1, DEDQA1, DEDTC1, HSTURB1, DHSDQA1, DHSDTC1,              &
         ETURB2, DEDQA2, DEDTC2, HSTURB2, DHSDQA2, DHSDTC2,              &
         ETURB4, DEDQA4, DEDTC4, HSTURB4, DHSDQA4, DHSDTC4,              &
         ETURBS, DEDQAS, DEDTCS, HSTURBS, DHSDQAS, DHSDTCS,              &
         met_force%Tair, met_force%Qair,                                 &
         RA1, RA2, RA4, RAS, sunang, PARDIR, PARDIF,                     &
         SWNETF, SWNETS, met_force%LWdown, met_force%Psurf/100.,         &
         veg_param%lai, veg_param%grn, Z2,                               &
         SQSCAT, RSOIL1, RSOIL2, RDC,                                    &
         QSAT1, DQS1, ALW1, BLW1,                                        &
         QSAT2, DQS2, ALW2, BLW2,                                        &
         QSAT4, DQS4, ALW4, BLW4,                                        &
         QSATS, DQSS, ALWS, BLWS,                                        &
         cat_param%bf1, cat_param%bf2, cat_param%bf3, cat_param%vgwmax,  &
         cat_param%cdcr1, cat_param%cdcr2,                               &
         cat_param%psis, cat_param%bee, cat_param%poros,                 &
         cat_param%wpwet, cat_param%cond, cat_param%gnu,                 &
         cat_param%ars1, cat_param%ars2, cat_param%ars3,                 &
         cat_param%ara1, cat_param%ara2, cat_param%ara3, cat_param%ara4, &
         cat_param%arw1, cat_param%arw2, cat_param%arw3, cat_param%arw4, &
         cat_param%tsa1, cat_param%tsa2, cat_param%tsb1, cat_param%tsb2, &
         cat_param%atau, cat_param%btau,                                 &
         .false.,                                                        &
         cat_progn%tc1, cat_progn%tc2, cat_progn%tc4,                    &
         cat_progn%qa1, cat_progn%qa2, cat_progn%qa4,                    &
         cat_progn%capac,                                                &
         cat_progn%catdef, cat_progn%rzexc, cat_progn%srfexc,            &
         GHTCNT, cat_diagS%tsurf,                                        &
         WESN, HTSN, SNDZ,                                               &
         cat_diagF%evap, cat_diagF%shflux, cat_diagF%runoff,             &
         cat_diagF%eint, cat_diagF%esoi, cat_diagF%eveg, cat_diagF%esno, &
         cat_diagF%bflow, cat_diagF%runsrf, cat_diagF%snmelt,            &
         cat_diagF%lwup, SWLAND, cat_diagF%lhflux, cat_diagF%qinfil,     &
         cat_diagS%ar1, cat_diagS%ar2, RZEQ,                             &
         cat_diagF%ghflux, GHFLXSNO, GHFLXTSKIN,                         &
         cat_diagS%tpsn(1), cat_diagS%asnow,                             &
         cat_diagS%tp(1), cat_diagS%tp(2), cat_diagS%tp(3),              &
         cat_diagS%tp(4), cat_diagS%tp(5), cat_diagS%tp(6),              &
         cat_diagS%sfmc, cat_diagS%rzmc, cat_diagS%prmc,                 &
         ENTOT, WTOT, WCHANGE, ECHANGE,                                  &
         cat_diagF%hsnacc, cat_diagF%evacc, cat_diagF%shacc,             &
         SHSNOW1, AVETSNOW1, WAT10CM1, WATSOI1, ICESOI1,                 &
         LHSNOW1, LWUPSNOW1, LWDNSNOW1, NETSWSNOW,                       &
         TCSORIG1, TPSN1IN1, TPSN1OUT1,                                  &
         LHACC=cat_diagF%lhacc,                                          &
         TC1_0=TC1_0, TC2_0=TC2_0, TC4_0=TC4_0,                          &
         QA1_0=QA1_0, QA2_0=QA2_0, QA4_0=QA4_0,                          &
         EACC_0=EACC_0                               )   
    
    
    ! deal with AGCM fix at the end of subroutine catchment()
    !
    ! reichle+qliu, 9 Oct 2008
    
    if (GEOS5_coupled_style) then
       
       ! keep "modified" tc's and qa's, also modify fluxes
       
       cat_diagF%evap   = cat_diagF%evap   - cat_diagF%evacc
       cat_diagF%lhflux = cat_diagF%lhflux - cat_diagF%lhacc
       cat_diagF%shflux = cat_diagF%shflux - cat_diagF%shacc              
       cat_diagF%eacc_0 = 0.

    else
       
       ! revert to "unmodified" tc's and qa's
       ! note that almost all diagnostic fluxes (including
       ! "evap", "shflux", "lhflux", "echange", "wchange")
       ! EXCEPT "evacc", "shacc", and "lhacc" are based on these 
       ! unmodified tc's and qa's
       
       cat_progn%tc1    = TC1_0
       cat_progn%tc2    = TC2_0
       cat_progn%tc4    = TC4_0
       cat_progn%qa1    = QA1_0
       cat_progn%qa2    = QA2_0
       cat_progn%qa4    = QA4_0
       cat_diagF%eacc_0 = EACC_0
       
    end if
       

    ! convert ght and snow arrays into structures needed by LDASsa
    
    do k=1,N_gt
       
       cat_progn%ght(k)=ghtcnt(k,:)
       
    end do
    
    do k=1,N_snow
       
       cat_progn%wesn(k) = WESN(k,:)
       cat_progn%htsn(k) = HTSN(k,:)
       cat_progn%sndz(k) = SNDZ(k,:)
       
    end do
    

    ! diagnose snow temperatures  in layer 2:N_snow 
    ! ( tpsn1 already provided by catchment() )
    
    do k=2,N_snow
       
       call StieglitzSnow_calc_tpsnow( N_cat, cat_progn%htsn(k), cat_progn%wesn(k), &
            cat_diagS%tpsn(k), tmprealvec )
       
       cat_diagS%tpsn(k) = cat_diagS%tpsn(k) + Tzero   ! convert to Kelvin
    
    end do
    

    ! diagnose T2m, Q2m for Louis scheme
    
    if (sfc_turb_scheme==0) then
       
       tmprealvec = log(1.+ 2./zol_sc) / log(1. + dze/zol_sc)

       ! reichle, 8 Feb 2011:
       !
       ! Not entirely sure about the following...
       !
       ! "t2m" from GEOS_SurfaceGridComp.F90:  
       !
       !   TTM = TH - (SH/MAPL_CP + DSH  *DTS)*TMP/CT
       !
       ! where DTS = grid-avg( TC(:,N)-THATM ), which can be guessed
       ! from GEOS_CatchGridComp.F90:
       !
       !   SHSBT (:,N) = (SH  + DSH  *(TC(:,N)-THATM))*CFT(:,N)
       !
       ! Assuming that in (SurfaceGridComp) grid space DSH*DTS is negligible,
       ! we have: 
       
       cat_diagF%t2m = cat_diagS%tsurf - tmprealvec*cat_diagF%shflux/ch_tile/cpair 

       ! Note that the "t2m" computed here approximates the t2m over *land*
       ! (as opposed to the T2m in GEOS_SurfaceGridComp, which is over the
       !  entire grid cell, including lake, permanent ice, and ocean if present)

       ! Similarly for q2m:

       allocate(qa(N_cat,4))

       qa(:,1) = cat_progn%qa1
       qa(:,2) = cat_progn%qa2
       qa(:,3) = cat_progn%qa4
       
       do n=1,N_cat
          qa(n,4) = MAPL_EQsat(cat_diagS(n)%tpsn(1),met_force(n)%Psurf,OverIce=.true.)
       end do

       allocate(qa_tile(N_cat))

       call subtile2tile(N_cat,cat_diagS%ar1,cat_diagS%ar2,cat_diagS%asnow,qa,qa_tile)
       
       cat_diagF%q2m = qa_tile - tmprealvec*cat_diagF%evap/cq_tile
       
       deallocate(qa)

       deallocate(ch_tile)
       deallocate(cq_tile)
       deallocate(qa_tile) 

    end if

  end subroutine propagate_cat
    
  ! *****************************************************************************
  
  subroutine sfc_turb(N_cat, met_force, cat_progn, cat_diagS, zol, dze, lai,  &
       RA1, ETURB1, DEDQA1, DEDTC1, HSTURB1, DHSDQA1, DHSDTC1,                &
       RA2, ETURB2, DEDQA2, DEDTC2, HSTURB2, DHSDQA2, DHSDTC2,                &
       RA4, ETURB4, DEDQA4, DEDTC4, HSTURB4, DHSDQA4, DHSDTC4,                &
       RAS, ETURBS, DEDQAS, DEDTCS, HSTURBS, DHSDQAS, DHSDTCS,                & 
       t2m_tile, q2m_tile, ch_tile, cq_tile )
    

    ! qliu+reichle: 14 Aug 2008 - land interface to subroutine louissurface() 
    ! qliu+reichle: 30 Oct 2008 - force derivatives to be non-negative in case
    !                             extra derivatives are included
    ! reichle     : 29 Nov 2010 - added Helfand Monin-Obukhov scheme
    ! reichle     :  2 Feb 2011 - added T2m, Q2m diagnostics

    implicit none
    
    integer,  intent(in) :: N_cat
    
    type(met_force_type), dimension(N_cat), intent(in) :: met_force
    type(cat_progn_type), dimension(N_cat), intent(in) :: cat_progn
    type(cat_diagS_type), dimension(N_cat), intent(in) :: cat_diagS
    
    real, dimension(N_cat), intent(in) :: zol, dze, lai
    
    real, dimension(N_cat), intent(out) :: RA1, RA2, RA4, RAS
    real, dimension(N_cat), intent(out) :: ETURB1, DEDQA1, DEDTC1
    real, dimension(N_cat), intent(out) :: ETURB2, DEDQA2, DEDTC2
    real, dimension(N_cat), intent(out) :: ETURB4, DEDQA4, DEDTC4
    real, dimension(N_cat), intent(out) :: ETURBS, DEDQAS, DEDTCS
    real, dimension(N_cat), intent(out) :: HSTURB1,DHSDQA1,DHSDTC1
    real, dimension(N_cat), intent(out) :: HSTURB2,DHSDQA2,DHSDTC2
    real, dimension(N_cat), intent(out) :: HSTURB4,DHSDQA4,DHSDTC4
    real, dimension(N_cat), intent(out) :: HSTURBS,DHSDQAS,DHSDTCS
    
    real, dimension(N_cat), intent(out), optional :: t2m_tile, q2m_tile   ! Helfand MO
    
    real, dimension(N_cat), intent(out), optional :: ch_tile, cq_tile     ! Louis
    
    
    ! local variables

    integer :: k, n
    
    real    :: rhoair

    real, dimension(N_cat,4) :: tc, qa, z0, ch, cq
    
    real, dimension(      4) :: dt, dq, dcqdq, dcqdt, dchdt, dchdq
    
    ! Louis

    real,    dimension(:),   allocatable :: cn, ri, zt, zq, uuu, ucn, re
    real,    dimension(:,:), allocatable :: cm, ww, dch, dcq 
    
    ! Helfand Monin-Obukhov

    integer, dimension(:),   allocatable :: tmpintvec
    
    real,    dimension(:),   allocatable :: tmprealvec, psmb, psl
    real,    dimension(:),   allocatable :: RHOH,VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS

    real,    dimension(:,:), allocatable :: t2m,q2m
    real,    dimension(:),   allocatable :: u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m
    
    character(len=*), parameter :: Iam = 'sfc_turb'
    character(len=400) :: err_msg

    ! ------------------------------------------------------------------------
    !
    ! initialize optional output variables
    
    if (present(ch_tile )) ch_tile  = nodata_generic
    if (present(cq_tile )) cq_tile  = nodata_generic
    if (present(t2m_tile)) t2m_tile = nodata_generic
    if (present(q2m_tile)) q2m_tile = nodata_generic


    ! prepare inputs that are common to all sfc_turb_scheme choices
    
    do k=1,N_cat
       
       tc(k,1) = cat_progn(k)%tc1
       tc(k,2) = cat_progn(k)%tc2
       tc(k,3) = cat_progn(k)%tc4
       tc(k,4) = cat_diagS(k)%tpsn(1)

       qa(k,1) = cat_progn(k)%qa1
       qa(k,2) = cat_progn(k)%qa2
       qa(k,3) = cat_progn(k)%qa4
       qa(k,4) = MAPL_EQsat(cat_diagS(k)%tpsn(1),met_force(k)%Psurf,OverIce=.true.)
       
       z0(k,1:4) = zol(k)
       
    end do

    ! allocate and initialize variables as needed
    
    select case (sfc_turb_scheme)
       
    case (0)     ! Louis
       
       allocate(cn( N_cat))
       allocate(ri( N_cat))
       allocate(zt( N_cat))
       allocate(zq( N_cat))
       allocate(uuu(N_cat))
       allocate(ucn(N_cat))
       allocate(re( N_cat))
       
       allocate(cm( N_cat,4))
       allocate(ww( N_cat,4))
       allocate(dch(N_cat,4))
       allocate(dcq(N_cat,4))

       ! note that cm and z0 are intent(inout) in louissurface()

       cm = 0.
       ww = 0.
       
    case (1)     ! Helfand Monin-Obukhov
       
       allocate(tmpintvec( N_cat))       

       allocate(tmprealvec(N_cat))
       allocate(psmb(      N_cat)) 
       allocate(psl(       N_cat)) 
       allocate(RHOH(      N_cat))
       allocate(VKM(       N_cat))
       allocate(USTAR(     N_cat))
       allocate(XX(        N_cat))
       allocate(YY(        N_cat))
       allocate(CU(        N_cat))
       allocate(CT(        N_cat))
       allocate(RIB(       N_cat))
       allocate(ZETA(      N_cat))
       allocate(WS(        N_cat))
                           
       allocate(t2m(       N_cat,4))
       allocate(q2m(       N_cat,4))

       allocate(u2m(       N_cat))
       allocate(v2m(       N_cat))
       allocate(t10m(      N_cat))
       allocate(q10m(      N_cat))
       allocate(u10m(      N_cat))
       allocate(v10m(      N_cat))
       allocate(u50m(      N_cat))
       allocate(v50m(      N_cat))
              
    case default
       
       err_msg = 'Unknown selection for sfc_turb_scheme.'
       stop
       
    end select
    
    
    ! ----------------------------------------------------------------------------
    !
    ! call surface layer subroutine
    
    do n=1,4
       
       select case (sfc_turb_scheme)
          
          ! ----------------------------------------------------------------------
          
       case (0)     ! Louis
          
          ! note that cm and z0 are intent(inout) in louissurface()
          
          call louissurface(SFCTYPE, n, met_force%Wind, ww, met_force%Psurf,    &
               met_force%Tair, tc, met_force%Qair, qa, met_force%Rainf_C,       &
               lai, z0, dze,                                                    &
               cm, cn, ri, zt, zq, ch, cq, uuu, ucn, re, dch, dcq )
          
          if (incl_Louis_ww) then
             
             ! simplistic estimate of ww for off-line integrations
             ! (for a better job and better consistency with coupled integrations,
             !  ww would probably have to be included in cat_progn)
             ! reichle+qliu,  9 Oct 2008
             
             ww(:,n) = ch(:,n)*(tc(:,n)-met_force(:)%Tair-                      & 
                  (GRAV/CPAIR)*dze(:))/met_force(:)%Tair + VIREPS*cq(:,n)*      &
                  (qa(:,n)-met_force(:)%Qair)
             ww(:,n) = max(ww(:,n),0.0)
             ww(:,n) = (HPBL*GRAV*ww(:,n))**(2./3.)  
             
             call louissurface(SFCTYPE, n, met_force%Wind, ww, met_force%Psurf, &
                  met_force%Tair, tc, met_force%Qair, qa, met_force%Rainf_C,    &
                  lai, z0, dze,                                                 &
                  cm, cn, ri, zt, zq, ch, cq, uuu, ucn, re, dch, dcq )
             
          end if

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

       case (1)     ! Helfand Monin-Obukhov
          
          psmb = met_force%Psurf * 0.01            ! convert to millibar
          
          ! Approximate pressure at top of surface layer: 
          ! hydrostatic, eqn of state using avg temp and press
          
          tmprealvec = (dze*GRAV) / (RGAS*(met_force%Tair + tc(:,n)))
          
          psl = psmb * (1. - tmprealvec) / (1. + tmprealvec)
          
          ! surface type
          
          tmpintvec(1:N_cat) = SFCTYPE
          
          ! for wind use u-component=wind-speed, v-component=0.
          
          tmprealvec = 0. 

          call helfsurface( met_force%Wind, tmprealvec, met_force%Tair, tc(:,n), &
               met_force%Qair, qa(:,n), psl, psmb, z0(:,n), lai,                 &
               tmpintvec, dze, HELF_NITER, N_cat,                                &
               RHOH,ch(:,n),VKM,USTAR,XX,YY,CU,CT,RIB,ZETA,WS,                   &
               t2m(:,n),q2m(:,n),u2m,v2m,t10m,q10m,u10m,v10m,u50m,v50m,HELF_Z0)
          
          cq(:,n)  = ch(:,n) 

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

       end select
       
    end do
    
    ! -------------------------------------------------------------------------------
    !
    ! finalize output arguments
    
    do k=1,N_cat
       
       ! rhoair from GEOS_CatchGridComp.F90
       
       rhoair     = RGAS*met_force(k)%Tair*(1.+VIREPS*met_force(k)%Qair)
       rhoair     = met_force(k)%Psurf/rhoair
       
       RA1(k)     = rhoair/ch(k,1)
       RA2(k)     = rhoair/ch(k,2)
       RA4(k)     = rhoair/ch(k,3)
       RAS(k)     = rhoair/ch(k,4)
       
       dq         = qa(k,:) - met_force(k)%Qair
       
       ETURB1(k)  = cq(k,1)*dq(1)
       ETURB2(k)  = cq(k,2)*dq(2)
       ETURB4(k)  = cq(k,3)*dq(3)
       ETURBS(k)  = cq(k,4)*dq(4)
       
       dt         = tc(k,:) - met_force(k)%Tair
       
       HSTURB1(k) = CPAIR*ch(k,1)*dt(1)
       HSTURB2(k) = CPAIR*ch(k,2)*dt(2)
       HSTURB4(k) = CPAIR*ch(k,3)*dt(3)
       HSTURBS(k) = CPAIR*ch(k,4)*dt(4)
       
       ! derivatives
       
       DEDQA1(k)  = cq(k,1)
       DEDQA2(k)  = cq(k,2)
       DEDQA4(k)  = cq(k,3)
       DEDQAS(k)  = cq(k,4)
       
       DHSDQA1(k) = 0.
       DHSDQA2(k) = 0.
       DHSDQA4(k) = 0.
       DHSDQAS(k) = 0.
       
       DHSDTC1(k) = CPAIR*ch(k,1)
       DHSDTC2(k) = CPAIR*ch(k,2)
       DHSDTC4(k) = CPAIR*ch(k,3)
       DHSDTCS(k) = CPAIR*ch(k,4)
       
       DEDTC1(k)  = 0.
       DEDTC2(k)  = 0.
       DEDTC4(k)  = 0.
       DEDTCS(k)  = 0.
       

       if ((sfc_turb_scheme==0) .and. incl_Louis_extra_derivs) then

          ! overwrite derivatives with cross-terms
          
          ! include derivatives of drag coefficients and cross-derivatives
          ! (evap w.r.t. temperature, sensible w.r.t. humidity)
          
          ! "dcq" from louissurface() is d(cq)/d(TVA) where TVA = T*(1+VIREPS*Q)
          !
          ! moreover, cq=cq(Ri) where Ri is proportional to deltaTVA=TVA-TVS
          ! (where TVS is the virtual surface temperature) -> this produces a 
          ! minus sign
          !
          ! taken together, we need d(cq)/dq = d(cq)/dTV * dTV/dQ = -dcq * VIREPS*T
          
          dcqdq      = -dcq(k,:)*VIREPS*tc(k,:)
          
          DEDQA1(k)  = cq(k,1) + max( 0., dcqdq(1)*dq(1) )
          DEDQA2(k)  = cq(k,2) + max( 0., dcqdq(2)*dq(2) )
          DEDQA4(k)  = cq(k,3) + max( 0., dcqdq(3)*dq(3) )
          DEDQAS(k)  = cq(k,4) + max( 0., dcqdq(4)*dq(4) )

          ! similarly, d(cq)/dt = d(cq)/dTV * dTV/dT = -dcq * (1+VIREPS*Q)
          
          dcqdt      = -dcq(k,:)*(1.+VIREPS*qa(k,:))
          
          DEDTC1(k)  = max( 0., dcqdt(1)*dq(1) )
          DEDTC2(k)  = max( 0., dcqdt(2)*dq(2) ) 
          DEDTC4(k)  = max( 0., dcqdt(3)*dq(3) )
          DEDTCS(k)  = max( 0., dcqdt(4)*dq(4) )
          
          ! "dch" from louissurface() is d(ch)/d(TVA) where TVA = T*(1+VIREPS*Q)
          !
          ! moreover, ch=ch(Ri) where Ri is proportional to deltaTVA=TVA-TVS
          ! (where TVS is the virtual surface temperature) -> this produces a 
          ! minus sign
          !
          ! taken together, we need d(ch)/dT = d(ch)/dTV * dTV/dT = -dch * (1+VIREPS*Q)
          
          dchdt      = -dch(k,:)*(1.+VIREPS*qa(k,:))
          
          DHSDTC1(k) = CPAIR*(ch(k,1) + max( 0., dchdt(1)*dt(1)) )
          DHSDTC2(k) = CPAIR*(ch(k,2) + max( 0., dchdt(2)*dt(2)) )
          DHSDTC4(k) = CPAIR*(ch(k,3) + max( 0., dchdt(3)*dt(3)) )
          DHSDTCS(k) = CPAIR*(ch(k,4) + max( 0., dchdt(4)*dt(4)) )

          ! similarly, d(ch)/dQ = d(ch)/dTV * dTV/dQ = -dch * VIREPS*T
          
          dchdq      = -dch(k,:)*VIREPS*tc(k,:)

          DHSDQA1(k) = max( 0., CPAIR*dchdq(1)*dt(1) )  
          DHSDQA2(k) = max( 0., CPAIR*dchdq(2)*dt(2) )
          DHSDQA4(k) = max( 0., CPAIR*dchdq(3)*dt(3) )
          DHSDQAS(k) = max( 0., CPAIR*dchdq(4)*dt(4) )
          
       end if
       
    end do
    
    ! T2m, Q2m diagnostics
    
    if (sfc_turb_scheme==0 .and. present(ch_tile)  .and. present(cq_tile )) then
       
       ! Louis: compute tile-avg ch and cq so that T2m, Q2m can be diagnosed later
       
       call subtile2tile(N_cat,cat_diagS%ar1,cat_diagS%ar2,cat_diagS%asnow,ch,ch_tile)
       call subtile2tile(N_cat,cat_diagS%ar1,cat_diagS%ar2,cat_diagS%asnow,cq,cq_tile)
       
    end if
    
    if (sfc_turb_scheme==1 .and. present(t2m_tile) .and. present(q2m_tile)) then
       
       ! Helfand-MO: compute tile-avg T2m, Q2m
       
       call subtile2tile(N_cat,cat_diagS%ar1,cat_diagS%ar2,cat_diagS%asnow,t2m,t2m_tile)
       call subtile2tile(N_cat,cat_diagS%ar1,cat_diagS%ar2,cat_diagS%asnow,q2m,q2m_tile)
       
    end if
    
    ! -----------------------------------------------------------------
    !
    ! deallocate 
    
    select case (sfc_turb_scheme)
       
    case (0)     ! Louis
       
       deallocate(cn)
       deallocate(ri)
       deallocate(zt)
       deallocate(zq)
       deallocate(uuu)
       deallocate(ucn)
       deallocate(re)
       
       deallocate(cm)
       deallocate(ww)
       deallocate(dch)
       deallocate(dcq)
       
    case (1)     ! Helfand Monin-Obukhov

       deallocate(tmpintvec)
              
       deallocate(tmprealvec)
       deallocate(psmb)   
       deallocate(psl)    
       deallocate(RHOH)  
       deallocate(VKM)   
       deallocate(USTAR) 
       deallocate(XX)    
       deallocate(YY)    
       deallocate(CU)    
       deallocate(CT)    
       deallocate(RIB)   
       deallocate(ZETA)  
       deallocate(WS)    

       deallocate(t2m)   
       deallocate(q2m)   

       deallocate(u2m)   
       deallocate(v2m)   
       deallocate(t10m)  
       deallocate(q10m)  
       deallocate(u10m)  
       deallocate(v10m)  
       deallocate(u50m)  
       deallocate(v50m)  
    
    end select
    
  end subroutine sfc_turb

  ! *****************************************************************************

  SUBROUTINE PMONTH (         &
       NCH,    ITYP,          &  
       GREEN, ZLT, Z0, D,     &
       SQSCAT, Z2, ZVG,       &
       RSOIL1, RSOIL2, RDC    &
       )
    
    ! reichle+qliu,  8 Oct 2008 - updated to match GEOS5-MERRA parameters
    ! reichle       10 Feb 2011 - moved here from file pmonth_stage4.F90
    ! reichle        8 Dec 2014 - major clean-up, roundoff differences for two reasons:  
    !                              1) use of MAPL_PI here and
    !                              2) associated Dec 5 changes in subroutine propagate_cat()

    IMPLICIT NONE
    
    !         This subroutine sets seasonally varying
    !      variables.  NOTE: LAI is stored as LAI (SIB) / VCOVER (SIB).
    !    
    INTEGER ntyps
    parameter (NTYPS = 10)
    !    
    !     IntegerTYPe:  Vegetation types (biomes).
    !         1:  BrdEvr:  Broadleaf evergreen.
    !         2:  BrdDcd:  Broadleaf deciduous.
    !         3:  Needle:  Needleleaf trees.
    !         4:  GndCvr:  Groundcover (grass).

    INTEGER,                 intent(in)  :: NCH

    INTEGER, dimension(NCH), intent(in)  :: ITYP

    REAL,    dimension(NCH), intent(in)  :: GREEN, ZLT
    
    REAL,    dimension(NCH), intent(out) :: Z0, D, SQSCAT, Z2, ZVG, RSOIL1, RSOIL2, RDC
    !    
    ! local variables
    ! 
    INTEGER VEG, ChNo
    !    
    REAL VGRF11(NTYPS),         VGRF12(NTYPS),&
         VGTR11(NTYPS),         VGTR12(NTYPS),&
         VGROCA(NTYPS),         VGROTD(NTYPS),&
         VGRDRS(NTYPS),         VGZ2(NTYPS),  &
         VGRT(NTYPS)
    
    !    
    REAL ALPHAF, DUM1, DUM2, ROOTL, SCAT, VROOT
    !    
    ! ----- Sarith 4/19/06 -- from GEOScatchGridComp.F90
    !
    REAL, DIMENSION (NTYPS-2) :: VGRDA, VGRDB
    
    ! Correction to RDC formulation -Randy Koster, 4/1/2011
    !
    !data VGRDA / 285.9, 294.9, 652.9,  25.8,  100.7,   &
    !     22.9,  23.8, 23.8/
    !data VGRDB / 5.1 ,  7.2, 10.8,  4.8,  1.8,  5.1,  .000, .000/
    !
    data VGRDA / 285.9, 355.18, 660.24,  30.06,  100.7,  24.36,  23.8, 23.8/
    data VGRDB / 5.1 ,  7.2, 10.5,  4.8,  1.8,  5.1,  .000, .000/

    DATA VGRF11 /0.10, 0.10, 0.07, 0.105 & 
         , 0.10, 0.10, .001, .001, .001, .001/
    
    DATA VGRF12 /0.16, 0.16, 0.16, 0.360&
         , 0.16, 0.16, .001, .001, .001, .001/
    
    DATA VGTR11 /0.05, 0.05, 0.05, 0.070&
         , 0.05, 0.05, .001, .001, .001, .001/
    
    DATA VGTR12 /.001, .001, .001,  .220&
         , .001, .001, .001, .001, .001, .001/
    
    DATA VGROCA / &
         0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6,&
         .1E-6, .1E-6, .1E-6, .1E-6  / 
    
    DATA VGROTD /1.00, 1.00, 0.50, 0.50, 0.50, 0.20, 0.10, 0.10,&
         0.10, 0.10/ 
    
    DATA VGRDRS  /&
         0.75E13, 0.75E13, 0.75E13, 0.40E13, 0.75E13, 0.75E13,&
         0.10E13, 0.10E13, 0.10E13, 0.10E13  /
    
    ! use GEOS5/MERRA values for VGZ2, reichle+qliu, 8 Oct 2008
    ! DATA VGZ2 /35.0, 20.0, 17.0, 0.6, 5.0, 0.6, 0.1, 0.1, 0.1, 0.1/
    DATA VGZ2 /35.0, 20.0, 17.0, 0.6, 0.5, 0.6, 0.01, 0.01, 0.01, 0.01/ ! Dorman and Sellers (1989)
    
    DATA VGRT  / 19700., 7000., 9400., 7000., 7000., 14000., 1., 1., 1., 1./

    ! -----------------------------------------------------------------------------
    
    DO ChNo = 1, NCH
       
       VEG = ITYP(ChNo)
       
       if(VEG.ne.9) then 

          !     COMPUTE SOME PARAMETERS DIRECTLY:
          !  LAI and type dependent parameters; RDC formulation can use veg frac in next version. 
          !  see the last line in the do loop - Sarith
          !         RDC (ChNo) = WLO * VGRDC (KLOW, ITYP (ChNo)) +&
          !              WHI * VGRDC (KLOW + 1, ITYP (ChNo))
          
          ! use GEOS5/MERRA formulation for ROOTL, reichle+qliu, 8 Oct 2008
          !ROOTL = WLO * VGROTL (KLOW, ITYP (ChNo)) +&
          !     WHI * VGROTL (KLOW + 1, ITYP (ChNo))
          ROOTL = VGRT(VEG)

          ! the equations below are identical to those in AGCM tag "Ganymed-4_0"
          ! (GEOS_CatchGridComp.F90)
          ! - reichle, 5 Dec 2014
          
          VROOT  = ROOTL * VGROCA (VEG)
          DUM1   = ALOG (VROOT / (1. - VROOT))
          DUM2   = 1. / (8. * MAPL_PI * ROOTL)
          ALPHAF = DUM2 * (VROOT - 3. -2. * DUM1)

          RSOIL1 (ChNo) =                              &
               VGRDRS (VEG) / (ROOTL * VGROTD (VEG))

          RSOIL2 (ChNo) = ALPHAF / VGROTD (VEG)
          
          SCAT = GREEN (ChNo) *                        &
               (VGTR11 (VEG) + VGRF11 (VEG)) +         &
               (1. - GREEN (ChNo)) *                   &
               (VGTR12 (VEG) + VGRF12 (VEG))
          SQSCAT (ChNo) = SQRT (1. - SCAT)
          
          Z2 (ChNo) = VGZ2 (VEG)
          
          ! Time varying roughness length as in the AGCM -Sarith 4/18/06

          ZVG(ChNo) = (Z2 (ChNo)- (Z2 (ChNo)-MIN_VEG_HEIGHT)*exp(-ZLT (ChNo)))

          Z0 (ChNo) =  Z0_BY_ZVEG * ZVG(ChNo)
          D  (ChNo) =  D0_BY_ZVEG * ZVG(ChNo)    ! zero-plane displacement height
          
          ! Correction to RDC formulation -Randy Koster, 4/1/2011
          ! bad:       RDC = max(VGRDA(VEG)*min(VGRDB(VEG),LAI           ),0.001)
          ! correct:   RDC = max(VGRDA(VEG)*min(1.,        LAI/VGRDB(VEG)),0.001)

          !RDC(ChNo)   = max(VGRDA(ITYP (ChNo))*min(VGRDB(ITYP (ChNo)),       &
          !     ZLT (ChNo)),0.001)
          
          RDC(ChNo) = max( VGRDA(VEG)*min( 1., ZLT(ChNo)/VGRDB(VEG) ), 0.001)
          
       endif
    ENDDO
    
    !    
    RETURN
  END SUBROUTINE PMONTH

  ! *****************************************************************************
    
end module process_cat

! ========================== EOF ====================================