program make_qflux

! Description:  Reads in a file containing climatological Q-fluxes from a
!               CSIRO AGCM control run, and generates a file containing the
!               quadratic coefficients that are required to run the model in
!               Q-flux mode.
!
! Notes:	Based on qfluxminvbragg.f, written by Ian Watterson
!
! Written by:   Steve Phipps  8 November 2002

  implicit none

  ! Define parameters
  integer, parameter  :: lon=64, lat=28, nmon=12, nn=36
  real, dimension(3), parameter :: fac=(/ 1.0/3.0, 1.0/2.0, 1.0 /)

  ! Declare variables
  integer :: kmonth, kmonth1, i3, ir, lg, ns, mg, km1
  real :: sum, dif
  real, dimension(nn, nn) :: ai
  real, dimension(lon, 2, lat, nmon) :: gm1cur
  real, dimension(lon, 2, lat) :: ochf, occur
  character(len=80) :: header

  ! Formats
  130 format('Quadratic fit to Q-fluxes for month ', i2, ', field ', i1)

  ! Read in matrix of coefficients required to generate quadratic
  ! coefficients from Q-fluxes
  open (unit=10, file="matrixinv", form="formatted", status="old", &
        action="read")
  read (10, '(4e21.14)') ai
  close (10)

  ! Read in Q-fluxes
  open (unit=11, file="qflux.new", form="unformatted", status="old", &
        action="read")
  do kmonth = 1, nmon
     read (11) header
     read (11) kmonth1
     read (11) (((gm1cur(mg, ns, lg, kmonth), mg=1, lon), lg=1, lat), ns=1, 2)
  end do
  close (11)

  ! Open file to which to write quadratic coefficients
  open (unit=12, file="qfluxqu.dat", form="unformatted", status="new", &
        action="write")

  ! Loop over months
  do kmonth = 1, nmon
     i3 = kmonth*3 - 3
     ochf = 0.0

     ! Generate quadratic coefficients
     do ir = 1, 3
        do lg = 1, lat
           do ns = 1, 2
              do mg = 1, lon
                 sum = 0.0
                 do km1 = 1, nmon
                    sum = sum + ai(i3+ir, (km1-1)*3+3) * &
                                gm1cur(mg, ns, lg, km1)
                 end do
                 occur(mg, ns, lg) = sum
                 ochf(mg, ns, lg) = ochf(mg, ns, lg) + sum * fac(ir)
              end do
           end do
        end do
        write (header, 130) kmonth, ir
        write (12) header
        write (12) kmonth, ir
        write (12) (((occur(mg, ns, lg), mg=1, lon), lg=1, lat), ns=1, 2)
     end do

     ! Check that coefficients are correct
     do lg = 1, lat
        do ns = 1, 2
           do mg = 1, lon
              dif = ochf(mg, ns, lg) - gm1cur(mg, ns, lg, kmonth)
              if (abs(dif) > 1.0e-8) write (*, *) kmonth, dif, mg, ns, lg
           end do
        end do
     end do

  end do

  ! Close output file
  close (12)

end program make_qflux