PRO solar_rad_climate

;;;;; 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


;l=0.001                  ;;;  fraction of eath that change from non glaciation
;q=0.002                  ;;;  Ratio of aveage radiaton in the same region
betta= 3.74              ;;;   W/M^2/dec C

a0=222.80                ;;;  W/M^2
a1=174.80                ;;;  Const
b0=2.228                 ;;;  Const
b1= 0.636                ;;;  Const
delb=b0-b1
dela=a0-a1

Q=  1338                  ;;;  W/M^2
Tlat0= -9.00              ;;;  surface temp deg C
Qp= 1338                 ;;;  W/M^2
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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;
;
;
;Read ceres NetCDF file
 infile='/home/kbah/code/aos801/budyko/input/aerosol-aux.1d.05806.nc'
 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))

Tsurf3=FLTARR(N_ELEMENTS(Tsurf1))
Tsurf4=FLTARR(N_ELEMENTS(Tsurf1))
Tsurf5=FLTARR(N_ELEMENTS(Tsurf1))

Slat(*)= PI*COS(lats(*)/180.0)  ;Slat(*)= PI*(COS(lats(*))/2.0) ;;Lat. dependent function.
Qp=(1.92*60*24*30)/(1000)*15.92
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)=(Qp*Slat(ii)*(1-alpha(ii))*(1+(delQ/Q))-  (a0- n*(dela)))/(betta +b0 -n*delb)   
                                
      ENDFOR
       
      Tsurf2(*)=Tsurf2(*)-120
      ;Tsurf1=Tsurf2
      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