!=======================================================================
!BOP
!
! !MODULE: ice_atmo - atm-ice interface: stability based flux calculations
!
! !DESCRIPTION:
!
! Atmospheric boundary interface (stability based flux calculations)
!
! !REVISION HISTORY:
!  SVN:$Id: ice_atmo.F90 49 2007-01-11 22:07:00Z eclare $
!
! author: Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
! 2004: Block structure added by William Lipscomb
! 2006: Converted to free source form (F90) by Elizabeth Hunke
!
! !INTERFACE:
!
      module ice_atmo
!
! !USES:
!
      use ice_kinds_mod
      use ice_constants
!
!EOP
!
      implicit none
      save

      character (len=char_len) :: &
         atmbndy ! atmo boundary method, 'default' ('ccsm3') or 'constant'

      logical (kind=log_kind) :: &
         calc_strair ! if true, calculate wind stress components

!=======================================================================

      contains

!=======================================================================
!BOP
!
! !IROUTINE: atmo_boundary_layer - compute coefficients for atm-ice fluxes, 
!                                  stress and Tref/Qref
!
! !INTERFACE:
!
      subroutine atmo_boundary_layer (nx_block, ny_block, &
                                      sfctype,  icells,   &
                                      indxi,    indxj,    & 
                                      Tsf,      potT,     &
                                      uatm,     vatm,     &  
                                      wsresp,   tau_est,  &
                                      ugust,              &
                                      uvel,     vvel,     &  
                                      wind,     zlvl,     &  
                                      Qa,       rhoa,     &
                                      strx,     stry,     &   
                                      Uref,               &
                                      Tref,     Qref,     &
                                      delt,     delq,     &
                                      lhcoef,   shcoef)

! !DESCRIPTION:
!
! Compute coefficients for atm/ice fluxes, stress, and reference
! temperature and humidity. NOTE: \\
! (1) all fluxes are positive downward,  \\
! (2) here, tstar = (WT)/U*, and qstar = (WQ)/U*,  \\
! (3a) wind speeds should all be above a minimum speed (eg. 1.0 m/s). \\
! (3b) (wind-ice) speeds should all be above a minimum speed (eg. 1.0 m/s). \\
!
! ASSUME:
!  The saturation humidity of air at T(K): qsat(T)  (kg/m**3)
!
! Code originally based on CSM1
!
! !REVISION HISTORY: same as module
!
! !USES:
!
#ifdef CCSMCOUPLED
        use shr_flux_mod, only: shr_flux_update_stress
#endif
! !INPUT/OUTPUT PARAMETERS:
!
      integer (kind=int_kind), intent(in) :: &
         nx_block, ny_block, & ! block dimensions
         icells                ! number of cells that require atmo fluxes

      integer (kind=int_kind), dimension(nx_block*ny_block), &
         intent(in) :: &
         indxi, indxj    ! compressed i and j indices

      character (len=3), intent(in) :: &
         sfctype      ! ice or ocean

      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
         Tsf      , & ! surface temperature of ice or ocean
         potT     , & ! air potential temperature  (K)
         uatm     , & ! x-direction wind speed (m/s)
         vatm     , & ! y-direction wind speed (m/s)
         wsresp   , & ! response of atmospheric boundary layer to stress (m/s/Pa)
         tau_est  , & ! estimated boundary layer equilibrium stress (Pa)
         ugust    , & ! gustiness from atmosphere (m/s)
         uvel     , & ! x-direction ice speed (m/s)
         vvel     , & ! y-direction ice speed (m/s)
         wind     , & ! wind speed (m/s)
         zlvl     , & ! atm level height (m)
         Qa       , & ! specific humidity (kg/kg)
         rhoa         ! air density (kg/m^3)

      real (kind=dbl_kind), dimension (nx_block,ny_block), &
         intent(inout) :: &
         strx     , & ! x surface stress (N)
         stry         ! y surface stress (N)

      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
         Uref     , & ! reference height wind speed (m/s)
         Tref     , & ! reference height temperature  (K)
         Qref     , & ! reference height specific humidity (kg/kg)
         delt     , & ! potential T difference   (K)
         delq     , & ! humidity difference      (kg/kg)
         shcoef   , & ! transfer coefficient for sensible heat
         lhcoef       ! transfer coefficient for latent heat
!
!EOP
!
       integer (kind=int_kind) :: &
         k     , & ! iteration index
         i, j  , & ! horizontal indices
         ij        ! combined ij index

      real (kind=dbl_kind) :: &
         TsfK  , & ! surface temperature in Kelvin (K)
         xqq   , & ! temporary variable
         psimh , & ! stability function at zlvl   (momentum)
         taufac, & ! factor limiting changes in tau during iteration
         fac   , & ! interpolation factor
         al2   , & ! ln(z10   /zTrf)
         psix2 , & ! stability function at zTrf   (heat and water)
         psimhs, & ! stable profile
         ssq   , & ! sat surface humidity     (kg/kg)
         qqq   , & ! for qsat, dqsfcdt
         TTT   , & ! for qsat, dqsfcdt
         qsat  , & ! the saturation humidity of air (kg/m^3)
         Lheat     ! Lvap or Lsub, depending on surface type

      real (kind=dbl_kind), dimension (icells) :: &
         ustar , & ! ustar (m/s)
         tstar , & ! tstar
         qstar , & ! qstar
         rdn   , & ! sqrt of neutral exchange coefficient (momentum)
         rhn   , & ! sqrt of neutral exchange coefficient (heat)
         ren   , & ! sqrt of neutral exchange coefficient (water)
         rd    , & ! sqrt of exchange coefficient (momentum)
         re    , & ! sqrt of exchange coefficient (water)
         rh    , & ! sqrt of exchange coefficient (heat)
         vmag  , & ! surface wind magnitude   (m/s)
         vmagit, & ! iteration loop surface wind magnitude   (m/s)
         wind0 , & ! original wind without limiter   (m/s)
         windit, & ! iteration loop wind without limiter   (m/s)
         tau   , & ! stress at zlvl (Pa)
         taupr , & ! stress from previous iteration (Pa)
         dtau  , & ! difference in stress vs previous iteration (Pa)
         dtaupr, & ! difference in stress over previous iteration (Pa)
         alz   , & ! ln(zlvl  /z10)
         thva  , & ! virtual temperature      (K)
         cp    , & ! specific heat of moist air
         hol   , & ! H (at zlvl  ) over L
         stable, & ! stability factor
         psixh     ! stability function at zlvl   (heat and water)

      real (kind=dbl_kind), parameter :: &
         cpvir = cp_wv/cp_air-c1, & ! defined as cp_wv/cp_air - 1.
         zTrf  = c2             , & ! reference height for air temp (m)
         umin  = c1                 ! minimum wind speed (m/s)

      ! local functions
      real (kind=dbl_kind) :: &
         xd    , & ! dummy argument
         psimhu, & ! unstable part of psimh
         psixhu    ! unstable part of psimx

      !------------------------------------------------------------
      ! Define functions
      !------------------------------------------------------------

      psimhu(xd)  = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
                  - c2*atan(xd) + pih
!ech                  - c2*atan(xd) + 1.571_dbl_kind

      psixhu(xd)  =  c2 * log((c1 + xd*xd)/c2)

      al2 = log(zref/zTrf)

      !------------------------------------------------------------
      ! Initialize
      !------------------------------------------------------------

      do j = 1, ny_block
      do i = 1, nx_block
         Uref(i,j) = c0
         Tref(i,j) = c0
         Qref(i,j) = c0
         delt(i,j) = c0
         delq(i,j) = c0
         shcoef(i,j) = c0
         lhcoef(i,j) = c0
      enddo
      enddo

      !------------------------------------------------------------
      ! Compute turbulent flux coefficients, wind stress, and
      ! reference temperature and humidity.
      !------------------------------------------------------------

      !------------------------------------------------------------
      ! define variables that depend on surface type
      !------------------------------------------------------------

      if (sfctype(1:3)=='ice') then

         qqq  = qqqice          ! for qsat
         TTT  = TTTice          ! for qsat
         Lheat = Lsub           ! ice to vapor
         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)
            wind0(ij) = max(wind(i,j), 0.01_dbl_kind)
            vmag(ij) = max(umin, wind0(ij) + ugust(i,j))
!---------- (3b) option by Andrew Roberts
!            wind0(ij)   = sqrt( (uatm(i,j)-uvel(i,j))**2 + (vatm(i,j)-vvel(i,j))**2)
!---------- (3b) option end
            rdn(ij)  = vonkar/log(zref/iceruf) ! neutral coefficient
         enddo   ! ij

      elseif (sfctype(1:3)=='ocn') then

         qqq  = qqqocn
         TTT  = TTTocn
         Lheat = Lvap           ! liquid to vapor
         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)
            wind0(ij) = max(wind(i,j), 0.01_dbl_kind)
            vmag(ij) = max(umin, wind0(ij) + ugust(i,j))
!---------- (3b) option by Andrew Roberts
!            wind0(ij)   = sqrt( (uatm(i,j)-uvel(i,j))**2 + (vatm(i,j)-vvel(i,j))**2)
!---------- (3b) option end
            rdn(ij)  = sqrt(0.0027_dbl_kind/vmag(ij) &
                    + .000142_dbl_kind + .0000764_dbl_kind*vmag(ij))
         enddo   ! ij

      endif   ! sfctype

      do ij = 1, icells
         i = indxi(ij)
         j = indxj(ij)

      !------------------------------------------------------------
      ! define some more needed variables
      !------------------------------------------------------------

         TsfK       = Tsf(i,j) + Tffresh     ! surface temp (K)
         qsat       = qqq * exp(-TTT/TsfK)   ! saturation humidity (kg/m^3)
         ssq        = qsat / rhoa(i,j)       ! sat surf hum (kg/kg)

         thva(ij)   = potT(i,j) * (c1 + zvir * Qa(i,j)) ! virtual pot temp (K)
         delt(i,j)  = potT(i,j) - TsfK       ! pot temp diff (K)
         delq(i,j)  = Qa(i,j) - ssq          ! spec hum dif (kg/kg)
         alz(ij)    = log(zlvl(i,j)/zref)
         cp(ij)     = cp_air*(c1 + cpvir*ssq)

      !------------------------------------------------------------
      ! first estimate of Z/L and ustar, tstar and qstar
      !------------------------------------------------------------

         ! neutral coefficients, z/L = 0.0
         rhn(ij) = rdn(ij)
         ren(ij) = rdn(ij)

         ! ustar,tstar,qstar
         ustar(ij) = rdn(ij) * vmag(ij)
         tstar(ij) = rhn(ij) * delt(i,j)
         qstar(ij) = ren(ij) * delq(i,j)

         ! Set up variables for velocity iteration.
         vmagit(ij) = vmag(ij)
         windit(ij) = wind0(ij)
         taupr(ij) = tau_est(i,j)
         dtau(ij) = 1.e100_dbl_kind

      enddo                     ! ij

      !------------------------------------------------------------
      ! iterate to converge on Z/L, ustar, tstar and qstar
      !------------------------------------------------------------

      do k=1,5

         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)

        ! compute stability & evaluate all stability functions
            hol(ij) = vonkar * gravit * zlvl(i,j) &
                   * (tstar(ij)/thva(ij) &
                    + qstar(ij)/(c1/zvir+Qa(i,j))) &
                   / ustar(ij)**2
            hol(ij)    = sign( min(abs(hol(ij)),c10), hol(ij) )
            stable(ij) = p5 + sign(p5 , hol(ij))
            xqq    = max(sqrt(abs(c1 - c16*hol(ij))) , c1)
            xqq    = sqrt(xqq)

            ! Jordan et al 1999
            psimhs = -(0.7_dbl_kind*hol(ij) &
                   + 0.75_dbl_kind*(hol(ij)-14.3_dbl_kind) &
                   * exp(-0.35_dbl_kind*hol(ij)) + 10.7_dbl_kind)
            psimh  = psimhs*stable(ij) &
                    + (c1 - stable(ij))*psimhu(xqq)
            psixh(ij)  = psimhs*stable(ij) &
                    + (c1 - stable(ij))*psixhu(xqq)

        ! shift all coeffs to measurement height and stability
            rd(ij) = rdn(ij) / (c1+rdn(ij)/vonkar*(alz(ij)-psimh))
            rh(ij) = rhn(ij) / (c1+rhn(ij)/vonkar*(alz(ij)-psixh(ij)))
            re(ij) = ren(ij) / (c1+ren(ij)/vonkar*(alz(ij)-psixh(ij)))

        ! update ustar, tstar, qstar using updated, shifted coeffs
            ustar(ij) = rd(ij) * vmagit(ij)
            tstar(ij) = rh(ij) * delt(i,j)
            qstar(ij) = re(ij) * delq(i,j)

#ifdef CCSMCOUPLED
            tau(ij) = rhoa(i,j) * ustar(ij) * rd(ij) * windit(ij)
            call shr_flux_update_stress(wind0(ij), wsresp(i,j), tau_est(i,j), &
                 tau(ij), taupr(ij), dtau(ij), dtaupr(ij), windit(ij))
            vmagit(ij) = max(umin, windit(ij) + ugust(i,j))
#endif

         enddo                  ! ij
      enddo                     ! end iteration

      if (calc_strair) then

      ! initialize
      do j = 1, ny_block
      do i = 1, nx_block
         strx(i,j) = c0
         stry(i,j) = c0
      enddo
      enddo

#ifdef CPRCRAY
!DIR$ CONCURRENT !Cray
#endif
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells
         i = indxi(ij)
         j = indxj(ij)

      !------------------------------------------------------------
      ! momentum flux
      !------------------------------------------------------------
      ! tau = rhoa(i,j) * ustar * ustar
      ! strx = tau * uatm(i,j) / vmag
      ! stry = tau * vatm(i,j) / vmag
      !------------------------------------------------------------

         tau(ij) = rhoa(i,j) * ustar(ij) * rd(ij) ! not the stress at zlvl(i,j)
         strx(i,j) = tau(ij) * (uatm(i,j)-uvel(i,j)) * (windit(ij) / wind0(ij))
         stry(i,j) = tau(ij) * (vatm(i,j)-vvel(i,j)) * (windit(ij) / wind0(ij))

      enddo                     ! ij

      endif                     ! calc_strair

#ifdef CPRCRAY
!DIR$ CONCURRENT !Cray
#endif
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells
         i = indxi(ij)
         j = indxj(ij)

      !------------------------------------------------------------
      ! coefficients for turbulent flux calculation
      !------------------------------------------------------------
      ! add windless coefficient for sensible heat flux
      ! as in Jordan et al (JGR, 1999)
      !------------------------------------------------------------

         shcoef(i,j) = rhoa(i,j) * ustar(ij) * cp(ij) * rh(ij) + c1
         lhcoef(i,j) = rhoa(i,j) * ustar(ij) * Lheat  * re(ij)

      !------------------------------------------------------------
      ! Compute diagnostics: 2m ref T & Q and 10m wind speed.
      !------------------------------------------------------------
         hol(ij)  = hol(ij)*zTrf/zlvl(i,j)
         xqq      = max( c1, sqrt(abs(c1-c16*hol(ij))) )
         xqq      = sqrt(xqq)
         psix2    = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psixhu(xqq)
         fac      = (rh(ij)/vonkar) &
                  * (alz(ij) + al2 - psixh(ij) + psix2)
         Tref(i,j)= potT(i,j) - delt(i,j)*fac
         Tref(i,j)= Tref(i,j) - p01*zTrf ! pot temp to temp correction
         fac      = (re(ij)/vonkar) &
                  * (alz(ij) + al2 - psixh(ij) + psix2)
         Qref(i,j)= Qa(i,j) - delq(i,j)*fac

         Uref(i,j)= vmag(ij) * rd(ij) / rdn(ij)
      enddo                     ! ij

      end subroutine atmo_boundary_layer

!=======================================================================
!BOP
!
! !IROUTINE: atmo_boundary_const - compute coeeficients for atm-ice fluxes
!
!
! !INTERFACE:
!
      subroutine atmo_boundary_const (nx_block, ny_block, &
                                      sfctype,  icells,   &
                                      indxi,    indxj,    & 
                                      uatm,     vatm,     &  
                                      wind,     rhoa,     &
                                      strx,     stry,     &   
                                      lhcoef,   shcoef)

! !DESCRIPTION:
!
! Compute coefficients for atm/ice fluxes, stress
! NOTE: \\
! (1) all fluxes are positive downward,  \\
! (2) reference temperature and humidity are NOT computed
!
! !REVISION HISTORY: same as module
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      integer (kind=int_kind), intent(in) :: &
         nx_block, ny_block, & ! block dimensions
         icells                ! number of cells that require atmo fluxes

      integer (kind=int_kind), dimension(nx_block*ny_block), &
         intent(in) :: &
         indxi, indxj    ! compressed i and j indices

      character (len=3), intent(in) :: &
         sfctype      ! ice or ocean

      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
         uatm     , & ! x-direction wind speed (m/s)
         vatm     , & ! y-direction wind speed (m/s)
         wind     , & ! wind speed (m/s)
         rhoa         ! air density (kg/m^3)

      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout):: &
         strx     , & ! x surface stress (N)
         stry         ! y surface stress (N)

      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
         shcoef   , & ! transfer coefficient for sensible heat
         lhcoef       ! transfer coefficient for latent heat
!
!EOP
!
       integer (kind=int_kind) :: &
         i, j, & ! horizontal indices
         ij      ! combined ij index

      real (kind=dbl_kind) :: &
         tau, &  ! stress at zlvl
         Lheat   ! Lvap or Lsub, depending on surface type

      !------------------------------------------------------------
      ! Initialize
      !------------------------------------------------------------

      do j = 1, ny_block
      do i = 1, nx_block
         shcoef(i,j) = c0
         lhcoef(i,j) = c0
      enddo
      enddo

      if (calc_strair) then

      do j = 1, ny_block
      do i = 1, nx_block
         strx(i,j) = c0
         stry(i,j) = c0
      enddo
      enddo

#ifdef CPRCRAY
!DIR$ CONCURRENT !Cray
#endif
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells
         i = indxi(ij)
         j = indxj(ij)

      !------------------------------------------------------------
      ! momentum flux
      !------------------------------------------------------------

         tau = rhoa(i,j) * 0.0012_dbl_kind * wind(i,j)
!AOMIP         tau = rhoa(i,j) * (1.10_dbl_kind + c4*p01*wind(i,j)) &
!AOMIP                         * wind(i,j) * p001
         strx(i,j) = tau * uatm(i,j)
         stry(i,j) = tau * vatm(i,j)

      enddo                     ! ij

      endif                     ! calc_strair

      !------------------------------------------------------------
      ! define variables that depend on surface type
      !------------------------------------------------------------

      if (sfctype(1:3)=='ice') then
         Lheat = Lsub           ! ice to vapor
      elseif (sfctype(1:3)=='ocn') then
         Lheat = Lvap           ! liquid to vapor
      endif   ! sfctype

#ifdef CPRCRAY
!DIR$ CONCURRENT !Cray
#endif
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells
         i = indxi(ij)
         j = indxj(ij)

      !------------------------------------------------------------
      ! coefficients for turbulent flux calculation
      !------------------------------------------------------------

         shcoef(i,j) = (1.20e-3_dbl_kind)*cp_air*rhoa(i,j)*wind(i,j)
         lhcoef(i,j) = (1.50e-3_dbl_kind)*Lheat *rhoa(i,j)*wind(i,j)

      enddo                     ! ij

      end subroutine atmo_boundary_const

!=======================================================================

      end module ice_atmo

!=======================================================================
