!### macros #####################################################
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#include "tm5.inc"
!#################################################################
 
module emission_co2_ocean

  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform air-sea co2 emissions
  !
  !       This version built for use with Takahashi et al. (2002)
  !       climatological delta pCO2, using Wanninkhof (1992)
  !       short-term quadratic gas transfer velocity and the Kettle et
  !       al. (2004) barometric correction.
  !
  !       EMISSIONS
  !          interface:     call emission_apply_ocean
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  ! Table of contents:
  !   module emission_co2_ocean
  !     subroutine declare_emission_co2_ocean
  !     end subroutine declare_emission_co2_ocean
  !     subroutine calc_emission_co2_ocean
  !     end subroutine calc_emission_co2_ocean
  !     subroutine emission_apply_co2_ocean(region)
  !     end subroutine emission_apply_co2_ocean
  !     subroutine free_emission_co2_ocean
  !     end subroutine free_emission_co2_ocean
  !   end module emission_co2_ocean
  !
  !-----------------------------------------------------------------------

  use GO,                  only: gol, goErr, goPr, goLabel
  use Dims,                only: nregions, nlon360,nlat180
  use Meteo,               only: u10m_dat,v10m_dat,ci_dat
  use Meteo,               only: Set
  use chem_param,          only: nmembersloc 
  use global_types,        only: emis_data
  use emission_common,     only: enkf_is_forward, enkf_is_inverse, EmisInputDir
  use emission_common,     only: co2_ocean, co2_ocean_m, co2_ocean_n, flux_means_are_priors, fluxmultiplier
  use emission_data,       only: do_add_2d, flux_to_gridbox_per_month

  implicit none

  ! --- const -------------------------
  character(len=*), parameter  ::  mname = 'emission_co2_ocean'

  private

  ! public routines:
  public   :: declare_emission_co2_ocean, free_emission_co2_ocean, calc_emission_co2_ocean, emission_apply_co2_ocean



  type(emis_data),dimension(:,:),allocatable,target       :: co2_ocean_regxmem

  character(len=200)              :: deltapco2_prefix


contains

  subroutine declare_emission_co2_ocean
    use GO,            only: TrcFile, Init, Done, ReadRc
    use global_data,   only: rcfile
    use dims,          only: im,jm,nregions, iglbsfc

    implicit none

    integer              :: region, imr, jmr, n
    integer              :: status
    type(TrcFile)        :: rcF

    character(len=*), parameter     :: rname = mname//'/declare_emission_co2_ocean'

    if(enkf_is_forward()) then
       allocate(co2_ocean_regxmem(nregions,nmembersloc))
       
       do region=1,nregions

          imr=im(region)
          jmr=jm(region)
          do n=1,nmembersloc
             allocate(co2_ocean_regxmem(region,n)%surf(imr,jmr))
          end do
       enddo
    end if

    call Set(   u10m_dat(iglbsfc), status, used=.true. )
    call Set(   v10m_dat(iglbsfc), status, used=.true. )
    call Set(   ci_dat(iglbsfc), status, used=.true. )


    call Init( rcF, rcfile, status )        
    call ReadRc( rcF, 'deltapco2.prefix', deltapco2_prefix ,status)
    call Done( rcF ,status)

  end subroutine declare_emission_co2_ocean

  subroutine calc_emission_co2_ocean

    ! purpose is to fill in the co2_ocean_regxmem structure, every
    ! three hours. This flux has to be coarsened to the relevant region.
    ! The application is done at a different time frequency by
    ! apply_emission_co2_ocean.

    use Dims,          only: dx,dy,nlon360, nlat180, im,jm, sec_month, newsrun,  &
         idate, itau, staggered, nread, newmonth, newday, iglbsfc
    use datetime,      only: calc_sm, tstamp, tau2date
    use emission_data, only: chardate, increase_itau
    use emission_common,       only: ReadFromFile
    use toolbox,       only: escape_tm
    use toolbox,       only: coarsen_emission_2d
    use go_string,     only: goNum2str
    use ParTools,      only: myid, root ! for debugging
    use chem_param,    only: xmc

    implicit none

    character(len=*), parameter     :: rname = mname//'/calc_emission_co2_ocean'

    integer                         :: i,j,ii,jj,status
    integer, parameter              :: annual_field = 0, rank2 = 2, level1 = 1
    integer, dimension(6)           :: idate_ice
    real,dimension(nlon360,nlat180),save :: ww,k92, schmidtno,solubility,deltapco2
    real,pointer,dimension(:,:)     :: u10m,v10m,ci
    integer, parameter              :: add_field=0
    real,dimension(nlon360,nlat180) :: co2_ocean_mem ! flux for a given ensemble member
    integer                         :: mem

    if(mod(itau,nread) /= 0) return  ! only every nread hours

    ! get necessary fields

    if(newmonth) then

       call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', &
            'DELTAPCO2', deltapco2, status)
       call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', &
            'SOLUBILITY', solubility, status)
       call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', &
            'SCHMIDTNO', schmidtno, status)

       if (status /= 0) call escape_tm('stop in '//rname)

    endif

    u10m=> u10m_dat(iglbsfc)%data(:,:,1)
    v10m=> v10m_dat(iglbsfc)%data(:,:,1)
    ci=> ci_dat(iglbsfc)%data(:,:,1)

    k92=0.0

    ww=u10m**2+v10m**2  ! note: sqrt has been taken away since we square it in the next formula!

    ! make diffusion in m/s according to Wanninkhof92 param on 1x1 degree

    where(schmidtno(:,:).ge.0) k92=(0.31*(ww)/sqrt(schmidtno/660.0))*1.e-2/3600.   ! only over oceans


    ! make fluxes of co2 in mol/m2/s with tm5 meteo.  1e-6 to change deltapco2 from uatm to atm.
    co2_ocean=0

    where((solubility(:,:).ge.0).and.(deltapco2(:,:).gt.-1e10)) co2_ocean=k92*solubility*deltapco2*1e-6


    co2_ocean=co2_ocean*(1.0-ci) ! multiply out the sea-ice for this day

    ! accumulate weekly average
    if(allocated(co2_ocean_m)) then
       if(flux_means_are_priors .eqv. .TRUE.) then
          co2_ocean_m=co2_ocean_m+co2_ocean ! mol/m2/s
       else
          co2_ocean_m=co2_ocean_m+co2_ocean*fluxmultiplier(:,:,1) ! mol/m2/s, tracet 0 on PE0 is mean
       endif
       co2_ocean_n=co2_ocean_n+1
    endif

    nullify(u10m)
    nullify(v10m)
    nullify(ci)

    ! The rest only useful to the forward code
    if(enkf_is_inverse()) return

    ! Apply net flux scaling by ensemble member, change units, and
    ! coarsen to the regions.  When parameters are used beyond simple
    ! scaling of the net flux, some piece of logic from above will
    ! have to be subsumed by this loop to get the parameter
    ! deviations.

    do mem=1,nmembersloc
       co2_ocean_mem=co2_ocean*fluxmultiplier(:,:,mem) ! parameter applied, mol/m2/s
       ! accumulate weekly average
       co2_ocean_mem=co2_ocean_mem*xmc/1.e3 ! from mol/m2/s to kgC/m2/s
       call flux_to_gridbox_per_month(co2_ocean_mem) ! to kgC/area/month
       call coarsen_emission_2d('co2_ocean',360,180,co2_ocean_mem,co2_ocean_regxmem(:,mem),add_field)
    enddo


  end subroutine calc_emission_co2_ocean

  subroutine emission_apply_co2_ocean(region)

    ! This routine should only be called by the forward enkf code.  It
    ! applies the ocean co2 flux by member and region to the ico2ocean
    ! tracer.  In the inverse code, the calculated flux from each
    ! module is applied to the (single) co2 tracer.

    use dims,          only: okdebug
    use ParTools,      only: tracer_active, tracer_active
    use chem_param, only         : ico2_oce, xmc, xmco2


    implicit none

    integer, intent(in)        :: region
    integer                    :: mem

    character(len=*), parameter     :: rname = mname//'/emission_apply_co2_ocean'

    if(enkf_is_forward()) then
       if(tracer_active(ico2_oce)) then
             call do_add_2d(region,ico2_oce,1,co2_ocean_regxmem(region,1)%surf,xmco2,xmc)
       endif
    endif

    if(okdebug) write(*,*) 'end of emission_apply_co2_ocean'

  end subroutine emission_apply_co2_ocean

  subroutine free_emission_co2_ocean
    use dims,       only: nregions
    implicit none
    integer                                          :: region, n

    character(len=*), parameter     :: rname = mname//'/free_emission_co2_ocean'

    if(enkf_is_forward()) then
       do region=1,nregions
          do n=1,nmembersloc
             deallocate(co2_ocean_regxmem(region,n)%surf)
          end do
       enddo
    
       deallocate(co2_ocean_regxmem)
    end if
  end subroutine free_emission_co2_ocean


end module emission_co2_ocean
