!
#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
!
#define PRISM_ERR call prism_error(status,gol); call goErr
#define IF_PRISM_NOTOK_RETURN(action) if (status/=0) then; PRISM_ERR; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:       PRISM_PUTGET
!
! !DESCRIPTION: "Get" and "Put" PRISM
!\\
!\\
! !INTERFACE: 
!
MODULE PRISM_PUTGET
  !
  ! !USES:
  !
  USE MOD_OASIS
  use GO,           only : gol, goPr, goErr
  use GO,           only : TDate, TIncrDate, IncrDate, Pretty, operator(+)
  use TM5_PRISM,    only : comp_id, coupled_to_lpj, coupled_to_pis
  use TM5_Prism,    only : ifs_cpl_freq, lpj_cpl_freq, pis_cpl_freq, SetPrismTime, CplVar, ncplvar
  use dims,         only : im, jm, lm, lme, ndyn_max
  use partools,     only : isroot, par_broadcast, myid
  use tm5_distgrid, only : dgrid, Get_DistGrid, GATHER, SCATTER, Dist_Arr_Stat
  use tracer_data,  only : mass_dat
  use chem_param,   only : xmc, xmco2, ico2
  use binas,        only : xmair
  !
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  TM5_Prism_Puts, TM5_Prism_Gets
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'prism_putget'
  character(len=256)           :: error_message
  !
  ! !REVISION HISTORY: 
  !   14 Dec 2010 - Ph. Le Sager - added vertical regridding in "put".
  !    8 Oct 2013 - Ph. Le Sager - dummy CO2 exchange with LPJ-Guess
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS  

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TM5_PRISM_PUTS
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TM5_PRISM_PUTS( t, status )
    !
    ! !USES:
    !
    use meteodata,    only : m_dat
#ifdef with_ecearth_optics
    use ecearth_optics, only : optics_dat
    use meteodata,      only : gph_dat
#endif
    ! for vertical regridding
    use Grid,         only : FillLevelsParents
    use meteodata,    only : levi, sp_dat
#ifdef with_m7    
    use chem_param,   only : isoanus, isoaais, isoaacs, isoacos, isoaaii
#endif
    !
    ! !INPUT PARAMETERS:
    !
    type(TDate), intent(in) :: t
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,    intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   14 Dec 2010 - P. Le Sager - added vertical regridding of fields that
    !                  are passed back to IFS
    !   18 Sep 2013 - P. Le Sager - TM5-MP
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: ivar, ilev, prism_t, region, imr, jmr, lmr
    integer :: i1, i2, j1, j2, k
    integer :: soa_itr

    real, allocatable, target :: mmr_tm5(:,:,:) ! local field on TM5 levels
    real, allocatable, target :: mmr_ifs(:,:,:) ! local field on IFS levels
    
    real, pointer :: field3D(:,:,:)  ! field being send
!    real(ip_realwp_p), allocatable :: field3D_ip(:,:,:)  ! global data to send, with oasis precision

    logical            :: subset  ! TM5 levels = subset of IFS levels ?

    character(len=10)  :: interp
    real, dimension(:,:),   pointer :: sp
    real, dimension(:,:,:), pointer :: mass, masssoa, airm, gph

    type(TDate)     :: lag_date
    type(TIncrDate) :: deltat

    real :: fscaleCO2

    fscaleCO2 = xmair / xmco2
    
    !-----------------------------------------------------------------
    ! Check if fields are effectively sent or not, to avoid extra work
    ! !! THIS ASSUME THAT COUPLED FREQUENCY FOR TM5-LPJG IS MULTIPLE OF THE ONE FOR TM5-IFS 
    !-----------------------------------------------------------------
    ! by design, the lag defined in the namcouple is exactly one dynamic time step, and positive
    deltat = IncrDate(sec=ndyn_max)
    lag_date=t+deltat

    call SetPrismTime( prism_t, lag_date, status )
    IF_NOTOK_RETURN(status=1)

    ! avoid unneeded work (only if coupled to IFS only)
    if ((.not. coupled_to_lpj).and.(.not. coupled_to_pis)) then
       if (modulo( prism_t, ifs_cpl_freq*3600) /= 0) return
    end if
    
    ! Convert from TM5 time structure to OASIS time structure
    call SetPrismTime( prism_t, t, status )
    IF_ERROR_RETURN(status=1)
!DBG    if ( isroot ) then
!DBG       write (gol,'("  prism time        : ",i6)') prism_t; call goPr
!DBG    end if
    
    ! Region dims
    region = 1
    imr = im(region) ; jmr = jm(region) ; lmr = lm(region)
    subset = lmr /= lme         ! check that if nb of levels exchanges ifs_nlev(=lme) read in tm5_prism is same as nb of TM5 levels 

#ifdef with_ecearth_optics
    gph => gph_dat(region)%data
#endif
    
    !------------------------------------------------------------------------
    ! COUPLED VARIABLES
    !------------------------------------------------------------------------

    ! storage (could go into init)
    !------------------
    call Get_DistGrid( dgrid(region), I_STRT=i1,  I_STOP=i2, J_STRT=j1,  J_STOP=j2 )

    allocate( mmr_tm5(i1:i2, j1:j2, lmr) )               ! local tracer on source level (TM)
    if (subset) allocate( mmr_ifs(i1:i2, j1:j2, lme) )   ! local tracer on target level (IFS)

#ifdef parallel_cplng
    field3D => mmr_tm5
#else
    if (isroot) then
       allocate( field3D(imr,jmr,lme) ) ! global gathered field
    else
       allocate( field3D(1,1,1) )
    endif
#endif
        
    ! Send to oasis3
    !----------------------------------------------------
    VAR: do ivar = 1, ncplvar

       if ( CplVar(ivar)%intent /= 'out' ) cycle

       select case ( CplVar(ivar)%name )
#ifdef with_ecearth_optics
       case ( 'AOD_01', 'AOD_02', 'AOD_03', 'AOD_04', 'AOD_05', 'AOD_06', 'AOD_07', &
            'AOD_08', 'AOD_09', 'AOD_10', 'AOD_11', 'AOD_12', 'AOD_13', 'AOD_14', &
            'SSA_01', 'SSA_02', 'SSA_03', 'SSA_04', 'SSA_05', 'SSA_06', 'SSA_07', &
            'SSA_08', 'SSA_09', 'SSA_10', 'SSA_11', 'SSA_12', 'SSA_13', 'SSA_14', &
            'ASF_01', 'ASF_02', 'ASF_03', 'ASF_04', 'ASF_05', 'ASF_06', 'ASF_07', &
            'ASF_08', 'ASF_09', 'ASF_10', 'ASF_11', 'ASF_12', 'ASF_13', 'ASF_14' )
         select case ( CplVar(ivar)%name )
         case ( 'AOD_01' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,1)
         case ( 'SSA_01' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,1)
         case ( 'ASF_01' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,1)
         case ( 'AOD_02' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,2)
         case ( 'SSA_02' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,2)
         case ( 'ASF_02' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,2)
         case ( 'AOD_03' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,3)
         case ( 'SSA_03' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,3)
         case ( 'ASF_03' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,3)
         case ( 'AOD_04' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,4)
         case ( 'SSA_04' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,4)
         case ( 'ASF_04' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,4)
         case ( 'AOD_05' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,5)
         case ( 'SSA_05' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,5)
         case ( 'ASF_05' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,5)
         case ( 'AOD_06' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,6)
         case ( 'SSA_06' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,6)
         case ( 'ASF_06' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,6)
         case ( 'AOD_07' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,7)
         case ( 'SSA_07' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,7)
         case ( 'ASF_07' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,7)
         case ( 'AOD_08' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,8)
         case ( 'SSA_08' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,8)
         case ( 'ASF_08' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,8)
         case ( 'AOD_09' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,9)
         case ( 'SSA_09' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,9)
         case ( 'ASF_09' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,9)
         case ( 'AOD_10' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,10)
         case ( 'SSA_10' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,10)
         case ( 'ASF_10' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,10)
         case ( 'AOD_11' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,11)
         case ( 'SSA_11' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,11)
         case ( 'ASF_11' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,11)
         case ( 'AOD_12' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,12)
         case ( 'SSA_12' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,12)
         case ( 'ASF_12' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,12)
         case ( 'AOD_13' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,13)
         case ( 'SSA_13' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,13)
         case ( 'ASF_13' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,13)
         case ( 'AOD_14' ) ; mmr_tm5 = optics_dat(region)%Ext(:,:,:,14)
         case ( 'SSA_14' ) ; mmr_tm5 = optics_dat(region)%a  (:,:,:,14)
         case ( 'ASF_14' ) ; mmr_tm5 = optics_dat(region)%g  (:,:,:,14)
         case default
           write (gol,'("unsupported optical variable name: ",a)') trim(CplVar(ivar)%name); call goErr
           TRACEBACK; status=1; return
         end select

         if (subset) then ! vertical remapping
           ! From email exchange b/w KNMI experts:
           ! How to regrid the optical properties (extinction coefficient
           ! (Ext) [1/m], single-scattering albedo (SSA) [-], and asymmetry
           ! factor (AF) [-])?
           !            
           ! For a temporal averaging, let us take the monthly mean AOD, and
           ! calculate the monthly means of the other optical parameters SSA
           ! and ASYM in the proper way, by weighting with the amount of
           ! scattering and absorption. The recipe is:
           !
           ! monthly mean AOD = <AOD> = 1/N Sum (AOD(i))   (i=1,...,N; i = day, N = number of days in the month)
           ! monthly mean SSA = <SSA> = 1/N Sum (SSA(i) x AOD(i)) / <AOD>
           ! monthly mean ASYM = <ASYM> = 1/N Sum (ASYM(i) x SSA(i) x AOD(i)) / (<SSA> x <AOD>)
           !
           ! BUT HERE WE REGRID VERTICALLY, NOT IN TIME:
           ! The SSA and AF values in the IFS sublayers should be the same
           ! as in the corresponding merged TM5 layer. Basically they are
           ! intensive variables w/r/t pressure levels.
           !
           ! For Ext we can do as follows:
           !
           ! - convert to AOD = Ext x layer height
           ! - assume that AOD is distributed according to the mass
           !    distribution in the sublayers. This is consistent with the
           !    assumption of constant (?) aerosol mixing ratios. 
           ! - convert back to Ext - could be done in IFS?
           !
           select case ( CplVar(ivar)%name(1:3) )
           case('AOD')
             ! convert to AOD
!PLS if new_units:           case('AOD', 'SSA', 'ASF' )
!PLS if new_units:             ! all optical properties are now proportional to extinction (1/m)
             do ilev=1,lmr
               mmr_tm5(:,:,ilev) = mmr_tm5(:,:,ilev) * (gph(i1:i2, j1:j2,ilev+1) - gph(i1:i2, j1:j2,ilev))
             end do
             interp='sum'
           case('SSA','ASF')       ! remove this line if new_units 
             interp='mass-aver'    ! remove this line if new_units 
           case default
             TRACEBACK; status=1; return
           end select

           sp => sp_dat(region)%data(i1:i2,j1:j2,1)

           call FillLevelsParents(levi,'n',mmr_tm5, trim(interp), mmr_ifs, status, sp)
           IF_NOTOK_RETURN(status=1)

           ! ! convert AOD back to Ext (dh_ifs unknown)
           ! if (trim(interp) == 'sum') then
           !    mmr_ifs = mmr_ifs * dh_ifs
           ! endif

#ifdef parallel_cplng
           field3D => mmr_ifs
#else
           call GATHER( dgrid(region), mmr_ifs, field3D, 0, status )
           IF_NOTOK_RETURN(status=1)
#endif           
         else
            
!PLS If using new_units in ./ecearth_optics.F90 then this should be:
!PLS            ! convert extinctions to optical depths
!PLS            if (CplVar(ivar)%name(1:3) == 'AOD' .or. &
!PLS                 CplVar(ivar)%name(1:3) == 'SSA' .or. &
!PLS                 CplVar(ivar)%name(1:3) == 'ASF') then
!PLS and IFS code needs changes too!
            
           ! convert Ext to AOD
           if (CplVar(ivar)%name(1:3) == 'AOD') then
              do ilev=1,lmr
                 mmr_tm5(:,:,ilev) = mmr_tm5(:,:,ilev) * (gph(i1:i2, j1:j2,ilev+1) - gph(i1:i2, j1:j2,ilev))
              end do
           endif
           
#ifdef parallel_cplng
           field3D => mmr_tm5
#else
           call GATHER( dgrid(region), mmr_tm5, field3D, 0, status )
           IF_NOTOK_RETURN(status=1)
#endif
         endif

#endif
         case ('ppmCO2')           ! CO2 concentration (ppmv) to LPJ-Guess/PISCES
            
            mass => mass_dat(region)%rm(i1:i2,j1:j2,:,CplVar(ivar)%itr)
            airm => m_dat(region)%data(i1:i2,j1:j2,:)

            mmr_tm5(:,:,1) = mass(:,:,1)*fscaleCO2*1e6 / airm(:,:,1)  ! Convert from mass to ppmv
            
#ifdef parallel_cplng
            field3D => mmr_tm5
#else
            call GATHER( dgrid(region), mmr_tm5(:,:,1), field3D(:,:,1), 0, status )
            IF_NOTOK_RETURN(status=1)
#endif
            
         case default             ! TRACER MASS

            mass => mass_dat(region)%rm(i1:i2,j1:j2,:,CplVar(ivar)%itr)
            airm => m_dat(region)%data(i1:i2,j1:j2,:)
            sp   => sp_dat(region)%data(i1:i2,j1:j2,1)

            select case ( CplVar(ivar)%name )
#ifdef with_m7               
            case ('OC1','OC2','OC3','OC4','OC5')

               ! Select itracer for SOA compound to be added to OCn field to be send to IFS.  This
               ! can be done because SOA and POA are assumed to have the same mass density,
               ! hygroscopicity, and LW absorption coefficient.  The nucleation mode can be ignored,
               ! because it is neglected in the cloud activation scheme as well as in the mass-based
               ! calculation of LW absorption.
               
               select case (CplVar(ivar)%name)
               case ('OC1')
                  soa_itr=isoanus
               case ('OC2')
                  soa_itr=isoaais
               case ('OC3')
                  soa_itr=isoaacs
               case ('OC4')
                  soa_itr=isoacos
               case ('OC5')
                  soa_itr=isoaaii
               end select
               masssoa => mass_dat(region)%rm(i1:i2,j1:j2,:,soa_itr)
               mmr_tm5 = (mass + masssoa) / airm  ! Convert from mass to mass-mixing-ratio
#endif
            case default
               mmr_tm5 = mass / airm  ! Convert from mass to mass-mixing-ratio
            end select

            if (subset) then ! vertical remapping
               call FillLevelsParents(levi,'n',mmr_tm5,'mass-aver', mmr_ifs, status, sp)
               IF_NOTOK_RETURN(status=1)
#ifdef parallel_cplng
               field3D => mmr_ifs
#else
               call GATHER( dgrid(region), mmr_ifs, field3D, 0, status )
               IF_NOTOK_RETURN(status=1)
#endif
            else
#ifdef parallel_cplng
               field3D => mmr_tm5
#else
               call GATHER( dgrid(region), mmr_tm5, field3D, 0, status )
               IF_NOTOK_RETURN(status=1)
#endif
            endif

        end select

       ! ---- SEND ----
#ifndef parallel_cplng
       if ( isroot ) then
#endif
          do ilev = 1, CplVar(ivar)%nlev

             ! reverse layer order if needed, but skip surface fields
             if ( CplVar(ivar)%nlev /= 1 ) then
                k=CplVar(ivar)%nlev+1-ilev
                if (subset) k=ilev+lme-CplVar(ivar)%nlev
             else
                k = ilev        ! i.e. 1
             endif
             
             call oasis_put( CplVar(ivar)%var_id(ilev), prism_t, field3D(:,:,k), status )
             
!DBG             if(isroot .and.(ilev == 1).and.(ivar==53).and.(status/=OASIS_Ok)) then
!DBG                write (gol,*)"  prism sent        : ", prism_t, status; call goPr
!DBG             endif

             select case ( status )
             case ( OASIS_Sent,    OASIS_LocTrans,  OASIS_ToRest,    OASIS_Output, &
                    OASIS_SentOut, OASIS_ToRestOut, OASIS_WaitGroup, OASIS_Ok )
               continue
             case default
               TRACEBACK
               write (error_message,'("from OASIS_PUT : ",i6)') status
               call oasis_abort( comp_id, rname, error_message )
             end select

           end do
#ifndef parallel_cplng
        end if
#endif
     end do VAR

    ! Done
    deallocate( mmr_tm5  )
#ifndef parallel_cplng
    deallocate( field3D )
#endif
    if (subset) deallocate(mmr_ifs)
    status = 0

  END SUBROUTINE TM5_PRISM_PUTS
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:      TM5_Prism_Gets
  !
  ! !DESCRIPTION:   To receive tracer fields from other models; Meteo from IFS 
  !                 is received in tmm_mf_prism.F90.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TM5_PRISM_GETS( t, status )
    !
    ! !USES:
    !
    use grid, only : TllGridInfo, AreaOper, done
    !
    ! !INPUT PARAMETERS:
    !
    type(TDate), intent(in)    ::  t
    !                         
    ! !OUTPUT PARAMETERS:     
    !                         
    integer,     intent(out)   ::  status
    !
    ! !REMARKS: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: ivar, prism_t, region, imr, jmr, istat
    integer :: i1, i2, j1, j2
    integer :: ii,jj,kk
    type(TllGridInfo) :: lli

    real, allocatable :: mmr_tm5(:,:)    ! local field
    real, allocatable :: tempo_fld(:,:)  ! global field to receive fields
    real, allocatable :: oasis_fld(:,:)  ! global field to accumulate received field
    logical :: has_data

    ! --- begin -----------------------------------
    
    if ((.not. coupled_to_lpj) .and. (.not. coupled_to_pis)) then
       status=0
       return
    end if
    
!DBG    write (gol,'("get C fluxes from LPJ-Guess and/or PISCES")'); call goPr
    !
    ! if ( (modulo(t%hour,ifs_cpl_freq)/=0) .or. any((/t%min,t%sec,t%mili/)/=0) ) then
    !    write (gol,'("  skip; only every ",i2," hour ...")') ifs_cpl_freq; call goPr
    !    status=0; return
    ! end if

    ! Storage
    !----------------------------------------------------
    region = 1
    call Get_DistGrid( dgrid(region), I_STRT=i1,  I_STOP=i2, J_STRT=j1,  J_STOP=j2, lli=lli )
    imr = im(region) ; jmr = jm(region)
    
#ifdef parallel_cplng
    allocate( oasis_fld(i1:i2, j1:j2) )
    allocate( tempo_fld(i1:i2, j1:j2) )
#else
    allocate( mmr_tm5(i1:i2, j1:j2) )
    if (isroot) then
       allocate( oasis_fld(imr,jmr) )
       allocate( tempo_fld(imr,jmr) )
    else
       allocate( oasis_fld(1,1) )
       allocate( tempo_fld(1,1) )
    endif
#endif    
    
    ! Prism time
    !----------------------------------------------------
    call SetPrismTime( prism_t, t, status )
    IF_ERROR_RETURN(status=1)
!DBG    if ( isroot ) then
!DBG       write (gol,'("  prism time        : ",i6)') prism_t; call goPr
!DBG    end if

    ! Receive variable from LPJ-Guess and/or PISCES
    !----------------------------------------------------
    has_data = .false.
    
#ifndef parallel_cplng
    IF (isroot) THEN
#endif       
       oasis_fld = 0.0
       
       VARIA: DO ivar = 1, ncplvar

          SELECT CASE ( CplVar(ivar)%cpl_name )
             
          CASE ( 'land_c_flux_nat', 'land_c_flux_ant', 'land_c_npp', 'oce_c_flux' )
             
             TEMPO_FLD = 0.0
             CALL OASIS_GET( CplVar(ivar)%var_id(1), prism_t, tempo_fld, status )
             
             SELECT CASE ( status )
             CASE (OASIS_Recvd, OASIS_FromRest, OASIS_Input, OASIS_RecvOut, OASIS_FromRestOut )
               
                ! Convert C-fluxes into kg(CO2)/m2 
                SELECT CASE ( CplVar(ivar)%cpl_name )
                 
                CASE ( 'land_c_flux_nat' ) ! LPJG: kg(C)/m2/day
                
                  TEMPO_FLD = TEMPO_FLD * xmco2/xmc

                CASE ( 'land_c_flux_ant' ) ! LPJG: kg(C)/m2/day
                
                  TEMPO_FLD = TEMPO_FLD * xmco2/xmc
                  
                CASE ( 'land_c_npp' ) ! LPJG: kg(C)/m2/day

                  TEMPO_FLD = TEMPO_FLD * xmco2/xmc
                  
                CASE ( 'oce_c_flux' ) ! PISCES: molC/m2/s
                 
                  TEMPO_FLD = TEMPO_FLD * pis_cpl_freq * 3600 * xmco2 / 1000
                   
                END SELECT
               
                OASIS_FLD = OASIS_FLD + TEMPO_FLD
                HAS_DATA = .TRUE.
                
             CASE ( OASIS_OK )

                CONTINUE

             CASE DEFAULT
                TRACEBACK
                WRITE (error_message,'("PRISM_PUTGET: Error in OASIS_GET:",i6)') status
                CALL OASIS_ABORT( comp_id, rname, error_message )

             END SELECT
          END SELECT
       END DO VARIA
       
#ifndef parallel_cplng
    END IF
    
    call par_broadcast(has_data, status)
    IF_NOTOK_RETURN(status=1)       
#endif
  
    ! add data if any
    if ( has_data ) then

      ! DBG EXAMPLE - this will print the min/max/mean of the oasis_fld array
      !call Dist_Arr_Stat(dgrid(region), 'oasisfld', oasis_fld, 0, status)
      !IF_NOTOK_RETURN(status=1)
      
       call AreaOper( lli, oasis_fld, '*', 'm2', status) ! into kg(CO2)/box
       IF_NOTOK_RETURN(status=1)
       
#ifndef parallel_cplng

       CALL SCATTER( dgrid(region), mmr_tm5, oasis_fld, 0, status )
       IF_NOTOK_RETURN(status=1)
       
       mass_dat(region)%rm(i1:i2,j1:j2,1,ico2) = mass_dat(region)%rm(i1:i2,j1:j2,1,ico2) + mmr_tm5
#else
       ! spread out flux over 2 lowest levels to avoid negative mass
       mass_dat(region)%rm(i1:i2,j1:j2,1,ico2) = mass_dat(region)%rm(i1:i2,j1:j2,1,ico2) + oasis_fld/2.0
       mass_dat(region)%rm(i1:i2,j1:j2,2,ico2) = mass_dat(region)%rm(i1:i2,j1:j2,2,ico2) + oasis_fld/2.0
#endif

       ! safety check for negative co2 concentrations, see issue #706
       if ( minval( mass_dat(region)%rm(i1:i2,j1:j2,1:2,ico2) ) .lt. 0 ) then
          do kk=1,2
             do jj=j1,j2
                do ii=i1,i2
                   if ( ( mass_dat(region)%rm(ii,jj,kk,ico2) ) .lt. 0 ) then
                      write (gol,'("negative co2 mass after OASIS_GET prism_t= ",i20," ii= ",i6," jj= ",i6," kk= ",i6,")")') prism_t,ii,jj,kk; call goErr
                      write (gol,'("mass b4= ",e15.6," mass af= ",e15.6," oasisfld= ",e15.6,")")') (mass_dat(region)%rm(ii,jj,kk,ico2)-oasis_fld(ii,jj)/2.0),(mass_dat(region)%rm(ii,jj,kk,ico2)),oasis_fld(ii,jj)/2.0; call goErr
                   end if
                end do
             end do
          end do
          status = 1
          write (gol,'("negative co2 mass after OASIS_GET")'); call goErr
          call par_broadcast(status, istat, myid)
          IF_NOTOK_RETURN(status=1)
       endif

    endif
    
    ! done
    call done(lli, status)
    IF_NOTOK_RETURN(status=1)
    
    DEALLOCATE( oasis_fld, tempo_fld )
    if(allocated(mmr_tm5)) deallocate(mmr_tm5)
    status = 0

  END SUBROUTINE TM5_PRISM_GETS
  !EOC

END MODULE PRISM_PUTGET