source: trunk/src/cor30a.pro @ 5

;----------------------------------------------------------------------------------------------
function cor30a, u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave
; 
; COARE v3 algorithm to compute fluxes
;       version with shortened iteration  modified Rt and Rq
;       uses wave information wave period in s and wave ht in m
;       no wave, standard coare 2.6 charnock:  jwave=0 
;       Oost et al.  zo=50/2/pi L (u*/c)^4.5 if jwave=1
;       taylor and yelland  zo=1200 h*(L/h)^4.5 jwave=2
;
; u         ;wind speed (m/s)  at height zu (m)
; us        ;surface current speed in the wind direction (m/s)
; ts        ;bulk water temperature (C) if jcool=1, interface water T if jcool=0
; t         ;bulk air temperature (C), height zt
; Qs        ;bulk water spec hum (g/kg) if jcool=1, ...
; Q         ;bulk air spec hum (g/kg), height zq
; Rs        ;downward solar flux (W/m^2)    (modified because of cool skin)
; Rl        ;downard IR flux (W/m^2)        (modified because of cool skin)
; rain      ;rain rate (mm/hr)
; zi        ;PBL depth (m)
; P         ;Atmos surface pressure (mb)
; zu        ;wind speed measurement height (m)
; zt        ;air T measurement height (m)
; zq        ;air q measurement height (m)
; lat       ;latitude (deg, N=+)
; jcool     ;implement cool calculation skin switch, 0=no, 1=yes
; jwave     ;implement wave dependent roughness model
; twave     ;wave period (s)
; hwave     ;wave height (m)
;
;Example of values:
;         u  us  ts  ta   qs   qa   Qsw IRd   r  pbl  Ps   zu   zt  zq lat    
;x=cor30a(5.5,0,28.7,27.2,24.2,18.5,141.,419.,0.,600.,1010.,15.,15.,15.,0.,1,1,5.,1.) ;sample data stream
;
;Result with these sample values with Matlab code:
;     8.64830      101.640    0.0352910  2.17780e-05  0.000115000  0.000115000 
;    -29.5800     0.175430   -0.0423670 -0.000205610     0.250950  0.000351300
; 0.000969740      0.00000  8.11390e-05  0.000997340   0.00121410   0.00121400 
; 0.000941350   0.00107910   0.00107910     0.780060
;
;Result obtained with this idl routine:
;     8.64829      101.640    0.0352913  2.17780e-05  0.000115000  0.000115000 
;    -29.5802     0.175432   -0.0423667 -0.000205610     0.250948  0.000351304
; 0.000969737      0.00000  8.11394e-05  0.000997343   0.00121407   0.00121400
; 0.000941351   0.00107908   0.00107908     0.780056
;
; Maximum error on any parameter: .002 %    validated!
;
;----------------------------------------------------------------------------------------------

Qs=Qs/1000.
Q=Q/1000.

;***********   set constants *************
pi=!pi
Beta=1.2 
von=0.4 
fdg=1.00 
tdk=273.16 
;grav=grv(lat) 
grav=9.8
;*************  air constants ************
Rgas=287.1 
LLe=(2.501-.00237*ts)*1e6 
cpa=1004.67 
cpv=cpa*(1+0.84*Q) 
rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q)) 
visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t) 
;************  cool skin constants  *******
Al=2.1e-5*(ts+3.2)^0.79 
be=0.026 
cpw=4000 
rhow=1022 
visw=1e-6 
tcw=0.6 
bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa) 
wetc=0.622*LLe*Qs/(Rgas*(ts+tdk)^2) 
     
;***************   wave parameters  *********
lwave=grav/2/pi*twave^2 
cwave=grav/2/pi*twave 
     
;**************  compute aux stuff *******
Rns=Rs*.945 
Rnl=0.97*(5.67e-8*(ts-0.3*jcool+tdk)^4-Rl) 



;***************   Begin bulk loop *******
     
;***************  first guess ************
du=u-us 
dt=ts-t-.0098*zt 
dq=Qs-Q 
ta=t+tdk 
ug=0.5 
dter=0.3  
dqer=wetc*dter 
ut=sqrt(du*du+ug*ug) 
u10=ut*alog(10./1e-4)/alog(zu/1e-4) 
usr=.035*u10 
zo10=0.011*usr*usr/grav+0.11*visa/usr 
Cd10=(von/alog(10./zo10))^2 
Ch10=0.00115 
Ct10=Ch10/sqrt(Cd10) 
zot10=10/exp(von/Ct10) 
Cd=(von/alog(zu/zo10))^2 
Ct=von/alog(zt/zot10) 
CC=von*Ct/Cd 
Ribcu=-zu/zi/.004/Beta^3 
Ribu=-grav*zu/ta*((dt-dter*jcool)+.61*ta*dq)/ut^2 
nits=3 
;;if (Ribu le 0.) then begin
;;  zetu=CC*Ribu/(1+Ribu/Ribcu) 
;;endif else begin
;;  zetu=CC*Ribu*(1+27./9*Ribu/CC) 
;;endelse 
sw=(Ribu le 0.)
zetu=sw*(CC*Ribu/(1+Ribu/Ribcu))+(1-sw)*(CC*Ribu*(1+27./9*Ribu/CC))
;;
L10=zu/zetu 
;;if (zetu gt 50 ) then nits=1 
usr=ut*von/(alog(zu/zo10)-psiu(zu/L10)) 
tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot10)-psit(zt/L10)) 
qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zot10)-psit(zq/L10)) 

tkt=.001 
        
;;charn=0.011 
;;if (ut gt 10.) then charn=0.011+(ut-10)/(18.-10)*(0.018-0.011) 
;;if (ut gt 18.) then charn=0.018 
charn=(((0.011+(ut-10)/(18.-10)*(0.018-0.011)) > .011) < .018)
;;
   
;***************  bulk loop ************
for i=1,nits do begin
  zet=von*grav*zu/ta*(tsr*(1+0.61*Q)+.61*ta*qsr)/(usr*usr)/(1+0.61*Q) 
  case jwave of
    0: zo=charn*usr*usr/grav+0.11*visa/usr 
    1: zo=50./2/pi*lwave*(usr/cwave)^4.5+0.11*visa/usr ;Oost et al
    2: zo=1200*hwave*(hwave/lwave)^4.5+0.11*visa/usr  ;Taylor and Yelland
  endcase
  rr=zo*usr/visa 
  L=zu/zet 
;;zoq=min([1.15e-4,5.5e-5/rr^.6]) 
  zoq=(5.5e-5/rr^.6 < 1.15e-4)
;;
  zot=zoq 
  usr=ut*von/(alog(zu/zo)-psiu(zu/L)) 
  tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot)-psit(zt/L)) 
  qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zoq)-psit(zq/L)) 
  Bf=-grav/ta*usr*(tsr+.61*ta*qsr) 
;;if (Bf gt 0) then begin
;;  ug=Beta*(Bf*zi)^.333 
;;endif else begin
;;  ug=.2 
;;endelse 
  sw=(Bf gt 0)
  ug=sw*(Beta*(Bf*zi)^.333)+(1-sw)*.2
;;
  ut=sqrt(du*du+ug*ug) 
  Rnl=0.97*(5.67e-8*(ts-dter*jcool+tdk)^4-Rl) 
  hsb=-rhoa*cpa*usr*tsr 
  hlb=-rhoa*LLe*usr*qsr 
  qout=Rnl+hsb+hlb 
  dels=Rns*(.065+11*tkt-6.6e-5/tkt*(1-exp(-tkt/8.0e-4)))        ; Eq.16 Shortwave
  qcol=qout-dels 
  alq=Al*qcol+be*hlb*cpw/LLe                                    ; Eq. 7 Buoy flux water

;;  if (alq gt 0) then begin
;;    xlamx=6./(1+(bigc*alq/usr^4)^.75)^.333                    ; Eq 13 Saunders
;;    tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr)                      ;Eq.11 Sub. thk
;;  endif else begin
;;    xlamx=6.0 
;;;;  tkt=min([.01,xlamx*visw/(sqrt(rhoa/rhow)*usr)])           ;Eq.11 Sub. thk
;;    tkt=(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01)
;;;;
;;  endelse 
  sw=(alq gt 0)
  xlamx=sw*(6./(1+(bigc*alq/usr^4)^.75)^.333)+(1-sw)*6.0
  tkt=sw*(xlamx*visw/(sqrt(rhoa/rhow)*usr))+(1-sw)*(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01)
;;
     
  dter=qcol*tkt/tcw ;  Eq.12 Cool skin
  dqer=wetc*dter 
     
endfor      ;bulk iter loop

tau=rhoa*usr*usr*du/ut                 ;stress
hsb=-rhoa*cpa*usr*tsr                  ;sens
hlb=-rhoa*LLe*usr*qsr                  ;lat
                                       ; net solar Rns
                                       ; net lw    Rnl

;****************   rain heat flux ********
dwat=2.11e-5*((t+tdk)/tdk)^1.94                         ;! water vapour diffusivity
dtmp=(1.+3.309e-3*t-1.44e-6*t*t)*0.02411/(rhoa*cpa)     ;!heat diffusivity
alfac= 1/(1+(wetc*LLe*dwat)/(cpa*dtmp))                 ;! wet bulb factor
RF= rain*alfac*cpw*((ts-t-dter*jcool)+(Qs-Q-dqer*jcool)*LLe/cpa)/3600.


;y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau]]

;;****************   Webb et al. correection  ************
;wbar=1.61*hlb/LLe/(1+1.61*Q)/rhoa+hsb/rhoa/cpa/ta      ;formulation in hlb already includes webb
;hl_webb=rhoa*wbar*Q*LLe 
;;**************   compute transfer coeffs relative to ut @meas. ht **********
;;Cd=tau/rhoa/ut/max([.1,du])
;Cd=tau/rhoa/ut/(du > .1) 
;;;
Ch=-usr*tsr/ut/(dt-dter*jcool) 
Ce=-usr*qsr/(dq-dqer*jcool)/ut 
;;************  10-m neutral coeff realtive to ut ********
;Cdn_10=von*von/alog(10./zo)/alog(10./zo) 
;Chn_10=von*von*fdg/alog(10./zo)/alog(10./zot) 
;Cen_10=von*von*fdg/alog(10./zo)/alog(10./zoq) 

y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau],[Ch],[Ce]]   
;;y=[hsb,hlb,tau,zo,zot,zoq,L,usr,tsr,qsr,dter,dqer,tkt,RF,wbar,Cd,Ch,Ce,Cdn_10,Chn_10,Cen_10,ug ] 
;   1   2   3   4  5   6  7  8   9  10   11   12  13  14  15  16 17 18    19      20    21  22
;       hsb=                    sensible heat flux (w/m^2)
;       hlb=                    latent heat flux (w/m^2)
;       RF=                     rain heat flux(w/m^2)
;       wbar=                   webb mean w (m/s)
;       tau=                    stress (nt/m^2)
;       zo=                     velocity roughness length (m)
;       zot                     temperature roughness length (m)
;       zoq=                    moisture roughness length (m)
;       L=                      Monin_Obukhov stability length
;       usr=                    turbulent friction velocity (m/s), including gustiness
;       tsr                     temperature scaling parameter (K)
;       qsr                     humidity scaling parameter (g/g)
;       dter=                   cool skin temperature depression (K)
;       dqer=                   cool skin humidity depression (g/g)
;       tkt=                    cool skin thickness (m)
;       Cd=                     velocity drag coefficient at zu, referenced to u
;       Ch=                     heat transfer coefficient at zt
;       Ce=                     moisture transfer coefficient at zq
;       Cdn_10=                 10-m velocity drag coeeficient, including gustiness
;       Chn_10=                 10-m heat transfer coeeficient, including gustiness
;       Cen_10=                 10-m humidity transfer coeeficient, including gustiness



return, y

end
;----------------------------------------------------------------------------------------------