ice_atmo.f90

Go to the documentation of this file.

!/===========================================================================/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!=======================================================================
!BOP
!
! !MODULE: ice_atmo - atm-ice interface: stability based flux calculations
!
! !DESCRIPTION:
!
! Atmospheric boundary interface (stability based flux calculations)
!
! !REVISION HISTORY:
!
! author: Elizabeth C. Hunke, LANL
!
! Vectorized by Clifford Chen (Fujitsu) and William H. Lipscomb (LANL)
!
! !INTERFACE:
!
      module ice_atmo
!
! !USES:
!
      use ice_domain
      use ice_constants
      use ice_flux
      use ice_state
!
!EOP
!
      implicit none

!=======================================================================
 
      contains

!=======================================================================
!BOP
!
! !IROUTINE: atmo_boundary_layer - compute coefficients for atm-ice fluxes, stress and Tref
!
! !INTERFACE:
!
      subroutine atmo_boundary_layer (ni, sfctype, Tsf,  &
                                     strx, stry, Trf, Qrf, delt, delq)
!
! !DESCRIPTION:
!
! Compute coefficients for atm/ice fluxes, stress, and reference 
! temperature. NOTE: \! (1) all fluxes are positive downward, \! (2) here, tstar = (WT)/U*, and qstar = (WQ)/U*, \! (3) wind 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:
!
! author: Elizabeth C. Hunke, LANL
!
! !USES:
!
      use ice_grid       ! for tmask
!
! !INPUT/OUTPUT PARAMETERS:
!
      integer (kind=int_kind), intent(in) :: &
        ni        ! thickness category index

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

      real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
        tsf      ! surface temperature of ice or ocean

      real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) ::&
         strx     &! x surface stress (N)
      ,  stry     &! y surface stress (N)
      ,  trf      &! reference height temperature  (K)
      ,  qrf      &! reference height specific humidity (kg/kg)
      ,  delt     &! potential T difference   (K)
      ,  delq     ! humidity difference      (kg/kg)
!
!EOP
!
       integer (kind=int_kind) ::&
         k       & ! iteration index
      ,  i, j     ! horizontal indices

      real (kind=dbl_kind) :: &
         tsfk     &! surface temperature in Kelvin (K)
      ,  xqq      &! temporary variable
      ,  psimh    &! stability function at zlvl   (momentum)
      ,  tau      &! stress at zlvl
      ,  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, dqsatdt
      ,  ttt      &! for qsat, dqsatdt
      ,  qsat     &! the saturation humidity of air (kg/m^3)
      ,  lheat    ! Lvap or Lsub, depending on surface type

      real (kind=dbl_kind), dimension (1:(ihi-ilo+1)*(jhi-jlo+1)) ::&
         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)
      ,  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)

      integer (kind=int_kind), dimension (1:(ihi-ilo+1)*(jhi-jlo+1)) ::&
         indxi   & ! compressed index in i-direction
      ,  indxj    ! compressed index in j-direction

      integer (kind=int_kind) ::&
         icells  & ! number of cells with aicen > puny or tmask true
      ,  ij       ! combined ij index

      real (kind=dbl_kind), parameter ::   &
         cpvir = cp_wv/cp_air - c1i &! defined as cp_wv/cp_air - 1.
      ,  ztrf  = c2i                &! reference height for air temp (m)
      ,  umin  = c1i                 ! 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((c1i+xd*(c2i+xd))*(c1i+xd*xd)/c8i)  &
                  - c2i*atan(xd) + pi*p5
!ech     &            - c2i*atan(xd) + 1.571_dbl_kind

      psixhu(xd)  =  c2i * log((c1i + xd*xd)/c2i)

      al2 = log(zref/ztrf)

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

      do j = jlo,jhi
      do i = ilo,ihi
         lhcoef(i,j) = c0i
         shcoef(i,j) = c0i
         strx(i,j)   = c0i
         stry(i,j)   = c0i
         trf(i,j)    = c0i
         qrf(i,j)    = c0i
         delt(i,j)   = c0i
         delq(i,j)   = c0i
      enddo
      enddo

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

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

      !------------------------------------------------------------
      ! identify grid cells with ice
      !------------------------------------------------------------
         icells = 0
         do j = jlo,jhi
         do i = ilo,ihi
            if ( aicen(i,j,ni) > puny) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo
         enddo

      !------------------------------------------------------------
      ! define some needed variables
      !------------------------------------------------------------
         qqq  = qqqice          ! for qsat
         ttt  = tttice          ! for qsat
         lheat = lsub           ! ice to vapor
         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)
            vmag(ij) = max(umin, wind(i,j))
            rdn(ij)  = vonkar/log(zref/iceruf) ! neutral coefficient
         enddo   ! ij

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

      !------------------------------------------------------------
      ! identify ocean cells
      !------------------------------------------------------------
         icells = 0
         do j = jlo,jhi
         do i = ilo,ihi
            if ( tmask(i,j) ) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo
         enddo

       !------------------------------------------------------------
       ! define some needed variables
       !------------------------------------------------------------
         qqq  = qqqocn
         ttt  = tttocn
         lheat = lvap       ! liquid to vapor
         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)
            vmag(ij) = max(umin, wind(i,j))
            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) * (c1i + 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*(c1i + 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)

      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)/(c1i/zvir+qa(i,j)))                    &
                   / ustar(ij)**2                                      
            hol(ij)    = sign( min(abs(hol(ij)),c10i), hol(ij) )        
            stable(ij) = p5 + sign(p5 , hol(ij))                       
            xqq    = max(sqrt(abs(c1i - c16*hol(ij))) , c1i)             
            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)                                &
                    + (c1i - stable(ij))*psimhu(xqq)                    
            psixh(ij)  = psimhs*stable(ij)                            &
                    + (c1i - stable(ij))*psixhu(xqq)

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

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

         enddo                  ! ij
      enddo                     ! end iteration


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

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

!ccsm      shcoef(i,j) = rhoa(i,j) * ustar(ij) * cp(ij) * rh(ij)
      ! 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) + c1i
         lhcoef(i,j) = rhoa(i,j) * ustar(ij) * lheat  * re(ij)
      !------------------------------------------------------------
      ! momentum flux
      !------------------------------------------------------------
      ! tau = rhoa(i,j) * ustar * ustar 
      ! strx = tau * uatm(i,j) / vmag 
      ! stry = tau * vatm(i,j) / vmag 
      !------------------------------------------------------------

         tau = rhoa(i,j) * ustar(ij) * rd(ij) ! not the stress at zlvl(i,j)
         strx(i,j) = tau * uatm(i,j)
         stry(i,j) = tau * vatm(i,j)

      !------------------------------------------------------------
      ! Compute diagnostics: 2m ref T & Q
      !------------------------------------------------------------
         hol(ij) = hol(ij)*ztrf/zlvl(i,j)
         xqq     = max( c1i, sqrt(abs(c1i-c16*hol(ij))) )
         xqq     = sqrt(xqq)
         psix2   = -c5i*hol(ij)*stable(ij) + (c1i-stable(ij))*psixhu(xqq)  
         fac     = (rh(ij)/vonkar) * (alz(ij) + al2 - psixh(ij) + psix2)
         trf(i,j)= pott(i,j) - delt(i,j)*fac
         trf(i,j)= trf(i,j) - p01*ztrf ! pot temp to temp correction
         fac     = (re(ij)/vonkar) * (alz(ij) + al2 - psixh(ij) + psix2)
         qrf(i,j)= qa(i,j) - delq(i,j)*fac

      enddo                     ! ij

      end subroutine atmo_boundary_layer

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

      end module ice_atmo

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