source: trunk/MESOSCALE_DEV/PLOT/RESERVE/sea_level.idl @ 603

if ((winds(0) eq 'tk') and (winds(1) eq 'HGT')) then begin

print, 'computing sea-level pressure'
title='Reduced surf. pressure (Pa)'
;title='Equiv. Temperature (K)'
pression=what_I_plot


;what_I_plot = what_I_plot + 30.
;overvector_x = overvector_x*0. + overvector_x(0,0) 

;;
;; TGCM
;;
;nx=n_elements(overvector_y(*,0))
;ny=n_elements(overvector_y(0,*))
;tgcm=(  overvector_x(0,0)+ $
;        overvector_x(nx-1,0)+ $
;        overvector_x(0,ny-1)+ $
;        overvector_x(nx-1,ny-1) ) / 4.
;print, overvector_x(0,0), overvector_x(nx-1,0), overvector_x(0,ny-1), overvector_x(nx-1,ny-1)
;print, tgcm
;;
;; USUAL SYSTEMATIC BIAS
;;
;what_I_plot = what_I_plot + 20.
;;
;; ROUGH T EFFECT ON RADIATIVE TRANSFER
;;
;print, mean(overvector_x)
;caca=tgcm - overvector_x
;;caca=tgcm - mean(overvector_x)
;what_I_plot = what_I_plot - caca
;;
;; TEMP FOR HYPSOMETRIC
;; - attenue un peu - peu d'effets
;overvector_x = overvector_x*0. + tgcm 
;;
;; NOISE OF 2.2 Pa
;;
;what_I_plot = what_I_plot + 2.2*(2.*(RANDOMU(seed,nx,ny)-0.5))

;;; NON
;;;ecart_t = 222. - overvector_x                ;; tgcm - tmeso
;;;ecart_p = - ecart_t                  ;; pgcm - pmeso
;;;yeahyeah = what_I_plot - ecart_p     ;; pmeso + pgcm - pmeso = pgcm
;;;what_I_plot = yeahyeah

;;;;marche bof
;;overvector_y=CONGRID(overvector_y,nx*50,ny*50)
;;overvector_y=shift(overvector_y,1,1)
;;overvector_y=CONGRID(overvector_y,nx,ny)



;;augmen=10./100.
;;augmen=0.
;;mmm=mean(overvector_y)
;;overvector_y=overvector_y-mmm
;;overvector_y=overvector_y*(1.+augmen)
;;print, max(overvector_y), min(overvector_y)
;;overvector_y=mmm+overvector_y
;;print, max(overvector_y), mean(overvector_y), min(overvector_y)


;;*******************************BRUIT
;;;;50m limite ...
;noise=25.
;print, max(overvector_y), mean(overvector_y), min(overvector_y)
;overvector_y=overvector_y + noise*(2.*(RANDOMU(seed,nx,ny)-0.5))
;print, max(overvector_y), mean(overvector_y), min(overvector_y)
;;*******************************BRUIT










;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
H=192.*overvector_x/3.72        ;;unit is m
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
zref=mean(overvector_y)         ;;topo (m)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


goto, notest
;;;****************************************************************************************
        ;;; TEMPERATURE
        ;temp_meso=overvector_x    ;;save if needed
        ;overvector_x=mean(overvector_x)
        ;H=192.*overvector_x/3.72
        ;tempbias=temp_meso-overvector_x

;; PRESSURE
;ptopo=610.*exp(-overvector_y/H) & what_I_plot=ptopo    ;; perfect case

                ;w1=where((overvector_y - median(overvector_y)) gt 200.) & w2=where((overvector_y - median(overvector_y)) lt -200.)
                ;tempbias=overvector_x*0. & tempbias[w1] = -20. & tempbias[w2] = +20.

;; TEMPBIAS IN SCALE HEIGHT 
;tempbias=mean(overvector_x)-overvector_x
;tempbias=-10.

;; ZREF
;zref=0.        

;; REFPRES
;refpres=610.


;;
;; perfect case
;;
;ptopo=610.*exp(-overvector_y/H) & what_I_plot=ptopo 
;tempbias=222.-overvector_x
;zref=0.
;refpres=610.

;;
;; comp OMEGA
;;
;tcst = 222. ;;230. 226. 206.
;tempbias=tcst-overvector_x
zref=-945.43469
;refpres=0.
;zref=-200.
;zref=-400.
;zref=-95.
;;;****************************************************************************************
notest:

;goto, caca
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                       ;;
;; HYPSOMETRIC EQUATION  ;;
;;                       ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; recalculate H with temperature bias
;
if (n_elements(tempbias) eq 0) then tempbias=0. ;; verif: donne champ constant quand perfect
H=192.*(overvector_x+tempbias)/3.72
print, 'zref', zref, ' meanH', mean(H)
;
; hypsometric
;
what_I_plot=what_I_plot*exp((overvector_y-zref)/H)
;
; anomaly
;
if (n_elements(refpres) eq 0) then refpres=median(what_I_plot)
if (refpres eq -9999) then refpres=median(what_I_plot)
what_I_plot=what_I_plot-refpres
refpres=-9999.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
caca:


goto, coucou


        ;;for coco=0,100 do begin
        ;;altchosen = 100. - float(coco*10.)
        ;;print, altchosen


;altchosen=-915.93644
;altchosen=200.
;altchosen=max(overvector_y)
;altchosen=min(overvector_y)
;altchosen=mean(overvector_y)
;altchosen=0.
;altchosen=-2000.
;altchosen=-200.
;altchosen=200.
;altchosen=-400.
;altchosen=1000.
;altchosen=200.  ;; le meilleur match entre H et RT1km/g ATTENTION, NON !!
;altchosen=-200.  ;; le meilleur pour 13h
;altchosen=-400. ;;Tequiv trop hautes
;altchosen=-90.
;altchosen=-93.
;altchosen=-95.
;altchosen=-400.
;altchosen=-890.

altchosen=-50.
;read, altchosen, prompt='alt?'

;;-----------------------;;
;; find nearest point !  ;;
;;-----------------------;;
dis=overvector_y-altchosen & w=where(abs(dis) eq min(abs(dis))) & zref=overvector_y[w(0)] & Pref=pression[w(0)]
print, zref, Pref


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;
;hache=192.* overvector_x / 3.72 & what_I_plot2=hache*alog(Pref/what_I_plot) - (overvector_y-altchosen)
;what_I_plot=what_I_plot2
;print, min(what_I_plot2), max(what_I_plot2), mean(what_I_plot2)
;goto, coucou
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;-------------------------;;
;; POSSIBLE 0/0 PROBLEMS   ;;
;;-------------------------;;
w=where(abs(alog(Pref/what_I_plot)) lt 0.02) & if (w(0) ne -1) then what_I_plot[w]=!Values.F_NAN
w=where(abs(overvector_y-zref) lt 100.) & if (w(0) ne -1) then what_I_plot[w]=!Values.F_NAN 

;;-------------------------;;
;; EQUIVALENT TEMPERATURE  ;;
;;-------------------------;;
;
; H
;
what_I_plot=(overvector_y-zref)/alog(Pref/what_I_plot)/1000.
;
; T
;
what_I_plot=what_I_plot*3.72*1000./192.                 
;what_I_plot=what_I_plot - overvector_x


;;-----------------------;;
;; MOYENNE OK AVEC NANs  ;;
;;-----------------------;;
;
;w=where(FINITE(what_I_plot) ne 0) & moyen=mean(what_I_plot[w]) & print, moyen
;what_I_plot=what_I_plot-moyen

coucou:

        ;;endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




;what_I_plot=what_I_plot-mean(what_I_plot)



;;M1
;what_I_plot=what_I_plot-ptopo-mean(what_I_plot-ptopo)


;what_I_plot=caca
;what_I_plot=tempbias
;what_I_plot=topo_ecart



;level_1km=6
;;
;; mean temperature for the whole period
;;
;temptmp=reform(u[*,*,level_1km,*])
;meantemp=total(temptmp,3)/n_elements(temptmp(0,0,*))
;Hmean=192.*meantemp/3.72
;;
;; mean surface pressure for the whole period
;;
;temppres=reform(var[*,*,theloop,*])
;meanpres=total(temppres,3)/n_elements(temppres(0,0,*))
;
;
;        ;** HYDROSTATIC REDUCTION TO LEVEL OF REFERENCE
;       
;       overvector_x=reform(u[*,*,level_1km,nt])        ;;T 1km
;       H=192.*overvector_x/3.72        ;unit is m
;
;        overvector_y=reform(v[*,*,0,nt])        ;;topo (m)
;        zref=mean(overvector_y)
;               ;zref=median(overvector_y)
;               ;zref=(max(overvector_y)-min(overvector_y))/2.
;               ;zref=3000.
;        print, 'zref', zref
;
;       ;
;       ; HYPSOMETRIC EQUATION
;       ;
;;**********
;H=Hmean
;;**********
;       what_I_plot=what_I_plot*exp((overvector_y-zref)/H)
;       meanpres=meanpres*exp((overvector_y-zref)/Hmean)
;
;
;;
;; pressure anomaly
;;
;what_I_plot=what_I_plot-meanpres
;
;
;;presoro=610.*exp(overvector_y/10000.)
;;what_I_plot=what_I_plot-presoro


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




overvector_x=0
overvector_y=0

endif