Index: LMDZ6/trunk/libf/phylmdiso/cltracrn.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/cltracrn.F90	(revision 3939)
+++ LMDZ6/trunk/libf/phylmdiso/cltracrn.F90	(revision 3940)
@@ -1,288 +1,1 @@
-!$Id $
-
-SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &
-     cdrag,coef,t,ftsol,pctsrf,               &
-     tr,trs,paprs,pplay,delp,                 &
-     masktr,fshtr,hsoltr,tautr,vdeptr,        &
-     lat,d_tr,d_trs )
-  
-  USE dimphy
-  USE traclmdz_mod, ONLY : id_rn, id_pb
-  USE indice_sol_mod
-
-  IMPLICIT NONE
-!======================================================================
-! Auteur(s): Alex/LMD) date:  fev 99
-!            inspire de clqh + clvent
-! Objet: diffusion verticale de traceurs avec quantite de traceur ds 
-!        le sol ( reservoir de sol de radon ) 
-!        
-! note : pour l'instant le traceur dans le sol et le flux sont
-!        calcules mais ils ne servent que de diagnostiques
-!        seule la tendance sur le traceur est sortie (d_tr)
-!---------------------------------------------------------------------
-! Arguments:
-! itr......input-R-  le type de traceur : id_rn(radon), id_pb(plomb)
-! dtime....input-R-  intervalle du temps (en secondes) ~ pdtphys
-! u1lay....input-R-  vent u de la premiere couche (m/s)
-! v1lay....input-R-  vent v de la premiere couche (m/s)
-! cdrag....input-R-  cdrag
-! coef.....input-R-  le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
-! t........input-R-  temperature (K)
-! paprs....input-R-  pression a inter-couche (Pa)
-! pplay....input-R-  pression au milieu de couche (Pa)
-! delp.....input-R-  epaisseur de couche (Pa)
-! ftsol....input-R-  temperature du sol (en Kelvin)
-! tr.......input-R-  traceurs
-! trs......input-R-  traceurs dans le sol
-! masktr...input-R-  Masque reservoir de sol traceur (1 = reservoir)
-! fshtr....input-R-  Flux surfacique de production dans le sol
-! tautr....input-R-  Constante de decroissance du traceur
-! vdeptr...input-R-  Vitesse de depot sec dans la couche brownienne
-! hsoltr...input-R-  Epaisseur equivalente du reservoir de sol
-! lat......input-R-  latitude en degree
-! d_tr.....output-R- le changement de "tr"
-! d_trs....output-R- le changement de "trs"
-!======================================================================
-  include "YOMCST.h"
-!
-!Entrees
-  INTEGER,INTENT(IN)                     :: itr
-  REAL,INTENT(IN)                        :: dtime
-  REAL,DIMENSION(klon),INTENT(IN)        :: u1lay, v1lay
-  REAL,DIMENSION(klon),INTENT(IN)        :: cdrag
-  REAL,DIMENSION(klon,klev),INTENT(IN)   :: coef, t
-  REAL,DIMENSION(klon,nbsrf),INTENT(IN)  :: ftsol, pctsrf 
-  REAL,DIMENSION(klon,klev),INTENT(IN)   :: tr 
-  REAL,DIMENSION(klon),INTENT(IN)        :: trs
-  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs 
-  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay, delp
-  REAL,DIMENSION(klon),INTENT(IN)        :: masktr 
-  REAL,DIMENSION(klon),INTENT(IN)        :: fshtr 
-  REAL,INTENT(IN)                        :: hsoltr
-  REAL,INTENT(IN)                        :: tautr
-  REAL,INTENT(IN)                        :: vdeptr
-  REAL,DIMENSION(klon),INTENT(IN)        :: lat   
-
-!Sorties
-  REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
-  REAL,DIMENSION(klon),INTENT(OUT) :: d_trs  ! (diagnostic) traceur ds le sol
-
-!Locales
-  REAL,DIMENSION(klon,klev) :: flux_tr  ! (diagnostic) flux de traceur
-  INTEGER                   :: i, k, n, l
-  REAL,DIMENSION(klon)      :: rotrhi
-  REAL,DIMENSION(klon,klev) :: zx_coef
-  REAL,DIMENSION(klon)      :: zx_buf
-  REAL,DIMENSION(klon,klev) :: zx_ctr
-  REAL,DIMENSION(klon,klev) :: zx_dtr
-  REAL,DIMENSION(klon)      :: zx_trs
-  REAL                      :: zx_a, zx_b
-  
-  REAL,DIMENSION(klon,klev) :: local_tr
-  REAL,DIMENSION(klon)      :: local_trs
-  REAL,DIMENSION(klon)      :: zts      ! champ de temperature du sol
-  REAL,DIMENSION(klon)      :: zx_alpha1, zx_alpha2
-!======================================================================
-!AA Pour l'instant les 4 types de surface ne sont pas pris en compte
-!AA On fabrique avec zts un champ de temperature de sol  
-!AA que le pondere par la fraction de nature de sol.
- 
-  DO i = 1,klon
-     zts(i) = 0. 
-  ENDDO
-
-  DO n=1,nbsrf
-     DO i = 1,klon
-        zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
-     ENDDO
-  ENDDO
-
-  DO i = 1,klon
-     rotrhi(i) = RD * zts(i) / hsoltr 
-  ENDDO
-
-  DO k = 1, klev
-     DO i = 1, klon
-        local_tr(i,k) = tr(i,k)
-     ENDDO
-  ENDDO
-
-  DO i = 1, klon
-     local_trs(i) = trs(i)
-  ENDDO
-!======================================================================
-!AA   Attention si dans clmain zx_alf1(i) = 1.0 
-!AA   Il doit y avoir coherence (dc la meme chose ici)
-
-  DO i = 1, klon
-!AA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
-     zx_alpha1(i) = 1.0
-     zx_alpha2(i) = 1.0 - zx_alpha1(i)
-  ENDDO
-!======================================================================
-  DO i = 1, klon
-     zx_coef(i,1) = cdrag(i)*(1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
-          *pplay(i,1)/(RD*t(i,1))
-     zx_coef(i,1) = zx_coef(i,1) * dtime*RG
-  ENDDO
-
-  DO k = 2, klev
-     DO i = 1, klon
-        zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
-             *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
-        zx_coef(i,k) = zx_coef(i,k) * dtime*RG
-     ENDDO
-  ENDDO
-!======================================================================
-  DO i = 1, klon
-     zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)
-     zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
-     zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
-  ENDDO
-
-  DO l = klev-1, 2 , -1
-     DO i = 1, klon
-        zx_buf(i) = delp(i,l)+zx_coef(i,l)      &
-             +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
- 
-        zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) &
-             + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
-        zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
-     ENDDO
-  ENDDO
-
-  DO i = 1, klon
-     zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2))  &
-          + masktr(i) * zx_coef(i,1)                        &
-          *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )
-
-     zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1)                &
-          + zx_ctr(i,2)                                     &
-          *(zx_coef(i,2)                                    &
-          - masktr(i) * zx_coef(i,1)                        &
-          *zx_alpha2(i) ) ) / zx_buf(i)
-     zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
-  ENDDO
-!======================================================================
-! Calculer d'abord local_trs nouvelle quantite dans le reservoir
-! de sol
-!=====================================================================
-
-  DO i = 1, klon
-!-------------------------
-! Au dessus des continents
-!--
-! Le pb peut se deposer partout : vdeptr = 10-3 m/s
-! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
-!-------------------------------------------------------------------
-     IF ( NINT(masktr(i)) .EQ. 1 ) THEN
-        zx_trs(i) = local_trs(i)
-        zx_a = zx_trs(i)                                           &
-             +fshtr(i)*dtime*rotrhi(i)                             &
-             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
-             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
-             +zx_alpha2(i)*zx_ctr(i,2))
-! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance
-        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG            &
-             * (1.-zx_dtr(i,1)                                     &
-             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))             &
-             + dtime / tautr                                       & 
-             + dtime * vdeptr / hsoltr 
-        zx_trs(i) = zx_a / zx_b
-        local_trs(i) = zx_trs(i)
-     ENDIF
-!--------------------------------------------------------
-! Si on est entre 60N et 70N on divise par 2 l'emanation
-!--------------------------------------------------------
-
-     IF ( (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR.      &
-          (itr.eq.id_pb.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN
-        zx_trs(i) = local_trs(i)
-        zx_a = zx_trs(i)                                           &
-             +(fshtr(i)/2.)*dtime*rotrhi(i)                        & 
-             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
-             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
-             +zx_alpha2(i)*zx_ctr(i,2))
-        !
-        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG  &
-             * (1.-zx_dtr(i,1)                           &
-             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))   &
-             + dtime / tautr                             &
-             + dtime * vdeptr / hsoltr
-        ! 
-        zx_trs(i) = zx_a / zx_b
-        local_trs(i) = zx_trs(i)
-     ENDIF
-
-!----------------------------------------------
-! Au dessus des oceans et aux hautes latitudes
-!--
-! au dessous de -60S  pas d'emission de radon au dessus 
-! des oceans et des continents
-!---------------------------------------------------------------
-
-     IF ( (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.0).OR.       &
-          (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
-        zx_trs(i) = 0.
-        local_trs(i) = 0.
-     END IF
-!--
-! au dessus de 70 N pas d'emission de radon au dessus 
-! des oceans et des continents
-!--------------------------------------------------------------
-     IF ( (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.0).OR.    &
-          (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
-        zx_trs(i) = 0.
-        local_trs(i) = 0.
-     END IF
-!---------------------------------------------
-! Au dessus des oceans la source est nulle
-!--------------------------------------------
-
-     IF (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.0) THEN
-        zx_trs(i) = 0.
-        local_trs(i) = 0.
-     END IF
-
-  ENDDO    ! sur le i=1,klon
-!
-!======================================================================
-! Une fois on a zx_trs, on peut faire l'iteration 
-!====================================================================== 
-
-  DO i = 1, klon
-     local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
-  ENDDO
-  DO l = 2, klev
-     DO i = 1, klon
-        local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
-     ENDDO
-  ENDDO
-!======================================================================
-! Calcul du flux de traceur (flux_tr): UA/(m**2 s)
-!======================================================================
-  DO i = 1, klon
-     flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG                      &
-          * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2) &
-          -zx_trs(i)) / dtime
-  ENDDO
-  DO l = 2, klev
-     DO i = 1, klon
-        flux_tr(i,l) = zx_coef(i,l)/RG                    &
-             * (local_tr(i,l)-local_tr(i,l-1)) / dtime
-     ENDDO
-  ENDDO
-!======================================================================
-! Calcul des tendances du traceur ds le sol et dans l'atmosphere
-!======================================================================
-  DO l = 1, klev
-     DO i = 1, klon
-        d_tr(i,l) = local_tr(i,l) - tr(i,l)
-     ENDDO
-  ENDDO
-  DO i = 1, klon
-     d_trs(i) = local_trs(i) - trs(i)
-  ENDDO
-
-END SUBROUTINE cltracrn
+link ../phylmd/cltracrn.F90
