#include <misc.h>#include <params.h>subroutine flxoce(indx    ,npts    ,pmidm1  ,ubot    ,vbot    , &                  tbot    ,qbot    ,thbot   ,zbot    ,ts      , &                  ltheat  ,shf     ,lhf     ,taux    ,tauy    , &                  lwup    ,tref    )!----------------------------------------------------------------------- ! ! Purpose: ! Compute ocean to atmosphere surface fluxes of sensible, latent heat! and stress components:!! Method: ! Assume:!   1) Neutral 10m drag coeff: !         cdn = .0027/U10N + .000142 + .0000764 U10N!   2) Neutral 10m stanton number: !         ctn = .0327 sqrt(cdn), unstable!         ctn = .0180 sqrt(cdn), stable!   3) Neutral 10m dalton number:  !         cen = .0346 sqrt(cdn)!   4) The saturation humidity of air at T(K): !         qsat(T)  (kg/m^3)! Note:!   1) here, tstar = <WT>/U*, and qstar = <WQ>/U*.!   2) wind speeds should all be above a minimum speed (umin)! ! Author: Bill Large/M.Vertenstein, Sep. 1995! !-----------------------------------------------------------------------  use precision  use ppgrid  use physconst, only: rair, cpair, cpvir, zvir, gravit, stebol  implicit none#include <parpbl.h>!------------------------------Arguments--------------------------------  integer , intent(in) :: indx(pcols)    ! column index array (land)  integer , intent(in) :: npts           ! Number of land points  real(r8), intent(in) :: pmidm1(pcols)  ! Bottom level pressure  real(r8), intent(in) :: ubot(pcols)    ! Bottom level u wind  real(r8), intent(in) :: vbot(pcols)    ! Bottom level v wind  real(r8), intent(in) :: tbot(pcols)    ! Bottom level temperature  real(r8), intent(in) :: qbot(pcols)    ! Bottom level specific humidity  real(r8), intent(in) :: thbot(pcols)   ! Bottom level potential temperature  real(r8), intent(in) :: zbot(pcols)    ! Bottom level height above surface  real(r8), intent(in) :: ts(pcols)      ! Surface temperature  real(r8), intent(in) :: ltheat(pcols)  ! Latent heat for given srf conditions  real(r8), intent(inout) :: shf(pcols)  ! Initial sensible heat flux (W/m2)  real(r8), intent(inout) :: lhf(pcols)  ! Initial latent heat flux (W/m2)  real(r8), intent(inout) :: taux(pcols) ! X surface stress (N/m2)  real(r8), intent(inout) :: tauy(pcols) ! Y surface stress (N/m2)  real(r8), intent(inout) :: lwup(pcols) ! Longwave up flux at surface (W/m2)  real(r8), intent(inout) :: tref(pcols) ! 2m reference temperature !-----------------------------------------------------------------------!---------------------------Local variables-----------------------------  integer i,ii            ! column indices  real(r8) ssq            ! Surface saturation specific humidity  real(r8) ustar          ! ustar  real(r8) tstar          ! tstar  real(r8) qstar          ! qstar  real(r8) u10n           ! neutral 10 m wind speed over ocean  real(r8) vmag           ! Surface wind magnitude  real(r8) thvbot         ! Bottom lev virtual potential temp  real(r8) delt           ! potential T difference (K)  real(r8) delq           ! specific humidity difference (kg/kg)  real(r8) rdn            ! sqrt of neutral exchange coeff (momentum)  real(r8) rhn            ! sqrt of neutral exchange coeff (heat)  real(r8) ren            ! sqrt of neutral exchange coeff (tracers)            real(r8) rd             ! sqrt of exchange coeff (momentum)  real(r8) rh             ! sqrt of exchange coeff (heat)  real(r8) re             ! sqrt of exchange coeff (tracers)  real(r8) hol            ! Ref hgt (10m) / monin-obukhov length  real(r8) xsq            ! Temporary variable  real(r8) xqq            ! Temporary variable  real(r8) alz            ! ln(zbot/z10)  real(r8) cp             ! Specific heat of moist air  real(r8) tau            ! Reference height stress  real(r8) psimh          ! Stability funct at ref lev (momentum)  real(r8) psixh          ! Stability funct at ref lev (heat & tracers)   real(r8) stable         ! Stability factor  real(r8) rbot(pcols)    ! Density at bottom model level  real(r8) bn             ! exchange coef funct for interpolation  real(r8) bh             ! exchange coef funct for interpolation  real(r8) fac            ! interpolation factor  real(r8) ln0            ! log factor for interpolation  real(r8) ln3            ! log factor for interpolation  real(r8), parameter :: ztref =2.0 ! reference height for air temperature!-----------------------------------------------------------------------!--------------------------Statement functions--------------------------  real(r8) psimhu         ! Unstable part of psimh  real(r8) psixhu         ! Unstable part of psixh  real(r8) qsat           ! Saturation specific humidty of air  real(r8) cdn            ! Neutral drag coeff at bottom model level  real(r8) xd             ! Dummy argument  real(r8) Tk             ! Temperature (K)  real(r8) Umps           ! Wind velocity (m/sec)!  qsat(Tk)   = 640380. / exp(5107.4/Tk)  cdn(Umps)  = 0.0027 / Umps + .000142 + .0000764 * Umps  psimhu(xd) = log((1.+xd*(2.+xd))*(1.+xd*xd)/8.) - 2.*atan(xd) + 1.571  psixhu(xd) = 2. * log((1. + xd*xd)/2.)!-----------------------------------------------------------------------!! Loop over ocean points!  do ii=1,npts     i = indx(ii)!!---------------------------------------------------------------! Set up necessary variables!---------------------------------------------------------------!     rbot(i)= pmidm1(i) / (rair*tbot(i))     vmag   = max(umin, sqrt(ubot(i)**2+vbot(i)**2))     thvbot = thbot(i) * (1.0 + zvir*qbot(i))     ssq    = 0.98 * qsat(ts(i)) / rbot(i)     delt   = thbot(i) - ts(i)     delq   = qbot(i) - ssq       alz    = log(zbot(i)/zref)      cp     = cpair*(1. + cpvir*ssq) !!---------------------------------------------------------------! First iteration to converge on Z/L and hence the fluxes!---------------------------------------------------------------!! Initial guess for roots of neutral exchange coefficients, ! assume z/L=0. and u10n is approximated by vmag.! Stable if (thbot > ts ).!     stable = 0.5 + sign(0.5_r8 , delt)     rdn  = sqrt(cdn(vmag))     rhn  = (1.-stable) * 0.0327 + stable * 0.018      ren  = 0.0346 !! Initial guess of ustar,tstar and qstar!     ustar = rdn*vmag     tstar = rhn*delt     qstar = ren*delq!! Compute stability and evaluate all stability functions! Stable if (thbot > ts or hol > 0 )!     hol = xkar *gravit *zbot(i) * (tstar/thvbot + qstar/(1./zvir+qbot(i))) / ustar**2     hol = sign( min(abs(hol),10._r8), hol )     stable = 0.5 + sign(0.5_r8 , hol)     xsq   = max(sqrt(abs(1. - 16.*hol)) , 1._r8)     xqq   = sqrt(xsq)     psimh = -5. * hol * stable + (1.-stable)*psimhu(xqq)     psixh = -5. * hol * stable + (1.-stable)*psixhu(xqq)!! Shift 10m neutral wind speed using old rdn coefficient!     rd   = rdn / (1.+rdn/xkar*(alz-psimh))     u10n = vmag * rd / rdn!! Update the neutral transfer coefficients at 10m and neutral stability!     rdn = sqrt(cdn(u10n))     ren = 0.0346     rhn = (1.-stable) * 0.0327 + stable * 0.018 !! Shift all coeffs to measurement height and stability!     rd = rdn / (1.+rdn/xkar*(alz-psimh))      rh = rhn / (1.+rhn/xkar*(alz-psixh))      re = ren / (1.+ren/xkar*(alz-psixh))!! Update ustar, tstar, qstar using updated, shifted coeffs !     ustar = rd * vmag      tstar = rh * delt      qstar = re * delq !!---------------------------------------------------------------! Second iteration to converge on Z/L and hence the fluxes!---------------------------------------------------------------!! Recompute stability & evaluate all stability functions  ! Stable if (thbot > ts or hol > 0 )!      hol = xkar * gravit * zbot(i) * (tstar/thvbot + qstar/(1./zvir+qbot(i))) / ustar**2     hol = sign( min(abs(hol),10._r8), hol )     stable = 0.5 + sign(0.5_r8 , hol)     xsq   = max(sqrt(abs(1. - 16.*hol)) , 1._r8)     xqq   = sqrt(xsq)     psimh = -5. * hol * stable + (1.-stable)*psimhu(xqq)     psixh = -5. * hol * stable + (1.-stable)*psixhu(xqq)!! Shift 10m neutral wind speed using old rdn coefficient!     rd   = rdn / (1.+rdn/xkar*(alz-psimh))     u10n = vmag * rd / rdn!! Update the neutral transfer coefficients at 10m and neutral stability!     rdn = sqrt(cdn(u10n))     ren = 0.0346     rhn = (1.-stable) * 0.0327 + stable * 0.018 !! Shift all coeffs to measurement height and stability     rd = rdn / (1.+rdn/xkar*(alz-psimh))      rh = rhn / (1.+rhn/xkar*(alz-psixh))      re = ren / (1.+ren/xkar*(alz-psixh))!!---------------------------------------------------------------! Compute the fluxes!---------------------------------------------------------------!! Update ustar, tstar, qstar using updated, shifted coeffs !     ustar = rd * vmag      tstar = rh * delt      qstar = re * delq !! Compute surface stress components!     tau     =  rbot(i) * ustar * ustar      taux(i) = -tau * ubot(i) / vmag      tauy(i) = -tau * vbot(i) / vmag !! Compute heat flux components at current surface temperature! (Define positive latent and sensible heat as upwards into the atm)!     shf(i) = -cp * tau * tstar / ustar      lhf(i) = -ltheat(i) * tau * qstar / ustar     lwup(i) = stebol * ts(i)**4 !!---------------------------------------------------------------! Following Geleyn(1988), interpolate ts to fixed height zref!---------------------------------------------------------------!! Compute function of exchange coefficients. Assume that ! cn = rdn*rdn, cm=rd*rd and ch=rh*rd, and therefore ! 1/sqrt(cn(i))=1/rdn and sqrt(cm(i))/ch(i)=1/rh !     bn = xkar/rdn     bh = xkar/rh!! Interpolation factor for stable and unstable cases!     ln0 = log(1.0 + (ztref/zbot(i))*(exp(bn) - 1.0))     ln3 = log(1.0 + (ztref/zbot(i))*(exp(bn - bh) - 1.0))     fac = (ln0 - ztref/zbot(i)*(bn - bh))/bh * stable + (ln0 - ln3)/bh * (1.-stable)     fac = min(max(fac,0._r8),1._r8)!! Actual interpolation!     tref(i) = ts(i) + (tbot(i) - ts(i))*fac  end do!  returnend subroutine flxoce