Global areal average including pole

!______________________________________________________________________________________
!
! Original Programmer : Taufiq Hassan
! Creation Date : Aug. 15, 2014
! Last Update : Sep. 1, 2014 (by Taufiq)
!======================================================================================!
! This Program will read a T42 gridded binary file and extract any variable within.    !
! It will calcuLATe the global average of any variable with LATitude weighting.        !           !
! To account for missing AREA around the poles, I had to give extra weighting to       !
! the first (LAT=-87.8638) and the last (LAT=87.8638) LATitude values.             !
!======================================================================================!
program convert
implicit none
!_______________________Variable Declaration_________________________________________!
real, parameter :: UNDEF = -9999.0, pi=3.14159265
integer,parameter::NX=128,NY=64,NZ=1,NT=12,NV=4,NS=0
character*120 :: name
integer:: j,k,i
real:: OTOA(NX,NY),TOA(NX,NY),SUMLON(NY),LON_AVG(NY),LAT(NY),AREA(NY),DIFF(NX,NY),SUMM(NX,NY),TIME_AVG(NX,NY)
real:: top,btm,true_avg,L,ratio
!____________________________________________________________________________________!
!==================Opening and retreiving variable===================================!
do i=75,100
  if (i<100) then
    write(name,'(A,I2,A)') 'SSA_',i,'.data'
    else
    write(name,'(A,I3,A)') 'SSA_',i,'.data'
  end if
  open(unit = 10, File=name,status='unknown', &
        form='unformatted',access = 'direct', recl=4*NX*NY)
  open(unit = 11, File='NOaerosol_allsky.data',status='unknown', &
        form='unformatted',access = 'direct', recl=4*NX*NY)
  SUMM(:,:)=0.0
  do k=0,11
  read (10,rec=2+NV*k) OTOA !for any month
  read (11,rec=1+3*k) TOA
  ! print*, TOA(4,5)
  DIFF=TOA-OTOA
  ! print*, DIFF(8,9)
  SUMM=SUMM+DIFF
  end do
  TIME_AVG=SUMM/12.
! print*, TIME_AVG(20,5)
!____________________________________________________________________________________!
!===========The mean variable at each LATitude for the selected month================!
SUMLON=0.0
do k=1,NY
  do j=1,NX
      SUMLON(k)=SUMLON(k)+TIME_AVG(j,k)
  enddo
enddo
LON_AVG=SUMLON/128.0
!------------------------------------------------------------------------------------!
! using LATitude as an array
LAT=(/-87.8638, -85.09653, -82.31291, -79.5256, -76.7369, -73.94752, &
-71.15775, -68.36776, -65.57761, -62.78735, -59.99702, -57.20663, &
-54.4162, -51.62573, -48.83524, -46.04473, -43.2542, -40.46365, &
-37.67309, -34.88252, -32.09195, -29.30136, -26.51077, -23.72017, &
-20.92957, -18.13897, -15.34836, -12.55776, -9.767145, -6.976533, &
-4.185921, -1.395307, 1.395307, 4.185921, 6.976533, 9.767145, 12.55776, &
15.34836, 18.13897, 20.92957, 23.72017, 26.51077, 29.30136, 32.09195, &
34.88252, 37.67309, 40.46365, 43.2542, 46.04473, 48.83524, 51.62573, &
54.4162, 57.20663, 59.99702, 62.78735, 65.57761, 68.36776, 71.15775, &
73.94752, 76.7369, 79.5256, 82.31291, 85.09653, 87.8638/)
! =====================calcuLATing reLATive AREA for each grid ========================!
do j=2,63
AREA(j)=cos(LAT(j)*pi/180.)*(LAT(j+1)-LAT(j-1))/2.0
enddo
AREA(1)=cos(LAT(1)*pi/180.)*((LAT(1+1)+LAT(1))/2.0+90.0)
AREA(NY)=cos(LAT(NY)*pi/180.)*(90.0-(LAT(NY-1)+LAT(NY))/2.0)
!=========Taking cosine weighted mean ================================================!
top=0.0; btm=0.0
do j=1,64
 top=top+AREA(j)*LON_AVG(j)
 btm=btm+AREA(j)
enddo
!======================calcuLATe true AREA average===================================!
true_avg=top/btm
print*, true_avg
close(10)
close(11)
enddo
end