source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/traclmdz_mod.F90 @ 5425

!$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.  
!

  IMPLICIT NONE

  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)

!IM ajout traceurs RR 
  INTEGER,SAVE :: id_dry !traceur dry intrusions
!$OMP THREADPRIVATE(id_dry)
  INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
!$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
  INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0 ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
!                                            ! qui ne sont pas transportes par la convection
!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)

  INTEGER, SAVE:: id_o3
  !$OMP THREADPRIVATE(id_o3)
  ! index of ozone tracer with Cariolle parameterization
  ! 0 means no ozone tracer

  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
    
    ! 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, t_seri, pplay, sh, aerosol, lessivage)
    ! This subroutine allocates and initialize module variables and control variables.
    USE dimphy
    USE infotrac
    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
    USE press_coefoz_m, ONLY: press_coefoz
    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl

    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)
!IM traceurs RR   REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri! Concentration Traceur [U/KgA]  
    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA]  
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique

! Output variables
    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
    LOGICAL,INTENT(OUT)                  :: lessivage
        
! Local variables    
    INTEGER :: ierr, it, iiq, i, k
    REAL,DIMENSION(klon,klev)      :: qsat   ! pression de la vapeur a saturation
    
! --------------------------------------------
! 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, O3, CO2,...
! --------------------------------------------
    id_be=0
    id_o3=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'
       ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
          ! Recherche de l'ozone : parametrization de la chimie par Cariolle
          id_o3=it
          CALL alloc_coefoz   ! allocate ozone coefficients
          CALL press_coefoz   ! read input pressure levels
       END IF    
    END DO

    id_dry=0

    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "dry" .OR. tname(iiq) == "Dry" ) THEN
        id_dry=it
       END IF    
    END DO  

    id_pcsat=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
        id_pcsat=it
      END IF
    END DO

    id_pcocsat=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
        id_pcocsat=it
       END IF
    END DO

    id_pcq=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
        id_pcq=it
       END IF
    END DO

    id_pcs0=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
        id_pcs0=it
      END IF
    END DO

    id_pcos0=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
        id_pcos0=it
       END IF
    END DO

    id_pcq0=0
    DO it=1,nbtr
       iiq=niadv(it+2)
       IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
        id_pcq0=it
       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

!IM initialisation traceurs pseudo-vapeurs
    call q_sat(klon*klev,t_seri,pplay,qsat)
    IF ( id_pcsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcsat) = qsat(i,k)
       ELSE
        tr_seri(i,k,id_pcsat) = 100.
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcocsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        IF ( pctsrf (i, is_oce) > 0. ) THEN
         tr_seri(i,k,id_pcocsat) = qsat(i,k)
        ELSE
         tr_seri(i,k,id_pcocsat) = 100.
        END IF
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcq) = sh(i,k)
       ELSE
        tr_seri(i,k,id_pcq) = 100.
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcs0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcs0) = qsat(i,k)
       ELSE
        tr_seri(i,k,id_pcs0) = 100.
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcos0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        IF ( pctsrf (i, is_oce) > 0. ) THEN
         tr_seri(i,k,id_pcos0) = qsat(i,k)
        ELSE
         tr_seri(i,k,id_pcos0) = 100.
        END IF
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcq0) = sh(i,k)
       ELSE
        tr_seri(i,k,id_pcq0) = 100.
       END IF
      END DO
     END DO
    END IF
 
  END SUBROUTINE traclmdz_init

  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
       tr_seri, source, solsym, d_tr_cl, zmasse)
    
    USE dimphy
    USE infotrac
    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
    USE o3_chem_m, ONLY: o3_chem
    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
    INCLUDE "YOMCST.h"
    INCLUDE "indicesol.h"

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

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

!
!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)
    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2


!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
!IM traceurs RR
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique

! 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
    INTEGER lmt_pas ! number of time steps of "physics" per day

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

!IM traceurs RR
    REAL,DIMENSION(klon,klev)      :: qsat   ! pression de la vapeur a saturation
    REAL :: amn, amx
!
!=================================================================
!  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
  
!IM ajout traceurs RR
    call q_sat(klon*klev,t_seri,pplay,qsat)
    
    IF ( id_pcsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcsat) = qsat(i,k)
       END IF
      END DO 
     END DO
    END IF

    IF ( id_pcocsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        IF ( pctsrf (i, is_oce) > 0. ) THEN
         tr_seri(i,k,id_pcocsat) = qsat(i,k)
        END IF
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcq) = sh(i,k)
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcs0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcs0) = qsat(i,k)
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcos0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        IF ( pctsrf (i, is_oce) > 0. ) THEN
         tr_seri(i,k,id_pcos0) = qsat(i,k)
        END IF
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).GE.85000.) THEN
        tr_seri(i,k,id_pcq0) = sh(i,k)
       END IF
      END DO
     END DO
    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 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, zmasse * rg, &
               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
          
!IM traceurs RR
    IF ( id_pcsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
       END IF 
      END DO
     END DO
    END IF

    IF ( id_pcocsat /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq /= 0 ) THEN
     DO k = 1, klev 
      DO i = 1, klon 
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
       END IF
      END DO  
     END DO  
    END IF  

    IF ( id_pcs0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcos0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
       END IF
      END DO
     END DO
    END IF

    IF ( id_pcq0 /= 0 ) THEN
     DO k = 1, klev
      DO i = 1, klon
       IF ( pplay(i,k).LT.85000.) THEN
        tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
       END IF
      END DO
     END DO
    END IF
!======================================================================
!   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

!======================================================================
!   Parameterization of ozone chemistry
!======================================================================

    IF (id_o3 /= 0) then
       lmt_pas = NINT(86400./pdtphys)
       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
          ! Once per day, update the coefficients for ozone chemistry:
          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
       END IF
       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
            xlon, tr_seri(:, :, id_o3))
    END IF

!======================================================================
!   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
    
    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