; =========================================================================
; create idealized WRF sounding
;
; Manuela Lehner
; April 2011
; =========================================================================

; output file
outfile = '~/Desktop/input_sounding'

; vertical sounding resolution and top of the sounding atmopshere
ztop = 20000.  ;[m]
dz = 20.       ;[m]

; surface values of pressure, potential temperature, mixing ratio, u, and v
psfc = 1000.
tsfc = 290.0
msfc = 0.   ; 4.
usfc = 0.00
vsfc = 5.00

; temperature profile: specify theta at the top of each layer (linear increase or decrease)
use_n = 0    ; N^2 in tlayer
tlayer = [340]
tztop = [ztop]

;use_n = 0
;tlayer = [330., 500]
;tztop = [11000., ztop]

; mixing ratio profile: specify mixing ratio at the top of each layer (linear increase or decrease)
mlayer = [0.]
mztop = [ztop]
;mlayer = [4., 0.]
;mztop = [1000, ztop]

; wind profile: specify u and v at the top of each layer (linear increase or decrease)
;ulayer = [0., 7., 1., 1.]
;uztop = [180., 230., 280., ztop]
;uintmet = ['lin', 'sin', 'sin', 'lin']
ulayer = [0.00]
uztop = [ztop]
uintmet = ['lin']

;ulayer = [0., 7., 7., 1., 1.]
;uztop = [130., 140., 180., 230., ztop]
;uintmet = ['lin', 'sin', 'lin', 'sin', 'lin']
;ulayer = [0., 3., 1., 1.]
;uztop = [130., 180., 230., ztop]
;uintmet = ['lin', 'sin', 'sin', 'lin']

vlayer = [5.]
vztop = [ztop]
vintmet = ['lin']

; add patch information for input_sounding_1?
add_patch_info = 0
npatches = 1
xpatch_beg = [  1, 139, 164]
xpatch_end = [138, 163, 301]
ypatch_beg = [  1,   1,   1]
ypatch_end = [  5,   5,   5]


; ==================================================================================================
; START PROGRAM HERE
; ==================================================================================================

; ----- some constants -----------------------------------------------------------------------------
p0 = 1000.
Rl = 287.
Rw = 461.5
cp = 1004.
lv = 2.5e6
T0 = 273.16
es0 = 6.112e-2

; ------ calculate profiles ------------------------------------------------------------------------
nlevs = fix( ztop / dz ) + 1

; ----- height above ground -----
zmsl = findgen( nlevs ) * dz

; ----- potential temperature -----
tprof = fltarr( nlevs ) * !values.f_nan
tprof[0] = tsfc
if ( use_n eq 0 ) then begin
  for ilayer = 0, n_elements( tztop ) - 1 do begin
    if ( ilayer eq 0 ) then begin
      tzbottom = 0. 
      tbottom = tsfc
    endif else begin
      tzbottom = tztop[ilayer-1]
      tbottom = tlayer[ilayer-1]
    endelse
    ind = where( (zmsl gt tzbottom) and (zmsl le tztop[ilayer]) )
    tprof[ind] = tbottom + ( tlayer[ilayer] - tbottom ) / ( tztop[ilayer] - tzbottom ) * ( zmsl[ind] - tzbottom )
  endfor
endif else begin
  for ilev = 1, nlevs - 1 do begin
    ind = where( tztop ge zmsl[ilev] )
    ind = ind[0]
    dtheta_dz = tlayer[ind] * tprof[ilev-1] / 9.81
    tprof[ilev] = tprof[ilev-1] + dtheta_dz * ( zmsl[ilev] - zmsl[ilev-1] )
  endfor
endelse

; ---- mixing ratio ------
mprof = fltarr( nlevs ) * !values.f_nan
mprof[0] = msfc
for ilayer = 0, n_elements( mztop ) - 1 do begin
  if ( ilayer eq 0 ) then begin
    mzbottom = 0. 
    mbottom = msfc
  endif else begin
    mzbottom = mztop[ilayer-1]
    mbottom = mlayer[ilayer-1]
  endelse
  ind = where( (zmsl gt mzbottom) and (zmsl le mztop[ilayer]) )
  mprof[ind] = mbottom + ( mlayer[ilayer] - mbottom ) / ( mztop[ilayer] - mzbottom ) * ( zmsl[ind] - mzbottom )
endfor

; ----- u wind comonent -----
uprof = fltarr( nlevs ) * !values.f_nan
uprof[0] = usfc
for ilayer = 0, n_elements( uztop ) - 1 do begin
  if ( ilayer eq 0 ) then begin
    uzbottom = 0. 
    ubottom = usfc
  endif else begin
    uzbottom = uztop[ilayer-1]
    ubottom = ulayer[ilayer-1]
  endelse
  ind = where( (zmsl gt uzbottom) and (zmsl le uztop[ilayer]) )
  case ( uintmet[ilayer] ) of
    'lin': uprof[ind] = ubottom + ( ulayer[ilayer] - ubottom ) /     $
           ( uztop[ilayer] - uzbottom ) * ( zmsl[ind] - uzbottom )
    'sin': begin
           if ( get_sign( [ulayer[ilayer]-ubottom] ) gt 0 ) then                     $
             uprof[ind] = ubottom + ( ulayer[ilayer] - ubottom ) *                   $
             sin( ( zmsl[ind] - uzbottom ) / ( uztop[ilayer] - uzbottom ) * !pi/2. ) $
           else                                                                      $
             uprof[ind] = ubottom + ( ulayer[ilayer] - ubottom ) *                   $
             (sin( ( zmsl[ind] - uzbottom ) / ( uztop[ilayer] - uzbottom ) * !pi/2. + 3./2.*!pi ) + 1.)
          end
    else: message, 'Unknown method'
  endcase
endfor

; ----- v wind comonent -----
vprof = fltarr( nlevs ) * !values.f_nan
vprof[0] = vsfc
for ilayer = 0, n_elements( vztop ) - 1 do begin
  if ( ilayer eq 0 ) then begin
    vzbottom = 0. 
    vbottom = vsfc
  endif else begin
    vzbottom = vztop[ilayer-1]
    vbottom = vlayer[ilayer-1]
  endelse
  ind = where( (zmsl gt vzbottom) and (zmsl le vztop[ilayer]) )
  case ( vintmet[ilayer] ) of
    'lin': vprof[ind] = vbottom + ( vlayer[ilayer] - vbottom ) /     $
           ( vztop[ilayer] - vzbottom ) * ( zmsl[ind] - vzbottom )
    'sin': begin
           if ( get_sign( [vlayer[ilayer]-vbottom] ) gt 0 ) then                     $
             vprof[ind] = vbottom + ( vlayer[ilayer] - vbottom ) *                   $
             sin( ( zmsl[ind] - vzbottom ) / ( vztop[ilayer] - vzbottom ) * !pi/2. ) $
           else                                                                      $
             vprof[ind] = vbottom + ( vlayer[ilayer] - vbottom ) *                   $
             (sin( ( zmsl[ind] - vzbottom ) / ( vztop[ilayer] - vzbottom ) * !pi/2. + 3./2.*!pi ) + 1.)
          end
    else: message, 'Unknown method'
  endcase
endfor



; ==================================================================================================
; WRITE SOUNDING FILE
; ==================================================================================================

openw, 1, outfile

if ( add_patch_info eq 1 ) then begin
  printf, 1, npatches
  for ipatch = 0, npatches - 1 do begin
    printf, 1, xpatch_beg[ipatch], xpatch_end[ipatch], ypatch_beg[ipatch], ypatch_end[ipatch]
  endfor
endif

printf, 1, psfc, tsfc, msfc
for iz = 0, nlevs - 1 do begin
  printf, 1, zmsl[iz], tprof[iz], mprof[iz], uprof[iz], vprof[iz]
endfor

close, 1



; ==================================================================================================
; PLOT SOUNDING 
; ==================================================================================================

pt = plot( tprof, zmsl,                          $
           xtitle = 'Potential temperature (K)', $
           ytitle = 'Height (m)',                $
           xrange = [min(tprof), max(tprof)],    $
           yrange = [min(zmsl), max(zmsl)],      $
           layout = [2,2,1] )

pm = plot( mprof, zmsl,                          $
           xtitle = 'Mixing ratio (g kg$^{-1}$)',$
           ytitle = 'Height (m)',                $
           xrange = [min(mprof), max(mprof)],    $
           yrange = [min(zmsl), max(zmsl)],      $
           layout = [2,2,2], /current )

pu = plot( uprof, zmsl,                          $
           xtitle = 'u (m s$^{-1}$)',            $
           ytitle = 'Height (m)',                $
           xrange = [min(uprof), max(uprof)],    $
           yrange = [min(zmsl), max(zmsl)],      $
           layout = [2,2,3], /current )

pv = plot( vprof, zmsl,                          $
           xtitle = 'v (m s$^{-1}$)',            $
           ytitle = 'Height (m)',                $
           xrange = [min(vprof), max(vprof)],    $
           yrange = [min(zmsl), max(zmsl)],      $
           layout = [2,2,4], /current )

END