begin
; provide path & filename
  f1 = addfile("../3Dyn_vitt/dist_solver/build/testout.nc","r")
  
  itime2get = 0    ; time index starts with 0
  hgt2get   = 120. ; in km
  eps       = 0.01 ; km allowed difference for finding right height
 
  hgt= f1->alt_s1
  ;hgt= hgt*1.e-3  ;  [km]
  
  ihgt2get = closest_val(hgt2get,hgt)  ; get fixed height closest to desired height
  print("hgt= " + hgt(ihgt2get))
  hfix2get = hgt(ihgt2get)
    
  sH_1	     = f1->sigma_hal_s1(itime2get,:,:)   
  sP_1       = f1->sigma_ped_s1(itime2get,:,:)   
  
  qdlat_fld1 = f1->lat_s1      
  hgt_fld1   = f1->alt_s1  ; already in [km]
  mlon1      = f1->lon_s1
  time       = f1->time(itime2get)
  nmlon1     = dimsizes(mlon1)
  
  ihgt_tmp = new((/dimsizes(hgt_fld1)/),integer)
  icount = 0
  do i = 0,dimsizes(hgt_fld1)-1
    if(abs(hgt_fld1(i)-hfix2get).lt.eps) then
      ihgt_tmp(icount) = i
      icount = icount + 1  
    end if
  end do
  nind = icount 
  ihgt_ind = ihgt_tmp(0:nind-1)  ; index of the points which are closest to hfix2get
 
  qdlat_wanted    = qdlat_fld1(ihgt_ind)
  ip = dim_pqsort(qdlat_wanted, 1)  ; do not sort
  qdlat_wanted    =qdlat_fld1(ihgt_ind(ip))
  
  ihgt_sort = ihgt_ind(ip)
  sH1_wanted      = sH_1(ihgt_ind(ip),:)
  sP1_wanted      = sP_1(ihgt_ind(ip),:)
  
  
  ;print(qdlat_fld1(ihgt_ind)+" " + hgt_fld1(ihgt_ind))
  
; test if 0 is double 
  nlat_wanted = dimsizes(qdlat_wanted) 
  nd = ind(qdlat_wanted.eq.0)
  nzero = 0
  if(dimsizes(nd).gt.1) then
    nzero = 1
  end if 
  
  nnew = nlat_wanted-nzero
  nnew_h = nnew/2
   
; assign coordinates
  if(nzero.eq.0) then
    sH1_wanted&pflpts1 = qdlat_wanted 
    sH1_wanted&lon_s1        = mlon1
    sP1_wanted&pflpts1 = qdlat_wanted 
    sP1_wanted&lon_s1        = mlon1
  end if
  
  if(nzero.eq.1) then
    qdlat_sorted          = qdlat_wanted(0:nnew-1)
    qdlat_sorted(nnew_h:) = qdlat_wanted(nnew_h+1:nlat_wanted-1)
    delete(qdlat_wanted)
    qdlat_wanted = qdlat_sorted
    
    H_sorted                 = sH1_wanted(0:nnew-1,:)
    H_sorted(nnew_h:,:)      = sH1_wanted(nnew_h+1:nlat_wanted-1,:)
    delete(sH1_wanted)
    sH1_wanted		     = H_sorted
    sH1_wanted &fldpts_qdlat1 = qdlat_wanted 
    sH1_wanted &mlon1	     = mlon1
    delete(H_sorted)
    
    H_sorted                 = sP1_wanted(0:nnew-1,:)
    H_sorted(nnew_h:,:)      = sP1_wanted(nnew_h+1:nlat_wanted-1,:)
    delete(sP1_wanted)
    sP1_wanted		     = H_sorted
    sP1_wanted &fldpts_qdlat1 = qdlat_wanted 
    sP1_wanted &mlon1	     = mlon1
    delete(H_sorted)
    
  end if

 ; print(qdlat_wanted)
;------------------------------------------------------------------------  
  wks = gsn_open_wks("png","ceg_fld")
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
  plot = new(10,graphic)                          ; create a plot array
  
  res                       = True            ; plot mods desired
  res@gsnDraw               = False           ; don't draw
  res@gsnFrame              = False           ; don't advance frame
  res@gsnMaximize           = True            ; largest plot possible

;---This resource not needed in V6.1.0
  res@gsnSpreadColors       = True            ; use full range of colormap

;---This resource defaults to True in NCL V6.1.0
  res@lbLabelAutoStride     = True            ; nice stride on label bar
  res@vpWidthF      = 0.8
  res@vpHeightF     = 0.4

  res@cnFillOn              = True            ; turn on color
  res@cnLinesOn             = False           ; no contour lines
  res@cnFillMode            = "RasterFill"       ; Raster Mode
  res@cnLineLabelsOn        = False           ; no line labels
  res@gsnYAxisIrregular2Linear = True
  res@tiXAxisString   = "qd. longitude"
  res@tiYAxisString   = "qd. latitude INDEX"
  res@trYTensionF    = 1.
  res@tiMainString          = "h= "+decimalPlaces(hfix2get ,1,True) + " km" +"; time= "+\
    decimalPlaces(time ,3,True)+" days"
  res@gsnRightString        = "[S/m]"

;  H_wanted(:,dimsizes(mlon1)-1) = (/H_wanted(:,0)/) 
  plot(0) = gsn_csm_contour(wks,sH1_wanted,res)
  plot(1) = gsn_csm_contour(wks,sP1_wanted,res)

  resP                = True 
  resP@gsnMaximize                = True                ; maximize plots
  resP@gsnPanelYWhiteSpacePercent = 5
  resP@gsnPanelXWhiteSpacePercent = 5
  gsn_panel(wks,plot(0:1),(/2,1/),resP)             ; now draw as one plot 
               
  
end