source: trunk/src/cor30a.pro @ 152

;+
;
; ==========
; cor30a.pro
; ==========
;
; .. function:: cor30a(u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave)
;
; DESCRIPTION
; ===========
;
; 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
;
;     :param u: wind speed (m/s)  at height zu (m)
;     :param us: surface current speed in the wind direction (m/s)
;     :param ts: bulk water temperature (C) if jcool=1, interface water T if jcool=0
;     :param t: bulk air temperature (C), height zt
;     :param Qs: bulk water spec hum (g/kg) if jcool=1, ...
;     :param Q: bulk air spec hum (g/kg), height zq
;     :param Rs: downward solar flux (W/m^2)    (modified because of cool skin)
;     :param Rl: downard IR flux (W/m^2)        (modified because of cool skin)
;     :param rain: rain rate (mm/hr)
;     :param zi: PBL depth (m)
;     :param P: Atmos surface pressure (mb)
;     :param zu: wind speed measurement height (m)
;     :param zt: air T measurement height (m)
;     :param zq: air q measurement height (m)
;     :param lat: latitude (deg, N=+)
;     :param jcool: implement cool calculation skin switch, 0=no, 1=yes
;     :param jwave: implement wave dependent roughness model
;     :param twave: wave period (s)
;     :param hwave: wave height (m)
;
; EXAMPLES
; ========
;
; ::
;
;                  u  us  ts  ta   qs   qa   Qsw IRd   r  pbl  Ps   zu   zt  zq lat
;
;    IDL> 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.)
;
; 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!
;
; SEE ALSO
; ========
;
; used by :ref:`TropFlux_19890101_20091231.pro`
;
; call :func:`psiu`, :func:`psit`.
;
; TODO
; ====
;
; coding rules
;
; check for the new module provide by pk 20110811 by mail to fp
; "You can use this coare algorithm in your depositary.  It is same as the one you
; have, but with some functions activated, which will be useful in the next stage."
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20110830T085907Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add ug
;
; - fplod 20110822T085404Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add ref to psiu and psit called her and provided by pk today
;
; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
; - pbk 2008
;
;   * creation
;
;-
function cor30a, u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave

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],[ug]]
;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