c swd -- bulkf formula used in bulkf and ice pkgs

c taken from exf package 
#include "CPP_OPTIONS.h"
      subroutine bulkf_formula_lanl(
     I                           uw, vw, us, Ta, Qa, nc, tsf_in,
     I                           flwup, flha, fsha, df0dT,
     I                           ust, vst, evp, ssq, iceornot, windread
     &                           )

       IMPLICIT NONE
c
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "FFIELDS.h"
#ifdef ALLOW_BULK_FORCE
#include "BULKF_ICE_CONSTANTS.h"
#include "BULKF.h"
#endif

c
c calculate bulk formula fluxes over open ocean or seaice
c
c input
      _RL us                 ! wind speed
      _RL uw                 ! zonal wind speed (at grid center)
      _RL vw                 ! meridional wind speed (at grid center)
      _RL Ta                 ! air temperature at ht
      _RL Qa                 ! specific humidity at ht
      _RL tsf_in             ! surface temperature (either ice or open water)
      _RL nc                 ! fraction cloud cover
      integer iceornot       ! 0=open water, 1=ice cover
      logical windread       !
c output
      _RL flwup              ! upward long wave radiation
      _RL flha                ! latent heat flux
      _RL fsha                ! sensible heat flux
      _RL df0dT              ! derivative of heat flux with respect to Tsf
      _RL ust                ! zonal wind stress (at grid center)
      _RL vst                ! meridional wind stress (at grid center)
      _RL evp               ! evaporation rate (over open water)
      _RL ssq                ! surface specific humidity (kg/kg)
c
c local variables  
      _RL tsf                ! surface temperature in K
      _RL ht                 ! reference height (m)
      _RL hq                 ! reference height for humidity (m)
      _RL hu                 ! reference height for wind speed (m)
      _RL zref               ! reference height
      _RL czol
      _RL usm                ! wind speed limited
      _RL t0                 ! virtual temperature (K)
      _RL deltap             ! potential temperature diff (K)
      _RL delq               ! specific humidity difference
      _RL stable
      _RL rdn   ,ren, rhn 
      _RL ustar
      _RL tstar
      _RL qstar
      _RL huol
      _RL htol
      _RL hqol
      _RL xsq
      _RL x
      _RL re
      _RL rh
      _RL tau
      _RL psimh
      _RL psixh
      _RL rd
      _RL uzn
      _RL shn
      _RL aln
      _RL cdalton
      _RL dflhdT             ! derivative of latent heat with respect to T
      _RL dfshdT             ! derivative of sensible heat with respect to T
      _RL dflwupdT           ! derivative of long wave with respect to T
      _RL mixratio
      _RL ea
      _RL psim_fac
      _RL umin
      _RL lath
      _RL csha
      _RL clha
      _RL zice
      integer niter_bulk, iter

#ifdef ALLOW_BULK_FORCE

c --- external functions ---
      _RL       exf_BulkCdn
      external  exf_BulkCdn
      _RL       exf_BulkqSat
      external  exf_BulkqSat
      _RL       exf_BulkRhn
      external  exf_BulkRhn

cQQQQQQQQQQQQQQQQQQQQq
c -- compute turbulent surface fluxes
              ht =  2.d0
              hq =  2.d0
              hu = 10.d0
              zref = 10.d0
              zice = 0.0005
              aln = log(ht/zref)
              niter_bulk = 5
              cdalton = 0.0346000
              czol = zref*xkar*gravity
              psim_fac=5.d0
              umin=1.d0
c 
              lath=Lvap
              if (iceornot.gt.0) lath=Lvap+Lfresh
              Tsf=Tsf_in+Tf0kel
c     Wind speed 
              if (us.eq.0) then
                us = sqrt(uw*uw + vw*vw)
              endif
              usm = max(us,umin)
cQQQ try to limit drag
cQQ           usm = min(usm,5.d0)
c
              t0     = Ta*(1.d0 + humid_fac*Qa)
cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
         ssq    =3.797915*exp(lath*(7.93252d-6-2.166847d-3/Tsf))/p0
cBB debugging
cBB             print*,'ice, ssq',  ssq
c
              deltap = ta  - tsf + gamma_blk*ht
              delq   = Qa - ssq
c
c initialize estimate exchange coefficients
              rdn=xkar/(log(zref/zice))
              rhn=rdn
              ren=rdn
c calculate turbulent scales
              ustar=rdn*usm
              tstar=rhn*deltap
              qstar=ren*delq
c
c interation with psi-functions to find transfer coefficients
              do iter=1,niter_bulk
                 huol   = czol/ustar**2 *(tstar/t0 +
     &                    qstar/(1.d0/humid_fac+Qa))
                 huol   = sign( min(abs(huol),10.d0), huol)
                 stable = 5. _d -1 + sign(5. _d -1 , huol)
                 xsq    = max(sqrt(abs(1.d0 - 16.d0*huol)),1.d0)
                 x      = sqrt(xsq)
                 psimh = -5.d0*huol*stable + (1.d0-stable)*
     &                    (2.d0*log(5.e-1*(1.d0+x)) +
     &                    2.d0*log(5.e-1*(1.d0+xsq)) -
     &                    2.d0*atan(x) + pi*.5d0)
                 psixh  = -5.d0*huol*stable + (1.d0-stable)*
     &                     (2.d0*log(5.e-1*(1.d0+xsq)))
   
c Update the transfer coefficients 
   
                 rd = rdn/(1.d0 + rdn*(aln-psimh)/xkar)
                 rh = rhn/(1.d0 + rhn*(aln-psixh)/xkar)
                 re = rh
c  Update ustar, tstar, qstar using updated, shifted coefficients.
                 ustar = rd*usm
                 qstar = re*delq
                 tstar = rh*deltap
                enddo
c
                tau   = rhoa*ustar**2
                tau   = tau*us/usm
                csha   = rhoa*cpair*us*rh*rd   
                clha   = rhoa*lath*us*re*rd   
c
                fsha  = csha*deltap
                flha  = clha*delq
                evp   = flha/lath/rhosw
c
c up long wave radiation
cQQ           mixratio=Qa/(1-Qa)
cQQ           ea=p0*mixratio/(0.62197+mixratio)
cQQ           flwup=-0.985*stefan*tsf**4
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
              if (iceornot.eq.0) then
               flwup=-ocean_emissivity*stefan*tsf**4
              else
               if (iceornot.eq.2) then
                 flwup=-snow_emissivity*stefan*tsf**4
               else 
                 flwup=-ice_emissivity*stefan*tsf**4
               endif
              endif
#ifdef ALLOW_THERM_SEAICE
cswdcou -- remove ---
cQQ    if (iceornot.gt.0) then
cswdcou --- end remove ---
c derivatives with respect to Tsf
cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
              dflhdT = -clha*ssq*Lath*2.166847E-3/(Tsf**2)
              dfshdT = -csha
              if (iceornot.eq.0) then
               dflwupdT=-4.*ocean_emissivity*stefan*tsf**3
              else
               if (iceornot.eq.2) then
                dflwupdT=-4.*snow_emissivity*stefan*tsf**3
               else
                dflwupdT=-4.*ice_emissivity*stefan*tsf**3
               endif
              endif 
cQQ           dflwupdT=-4.*0.985*stefan*tsf**3
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
c total derivative with respect to surface temperature
              df0dT=dflwupdT+dfshdT+dflhdT
cBB
cBB     print*,'derivatives:',df0dT,dflwupdT,dfshdT,dflhdT
cswdcou -- remove ---
cQQ     endif
cswdcou
#endif /*ALLOW_THERM_SEAICE*/
c
c wind stress at center points
              if (.NOT.windread) then
                 ust = rhoa*exf_BulkCdn(usm)*us*uw
                 vst = rhoa*exf_BulkCdn(usm)*us*vw
               endif
#endif /*ALLOW_BULK_FORCE*/

      return
      end