;
pro osc1
;
phycon
;
sav_filename = 'osc1.sav'
;
defsysv,   '!gamma_d', !g / !Cpd
defsysv,   '!gamma_A', 4.d-3
defsysv,      '!TkAo', 283.15d
defsysv,        '!dt', 10.d
defsysv,        '!ho', 1000.d
;
print, sqrt((!g/!TkAo)*(!gamma_d - !gamma_A))
;
;..initialization of atmosphere
;
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
;
;..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, format='(a, 5f15.5)'
;
;..store
;
 if flag eq 0 then begin
  aTkP   = [TkP]
  aUP    = [UP]
  aTkA   = [TkA]
  azP    = [zP]
  asec   = [t]
  dydt   = osc1_derivative(t,[dTemp, UP, TkA, zP])
  adUpdt = [dydt[1]]
  flag   = 1
 endif else begin
  aTkP   = [aTkP, dTemp + TkA]
  aUP    = [aUP, UP]
  aTkA   = [aTkA, TkA]
  azP    = [azP, zP]
  asec   = [asec, t]
  dydt   = osc1_derivative(t,[dTemp, UP, TkA, zP])
  adUpdt = [adUpdt, dydt[1]]
 endelse
;
;..integrate
;
 result   = RK4([dTemp, UP, TkA, zP], osc1_derivative(t,[dTemp, UP, TkA, zP]), t, !dt, 'osc1_derivative', /double)
;
;..update
;
 dtemp    = result[0]
 TkP      = dTemp + result[2]
 UP       = result[1]
 TkA      = result[2]
 zP       = result[3]
 t        = t + !dt
;
; print, 'here 2 ', t, dTemp, UP, TkA, zP, format='(a, 6f15.4)'
; wait, 5
;
 if zp le 0 then goto, done
;
endwhile
;
done:
;
gamma_A = !gamma_A
TkAo    = !TkAo
ho      = !ho
save, filename = sav_filename, aTkP, aUP, aTkA, azP, asec, adUpdt, gamma_A, TkAo, ho
stop
;
end