program make_ano

 implicit none

!Enver start
 integer  lrecIn1,lrecIn2,lrecIn3
!Enver end
 integer, parameter :: imax = 144, jmax = 73
 integer, parameter :: dmax = 365, tmax = 9862

! sday : number(in julian) of the first day in the first year
! eday : number(in julian) of the last day in the last year
 integer, parameter :: sday = 1, eday = 365

! year1 : 1st year, year2 : last year
 integer, parameter :: year1 = 1979, year2 = 2005

! leap = 1 : leap year
! leap = 0 : no leap year
 integer, parameter :: leap = 1

! Enver linux depend on machine. However a better options is to use
!       inquiere functions

! linux = 4 : linux machine (recl*4)
! linux = 1 : other machine (recl)
 integer, parameter :: linux = 1

 real, dimension(imax,jmax) :: ano_var
 real, dimension(imax,jmax,dmax+1) :: ave_var, n_ave_var
 real, dimension(imax,jmax,dmax+1) :: cli_var
 real, dimension(imax,jmax,tmax) :: var
 real :: dmiss
 integer, dimension(tmax) :: julian
 integer :: i, j, iy, id, it, dd, ddmax, nd, yy
 data dmiss /missing/

!Enver start 
 inquire ( IOLENGTH=lrecIn1 ) var(:,:,1)
 inquire ( IOLENGTH=lrecIn2 ) cli_var(:,:,1)
 inquire ( IOLENGTH=lrecIn3 ) ano_var(:,:)
!Enver end
 open (1, file='/home/enver/work/programs/MJOWG/msd/level_1/u200_n1/data/daily.period.gdat', &
!Enver access='direct', recl=imax*jmax*linux, status='old')
 access='direct', recl=lrecIn1, status='old')
 open (2, file='/home/enver/work/programs/MJOWG/msd/level_1/u200_n1/data/daily.clim.period.gdat', &
 access='direct', recl=lrecIn2, status='unknown')
 open (3, file='/home/enver/work/programs/MJOWG/msd/level_1/u200_n1/data/daily.anom.period.gdat', &
 access='direct', recl=lrecIn3, status='unknown')


 if (leap.eq.1) then

 it = 0

! first year
 yy = year1
 nd = dmax
 if (mod(yy,4).eq.0.and.mod(yy,100).ne.0) nd = dmax + 1
 if (mod(yy,400).eq.0) nd = dmax + 1

 do id = sday, nd

 it = it + 1
 read (1, rec = it ) var(:,:,it)

 if (nd.eq.dmax+1) then
 julian(it) = id 
 else
  if (id.le.59) then
   julian(it) = id 
  else
   julian(it) = id + 1
  endif
 endif
   
 enddo ! id

! years other than first and last
 do yy = year1+1, year2-1
 nd = dmax
 if (mod(yy,4).eq.0.and.mod(yy,100).ne.0) nd = dmax + 1
 if (mod(yy,400).eq.0) nd = dmax + 1

 do id = 1, nd

 it = it + 1
 read (1, rec = it ) var(:,:,it)

 if (nd.eq.dmax+1) then
 julian(it) = id 
 else
  if (id.le.59) then
   julian(it) = id 
  else
   julian(it) = id + 1
  endif
 endif
   
 enddo ! id

 enddo ! yy

! last year
 yy = year2
 nd = dmax
 if (mod(yy,4).eq.0.and.mod(yy,100).ne.0) nd = dmax + 1
 if (mod(yy,400).eq.0) nd = dmax + 1

 do id = 1, eday

 it = it + 1
 read (1, rec = it) var(:,:,it)

 if (nd.eq.dmax+1) then
 julian(it) = id 
 else
  if (id.le.59) then
   julian(it) = id 
  else
   julian(it) = id + 1
  endif
 endif
   
 enddo ! id

 else

  do it = 1, tmax
   read (1, rec = it ) var(:,:,it)
   julian(it) = mod(it+sday-1,dmax)
   if (julian(it).eq.0) julian(it) = dmax
  enddo ! it

 endif

  ave_var = 0.
  n_ave_var = 0.

 do it = 1, tmax
  dd = julian(it)

   do j = 1, jmax
   do i = 1, imax
     if (var(i,j,it).ne.dmiss) then  
       ave_var(i,j,dd) = ave_var(i,j,dd) + var(i,j,it)
       n_ave_var(i,j,dd) = n_ave_var(i,j,dd) + 1
     endif
   enddo
   enddo

 enddo ! it

 if (leap.eq.1) then
  ddmax = dmax + 1
 else
  ddmax = dmax
 endif

  do id = 1, ddmax

   do j = 1, jmax
   do i = 1, imax
     if (n_ave_var(i,j,id).ne.0) then  
       cli_var(i,j,id) = ave_var(i,j,id)/n_ave_var(i,j,id)
     else
       cli_var(i,j,id) = dmiss
     endif
   enddo
   enddo

   write (2, rec=id) cli_var(:,:,id)

 enddo ! id 

 do it = 1, tmax
  dd = julian(it)
  
  do j = 1, jmax
  do i = 1, imax
     if (var(i,j,it).ne.dmiss.and.cli_var(i,j,dd).ne.dmiss) then
      ano_var(i,j) = var(i,j,it) - cli_var(i,j,dd)
     else
      ano_var(i,j) = dmiss
     endif
  enddo 
  enddo 

   write (3, rec=it) ano_var
 enddo ! it

end program make_ano