;
pro osc3
;
phycon
;
sav_filename = 'osc3.sav'
;
defsysv,   '!gamma_l', 2d-6
defsysv,         '!E', 1.5d-3
defsysv,   '!gamma_d', !g / !Cpd
defsysv,   '!gamma_A', 4.d-3
defsysv,      '!TkAo', 283.15d
defsysv,    '!wtot_A', 8d-3
defsysv,        '!dt', 10.d
defsysv,        '!ho', 1000.d
;
;..initialization of cloud properties
;
wl       = !gamma_l*!ho
TkA      = !TkAo
;..initialization of parcel
;
TkP      = TkA - 1.0d ; parcel negatively buoyant at start
dTemp    = TkP - TkA
UP       = 0.d        ; initially, parcel has no vertical velocity
zP       = !ho
wP       = 0.5 * (!wtot_A - wl)
;
;..initialize time
;
t        = 0d
;
;.."flag" is used to control the storage of calculation results
;
flag     = 0
;
while t le 3600 do begin
 print, 'here 1 ', t, dTemp, UP, TkA, zP, wl, wP, !wtot_A - wl, format='(a, 8f15.5)'
;
;..store
;
 if flag eq 0 then begin
  aTkP   = [TkP]
  aUP    = [UP]
  aTkA   = [TkA]
  azP    = [zP]
  awl    = [wl]
  awP    = [wP]
  asec   = [t]
  dydt   = osc3_derivative(t,[dTemp, UP, TkA, zP, wl, wP])
  adUpdt = [dydt[1]]
  flag   = 1
 endif else begin
  aTkP   = [aTkP, dTemp + TkA]
  aUP    = [aUP, UP]
  aTkA   = [aTkA, TkA]
  azP    = [azP, zP]
  awl    = [awl, wl]
  awP    = [awP, wP]
  asec   = [asec, t]
  dydt   = osc3_derivative(t,[dTemp, UP, TkA, zP, wl, wP])
  adUpdt = [adUpdt, dydt[1]]
 endelse
;
;..integrate
;
 result   = RK4([dTemp, UP, TkA, zP, wl, wP], osc3_derivative(t,[dTemp, UP, TkA, zP, wl, wP]), t, !dt, 'osc3_derivative', /double)
;
;..update
;
 dtemp    = result[0]
 TkP      = dTemp + result[2]
 UP       = result[1]
 TkA      = result[2]
 zP       = result[3]
 wl       = result[4]
 wP       = result[5]
 t        = t + !dt
;
; print, 'here 2 ', t, dTemp, UP, TkA, zP, wl, format='(a, 6f15.4)'
; wait, 5
;
 if zp le 0 then goto, done
;
endwhile
;
done:
;
gamma_A = !gamma_A
TkAo    = !TkAo
gamma_l = !gamma_l
wtot_A  = !wtot_A
ho      = !ho
save, filename = sav_filename, aTkP, aUP, aTkA, azP, awl, awP, asec, adUpdt, gamma_A, TkAo, gamma_l, wtot_A, ho
stop
;
end