! ==============================================================================
! {%DBL}_box
! generated: {%TIMEDATE}
!
! this module is automaticly generated by imdouble utility
! contains: some maintenance routines for budgeting configurations (frac)
! level: boxmodel
!
! {$DBL_INFO} ! this is a template file for imdouble utility
!
! [Gromov, MPIC, 2007-2008]
! ==============================================================================

! - general doubling parameters (as conditional defines) -----------------------

#include "{%CMODEL}_dbl_parameters.inc"

! ------------------------------------------------------------------------------

! {$CONF_PARAM}

MODULE messy_mecca_{%DBL}_box

  USE messy_mecca_kpp ! dp, nreact, nspec, ind_*, SPC_NAMES, EQN_TAGS
  USE caaba_io,       ONLY: open_output_file, write_output_file, close_file
  USE caaba_mem,      ONLY: C, cair, press

  USE {%CMODEL}_dbl_common_box

  IMPLICIT NONE

! netcdf handle for deltas, conc., etc. output
  INTEGER :: ncid_{%DBL}

! treshold value: below it, species might stop to sink to the others 
! (but can receive still)
  REAL(dp), PARAMETER :: THRES = 1.0E-40_dp * 2.5047E+19_dp  
!                                ?          * mean cair

! -----------------------------------------------------------------------------

! here constants and doubled species indices are to be defined
! {$TRAC_DECL} [%ind_@%] <-- boxmodel syntax  (%{%TAG}_@%) <-- isotracers syntax

! -----------------------------------------------------------------------------
  
! no. of "rejected" species (under threshold)
  INTEGER            :: {%DBL}_NREJCT

! output array: minor fractions (+total's fractions),
!               total concentrations (+regular) + NREJCT
  REAL(dp)           :: DOUT(({%NSPEC}+1)*({%NCLASS}-1)+ &
                             {%NCLASS}+1+1)

! classes fractions
  REAL(dp)           :: CF({%QSPEC},{%NCLASS})
  
! totals: concentration & budget fractions
  REAL(dp)           :: TCC({%NCLASS}), TCF({%NCLASS})
  
! -----------------------------------------------------------------------------

  PUBLIC {%DBL}_x0
  PUBLIC {%DBL}_emis
  PUBLIC {%DBL}_depos
  PUBLIC {%DBL}_pmix
  PUBLIC {%DBL}_process
  PUBLIC {%DBL}_calctotals
  PUBLIC {%DBL}_calcfractions
  PUBLIC {%DBL}_correct2reg
  PUBLIC {%DBL}_correct2dbl
! PUBLIC {%DBL}_fudge
  PUBLIC {%DBL}_resetPTs
  PUBLIC {%DBL}_init
  PUBLIC {%DBL}_result
  PUBLIC {%DBL}_finish

! ==============================================================================

CONTAINS

! ==============================================================================

  SUBROUTINE {%DBL}_x0
  
    IMPLICIT NONE

  ! tracers mixing ratios initialization (x0)

    INTEGER :: i

#ifndef INIUNIT_FRACMIN
 FATAL: initialization unit is not fracmin, please check configuration and former
#endif

! {$x0} [%#%] (%    CF({%TAG}_@,#) = $%)

#ifdef ZERO_TEST
    CF(:,1) = 1.0_dp
    DO i = 2, {%NCLASS}
      CF(:,i) = 0.0_dp
    ENDDO
#endif

  ! initializing isotopologues concentration according to "regular", then

  ! setting all classes fractions
    DO i = 1, {%NDSPEC}
      CF(i,1) = 1.0_dp - SUM(CF(i,2:{%NCLASS}))
    ENDDO

  ! initializing doubled tracers according to the fractions
    DO i = 1, {%NCLASS}
      C({%RDIND}(:,i)) = C({%RDIND}(:,0)) * CF(:,i)
    ENDDO

  ! updating totals in the system
    CALL {%DBL}_calctotals
    CALL {%DBL}_calcfractions
    
  END SUBROUTINE {%DBL}_x0



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_emis(ind_d, amount, fracs)
 
    IMPLICIT NONE
    
    INTEGER,  INTENT(IN)    :: ind_d
    REAL(dp), INTENT(IN)    :: amount
    REAL(dp), INTENT(IN)    :: fracs(:)

  ! filtering possible dummies
    IF ((ind_d .LT. 1) .OR. (ind_d .GT. {%NDSPEC})) RETURN

! uncomment to manage emission only through {%DBL}
!    C({%RDIND}(ind_d,0)) = C({%RDIND(ind_d,0)) + amount 

! emission of corresponding amount fractions into the box
    C({%RDIND}(ind_d,1:{%NCLASS})) = C({%RDIND}(ind_d,1:{%NCLASS})) + amount * fracs(:)

  END SUBROUTINE {%DBL}_emis



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_depos(ind_d, factor)

    IMPLICIT NONE

    INTEGER,  INTENT(IN)    :: ind_d
    REAL(dp), INTENT(IN)    :: factor

  ! filtering possible dummies
    IF ((ind_d .LT. 1) .OR. (ind_d .GT. {%NDSPEC})) RETURN

  ! simple deposition routine, introduces no selective deposition
    
    C({%RDIND}(ind_d,1:{%NCLASS})) = C({%RDIND}(ind_d,1:{%NCLASS})) * factor
    
  END SUBROUTINE {%DBL}_depos



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_pmix(TSL, dilF, ind_d, mix_amount, mix_fracs)
 
    IMPLICIT NONE

  ! pseudo-mixing of species ind_d with background concentration mix_amount of
  ! mix_deltas composition within TSL timestep with dilF dilution factor [1/s]
    
    REAL(dp), INTENT(IN)    :: TSL, dilF      ! timestep length, dilution factor 
    INTEGER,  INTENT(IN)    :: ind_d          ! spec. index
    REAL(dp), INTENT(IN)    :: mix_amount     ! backgr. concentration
    REAL(dp), INTENT(IN)    :: mix_fracs(:)   ! backgr. deltas
    REAL(dp)                :: corr, tot

  ! filtering possible dummies
    IF ((ind_d .LT. 1) .OR. (ind_d .GT. {%NDSPEC})) RETURN

  ! buget to correct to
    tot = SUM(C(RD{%A}IND(ind_d,1:{%NCLASS})))
    corr = tot + ( mix_amount - tot ) * TSL * dilF

  ! emission of background iso-composition
    CALL {%DBL}_emis(ind_d, mix_amount * TSL * dilF, mix_fracs)
    tot = SUM(C(RD{%A}IND(ind_d,1:{%NCLASS})))

  ! removal preserving current composition
    C(RD{%A}IND(ind_d,1:{%NCLASS})) = C(RD{%A}IND(ind_d,1:{%NCLASS})) / tot * corr

  END SUBROUTINE {%DBL}_pmix



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_process

    IMPLICIT NONE

    INTEGER  :: i, s
    REAL(dp) :: chkamnt

  ! calculating the number of specs falling below THRES
    {%DBL}_NREJCT = 0
    DO i = 1, {%NDSPEC}
      chkamnt = SUM(C({%RDIND}(i,1:{%NCLASS})))
      IF (chkamnt .LT. THRES) THEN
        {%DBL}_NREJCT = {%DBL}_NREJCT + 1
      ENDIF
    ENDDO           ! ndspec cycle

  ! every-step fractions/totals update
    CALL {%DBL}_calctotals
    CALL {%DBL}_calcfractions

  END SUBROUTINE {%DBL}_process



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_calctotals

    IMPLICIT NONE

    INTEGER  :: i

  ! here the number of total molecules is calculated from each species composition

  ! careful, {%ABBT}_T is calculated from regular!
    C(ind_{%CONF}T{%A}) = SUM( C({%RDIND}(:,0)) )

  ! classes concentrations
    DO i = 1, NDCLASS
      TCC(i) = SUM( C(RDIND(:,i)) )
    ENDDO

->>- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {>CONF:O3F}
    C(ind_O3F_N_T) = TCC(1)
    C(ind_O3F_Z_T) = TCC(2)
-<<- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {<CONF:O3F}

  END SUBROUTINE {%DBL}_calctotals



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_calcfractions
 
    IMPLICIT NONE

    INTEGER  :: i
    REAL(dp) :: tot

  ! calculating new delta-13C values
    DO i = 1, {%NDSPEC}
      tot = SUM(C({%RDIND}(i,1:{%NCLASS})))
      IF (tot .GT. 0.0_dp) THEN
        CF(i,:) = C({%RDIND}(i,1:{%NCLASS})) / tot 
      ELSE
        CF(i,:) = UNDEF
      ENDIF
    ENDDO        ! NISPEC-cycle

    ! totals   
    tot = SUM(TCC(:))
    IF (tot .GT. 0.0_dp) THEN
      TCF(:) = TCC(:) / tot
    ELSE
      TCF(:) = UNDEF
    ENDIF

  END SUBROUTINE {%DBL}_calcfractions
  


! -----------------------------------------------------------------------------

! correction of total isotopomers budget to "regular" species budget  

  SUBROUTINE {%DBL}_correct2reg

    IMPLICIT NONE

    INTEGER  :: i
    REAL(dp) :: tot
   
#ifdef CLASSES_1
  ! in case one class is defined, quitting
#ifdef DEBUG
    print *,'{%DBL}_correct2reg: no correction performed (one class)'
#endif
    return
#endif

  ! here is the ver. with corr. of ALL species to regular

    DO i = 1, {%NDSPEC}
      tot = SUM(C({%RDIND}(i,1:{%NCLASS})))
      IF (tot .LE. 0.0_dp) THEN
        C({%RDIND}(i,1:{%NCLASS})) = 0.0_dp
      ELSE
        C({%RDIND}(i,1:{%NCLASS})) = ( C({%RDIND}(i,1:{%NCLASS})) * C({%RDIND}(i,0)) ) / tot
      ENDIF
    ENDDO

  END SUBROUTINE {%DBL}_correct2reg



! -----------------------------------------------------------------------------
  
! correction of "regular" species budget to the total isotopologues budget

  SUBROUTINE {%DBL}_correct2dbl

    IMPLICIT NONE

    INTEGER  :: i

#ifdef CLASSES_1
  ! in case one class is defined, quitting
#ifdef DEBUG
    print *,'{%DBL}_correct2dbl: no correction performed (one class)'
#endif
    return
#endif

  ! here is the ver. with corr. of ALL species to regular
    DO i = 1, NDSPEC
      C({%RDIND}(i,0)) = SUM(C({%RDIND}(i,1:{%NCLASS})))
    ENDDO

!    C({%RDIND}(:,0)) = SUM(C({%RDIND}(:,1:{%NCLASS})),DIM=2)

  END SUBROUTINE {%DBL}_correct2dbl



! -----------------------------------------------------------------------------
  
  SUBROUTINE {%DBL}_resetPTs

  ! production tracers initialization (reset) routine

    IMPLICIT NONE
    
! {x$RESET_PTs}
! - currently disabled with use of DRPT{%ATOM}IND()

#ifdef USE_PT
    C(DRPT{%ATOM}IND(:)) = 0.0_dp    ! <-- boxmodel syntax
#endif

  END SUBROUTINE {%DBL}_resetPTs



! -----------------------------------------------------------------------------

! output file for doubled species info
  SUBROUTINE {%DBL}_init

    IMPLICIT NONE

! TODO: put additional tracers/variables+units after INIT_TRAC, INIT_UNIT

->>- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {>CONF:O3F}
    CALL open_output_file(ncid_{%DBL}, 'caaba_mecca_{%DBL}', &
      (/   &
! {$TAG_SPECS} [%fO3_@%]
       , &
{$ELSA}       'fON_T', 'fO3_T', &
{$ELSA}       'TON', 'TO3', 'TOR' &
       , &
{$ELSA}       'NREJCT' &
       /), (/   &
! {$TAG_SPECS} [%frac%]
       , &
{$ELSA}       'frac', 'frac', &
{$ELSA}       'mol/mol', 'mol/mol', 'mol/mol' &
       , &
{$ELSA}       'specs' &
       /), (/   &
! {$TAG_SPECS} [%@SRf_O_3(@)%]
       , &
{$ELSA}       '@SRf_N_O_N_-_O_3(TO)', '@SRf_O_3(TO)', &
{$ELSA}       '@SRTO_N_O_N_-_O_3', '@SRT_O_3', '@SRTO (regular)' &
       , &
{$ELSA}       '@SRnumber of rejected species' &
       /) )
-<<- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {<CONF:O3F}

  END SUBROUTINE {%DBL}_init



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_result(model_time)
  
    IMPLICIT NONE
    
    REAL(dp), INTENT(IN) :: model_time
    INTEGER              :: i

  ! last value is a common parameter
    DOUT(UBOUND(DOUT)) = REAL({%DBL}_NREJCT)
    
->>- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {>CONF:O3F}
! output array: minor fractions (+total's fractions),
!               total concentrations (+regular) + NREJCT
    DO i = 2, {%NCLASS}
      DOUT((i-2)*{%NSPEC}+1:(i-1)*{%NSPEC}) = CF(1:{%NSPEC},i)
    ENDDO
    
    DOUT(({%NSPEC})*({%NCLASS}-1)+1: &
         ({%NSPEC})*({%NCLASS}-1)+{%NCLASS}) = TCF(:)

    DOUT(({%NSPEC})*({%NCLASS}-1)+{%NCLASS}+1: &
         ({%NSPEC})*({%NCLASS}-1)+{%NCLASS}+{%NCLASS}) = TCC(:)
         
    DOUT(({%NSPEC})*({%NCLASS}-1)+{%NCLASS}+{%NCLASS}+1) = C(ind_{%CONF}T{%A})
    
    CALL write_output_file(ncid_{%DBL}, model_time, D{%A}OUT)
-<<- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {<CONF:O3F}

  END SUBROUTINE {%DBL}_result



! -----------------------------------------------------------------------------

  SUBROUTINE {%DBL}_finish

->>- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {>CONF:O3F}
    CALL close_file(ncid_{%DBL})
-<<- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ {<CONF:O3F}

  END SUBROUTINE {%DBL}_finish



! -----------------------------------------------------------------------------
  
END MODULE messy_mecca_{%DBL}_box

! *****************************************************************************

