PRO solar_rad_climate_budyko

;;;;; Following Budyko et. al 1968, utiliza satelite obs to studdy the effect of effect
;;;;   of clouds and surface albedo and its effect on climate.
;;; contants

alphap=0.33             ;;;  Albedo & planetary albedo
alpha1=0.32             ;;;  for Ts > -10 deg C, 0 to 60 lat 
alpha2=0.50             ;;;  Transition between these two, 60 to 70 lat
alpha3=0.62             ;;;  for Ts < -10 dGc C,  70 to 90 lat


betta= 0.235*1000              ;;;   W/M^2/dec C

a=14.0                ;;;  W/M^2
a1=3.0                ;;;  Const
b=0.14                 ;;;  Const
b1= 0.10                ;;;  Const


Q=  1.92   ;1338                  ;;;  W/M^2
Tlat0= -9.00              ;;;  surface temp deg C
Qp= 1.92    ;1338                  ;;;  W/M^2
delQp=0.05*Qp
Tp=15.0                  ;;;  Earth earth mean temperature 
delTp=0.5                ;;;  Change in Earth mean planetary temp.
delQ=0.05*Q                 ;;;  percent
n=0.50                   ;;;  CHANGE TO F(x)   ;cloud fraction, change to a function n(time,lat)
PI=22.0/7.0
;Qp=(1.92*60*24*30)/(1000)*15.92
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;

loadct, 39
;
;
;Read ceres NetCDF file
 infile='/home/kbah/code/aos801/budyko/input/aerosol-aux.1d.05807.nc'
 outfile='/home/kbah/code/aos801/budyko/output/'
 dname='npix'
 
 ncdfid= NCDF_OPEN(infile)
 dims = NCDF_DIMID( ncdfid, dname)
 
 fileinfo= NCDF_INQUIRE(ncdfid)
 NCDF_DIMINQ, ncdfid, dims, Name1, Size1

 varid0=NCDF_VARID(ncdfid, 'latitude')
 NCDF_VARGET,ncdfid, varid0, lats0
 
 varid1=NCDF_VARID(ncdfid, 'stceresm')
 NCDF_VARGET,ncdfid, varid1, Tsurf0
 
 ;goodpixels=WHERE(Tsurf0 Gt -999)
 goodpixels=WHERE(lats0 GE 0) 
 Tsurf1=Tsurf0(goodpixels)
 Tsurf1(*)=Tsurf1(*)-273.15  ;;Convert to Celsius
 lats=lats0(goodpixels)
 
 
 pp=WHERE((lats0 GT -999)and (tsurf0 GT -999))
 lats0=lats0(pp)
 tsurf0=tsurf0(pp)
 

;;define arrays
;lats=FINDGEN(91)
Tsurf2=FLTARR(N_ELEMENTS(Tsurf1))
Slat=FLTARR(N_ELEMENTS(Tsurf1))
alpha=FLTARR(N_ELEMENTS(Tsurf1))



Slat(*)= PI*COS(lats(*)/180.*(7.0))  ;Slat(*)= PI*(COS(lats(*))/2.0) ;;Lat. dependent function.
;Slat(*)=Slat(*)/max(slat)

WHILE (ARRAY_EQUAL(TSURF1, TSURF2) eq 0)DO BEGIN 
       
      FOR ii=0, N_ELEMENTS(lats)-1 DO BEGIN         
           IF (Tsurf1(ii) GT -10.0) THEN alpha(ii)=alpha1
           IF (Tsurf1(ii) EQ -10.0) THEN alpha(ii)=alpha2
           IF (Tsurf1(ii) LT -10.0) THEN alpha(ii)=alpha3        
      ENDFOR
      
      alphap=MEAN(alpha)       ;;Planetary albedo
      Tp=mean(Tsurf1)     
      ;Tp1=((Qp*(1-alphap)-(a0-n*dela))/(b0- n*delb))  ;-273.15   ; 15
      
      FOR ii=0, N_ELEMENTS(lats)-1 DO BEGIN          
           ;Tsurf2(ii)=(Qp*Slat(ii)*(1-alpha(ii))*(1+(delQ/Q)) - (a0- n*(dela)) + betta*(Tp+delTp))/(betta +b0 -n*delb)  
            Tsurf2(ii)=(Q*(1-alpha(ii))*(1+(delQp/Qp)) -  a+a1*n + betta*Tp  + (betta*Qp/b-b1*n)*((delQp/Qp)*(1-alphap- 0.30*Slat(ii)) -0.30*Slat(ii)))/(betta +b -b1*n)   
                                
      ENDFOR
       
      Tsurf2(*)=Tsurf2(*)*(-1.0) +32
      ;Tsurf1=Tsurf2
      print, 'stop'
      
      
      ;;;pLOT sub routines  ;;;;
      outfile = outfile+'test.ps'
      ttl = '!5testing tSurf plots'
      
      
      PLOT,lats[9000:18152], Tsurf1[9000:18152],THICK=2,$
        xs=3,xrange=[0,90],ys=1,yrange=[-10,40],LINESTYLE=0,$
        xtitle = '!5Latitude',ytitle = '!5Surface temperature [C]',$
        title = '!5Avg. Lat. -Temp. Distribution ', POSITION=[0.15,0.15,0.85,0.85],CHARSIZE=2.2
        oplot, lats, tsurf2,LINESTYLE=2,THICK=2 ;,COLOR=100

  
     
     outputdata=outfile + 'surfacetemp.dat'
     close,20
     openw, 20,outputdata
     printf,20, 'Tsurf1'
     printf,20, Tsurf1
     printf,20, '       '
     printf,20, '       '
     printf,20, 'Latitudes'
     printf,20, lats
     close, 20
      
      
     PRINT, 'STOP'
      
   ENDWHILE
   
print, 'stop'

END




; Tsurf(ii)=(Q*Slat(ii)*(1-alpha) - (a0 - n*dela) + betta*Tp)/(betta + (b0- n*delb))
; delTp= (Qp/b0-n*delb)*(delQp/Qp*(1-alpha-delalpha) -delalpha)   ;;changes in mean planetry temp.
; Tp=  ((Qp*(1-alpha) - (a0-n*dela))/(b0-n*delb))                 ;;Planetary Temperature
; T= (Q*(1-alpha)*(1+(delQp/Qp)) -  a+a1*n(ii) + betta*Tpm  + (betta*Qp/b-b1*n(ii))*((delQp/Qp)*(1-alphap- 0.30*S) -0.30*S))/(betta +b -b1*n(ii))
; Tlat0=Tlat(ii)
; alpha0=alpha      ;;after first loop alpha0=alpha