!        Generated by TAPENADE     (INRIA, Tropics team)
!  Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!  
! $Id: ice_atmo.F,v 1.10 2004/02/09 17:58:46 lipscomb Exp $
!=======================================================================
!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
  USE ice_constants
  USE ice_domain
  USE ice_flux
  USE ice_grid
  USE ice_kinds_mod
  USE ice_model_size
  USE ice_state
  REAL(KIND=DBL_KIND) :: umin

CONTAINS
  SUBROUTINE ATMO_BOUNDARY_LAYER(n, sfctype, tsf, strx, stry, trf, qrf, &
&    delt, delq)
    IMPLICIT NONE
!ech     &            - c2*atan(xd) + 1.571_dbl_kind
!
!
! !DESCRIPTION:
!
! Compute coefficients for atm/ice fluxes, stress, and reference 
! temperature. NOTE: \\
! (1) all fluxes are positive downward, \\
! (2) net heat flux = fswabs + flwup + (1-emissivity)flwdn + fsh + flh, \\
! (3) here, tstar = (WT)/U*, and qstar = (WQ)/U*, \\
! (4) 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:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
! thickness category index
!
! ice or ocean
!
! surface temperature of ice or ocean
!
! x surface stress (N)
! y surface stress (N)
! reference height temperature  (K)
! reference height specific humidity (kg/kg)
! potential T difference   (K)
! humidity difference      (kg/kg)
!
!EOP
!
! iteration index
! horizontal indices
!
! surface temperature in Kelvin (K)
! temporary variable
! stability function at zlvl   (momentum)
! stress at zlvl
! interpolation factor
! ln(z10   /zTrf)
! stability function at zTrf   (heat and water)
! stable profile 
! sat surface humidity     (kg/kg)
! for qsat, dqsatdt
! for qsat, dqsatdt
! the saturation humidity of air (kg/m^3)
! Lvap or Lsub, depending on surface type
!
! ustar (m/s)
! tstar
! qstar
! sqrt of neutral exchange coefficient (momentum)
! sqrt of neutral exchange coefficient (heat)
! sqrt of neutral exchange coefficient (water)
! sqrt of exchange coefficient (momentum)
! sqrt of exchange coefficient (water)
! sqrt of exchange coefficient (heat)
! surface wind magnitude   (m/s)
! ln(zlvl  /z10)
! virtual temperature      (K)
! specific heat of moist air
! H (at zlvl  ) over L
! stability factor
! stability function at zlvl   (heat and water)
!
! compressed index in i-direction
! compressed index in j-direction
!
! number of cells with aicen > puny or tmask true
! combined ij index
!
! defined as cp_wv/cp_air - 1.
! reference height for air temp (m)
! local functions
!jkim     &,  umin  = c1                ! minimum wind speed (m/s)
!
! dummy argument  
! unstable part of psimh
! unstable part of psimx
!------------------------------------------------------------
! Define functions
!------------------------------------------------------------
!
!
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: delq
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: delt
    INTEGER(KIND=INT_KIND),INTENT(IN) :: n
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: qrf
    CHARACTER(len=3),INTENT(IN) :: sfctype
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: strx
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: stry
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT) :: trf
    REAL(KIND=DBL_KIND), DIMENSION(ilo:ihi, jlo:jhi),INTENT(IN) :: tsf
    REAL(KIND=DBL_KIND),PARAMETER :: cpvir=cp_wv/cp_air-c1
    REAL(KIND=DBL_KIND),PARAMETER :: ztrf=c2
    REAL :: abs1
    REAL :: abs2
    REAL(KIND=DBL_KIND) :: arg1
    INTEGER(KIND=INT_KIND) :: icells, ij
    INTEGER(KIND=INT_KIND) :: indxi((ihi-ilo+1)*(jhi-jlo+1)), indxj((ihi&
&    -ilo+1)*(jhi-jlo+1))
    INTEGER(KIND=INT_KIND) :: i, j, k
    REAL :: min1, x1
    REAL(KIND=DBL_KIND) :: psimhu, psixhu, xd
    REAL(KIND=DBL_KIND) :: al2, fac, lheat, psimh, psimhs, psix2, qqq, &
&    qsat, ssq, tau, tsfk, ttt, xqq
    REAL(KIND=DBL_KIND) :: alz((ihi-ilo+1)*(jhi-jlo+1)), cp((ihi-ilo+1)*&
&    (jhi-jlo+1)), hol((ihi-ilo+1)*(jhi-jlo+1)), psixh((ihi-ilo+1)*(jhi-&
&    jlo+1)), qstar((ihi-ilo+1)*(jhi-jlo+1)), rd((ihi-ilo+1)*(jhi-jlo+1)&
&    ), rdn((ihi-ilo+1)*(jhi-jlo+1)), re((ihi-ilo+1)*(jhi-jlo+1)), ren((&
&    ihi-ilo+1)*(jhi-jlo+1)), rh((ihi-ilo+1)*(jhi-jlo+1)), rhn((ihi-ilo+&
&    1)*(jhi-jlo+1)), stable((ihi-ilo+1)*(jhi-jlo+1)), thva((ihi-ilo+1)*&
&    (jhi-jlo+1)), tstar((ihi-ilo+1)*(jhi-jlo+1)), ustar((ihi-ilo+1)*(&
&    jhi-jlo+1)), vmag((ihi-ilo+1)*(jhi-jlo+1))
    REAL :: x2
    REAL :: y1
    INTRINSIC EXP, MAX, SIGN, ABS, LOG, ATAN, MIN, SQRT


!
!------------------------------------------------------------
! Initialize
!------------------------------------------------------------
    arg1 = zref/ztrf
    al2 = LOG(arg1)
!
!
    DO j=jlo,jhi
      DO i=ilo,ihi
        lhcoef(i, j) = c0
        shcoef(i, j) = c0
        strx(i, j) = c0
        stry(i, j) = c0
        trf(i, j) = c0
        qrf(i, j) = c0
        delt(i, j) = c0
        delq(i, j) = c0
      END DO
    END DO
!------------------------------------------------------------
! Compute turbulent flux coefficients, wind stress, and 
! reference temperature and humidity.
!------------------------------------------------------------
! sfctype
!
!
    IF (sfctype(1:3) .EQ. 'ice') THEN
!------------------------------------------------------------
! identify grid cells with ice
!------------------------------------------------------------
!
      icells = 0
      DO j=jlo,jhi
        DO i=ilo,ihi
          IF (aicen(i, j, n) .GT. puny) THEN
            icells = icells + 1
            indxi(icells) = i
            indxj(icells) = j
          END IF
        END DO
      END DO
!------------------------------------------------------------
! define some needed variables
!------------------------------------------------------------
!
! for qsat
      qqq = qqqice
! for qsat
      ttt = tttice
! ice to vapor
      lheat = lsub
      DO ij=1,icells
        i = indxi(ij)
        j = indxj(ij)
        IF (umin .LT. wind(i, j)) THEN
          vmag(ij) = wind(i, j)
        ELSE
          vmag(ij) = umin
        END IF
! neutral coefficient
        arg1 = zref/iceruf
        rdn(ij) = vonkar/LOG(arg1)
      END DO
    ELSE IF (sfctype(1:3) .EQ. 'ocn') THEN
! ij
!
!------------------------------------------------------------
! 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
          END IF
        END DO
      END DO
!------------------------------------------------------------
! define some needed variables
!------------------------------------------------------------
!
      qqq = qqqocn
      ttt = tttocn
! liquid to vapor
      lheat = lvap
      DO ij=1,icells
        i = indxi(ij)
        j = indxj(ij)
        IF (umin .LT. wind(i, j)) THEN
          vmag(ij) = wind(i, j)
        ELSE
          vmag(ij) = umin
        END IF
        arg1 = 0.0027_dbl_kind/vmag(ij) + .000142_dbl_kind + &
&          .0000764_dbl_kind*vmag(ij)
        rdn(ij) = SQRT(arg1)
      END DO
    END IF
! ij
!
!
!
    DO ij=1,icells
      i = indxi(ij)
!------------------------------------------------------------
! define some more needed variables
!------------------------------------------------------------
      j = indxj(ij)
!
!
! surface temp (K)
      tsfk = tsf(i, j) + tffresh
! saturation humidity (kg/m^3)
      arg1 = -(ttt/tsfk)
      qsat = qqq*EXP(arg1)
! sat surf hum (kg/kg)
      ssq = qsat/rhoa(i, j)
!
! virtual pot temp (K)
      thva(ij) = pott(i, j)*(c1+zvir*qa(i, j))
! pot temp diff (K)
      delt(i, j) = pott(i, j) - tsfk
! spec hum dif (kg/kg)
      delq(i, j) = qa(i, j) - ssq
      arg1 = zlvl(i, j)/zref
      alz(ij) = LOG(arg1)
!------------------------------------------------------------
! first estimate of Z/L and ustar, tstar and qstar
!------------------------------------------------------------
! neutral coefficients, z/L = 0.0 
      cp(ij) = cp_air*(c1+cpvir*ssq)
!
!
      rhn(ij) = rdn(ij)
! ustar,tstar,qstar
      ren(ij) = rdn(ij)
!
      ustar(ij) = rdn(ij)*vmag(ij)
      tstar(ij) = rhn(ij)*delt(i, j)
      qstar(ij) = ren(ij)*delq(i, j)
    END DO
!
! ij
!------------------------------------------------------------
! iterate to converge on Z/L, ustar, tstar and qstar
!------------------------------------------------------------
!
!
!         do k=1,2             !!! ech thinks this should be at least 5
    DO k=1,5
!
      DO ij=1,icells
        i = indxi(ij)
! compute stability & evaluate all stability functions 
        j = indxj(ij)
!
        hol(ij) = vonkar*gravit*zlvl(i, j)*(tstar(ij)/thva(ij)+qstar(ij)&
&          /(c1/zvir+qa(i, j)))/ustar(ij)**2
        IF (hol(ij) .GE. 0.) THEN
          x1 = hol(ij)
        ELSE
          x1 = -hol(ij)
        END IF
        IF (x1 .GT. c10) THEN
          min1 = c10
        ELSE
          min1 = x1
        END IF
        hol(ij) = SIGN(min1, hol(ij))
        stable(ij) = p5 + SIGN(p5, hol(ij))
        IF (c1 - c16*hol(ij) .GE. 0.) THEN
          abs1 = c1 - c16*hol(ij)
        ELSE
          abs1 = -(c1-c16*hol(ij))
        END IF
        x2 = SQRT(abs1)
        IF (x2 .LT. c1) THEN
          xqq = c1
        ELSE
          xqq = x2
        END IF
        xqq = SQRT(xqq)
!ccsm       psimh  = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psimhu(xqq)
!ccsm       psixh(ij)  = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psixhu(xqq)
!
! Jordan et al 1999
        arg1 = -(0.35_dbl_kind*hol(ij))
        psimhs = -(0.7_dbl_kind*hol(ij)+0.75_dbl_kind*(hol(ij)-&
&          14.3_dbl_kind)*EXP(arg1)+10.7_dbl_kind)
!
        arg1 = (c1+xqq*(c2+xqq))*(c1+xqq*xqq)/c8
        psimh = psimhs*stable(ij) + (c1-stable(ij))*(LOG(arg1)-c2*ATAN(&
&          xqq)+pi*p5)
! shift all coeffs to measurement height and stability
        arg1 = (c1+xqq*xqq)/c2
        psixh(ij) = psimhs*stable(ij) + (c1-stable(ij))*(c2*LOG(arg1))
!
        rd(ij) = rdn(ij)/(c1+rdn(ij)/vonkar*(alz(ij)-psimh))
        rh(ij) = rhn(ij)/(c1+rhn(ij)/vonkar*(alz(ij)-psixh(ij)))
! update ustar, tstar, qstar using updated, shifted coeffs 
        re(ij) = ren(ij)/(c1+ren(ij)/vonkar*(alz(ij)-psixh(ij)))
!
        ustar(ij) = rd(ij)*vmag(ij)
        tstar(ij) = rh(ij)*delt(i, j)
        qstar(ij) = re(ij)*delq(i, j)
      END DO
    END DO
!
! ij
! end iteration
!
!
    DO ij=1,icells
      i = indxi(ij)
!------------------------------------------------------------
! coefficients for turbulent flux calculation
!------------------------------------------------------------
! add windless coefficient for sensible heat flux
! as in Jordan et al (JGR, 1999)
      j = indxj(ij)
!
!
!ccsm      shcoef(i,j) = rhoa(i,j) * ustar(ij) * cp(ij) * rh(ij)
      shcoef(i, j) = rhoa(i, j)*ustar(ij)*cp(ij)*rh(ij) + c1
!------------------------------------------------------------
! momentum flux
!------------------------------------------------------------
! tau = rhoa(i,j) * ustar * ustar 
! strx = tau * uatm(i,j) / vmag 
! stry = tau * vatm(i,j) / vmag 
!------------------------------------------------------------
      lhcoef(i, j) = rhoa(i, j)*ustar(ij)*lheat*re(ij)
!
!
! not the stress at zlvl(i,j)
      tau = rhoa(i, j)*ustar(ij)*rd(ij)
      strx(i, j) = tau*uatm(i, j)
!------------------------------------------------------------
! Compute diagnostics: 2m ref T & Q
!------------------------------------------------------------
      stry(i, j) = tau*vatm(i, j)
!
      hol(ij) = hol(ij)*ztrf/zlvl(i, j)
      IF (c1 - c16*hol(ij) .GE. 0.) THEN
        abs2 = c1 - c16*hol(ij)
      ELSE
        abs2 = -(c1-c16*hol(ij))
      END IF
      y1 = SQRT(abs2)
      IF (c1 .LT. y1) THEN
        xqq = y1
      ELSE
        xqq = c1
      END IF
      xqq = SQRT(xqq)
      arg1 = (c1+xqq*xqq)/c2
      psix2 = -(c5*hol(ij)*stable(ij)) + (c1-stable(ij))*(c2*LOG(arg1))
      fac = rh(ij)/vonkar*(alz(ij)+al2-psixh(ij)+psix2)
      trf(i, j) = pott(i, j) - delt(i, j)*fac
! pot temp to temp correction
      trf(i, j) = trf(i, j) - 0.01*ztrf
      fac = re(ij)/vonkar*(alz(ij)+al2-psixh(ij)+psix2)
      qrf(i, j) = qa(i, j) - delq(i, j)*fac
    END DO
  END SUBROUTINE ATMO_BOUNDARY_LAYER
END MODULE ICE_ATMO