!$Id $
!
MODULE traclmdz_mod
! 
! In this module all tracers specific to LMDZ are treated. This module is used 
! only if running without any other chemestry model as INCA or REPROBUS.  
!

  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr   ! Masque reservoir de sol traceur
!$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
!$OMP THREADPRIVATE(fshtr)
  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: hsoltr   ! Epaisseur equivalente du reservoir de sol
!$OMP THREADPRIVATE(hsoltr)
!
!Radioelements:
!--------------
!
  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: tautr    ! Constante de decroissance radioactive
!$OMP THREADPRIVATE(tautr)
  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
!$OMP THREADPRIVATE(vdeptr)
  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: scavtr   ! Coefficient de lessivage
!$OMP THREADPRIVATE(scavtr)
  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
!$OMP THREADPRIVATE(srcbe)

  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
!$OMP THREADPRIVATE(radio)  

  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs     ! Conc. radon ds le sol
!$OMP THREADPRIVATE(trs)

  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ] 
!$OMP THREADPRIVATE(id_be)

  LOGICAL,SAVE :: rnpb=.TRUE. ! Presence du couple Rn222, Pb210
!$OMP THREADPRIVATE(rnpb)


CONTAINS

  SUBROUTINE traclmdz_from_restart(trs_in)
    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc). 
    ! This subroutine is called from phyetat0 after the field trs_in has been read.
    
    USE dimphy
    USE infotrac
    IMPLICIT NONE
    
    ! Input argument
    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in 
    
    ! Local variables
    INTEGER :: ierr
    
    ! Allocate restart variables trs
    ALLOCATE( trs(klon,nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_from_restart', 'pb in allocation 1',1)
    
    ! Initialize trs with values read from restart file 
    trs(:,:) = trs_in(:,:)
    
  END SUBROUTINE traclmdz_from_restart


  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
    ! This subroutine allocates and initialize module variables and control variables.
    USE dimphy
    USE infotrac
    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl

    IMPLICIT NONE

    INCLUDE "indicesol.h"

! Input variables
    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri! Concentration Traceur [U/KgA]  

! Output variables
    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
    LOGICAL,INTENT(OUT)                  :: lessivage
        
! Local variables    
    INTEGER :: ierr, it, iiq
    
! --------------------------------------------
! Allocation
! --------------------------------------------

    ALLOCATE( scavtr(nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 9',1)
    scavtr(:)=1.
    
    ALLOCATE( radio(nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 11',1)
    radio(:) = .false.    ! Par defaut pas decroissance radioactive
    
    ALLOCATE( masktr(klon,nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 2',1)
    
    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 3',1)
    
    ALLOCATE( hsoltr(nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 4',1)
    
    ALLOCATE( tautr(nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 5',1)
    tautr(:)  = 0.
    
    ALLOCATE( vdeptr(nbtr), stat=ierr)
    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 6',1)
    vdeptr(:) = 0.


    lessivage  = .TRUE.
    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
    
!
! Recherche des traceurs connus : Be7, CO2,...
! --------------------------------------------
    id_be=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN  
          ! Recherche du Beryllium 7
          id_be=it
          ALLOCATE( srcbe(klon,klev) )
          radio(id_be) = .TRUE.
          aerosol(id_be) = .TRUE. ! le Be est un aerosol
          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
          WRITE(*,*) 'Initialisation srcBe: OK'
       END IF    
    END DO
!
! Valeurs specifiques pour les traceurs Rn222 et Pb210
! ----------------------------------------------
    IF (rnpb) THEN
        
       radio(1)= .TRUE.
       radio(2)= .TRUE.
       pbl_flg(1) = 0 ! au lieu de clsol=true ! CL au sol calcule
       pbl_flg(2) = 0 ! au lieu de clsol=true
       
       aerosol(2) = .TRUE. ! le Pb est un aerosol
       
       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
    END IF

!
! Initialisation de module carbon_cycle_mod
! ----------------------------------------------
    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
       CALL carbon_cycle_init(tr_seri, aerosol, radio)
    END IF

  END SUBROUTINE traclmdz_init

  SUBROUTINE traclmdz(                           &
       nstep,    pdtphys,      t_seri,           &
       paprs,    pplay,        cdragh,  coefh,   &
       yu1,      yv1,          ftsol,   pctsrf,  &
       xlat,     couchelimite,                   &
       tr_seri,  source,       solsym,  d_tr_cl)
    
    USE dimphy
    USE infotrac
    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
    
    IMPLICIT NONE
    
    INCLUDE "YOMCST.h"
    INCLUDE "indicesol.h"

!==========================================================================
!                   -- DESCRIPTION DES ARGUMENTS --
!==========================================================================

! Input arguments
!
!Configuration grille,temps:
    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)  
    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point 

!
!Physique: 
!--------
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)


!Couche limite:
!--------------
!
    REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
    LOGICAL,INTENT(IN)                   :: couchelimite

! Arguments necessaires pour les sources et puits de traceur:
    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)


! InOutput argument
    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]  

! Output argument
    CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit 
    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur

!=======================================================================================
!                        -- VARIABLES LOCALES TRACEURS --
!=======================================================================================

    INTEGER :: i, k, it

    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
    REAL,DIMENSION(klon,klev)      :: delp     ! epaisseur de couche (Pa)
    
    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
    REAL                           :: zrho      ! Masse Volumique de l'air KgA/m3

!
!
!=================================================================
!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
!=================================================================
!
    IF ( id_be /= 0 ) THEN
       DO k = 1, klev
          DO i = 1, klon
             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
          END DO
       END DO
       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
    END IF
  

    DO it=1,nbtr
       WRITE(solsym(it),'(i2)') it
    END DO
!======================================================================
!     -- Calcul de l'effet de la couche limite --
!======================================================================

    IF (couchelimite) THEN             
       source(:,:) = 0.0

       IF (id_be /=0) THEN
          DO i=1, klon
             zrho = pplay(i,1)/t_seri(i,1)/RD
             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
          END DO
       END IF

    END IF
    

    DO k = 1, klev
       DO i = 1, klon
          delp(i,k) = paprs(i,k)-paprs(i,k+1)
       END DO
    END DO
    
    DO it=1, nbtr
       IF (couchelimite .AND. pbl_flg(it) == 0 ) THEN ! couche limite avec quantite dans le sol calculee
          CALL cltracrn(it, pdtphys, yu1, yv1,     &
               cdragh, coefh,t_seri,ftsol,pctsrf,  &
               tr_seri(:,:,it),trs(:,it),          &
               paprs, pplay, delp,                 &
               masktr(:,it),fshtr(:,it),hsoltr(it),&
               tautr(it),vdeptr(it),               &
               xlat,d_tr_cl(:,:,it),d_trs)
          
          DO k = 1, klev
             DO i = 1, klon
                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
             END DO
          END DO
        
          ! Traceur dans le reservoir sol
          DO i = 1, klon
             trs(i,it) = trs(i,it) + d_trs(i)
          END DO
       END IF
    END DO
           
!======================================================================
!   Calcul de l'effet du puits radioactif
!======================================================================
    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
  
    DO it=1,nbtr
       IF(radio(it)) then     
          DO k = 1, klev
             DO i = 1, klon
                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
             END DO
          END DO
          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
       END IF
    END DO

!======================================================================
!   Calcul de cycle de carbon
!======================================================================
    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
    END IF

  END SUBROUTINE traclmdz


  SUBROUTINE traclmdz_to_restart(trs_out)
    ! This subroutine is called from phyredem.F where the module 
    ! variable trs is written to restart file (restartphy.nc)
    USE dimphy
    USE infotrac
    
    IMPLICIT NONE
    
    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
    INTEGER :: ierr

    IF ( ALLOCATED(trs) ) THEN
       trs_out(:,:) = trs(:,:)
    ELSE
       ! No previous allocate of trs. This is the case for create_etat0_limit.
       trs_out(:,:) = 0.0
    END IF
    
  END SUBROUTINE traclmdz_to_restart
  

END MODULE traclmdz_mod