c Last Update to Flux Routines: 29 Apr 1994
c
      subroutine inertial(cu,ct,cq,zetb,tmast,zb,sonix,
     *              us_ib,f_hs_ib,f_ls_ib,us_i,f_hs_i,f_ls_i,zet,zmol_i)
     
c subroutine to calculate intertial dissipation fluxes from turbulence
c     structure functions
c
c Chris W. Fairall, Peter A. Coppin: Fortran version 1.0 29-Apr-94
c 
c INPUTS
c
c cu - Cu^2, streamwise velocity structure function (m/s)^2/m^.667
c ct - Ct^2, temperature structure function k^2/m^.667
c cq - Cq^2, absolute humidity structure function (g/m^3)^2/m^.667
c these are usually obtained from spectral analysis 
c (see Fairall and Larsen (1986) BLM 34:287-301)
c
c using the inertial subrange of the 1-dimensional variance spectrum Sx(f) where f is the frequency
c
c      Cx^2 = 4.0 * (2*Pi/U)^0.667 * f^1.667 * Sx(f)
c
c Note: U is the RELATIVE wind speed (not corrected for platform speed!)
c       if vertical velocity is used let Cu^2=3/4 * Cw^2
c       the air density has been hard wired to 1.15 g/m^3, a later version may include 
c           calculations requiring input of wet bulb temp and atmospheric pressure,
c           variations in this value have a linear effect on sensible and latent heat and
c           a minor influence on z/L, the former can obviously be corrected.
c
c zetb - the Monin-Obukov parameter z/L obtained via the bulk algorithm
c tmast - the air temperature at the height of measurement
c zb - the height of measurement in m
c sonix - a parameter: 1 if temperature is speed of sound derived from a sonic, otherwise=0
c 
c OUTPUTS
c
c us_ib - ustar using bulk z/L
c f_hs_ib - sensible heat flux (w/m^2) using bulk z/L
c f_ls_ib - latent heat flux (w/m^2) using bulk z/L
c
c us_i - ustar using inertial calculated z/L
c f_hs_i - sensible heat flux (w/m^2) using inertial calculated z/L
c f_ls_i - latent heat flux (w/m^2) using inertial calculated z/L
c
c zet - inertial calculated z/L
c zmol_i - inertial calculated L
c




      tair=tmast+273.16
     
c constants

      vkon=0.4                 !Von Karman
      cp=1004.67              !spec heat dry air
      c_le=2500.-2.274*tmast     !Latent heat vaporization of H2O
      g=9.72                  !gravitational constant
      rho_air=1.15            !COARE air density
      
      IF (zetb .eq. -9999.) then
        zet=0                 !compute neutral fluxes
      ELSE
        zet=zetb
      END IF
      
      call stabcalc(zet,psi_fu,psi_ft)
          
      xfu=psi_fu
      xft=psi_ft

      us_ib=SQRT(cu*zb**0.667/psi_fu)
      ts_ib=SIGN(1.,zetb)*SQRT(ct*zb**0.667/psi_ft)
      qs_ib=-SQRT(cq*zb**0.667/psi_ft)         !assumes Qsea>Qair
      ts_ib=ts_ib-sonix*0.00051*tair*qs_ib/rho_air !sonic T moisture correction
c
      f_hs_ib=-rho_air*cp*us_ib*ts_ib
      f_ls_ib=-c_le*us_ib*qs_ib
c
c use bulk estimates to seed iteration for measured zet            
c
      us0=us_ib
      ts0=ts_ib
      qs0=qs_ib
c
c iterate 3 times
c
      DO 20 i=1,3
        zet=vkon*zb*g/tair*(ts0+0.00061*tair*qs0/rho_air)/(us0*us0)
        call stabcalc(zet,psi_fu,psi_ft)              
c
        us0=SQRT(cu*zb**0.667/psi_fu)
        ts0=SIGN(1.,zetb)*SQRT(ct*zb**0.667/psi_ft)
        qs0=-SQRT(cq*zb**0.667/psi_ft)         !assumes Qsea>Qair
        ts0=ts0-sonix*0.00051*tair*qs0/rho_air !sonic T moisture correction
      
20    CONTINUE      

      qs_i=qs0
      ts_i=ts0
      us_i=us0
      f_hs_i=-rho_air*cp*us_i*ts_i
      f_ls_i=-c_le*us_i*qs_i
      zmol_i=zb/zet
      
      
      RETURN
      END
      
c -----------------------------------------------------------------------
      
      SUBROUTINE stabcalc(zet,psi_fu,psi_ft)

c subroutine to calculate psi functions
      
      IF (zet .lt. 0.0) then
      
c  psi_fu from J.Edson 3/94
        psi_fu=4.0*(1./(1.-20.*zet)**0.33-1.0*zet-1/(7.-zet))**0.667
        psi_ft=5.0/(1.0-6.2*zet)**0.667
        
      ELSE
      
        psi_fu=4.0*(1+2.5*zet**0.667)
        psi_ft=5.0*(1+2.5*zet**0.667)

      END IF
      
      RETURN
      END

c -------------------------------------------------------------------------