#include "options.inc"


       module npzd_averages_module
c
c=======================================================================
c      perform averages of NPZD model
c      and write to NetCDF 
c      SPFLAME version:   c.eden
c=======================================================================
c
c      local cpp directives:
c
       implicit none
       real, allocatable, dimension(:,:,:) :: dnutmean
       real, allocatable, dimension(:,:,:) :: dco2mean
       real, allocatable, dimension(:,:) :: pco2_mean,o2_mean
       real :: airco2_mean
       integer :: counter 
       end module npzd_averages_module


       subroutine init_npzd_averages
       use spflame_module
       use npzd_averages_module
       implicit none
       if (my_pe==0) print*,' initialising NPZD diagnostic averages'
       allocate( dnutmean(is_pe:ie_pe,km,js_pe:je_pe) )
       allocate( dco2mean(is_pe:ie_pe,km,js_pe:je_pe) )
       allocate( pco2_mean(is_pe:ie_pe,js_pe:je_pe) ); pco2_mean=0.
       allocate( o2_mean(is_pe:ie_pe,js_pe:je_pe) ); o2_mean=0.
       airco2_mean=0.0; dnutmean=0.; dco2mean=0.; pco2_mean=0.
       o2_mean=0.; counter=0
       call read_unfinished_npzd_averages()
       if (my_pe==0) 
     &         print*,' done initialising NPZD diagnostic averages'
       end subroutine init_npzd_averages


       subroutine npzd_averages
       use spflame_module
       use npzd_averages_module
       use passive_tracer_module
       use npzd_module
       implicit none
       integer :: i,k,j

       if (.not. eulerback1 ) then
         counter=counter+1
         dnutmean=dnutmean+
     &         dnut_diag(is_pe:ie_pe,1:km,js_pe:je_pe)*redfield_CN
         if (enable_co2) then
          do j=js_pe,je_pe
           do k=1,km
            do i=is_pe,ie_pe
            dco2mean(i,k,j)=dco2mean(i,k,j)+
     &            (t(i,k,j,n_co2,taup1)-t(i,k,j,n_co2,taum1) )/c2dt
            enddo
           enddo
          enddo
         endif
         pco2_mean(is_pe:ie_pe,js_pe:je_pe)= 
     &       pco2_mean(is_pe:ie_pe,js_pe:je_pe)+
     &       pco2_water(is_pe:ie_pe,js_pe:je_pe)
         airco2_mean = airco2_mean + pco2_air
         o2_mean(is_pe:ie_pe,js_pe:je_pe)= 
     &       o2_mean(is_pe:ie_pe,js_pe:je_pe)
     &      +o2_sat(is_pe:ie_pe,js_pe:je_pe)
       endif ! eulerback1
       if (end_of_month) then
         call write_npzd_averages_cdf()
         counter=0 ; dnutmean=0.; dco2mean=0.
         pco2_mean=0.0; airco2_mean = 0.0; o2_mean=0.0
       endif
       if (last_time_step) call write_unfinished_npzd_averages()
       end subroutine npzd_averages


       subroutine write_unfinished_npzd_averages()
       use spflame_module
       use npzd_averages_module
       implicit none
       character (len=80) :: name,ident
       integer i,j,k,io,n,is,ie,js,je
       real, allocatable ::  buf(:,:)

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe
       write(name,'("npzd_averages_",i3,".dta")') sub_domain
       do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo

       if (my_pe==0) then
        print*,''
        print*,' writing  NPZD unfinished averages to file ',
     &          name(1:len_trim(name))
        call get_free_iounit(io)
        open(io,file=name,form='unformatted',status='unknown')
        write(io) counter
       endif

       allocate( buf(imt,jmt) )

       ident = 'dnut'
       if (my_pe==0) write(io) ident
       do k=1,km
        buf(is:ie,js:je)=dnutmean(is:ie,k,js:je)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo

       ident = 'dco2'
       if (my_pe==0) write(io) ident
       do k=1,km
        buf(is:ie,js:je)=dco2mean(is:ie,k,js:je)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo

       ident = 'pco2_water'
       if (my_pe==0) write(io) ident
       buf(is:ie,js:je)=pco2_mean(is:ie,js:je)
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
       if (my_pe==0) write(io) airco2_mean

       ident = 'o2_sat'
       if (my_pe==0) write(io) ident
       buf(is:ie,js:je)=o2_mean(is:ie,js:je)
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf

       if (my_pe==0) close(io)
       deallocate( buf )
       end subroutine write_unfinished_npzd_averages


       subroutine read_unfinished_npzd_averages()
       use spflame_module
       use npzd_averages_module
       implicit none
       character (len=80) :: name,ident
       integer i,j,k,io,n,is,ie,js,je
       real, allocatable ::  buf(:,:)
       logical :: ok

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe
       write(name,'("npzd_averages_",i3,".dta")') sub_domain
       do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo

       if (my_pe==0) then
        print*,''
        print*,' reading NPZD averages from file ',
     &          name(1:len_trim(name))
       endif

       call get_free_iounit(io)
       open(io,file=name,form='unformatted',status='old',err=200)
       read(io) counter
       allocate( buf(imt,jmt) )

 10    ident=''; ok = .true.
       if (my_pe==0) then 
        read(io,end=20) ident
        goto 30
 20     ok = .false.
 30     continue
       endif
       call bcast_logical(ok,1,0)
       call bcast_char(ident,80,0)

       if (ident(1:4) == 'dnut') then

        do k=1,km
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         dnutmean(is:ie,k,js:je)=buf(is:ie,js:je)
        enddo

       elseif (ident(1:4) == 'dco2' ) then

        do k=1,km
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         dco2mean(is:ie,k,js:je)=buf(is:ie,js:je)
        enddo

       elseif (ident(1:10) == 'pco2_water' ) then

        if (my_pe==0) read(io) buf; call pe0_send_2D(buf)
        pco2_mean(is:ie,js:je)=buf(is:ie,js:je)
        if (my_pe==0) read(io) airco2_mean
        call bcast_real(airco2_mean,1,0)

       elseif (ident(1:6) == 'o2_sat' ) then
        if (my_pe==0) read(io) buf; call pe0_send_2D(buf)
        o2_mean(is:ie,js:je)=buf(is:ie,js:je)

       else
         if (my_pe==0) print*,' cannot read ',ident,' in ',name
       endif
       if (ok) goto 10  ! read next item til end of file

       close(io)
       if (my_pe==0) then
        print*,' done '
        print*,''
       endif
       deallocate( buf )
       return

 200   if (my_pe==0) then
        print*,''
        print*,'---------------------------------------------'
        print*,' WARNING : cannot read file ',name(1:len_trim(name))
        print*,'---------------------------------------------'
        print*,''
       endif

       end subroutine read_unfinished_npzd_averages


       subroutine write_npzd_averages_cdf()
       use spflame_module
       use npzd_averages_module
       implicit none
#ifdef netcdf_diagnostics
c
c      write averaged quantities to a nice NetCDF file
c
       character (len=80) :: fname,name,unit
       integer :: year,month,day
#include "netcdf.inc"
       integer :: ncid,iret,i,j,k,dims(4),start(4),count(4)
       real    :: spval=-9.9e12, tt
       integer :: lon_tdim,lon_tid,lon_udim,lon_uid
       integer :: lat_tdim,lat_tid,lat_udim,lat_uid
       integer :: depth_tdim,depth_tid,depth_wdim,depth_wid
       integer :: itimedim,itimeid,dnutid,npe,n,dco2id
       integer :: pco2_id,o2_id,airco2_id
       integer :: is,ie,js,je
#ifdef netcdf_real4
      real (kind=4) :: var(is_pe:ie_pe,js_pe:je_pe)
      real (kind=4) , allocatable :: v2(:)
#else
      real          :: var(is_pe:ie_pe,js_pe:je_pe)
      real          , allocatable :: v2(:)
#endif

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(fname,'("npzd_averages_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            sub_domain,year,month,day
       do i=1,len_trim(fname);if (fname(i:i)==' ')fname(i:i)='0';enddo

       if (my_pe == 0) then
        print*,' Writing NPZD averages to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
c       dimensions
        lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
        Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
        Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
        Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
        depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
        depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
        iTimedim  = ncddef(ncid, 'Time', 1, iret)
c       grid variables
        dims(1)  = Lon_tdim
        Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
        dims(1)  = Lon_udim
        Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
        dims(1)  = Lat_tdim
        Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
        dims(1)  = Lat_udim
        Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
        dims(1)  = iTimedim
        iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
        dims(1)  = depth_wdim
        depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
        dims(1)  = depth_tdim
        depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c       attributes of the grid
        name = 'Longitude on T grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Longitude on U grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on T grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on U grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Depth of T Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
        call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
        name = 'Depth of W Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
        call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
        name = 'Time                    '; unit = 'days            '
        call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
        call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c       variables
        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        dnutid=ncvdef (ncid,'dnut', NCFLOAT,4,dims,iret)
        dco2id=ncvdef (ncid,'dco2', NCFLOAT,4,dims,iret)

        dims=(/lon_tdim,lat_tdim,itimedim,1/)
        pco2_id=ncvdef (ncid,'pco2_water', NCFLOAT,3,dims,iret)
        dims=(/itimedim,1,1,1/)
        airco2_id=ncvdef (ncid,'pco2_air', NCFLOAT,1,dims,iret)
        dims=(/lon_tdim,lat_tdim,itimedim,1/)
        o2_id=ncvdef (ncid,'o2_sat', NCFLOAT,3,dims,iret)
c       attributes of variables
        name = 'Nutrient consumption of DIC'; unit = 'mmol/m^3/s'
        call dvcdf(ncid,dnutid,name,len_trim(name),
     &                 unit,len_trim(unit),spval)
        name = 'Total change of DIC'; unit = 'mmol/m^3/s'
        call dvcdf(ncid,dco2id,name,len_trim(name),
     &                 unit,len_trim(unit),spval)
        name = 'Sea surface CO2 partial pressure'; unit = 'ppmv'
        call dvcdf(ncid,pco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        name = 'Atmospheric CO2 partial pressure'; unit = 'ppmv'
        call dvcdf(ncid,airco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        name = 'Sea surface O2 saturation concentration'; 
        unit = 'mol/m^3'
        call dvcdf(ncid,o2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)

c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)
c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif ! my_pe == 0

c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'dnut',dnutid)
         iret=nf_inq_varid(ncid,'dco2',dco2id)
         iret=nf_inq_varid(ncid,'pco2_water',pco2_id)
         iret=nf_inq_varid(ncid,'pco2_air',airco2_id)
         iret=nf_inq_varid(ncid,'o2_sat',o2_id)

c        divide with number of averaged time steps 
         dnutmean=dnutmean/counter
         dco2mean=dco2mean/counter
c        convert to SI units
c        apply land mask
         where( tmask(is:ie,1:km,js:je) == 0. ) dnutmean(:,:,:) = spval
         where( tmask(is:ie,1:km,js:je) == 0. ) dco2mean(:,:,:) = spval
         pco2_mean=pco2_mean/counter
         airco2_mean=airco2_mean/counter
         where(tmask(is:ie,1,js:je)==0.) pco2_mean(is:ie,js:je) = spval
         o2_mean=o2_mean/counter
         where(tmask(is:ie,1,js:je)==0.) o2_mean(is:ie,js:je) = spval
         ! write PE domain to file
         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var=dnutmean(:,k,:)
          iret= nf_put_vara_real (ncid,dnutid,start,count,var )
          var=dco2mean(:,k,:)
          iret= nf_put_vara_real (ncid,dco2id,start,count,var )
         enddo
         var(is:ie,js:je)=pco2_mean(is:ie,js:je)
         start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
         iret= nf_put_vara_real (ncid,pco2_id,start,count,var )
         var(is,js)=airco2_mean
         iret= nf_put_vara_real(ncid,airco2_id,1,1,var(is,js) )
         var=o2_mean(:,:)
         start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
         iret= nf_put_vara_real (ncid,o2_id,start,count,var )
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo
       if (my_pe==0) then
        print*,'done'
       endif
#endif
      end subroutine write_npzd_averages_cdf




