Changeset 1191


Ignore:
Timestamp:
Jun 24, 2009, 11:56:13 AM (15 years ago)
Author:
jghattas
Message:

Reecriture de phytrac et les routines concernes (Anthony Jamelot)

  • les suffix change de F -> F90 (nflxtr.F90,cltracrn.F90,initrrnpb.F90,cvltr.F90,minmaxqfi.F90,cltrac.F90,phytrac.F90)

Traitement d'un nouveau traceur berelium (optionel, toujours pour des
tests)(Anthony Jamelot)

  • radiornpb.F change du nom pour radio_decay.F90 car il traite maintenant tout les traceurs radioactives
  • ajoute init_be.F90

Nouveau interface dans phytrac pour serparer les calculs et appels
specifique a INCA avec les traitements des traceurs specifiques au LMDZ
(JG)

  • ajoute tracinca_mod.F90 pour les appeles a INCA
  • ajoute traclmdz_mod.F90 pour les calculs des traceurs specifiques a LMDZ
  • enleve fichier restartrac et ajoute la variable trs dans restartphy.nc

La convergence numerique a etait rompue uniquement pour les traceurs
LMDZ RN et PB.

Location:
LMDZ4/branches/LMDZ4-dev/libf
Files:
3 added
8 edited
8 moved

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90

    r1180 r1191  
    2121  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
    2222
    23 ! Variables for INCA
     23! conv_flg(it)=0 : convection desactivated for tracer number it
    2424  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
     25! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it
    2526  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
    2627
     28  CHARACTER(len=4),SAVE :: type_trac
     29 
    2730CONTAINS
    2831
    … …  
    8083    descrq(20)='SLP'
    8184    descrq(30)='PRA'
     85   
     86
     87    IF (config_inca=='none') THEN
     88       type_trac='lmdz'
     89    ELSE
     90       type_trac='inca'
     91    END IF
    8292
    8393!-----------------------------------------------------------------------
    … …  
    8797!
    8898!-----------------------------------------------------------------------
    89     IF (config_inca == 'none') THEN
     99    IF (type_trac == 'lmdz') THEN
    90100       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
    91101       IF(ierr.EQ.0) THEN
    … …  
    109119    END IF
    110120!
    111 ! Allocate variables depending on nqtrue
     121! Allocate variables depending on nqtrue and nbtr
    112122!
    113123    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
    114 
    115     IF (config_inca /= 'none') THEN
    116        ! Varaibles only needed in case of INCA
    117        ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
    118     END IF
    119        
     124    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
     125    conv_flg(:) = 1 ! convection activated for all tracers
     126    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
     127
    120128!-----------------------------------------------------------------------
    121129! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
    … …  
    144152!    Get choice of advection schema from file tracer.def or from INCA
    145153!---------------------------------------------------------------------
    146     IF (config_inca == 'none') THEN
     154    IF (type_trac == 'lmdz') THEN
    147155       IF(ierr.EQ.0) THEN
    148156          ! Continue to read tracer.def
    … …  
    172180       END DO
    173181
    174     ELSE  ! config_inca='aero' ou 'chem'
     182    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
    175183! le module de chimie fournit les noms des traceurs
    176184! et les schemas d'advection associes.
    … …  
    191199       END DO
    192200
    193     END IF ! config_inca
     201    END IF ! type_trac
    194202
    195203!-----------------------------------------------------------------------
    … …  
    319327!
    320328    DEALLOCATE(tnom_0, hadv, vadv)
    321     IF (config_inca /= 'none') DEALLOCATE(tracnam)
     329    DEALLOCATE(tracnam)
    322330
    323331999 FORMAT (i2,1x,i2,1x,a8)

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90

    r1180 r1191  
    2121  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
    2222
    23 ! Variables for INCA
     23! conv_flg(it)=0 : convection desactivated for tracer number it
    2424  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
     25! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it
    2526  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
    2627
     28  CHARACTER(len=4),SAVE :: type_trac
     29 
    2730CONTAINS
    2831
    … …  
    8083    descrq(20)='SLP'
    8184    descrq(30)='PRA'
     85   
     86
     87    IF (config_inca=='none') THEN
     88       type_trac='lmdz'
     89    ELSE
     90       type_trac='inca'
     91    END IF
    8292
    8393!-----------------------------------------------------------------------
    … …  
    8797!
    8898!-----------------------------------------------------------------------
    89     IF (config_inca == 'none') THEN
     99    IF (type_trac == 'lmdz') THEN
    90100       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
    91101       IF(ierr.EQ.0) THEN
    … …  
    109119    END IF
    110120!
    111 ! Allocate variables depending on nqtrue
     121! Allocate variables depending on nqtrue and nbtr
    112122!
    113123    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
    114 
    115     IF (config_inca /= 'none') THEN
    116        ! Varaibles only needed in case of INCA
    117        ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
    118     END IF
    119        
     124    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
     125    conv_flg(:) = 1 ! convection activated for all tracers
     126    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
     127
    120128!-----------------------------------------------------------------------
    121129! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
    … …  
    144152!    Get choice of advection schema from file tracer.def or from INCA
    145153!---------------------------------------------------------------------
    146     IF (config_inca == 'none') THEN
     154    IF (type_trac == 'lmdz') THEN
    147155       IF(ierr.EQ.0) THEN
    148156          ! Continue to read tracer.def
    … …  
    172180       END DO
    173181
    174     ELSE  ! config_inca='aero' ou 'chem'
     182    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
    175183! le module de chimie fournit les noms des traceurs
    176184! et les schemas d'advection associes.
    … …  
    191199       END DO
    192200
    193     END IF ! config_inca
     201    END IF ! type_trac
    194202
    195203!-----------------------------------------------------------------------
    … …  
    319327!
    320328    DEALLOCATE(tnom_0, hadv, vadv)
    321     IF (config_inca /= 'none') DEALLOCATE(tracnam)
     329    DEALLOCATE(tracnam)
    322330
    323331999 FORMAT (i2,1x,i2,1x,a8)

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h

    r776 r1191  
    11!
    2 ! $Header$
     2! $Id$
    33!
    4 C     ----------------------------------------------------------------
    5 C*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
    6 C     ----------------------------------------------------------------
    7 C
    8       COMMON /YOECUMF/
    9      L                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV,
    10      R                 ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP,
    11      R                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON
    12 C
     4!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
     5!                 veillez n'utiliser que des ! pour les commentaires
     6!                 et bien positionner les & des lignes de continuation
     7!                 (les placer en colonne 6 et en colonne 73)
     8!
     9!     ----------------------------------------------------------------
     10!*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
     11!     ----------------------------------------------------------------
     12!
     13      COMMON /YOECUMF/                                                  &
     14     &                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV,              &
     15     &                 ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP,          &
     16     &                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON
     17
    1318      LOGICAL          LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV
    1419      REAL ENTRPEN, ENTRSCV, ENTRMID, ENTRDD
    1520      REAL CMFCTOP, CMFCMAX, CMFCMIN, CMFDEPS, RHCDD, CPRCON
    16 c$OMP THREADPRIVATE(/YOECUMF/)
    17 C
    18 *if (DOC,declared) <> 'UNKNOWN'
    19 C*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
    20 C
    21 C     M.TIEDTKE       E. C. M. W. F.      18/1/89
    22 C
    23 C     NAME      TYPE      PURPOSE
    24 C     ----      ----      -------
    25 C
    26 C     LMFPEN    LOGICAL  TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON
    27 C     LMFSCV    LOGICAL  TRUE IF SHALLOW     CONVECTION IS SWITCHED ON
    28 C     LMFMID    LOGICAL  TRUE IF MIDLEVEL    CONVECTION IS SWITCHED ON
    29 C     LMFDD     LOGICAL  TRUE IF CUMULUS DOWNDRAFT      IS SWITCHED ON
    30 C     LMFDUDV   LOGICAL  TRUE IF CUMULUS FRICTION       IS SWITCHED ON
    31 C     ENTRPEN   REAL     ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
    32 C     ENTRSCV   REAL     ENTRAINMENT RATE FOR SHALLOW CONVECTION
    33 C     ENTRMID   REAL     ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
    34 C     ENTRDD    REAL     ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS
    35 C     CMFCTOP   REAL     RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC
    36 C     CMFCMAX   REAL     MAXIMUM MASSFLUX VALUE ALLOWED FOR
    37 C     CMFCMIN   REAL     MINIMUM MASSFLUX VALUE (FOR SAFETY)
    38 C     CMFDEPS   REAL     FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
    39 C     RHCDD     REAL     RELATIVE SATURATION IN DOWNDRAFTS
    40 C     CPRCON    REAL     COEFFICIENTS FOR DETERMINING CONVERSION
    41 C                        FROM CLOUD WATER TO RAIN
    42 *ifend
    43 C     ----------------------------------------------------------------
     21!$OMP THREADPRIVATE(/YOECUMF/)
     22!
     23!*if (DOC,declared) <> 'UNKNOWN'
     24!*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
     25!
     26!     M.TIEDTKE       E. C. M. W. F.      18/1/89
     27!
     28!     NAME      TYPE      PURPOSE
     29!     ----      ----      -------
     30!
     31!     LMFPEN    LOGICAL  TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON
     32!     LMFSCV    LOGICAL  TRUE IF SHALLOW     CONVECTION IS SWITCHED ON
     33!     LMFMID    LOGICAL  TRUE IF MIDLEVEL    CONVECTION IS SWITCHED ON
     34!     LMFDD     LOGICAL  TRUE IF CUMULUS DOWNDRAFT      IS SWITCHED ON
     35!     LMFDUDV   LOGICAL  TRUE IF CUMULUS FRICTION       IS SWITCHED ON
     36!     ENTRPEN   REAL     ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
     37!     ENTRSCV   REAL     ENTRAINMENT RATE FOR SHALLOW CONVECTION
     38!     ENTRMID   REAL     ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
     39!     ENTRDD    REAL     ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS
     40!     CMFCTOP   REAL     RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC
     41!     CMFCMAX   REAL     MAXIMUM MASSFLUX VALUE ALLOWED FOR
     42!     CMFCMIN   REAL     MINIMUM MASSFLUX VALUE (FOR SAFETY)
     43!     CMFDEPS   REAL     FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
     44!     RHCDD     REAL     RELATIVE SATURATION IN DOWNDRAFTS
     45!     CPRCON    REAL     COEFFICIENTS FOR DETERMINING CONVERSION
     46!                        FROM CLOUD WATER TO RAIN
     47!*ifend
     48!     ----------------------------------------------------------------

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90

    r1188 r1191  
    11!
    2 ! $Header$
     2! $Id $
    33!
    4       SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,
    5      s                  d_tr)
    6       USE dimphy
    7       IMPLICIT none
    8 c======================================================================
    9 c Auteur(s): O. Boucher (LOA/LMD) date: 19961127
    10 c            inspire de clvent
    11 c Objet: diffusion verticale de traceurs avec flux fixe a la surface
    12 c        ou/et flux du type c-drag
    13 c======================================================================
    14 c Arguments:
    15 c dtime----input-R- intervalle du temps (en second)
    16 c coef-----input-R- le coefficient d'echange (m**2/s) l>1
    17 c t--------input-R- temperature (K)
    18 c tr-------input-R- la q. de traceurs
    19 c flux-----input-R- le flux de traceurs a la surface
    20 c paprs----input-R- pression a inter-couche (Pa)
    21 c pplay----input-R- pression au milieu de couche (Pa)
    22 c delp-----input-R- epaisseur de couche (Pa)
    23 c cdrag----input-R- cdrag pour le flux de surface (non active)
    24 c tr0------input-R- traceurs a la surface ou dans l'ocean (non active)
    25 c d_tr-----output-R- le changement de tr
    26 c flux_tr--output-R- flux de tr
    27 c======================================================================
    28 cym#include "dimensions.h"
    29 cym#include "dimphy.h"
    30       REAL dtime
    31       REAL coef(klon,klev)
    32       REAL t(klon,klev), tr(klon,klev)
    33       REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
    34       REAL d_tr(klon,klev)
    35       REAL flux(klon), cdrag(klon), tr0(klon)
    36 c      REAL flux_tr(klon,klev)
    37 c======================================================================
    38 #include "YOMCST.h"
    39 c======================================================================
    40       INTEGER i, k
    41       REAL zx_ctr(klon,2:klev)
    42       REAL zx_dtr(klon,2:klev)
    43       REAL zx_buf(klon)
    44       REAL zx_coef(klon,klev)
    45       REAL local_tr(klon,klev)
    46       REAL zx_alf1(klon), zx_alf2(klon), zx_flux(klon)
    47 c======================================================================
    48       DO k = 1, klev
    49       DO i = 1, klon
    50          local_tr(i,k) = tr(i,k)
    51       ENDDO
    52       ENDDO
    53 c
     4SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,d_tr)
     5  USE dimphy
     6  IMPLICIT NONE
     7!======================================================================
     8! Auteur(s): O. Boucher (LOA/LMD) date: 19961127
     9!            inspire de clvent
     10! Objet: diffusion verticale de traceurs avec flux fixe a la surface
     11!        ou/et flux du type c-drag
     12!
     13! Arguments:
     14!-----------
     15! dtime....input-R- intervalle du temps (en secondes)
     16! coef.....input-R- le coefficient d'echange (m**2/s) l>1
     17! t........input-R- temperature (K)
     18! tr.......input-R- la q. de traceurs
     19! flux.....input-R- le flux de traceurs a la surface
     20! paprs....input-R- pression a inter-couche (Pa)
     21! pplay....input-R- pression au milieu de couche (Pa)
     22! delp.....input-R- epaisseur de couche (Pa)
     23! cdrag....input-R- cdrag pour le flux de surface (non active)
     24! tr0......input-R- traceurs a la surface ou dans l'ocean (non active)
     25! d_tr.....output-R- le changement de tr
     26! flux_tr..output-R- flux de tr
     27!======================================================================
     28  include "YOMCST.h"
     29!
     30! Entree
     31!
     32  REAL,INTENT(IN)                        :: dtime
     33  REAL,DIMENSION(klon,klev),INTENT(IN)   :: coef
     34  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t, tr
     35  REAL,DIMENSION(klon),INTENT(IN)        :: flux !(at/s/m2)
     36  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
     37  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay, delp
     38!
     39! Sorties
     40!
     41  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
     42!  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: flux_tr
     43!
     44! Local
     45!
     46  INTEGER                   :: i, k
     47  REAL,DIMENSION(klon)      :: cdrag, tr0
     48  REAL,DIMENSION(klon,klev) :: zx_ctr
     49  REAL,DIMENSION(klon,klev) :: zx_dtr
     50  REAL,DIMENSION(klon)      :: zx_buf
     51  REAL,DIMENSION(klon,klev) :: zx_coef
     52  REAL,DIMENSION(klon,klev) :: local_tr
     53  REAL,DIMENSION(klon)      :: zx_alf1,zx_alf2,zx_flux
    5454
    55 c======================================================================
    56       DO i = 1, klon
    57          zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
    58          zx_alf2(i) = 1.0 - zx_alf1(i)
    59          zx_flux(i) =  -flux(i)*dtime*RG
    60 c--pour le moment le flux est prescrit
    61 c--cdrag et zx_coef(1) vaut 0
    62          cdrag(i) = 0.0
    63          tr0(i) = 0.0
    64          zx_coef(i,1) = cdrag(i)*dtime*RG
    65       ENDDO
    66 c======================================================================
    67       DO k = 2, klev
    68       DO i = 1, klon
    69          zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
    70      .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
    71          zx_coef(i,k) = zx_coef(i,k)*dtime*RG
    72       ENDDO
    73       ENDDO
    74 c======================================================================
    75       DO i = 1, klon
    76          zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2)
    77          zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+
    78      .                  zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i)
    79          zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) /
    80      .                  zx_buf(i)
    81       ENDDO
    82 c
    83       DO k = 3, klev
    84       DO i = 1, klon
    85          zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
    86      .                  + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1))
    87          zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1)
    88      .                  +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i)
    89          zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i)
    90       ENDDO
    91       ENDDO
    92       DO i = 1, klon
    93          local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev)
    94      .                        +zx_coef(i,klev)*zx_ctr(i,klev) )
    95      .                   / ( delp(i,klev) + zx_coef(i,klev)
    96      .                       -zx_coef(i,klev)*zx_dtr(i,klev) )
    97       ENDDO
    98       DO k = klev-1, 1, -1
    99       DO i = 1, klon
    100          local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1)
    101       ENDDO
    102       ENDDO
    103 c======================================================================
    104 c== flux_tr est le flux de traceur (positif vers bas)
    105 c      DO i = 1, klon
    106 c         flux_tr(i,1) = zx_coef(i,1)/(RG*dtime)
    107 c      ENDDO
    108 c      DO k = 2, klev
    109 c      DO i = 1, klon
    110 c         flux_tr(i,k) = zx_coef(i,k)/(RG*dtime)
    111 c     .               * (local_tr(i,k)-local_tr(i,k-1))
    112 c      ENDDO
    113 c      ENDDO
    114 c======================================================================
    115       DO k = 1, klev
    116       DO i = 1, klon
    117          d_tr(i,k) = local_tr(i,k) - tr(i,k)
    118       ENDDO
    119       ENDDO
    120 c
    121       RETURN
    122       END
     55!======================================================================
     56
     57  DO k = 1, klev
     58     DO i = 1, klon
     59        local_tr(i,k) = tr(i,k)
     60     ENDDO
     61  ENDDO
     62
     63!======================================================================
     64
     65  DO i = 1, klon
     66     zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
     67     zx_alf2(i) = 1.0 - zx_alf1(i)
     68     zx_flux(i) =  -flux(i)*dtime*RG
     69! Pour le moment le flux est prescrit cdrag et zx_coef(1) vaut 0
     70     cdrag(i) = 0.0
     71     tr0(i) = 0.0
     72     zx_coef(i,1) = cdrag(i)*dtime*RG
     73     zx_ctr(i,1)=0.
     74     zx_dtr(i,1)=0.
     75  ENDDO
     76
     77!======================================================================
     78
     79  DO k = 2, klev
     80     DO i = 1, klon
     81        zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))   &
     82             *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
     83        zx_coef(i,k) = zx_coef(i,k)*dtime*RG 
     84     ENDDO
     85  ENDDO
     86
     87!======================================================================
     88
     89  DO i = 1, klon
     90     zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2)
     91     !
     92     zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+                  &
     93          zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i)
     94     !
     95     zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) /   &
     96          zx_buf(i)
     97  ENDDO
     98
     99  DO k = 3, klev
     100     DO i = 1, klon
     101        zx_buf(i) = delp(i,k-1) + zx_coef(i,k)      &
     102             + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1))
     103        zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1)  &
     104             +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i)
     105        zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i)
     106     ENDDO
     107  ENDDO
     108
     109  DO i = 1, klon
     110     local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev) &
     111          +zx_coef(i,klev)*zx_ctr(i,klev) )             &
     112          / ( delp(i,klev) + zx_coef(i,klev)            &
     113          -zx_coef(i,klev)*zx_dtr(i,klev) )
     114  ENDDO
     115
     116  DO k = klev-1, 1, -1
     117     DO i = 1, klon
     118        local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1)
     119     ENDDO
     120  ENDDO
     121
     122!======================================================================
     123!== flux_tr est le flux de traceur (positif vers bas)
     124!      DO i = 1, klon
     125!         flux_tr(i,1) = zx_coef(i,1)/(RG*dtime)
     126!      ENDDO
     127!      DO k = 2, klev
     128!      DO i = 1, klon
     129!         flux_tr(i,k) = zx_coef(i,k)/(RG*dtime)
     130!     .               * (local_tr(i,k)-local_tr(i,k-1))
     131!      ENDDO
     132!      ENDDO
     133!======================================================================
     134  DO k = 1, klev
     135     DO i = 1, klon
     136        d_tr(i,k) = local_tr(i,k) - tr(i,k)
     137     ENDDO
     138  ENDDO
     139 
     140END SUBROUTINE cltrac

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90

    r1188 r1191  
     1!$Id $
     2
     3SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &
     4     cdrag,coef,t,ftsol,pctsrf,               &
     5     tr,trs,paprs,pplay,delp,                 &
     6     masktr,fshtr,hsoltr,tautr,vdeptr,        &
     7     lat,d_tr,d_trs )
     8 
     9  USE dimphy
     10  IMPLICIT NONE
     11!======================================================================
     12! Auteur(s): Alex/LMD) date:  fev 99
     13!            inspire de clqh + clvent
     14! Objet: diffusion verticale de traceurs avec quantite de traceur ds
     15!        le sol ( reservoir de sol de radon )
     16!       
     17! note : pour l'instant le traceur dans le sol et le flux sont
     18!        calcules mais ils ne servent que de diagnostiques
     19!        seule la tendance sur le traceur est sortie (d_tr)
     20!---------------------------------------------------------------------
     21! Arguments:
     22! itr......input-R-  le type de traceur 1- Rn 2 - Pb
     23! dtime....input-R-  intervalle du temps (en secondes) ~ pdtphys
     24! u1lay....input-R-  vent u de la premiere couche (m/s)
     25! v1lay....input-R-  vent v de la premiere couche (m/s)
     26! cdrag....input-R-  cdrag
     27! coef.....input-R-  le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
     28! t........input-R-  temperature (K)
     29! paprs....input-R-  pression a inter-couche (Pa)
     30! pplay....input-R-  pression au milieu de couche (Pa)
     31! delp.....input-R-  epaisseur de couche (Pa)
     32! ftsol....input-R-  temperature du sol (en Kelvin)
     33! tr.......input-R-  traceurs
     34! trs......input-R-  traceurs dans le sol
     35! masktr...input-R-  Masque reservoir de sol traceur (1 = reservoir)
     36! fshtr....input-R-  Flux surfacique de production dans le sol
     37! tautr....input-R-  Constante de decroissance du traceur
     38! vdeptr...input-R-  Vitesse de depot sec dans la couche brownienne
     39! hsoltr...input-R-  Epaisseur equivalente du reservoir de sol
     40! lat......input-R-  latitude en degree
     41! d_tr.....output-R- le changement de "tr"
     42! d_trs....output-R- le changement de "trs"
     43!======================================================================
     44  include "YOMCST.h"
     45  include "indicesol.h"
    146!
    2       SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay,
    3      e              cdrag,coef,t,ftsol,pctsrf,
    4      e              tr,trs,paprs,pplay,delp,
    5      e              masktr,fshtr,hsoltr,tautr,vdeptr,
    6      e              lat,
    7      s              d_tr,d_trs )
    8 
    9       USE dimphy
    10       IMPLICIT none
    11 c======================================================================
    12 c Auteur(s): Alex/LMD) date:  fev 99
    13 c            inspire de clqh + clvent
    14 c Objet: diffusion verticale de traceurs avec quantite de traceur ds
    15 c        le sol ( reservoir de sol de radon )
    16 c       
    17 c note : pour l'instant le traceur dans le sol et le flux sont
    18 c        calcules mais ils ne servent que de diagnostiques
    19 c        seule la tendance sur le traceur est sortie (d_tr)
    20 c======================================================================
    21 c Arguments:
    22 c itr---  -input-R- le type de traceur 1- Rn 2 - Pb
    23 c dtime----input-R- intervalle du temps (en second)
    24 c u1lay----input-R- vent u de la premiere couche (m/s)
    25 c v1lay----input-R- vent v de la premiere couche (m/s)
    26 c cdrag----input-R- cdrag
    27 c coef-----input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
    28 c t--------input-R- temperature (K)
    29 c paprs----input-R- pression a inter-couche (Pa)
    30 c pplay----input-R- pression au milieu de couche (Pa)
    31 c delp-----input-R- epaisseur de couche (Pa)
    32 c ftsol----input-R- temperature du sol (en Kelvin)
    33 c tr-------input-R- traceurs
    34 c trs------input-R- traceurs dans le sol
    35 c masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
    36 c fshtr----input-R- Flux surfacique de production dans le sol
    37 c tautr----input-R- Constante de decroissance du traceur
    38 c vdeptr---input-R- Vitesse de depot sec dans la couche brownienne
    39 c hsoltr---input-R- Epaisseur equivalente du reservoir de sol
    40 c lat-----input-R- latitude en degree
    41 c d_tr-----output-R- le changement de "tr"
    42 c d_trs----output-R- le changement de "trs"
    43 c======================================================================
    44 cym#include "dimensions.h"
    45 cym#include "dimphy.h"
    46 #include "YOMCST.h"
    47 #include "indicesol.h"
    48 c======================================================================
    49       REAL dtime
    50       REAL u1lay(klon), v1lay(klon)
    51       REAL cdrag(klon)
    52       REAL coef(klon,klev)
    53       REAL t(klon,klev), ftsol(klon,nbsrf), pctsrf(klon,nbsrf)
    54       REAL tr(klon,klev), trs(klon)
    55       REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
    56       REAL masktr(klon)
    57       REAL fshtr(klon)
    58       REAL hsoltr
    59       REAL tautr
    60       REAL vdeptr
    61       REAL lat(klon)   
    62       REAL d_tr(klon,klev)
    63 c======================================================================
    64       REAL flux_tr(klon,klev)  ! (diagnostic) flux de traceur
    65       REAL d_trs(klon)         ! (diagnostic) traceur ds le sol
    66 c======================================================================
    67       INTEGER i, k, itr, n, l
    68       REAL rotrhi(klon)
    69       REAL zx_coef(klon,klev)
    70       REAL zx_buf(klon)
    71       REAL zx_ctr(klon,klev)
    72       REAL zx_dtr(klon,klev)
    73       REAL zx_trs(klon)
    74       REAL zx_a, zx_b
    75 
    76       REAL local_tr(klon,klev)
    77       REAL local_trs(klon)
    78       REAL zts(klon)
    79       REAL zx_alpha1(klon), zx_alpha2(klon)
    80 c======================================================================
    81 cAA Pour l'instant les 4 types de surface ne sont pas pris en compte
    82 cAA On fabrique avec zts un champ de temperature de sol 
    83 cAA que le pondere par la fraction de nature de sol.
    84 c
    85 c      print*,'PASSAGE DANS CLTRACRN'
    86 
    87       DO i = 1,klon
    88             zts(i) = 0.
    89       ENDDO
    90 c
    91       DO n=1,nbsrf
    92          DO i = 1,klon
    93             zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
    94          ENDDO
    95       ENDDO
    96 c
    97       DO i = 1,klon
    98           rotrhi(i) = RD * zts(i) / hsoltr
    99       END DO
    100 c
    101       DO k = 1, klev
    102       DO i = 1, klon
    103          local_tr(i,k) = tr(i,k)
    104       ENDDO
    105       ENDDO
    106 c
    107       DO i = 1, klon
    108          local_trs(i) = trs(i)
    109       ENDDO
    110 c======================================================================
    111 cAA   Attention si dans clmain zx_alf1(i) = 1.0
    112 cAA   Il doit y avoir coherence (dc la meme chose ici)
    113 
    114       DO i = 1, klon
    115 cAA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
    116          zx_alpha1(i) = 1.0
    117          zx_alpha2(i) = 1.0 - zx_alpha1(i)
    118       ENDDO
    119 c======================================================================
    120       DO i = 1, klon
    121          zx_coef(i,1) = cdrag(i)
    122      .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
    123      .                 * pplay(i,1)/(RD*t(i,1))
    124          zx_coef(i,1) = zx_coef(i,1) * dtime*RG
    125       ENDDO
    126 c
    127       DO k = 2, klev
    128       DO i = 1, klon
    129          zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
    130      .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
    131          zx_coef(i,k) = zx_coef(i,k) * dtime*RG
    132       ENDDO
    133       ENDDO
    134 c======================================================================
    135       DO i = 1, klon
    136          zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)
    137          zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
    138          zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
    139       ENDDO
    140 c
    141       DO l = klev-1, 2 , -1
    142       DO i = 1, klon
    143          zx_buf(i) = delp(i,l)+zx_coef(i,l)
    144      .              +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
    145          zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l)
    146      .                   + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
    147          zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
    148       ENDDO
    149       ENDDO
    150 c
    151       DO i = 1, klon
    152          zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2))
    153      .               + masktr(i) * zx_coef(i,1)
    154      .               *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )
    155          zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1)
    156      .                  + zx_ctr(i,2)
    157      .                  *(zx_coef(i,2)
    158      .                    - masktr(i) * zx_coef(i,1)
    159      .                    *zx_alpha2(i) ) ) / zx_buf(i)
    160          zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
    161       ENDDO
    162 c======================================================================
    163 c Calculer d'abord local_trs nouvelle quantite dans le reservoir
    164 c de sol
    165 c
    166 c-------------------------
    167 c Au dessus des continents
    168 c-------------------------
    169 c Le pb peut se deposer partout : vdeptr = 10-3 m/s
    170 c Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
    171 c
    172       DO i = 1, klon
    173 c
    174         IF ( NINT(masktr(i)) .EQ. 1 ) THEN
    175          zx_trs(i) = local_trs(i)
    176          zx_a = zx_trs(i)
    177      .         +fshtr(i)*dtime*rotrhi(i)
    178      .         +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
    179      .         *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))
    180      .         +zx_alpha2(i)*zx_ctr(i,2))
    181          zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
    182      .              * (1.-zx_dtr(i,1)
    183      .              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))
    184      .              + dtime / tautr
    185 cAA: Pour l'instant, pour aller vite, le depot sec est traite
    186 C comme une decroissance
    187      .              + dtime * vdeptr / hsoltr
    188          zx_trs(i) = zx_a / zx_b
    189          local_trs(i) = zx_trs(i)
    190         ENDIF
    191 c
    192 c Si on est entre 60N et 70N on divise par 2 l'emanation
    193 c--------------------------------------------------------
    194 c
    195         IF
    196      .  ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60.
    197      .     .AND.lat(i).LE.70.)
    198      .                      .OR.
    199      .     (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60.
    200      .     .AND.lat(i).LE.70.) )
    201      .  THEN
    202          zx_trs(i) = local_trs(i)
    203          zx_a = zx_trs(i)
    204      .         +(fshtr(i)/2.)*dtime*rotrhi(i)
    205      .         +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
    206      .         *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))
    207      .         +zx_alpha2(i)*zx_ctr(i,2))
    208          zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
    209      .              * (1.-zx_dtr(i,1)
    210      .              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))
    211      .              + dtime / tautr
    212      .              + dtime * vdeptr / hsoltr
    213          zx_trs(i) = zx_a / zx_b
    214          local_trs(i) = zx_trs(i)
    215         ENDIF
    216 c
    217 c----------------------------------------------
    218 c Au dessus des oceans et aux hautes latitudes
    219 c----------------------------------------------
    220 c
    221 c au dessous de -60S  pas d'emission de radon au dessus
    222 c des oceans et des continents
    223 c---------------------------------------------------------------
    224 
    225        IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0)
    226      .                 .OR.
    227      . (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.))
    228      . THEN
    229          zx_trs(i) = 0.
    230          local_trs(i) = 0.
    231         END IF
    232 
    233 c au dessus de 70 N pas d'emission de radon au dessus
    234 c des oceans et des continents
    235 c--------------------------------------------------------------
    236        IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0)
    237      .                 .OR.
    238      . (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.))
    239      . THEN
    240          zx_trs(i) = 0.
    241          local_trs(i) = 0.
    242         END IF
    243 
    244 c Au dessus des oceans la source est nulle
    245 c-----------------------------------------
    246 c
    247         IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
    248          zx_trs(i) = 0.
    249          local_trs(i) = 0.
    250         END IF
    251 c
    252       ENDDO    ! sur le i=1,klon
    253 c
    254 c======================================================================
    255 c==== une fois on a zx_trs, on peut faire l'iteration ========
    256 c
    257       DO i = 1, klon
    258          local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
    259       ENDDO
    260       DO l = 2, klev
    261       DO i = 1, klon
    262         local_tr(i,l)
    263      .    = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
    264       ENDDO
    265       ENDDO
    266 c======================================================================
    267 c== Calcul du flux de traceur (flux_tr): UA/(m**2 s)
    268 c
    269       DO i = 1, klon
    270         flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG
    271      .      * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2)
    272      .                   -zx_trs(i)) / dtime
    273       ENDDO
    274       DO l = 2, klev
    275       DO i = 1, klon
    276         flux_tr(i,l) = zx_coef(i,l)/RG
    277      .       * (local_tr(i,l)-local_tr(i,l-1)) / dtime
    278       ENDDO
    279       ENDDO
    280 c======================================================================
    281 c== Calcul des tendances du traceur ds le sol et dans l'atmosphere
    282 c
    283       DO l = 1, klev
    284       DO i = 1, klon
    285          d_tr(i,l) = local_tr(i,l) - tr(i,l)
    286       ENDDO
    287       ENDDO
    288       DO i = 1, klon
    289          d_trs(i) = local_trs(i) - trs(i)
    290       ENDDO
    291 c======================================================================
    292 c
    293       RETURN
    294       END
     47!Entrees
     48  INTEGER,INTENT(IN)                     :: itr
     49  REAL,INTENT(IN)                        :: dtime
     50  REAL,DIMENSION(klon),INTENT(IN)        :: u1lay, v1lay
     51  REAL,DIMENSION(klon),INTENT(IN)        :: cdrag
     52  REAL,DIMENSION(klon,klev),INTENT(IN)   :: coef, t
     53  REAL,DIMENSION(klon,nbsrf),INTENT(IN)  :: ftsol, pctsrf
     54  REAL,DIMENSION(klon,klev),INTENT(IN)   :: tr
     55  REAL,DIMENSION(klon),INTENT(IN)        :: trs
     56  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
     57  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay, delp
     58  REAL,DIMENSION(klon),INTENT(IN)        :: masktr
     59  REAL,DIMENSION(klon),INTENT(IN)        :: fshtr
     60  REAL,INTENT(IN)                        :: hsoltr
     61  REAL,INTENT(IN)                        :: tautr
     62  REAL,INTENT(IN)                        :: vdeptr
     63  REAL,DIMENSION(klon),INTENT(IN)        :: lat   
     64
     65!Sorties
     66  REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
     67  REAL,DIMENSION(klon),INTENT(OUT) :: d_trs  ! (diagnostic) traceur ds le sol
     68
     69!Locales
     70  REAL,DIMENSION(klon,klev) :: flux_tr  ! (diagnostic) flux de traceur
     71  INTEGER                   :: i, k, n, l
     72  REAL,DIMENSION(klon)      :: rotrhi
     73  REAL,DIMENSION(klon,klev) :: zx_coef
     74  REAL,DIMENSION(klon)      :: zx_buf
     75  REAL,DIMENSION(klon,klev) :: zx_ctr
     76  REAL,DIMENSION(klon,klev) :: zx_dtr
     77  REAL,DIMENSION(klon)      :: zx_trs
     78  REAL                      :: zx_a, zx_b
     79 
     80  REAL,DIMENSION(klon,klev) :: local_tr
     81  REAL,DIMENSION(klon)      :: local_trs
     82  REAL,DIMENSION(klon)      :: zts      ! champ de temperature du sol
     83  REAL,DIMENSION(klon)      :: zx_alpha1, zx_alpha2
     84!======================================================================
     85!AA Pour l'instant les 4 types de surface ne sont pas pris en compte
     86!AA On fabrique avec zts un champ de temperature de sol 
     87!AA que le pondere par la fraction de nature de sol.
     88 
     89  DO i = 1,klon
     90     zts(i) = 0.
     91  ENDDO
     92
     93  DO n=1,nbsrf
     94     DO i = 1,klon
     95        zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
     96     ENDDO
     97  ENDDO
     98
     99  DO i = 1,klon
     100     rotrhi(i) = RD * zts(i) / hsoltr
     101  ENDDO
     102
     103  DO k = 1, klev
     104     DO i = 1, klon
     105        local_tr(i,k) = tr(i,k)
     106     ENDDO
     107  ENDDO
     108
     109  DO i = 1, klon
     110     local_trs(i) = trs(i)
     111  ENDDO
     112!======================================================================
     113!AA   Attention si dans clmain zx_alf1(i) = 1.0
     114!AA   Il doit y avoir coherence (dc la meme chose ici)
     115
     116  DO i = 1, klon
     117!AA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
     118     zx_alpha1(i) = 1.0
     119     zx_alpha2(i) = 1.0 - zx_alpha1(i)
     120  ENDDO
     121!======================================================================
     122  DO i = 1, klon
     123     zx_coef(i,1) = cdrag(i)*(1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
     124          *pplay(i,1)/(RD*t(i,1))
     125     zx_coef(i,1) = zx_coef(i,1) * dtime*RG
     126  ENDDO
     127
     128  DO k = 2, klev
     129     DO i = 1, klon
     130        zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
     131             *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
     132        zx_coef(i,k) = zx_coef(i,k) * dtime*RG
     133     ENDDO
     134  ENDDO
     135!======================================================================
     136  DO i = 1, klon
     137     zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)
     138     zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
     139     zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
     140  ENDDO
     141
     142  DO l = klev-1, 2 , -1
     143     DO i = 1, klon
     144        zx_buf(i) = delp(i,l)+zx_coef(i,l)      &
     145             +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
     146 
     147        zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) &
     148             + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
     149        zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
     150     ENDDO
     151  ENDDO
     152
     153  DO i = 1, klon
     154     zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2))  &
     155          + masktr(i) * zx_coef(i,1)                        &
     156          *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )
     157
     158     zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1)                &
     159          + zx_ctr(i,2)                                     &
     160          *(zx_coef(i,2)                                    &
     161          - masktr(i) * zx_coef(i,1)                        &
     162          *zx_alpha2(i) ) ) / zx_buf(i)
     163     zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
     164  ENDDO
     165!======================================================================
     166! Calculer d'abord local_trs nouvelle quantite dans le reservoir
     167! de sol
     168!=====================================================================
     169
     170  DO i = 1, klon
     171!-------------------------
     172! Au dessus des continents
     173!--
     174! Le pb peut se deposer partout : vdeptr = 10-3 m/s
     175! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
     176!-------------------------------------------------------------------
     177     IF ( NINT(masktr(i)) .EQ. 1 ) THEN
     178        zx_trs(i) = local_trs(i)
     179        zx_a = zx_trs(i)                                           &
     180             +fshtr(i)*dtime*rotrhi(i)                             &
     181             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
     182             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
     183             +zx_alpha2(i)*zx_ctr(i,2))
     184! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance
     185        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG            &
     186             * (1.-zx_dtr(i,1)                                     &
     187             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))             &
     188             + dtime / tautr                                       &
     189             + dtime * vdeptr / hsoltr
     190        zx_trs(i) = zx_a / zx_b
     191        local_trs(i) = zx_trs(i)
     192     ENDIF
     193!--------------------------------------------------------
     194! Si on est entre 60N et 70N on divise par 2 l'emanation
     195!--------------------------------------------------------
     196
     197     IF ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR.      &
     198          (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN
     199        zx_trs(i) = local_trs(i)
     200        zx_a = zx_trs(i)                                           &
     201             +(fshtr(i)/2.)*dtime*rotrhi(i)                        &
     202             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
     203             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
     204             +zx_alpha2(i)*zx_ctr(i,2))
     205        !
     206        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG  &
     207             * (1.-zx_dtr(i,1)                           &
     208             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))   &
     209             + dtime / tautr                             &
     210             + dtime * vdeptr / hsoltr
     211        !
     212        zx_trs(i) = zx_a / zx_b
     213        local_trs(i) = zx_trs(i)
     214     ENDIF
     215
     216!----------------------------------------------
     217! Au dessus des oceans et aux hautes latitudes
     218!--
     219! au dessous de -60S  pas d'emission de radon au dessus
     220! des oceans et des continents
     221!---------------------------------------------------------------
     222
     223     IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR.       &
     224          (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
     225        zx_trs(i) = 0.
     226        local_trs(i) = 0.
     227     END IF
     228!--
     229! au dessus de 70 N pas d'emission de radon au dessus
     230! des oceans et des continents
     231!--------------------------------------------------------------
     232     IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR.    &
     233          (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
     234        zx_trs(i) = 0.
     235        local_trs(i) = 0.
     236     END IF
     237!---------------------------------------------
     238! Au dessus des oceans la source est nulle
     239!--------------------------------------------
     240
     241     IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
     242        zx_trs(i) = 0.
     243        local_trs(i) = 0.
     244     END IF
     245
     246  ENDDO    ! sur le i=1,klon
     247!
     248!======================================================================
     249! Une fois on a zx_trs, on peut faire l'iteration
     250!======================================================================
     251
     252  DO i = 1, klon
     253     local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
     254  ENDDO
     255  DO l = 2, klev
     256     DO i = 1, klon
     257        local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
     258     ENDDO
     259  ENDDO
     260!======================================================================
     261! Calcul du flux de traceur (flux_tr): UA/(m**2 s)
     262!======================================================================
     263  DO i = 1, klon
     264     flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG                      &
     265          * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2) &
     266          -zx_trs(i)) / dtime
     267  ENDDO
     268  DO l = 2, klev
     269     DO i = 1, klon
     270        flux_tr(i,l) = zx_coef(i,l)/RG                    &
     271             * (local_tr(i,l)-local_tr(i,l-1)) / dtime
     272     ENDDO
     273  ENDDO
     274!======================================================================
     275! Calcul des tendances du traceur ds le sol et dans l'atmosphere
     276!======================================================================
     277  DO l = 1, klev
     278     DO i = 1, klon
     279        d_tr(i,l) = local_tr(i,l) - tr(i,l)
     280     ENDDO
     281  ENDDO
     282  DO i = 1, klon
     283     d_trs(i) = local_trs(i) - trs(i)
     284  ENDDO
     285
     286END SUBROUTINE cltracrn

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90

    r1188 r1191  
    1 c
    2 c $Header$
    3 c
    4       SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx)
    5       USE dimphy
    6       IMPLICIT NONE
    7 c=====================================================================
    8 c Objet : convection des traceurs / KE
    9 c Auteurs: M-A Filiberti and J-Y Grandpeix
    10 c=====================================================================
    11 c
    12 cym#include "dimensions.h"
    13 cym#include "dimphy.h"
    14 #include "YOMCST.h"
    15 #include "YOECUMF.h"
    16 c
    17       REAL pdtime
    18       REAL paprs(klon,klev+1) ! pression aux 1/2 couches (bas en haut)
    19       REAL pplay(klon,klev)  ! pression pour le milieu de chaque couche
    20       REAL x(klon,klev)        ! q de traceur (bas en haut)
    21       REAL dx(klon,klev)     ! tendance de traceur  (bas en haut)
    22       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
    23       REAL upd(klon,klev)      ! saturated updraft mass flux
    24       REAL dnd(klon,klev)      ! saturated downdraft mass flux
    25 c
    26 c--variables locales     
    27       real zed(klon,klev),zmd(klon,klev,klev)
    28       real za(klon,klev,klev)
    29       real zmfd(klon,klev),zmfa(klon,klev)
    30       real zmfp(klon,klev),zmfu(klon,klev)
    31       integer i,k,j
    32 c test conservation
    33 c      real conserv
    34 c =========================================
    35 c calcul des tendances liees au downdraft
    36 c =========================================
    37       zed(:,:)=0.
    38       zmfd(:,:)=0.
    39       zmfa(:,:)=0.
    40       zmfu(:,:)=0.
    41       zmfp(:,:)=0.
    42       zmd(:,:,:)=0.
    43       za(:,:,:)=0.
    44 c entrainement
    45       do k=1,klev-1
    46         do i=1,klon
    47           zed(i,k)=max(0.,mp(i,k)-mp(i,k+1))
    48         end do
    49       end do
    50 c
    51 c calcul de la matrice d echange
    52 c matrice de distribution de la masse entrainee en k
    53 c
    54       do k=1,klev
    55         do i=1,klon
    56           zmd(i,k,k)=zed(i,k)
    57         end do
    58       end do
    59       do k=2,klev
    60         do j=k-1,1,-1
    61           do i=1,klon
    62           if(mp(i,j+1).ne.0) then
    63           zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1))
    64           endif
    65           end do
    66         end do
    67       end do
    68       do k=1,klev
    69         do j=1,klev-1
    70           do i=1,klon
    71           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
    72           end do
    73         end do
    74       end do
    75 c
    76 c rajout du terme lie a l ascendance induite
    77 c
    78         do j=2,klev
    79          do i=1,klon
    80           za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
    81          end do
    82         end do
    83 C
    84 c tendances
    85 c           
    86       do k=1,klev
    87         do j=1,klev
    88           do i=1,klon
    89           zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j))
    90           end do
    91         end do
    92       end do
    93 c
    94 c =========================================
    95 c calcul des tendances liees aux flux satures
    96 c =========================================
    97       do j=1,klev
    98         do i=1,klon
    99           zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j))
    100         end do
    101       end do
    102       do k=1,klev
    103         do j=1,klev
    104           do i=1,klon
    105           zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j))
    106           end do
    107         end do
    108       end do
    109       do j=1,klev-1
    110         do i=1,klon
    111           zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j))
    112         end do
    113       end do
    114       do j=2,klev
    115         do i=1,klon
    116           zmfu(i,j)=zmfu(i,j)
    117      .             +min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1))
    118         end do
    119       end do
     1!
     2! $Id $
     3!
     4SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx)
     5  USE dimphy
     6  IMPLICIT NONE
     7!=====================================================================
     8! Objet : convection des traceurs / KE
     9! Auteurs: M-A Filiberti and J-Y Grandpeix
     10!=====================================================================
    12011
    121 c =========================================
    122 c--calcul final des tendances
    123 c =========================================
    124       do k=1, klev
    125         do i=1, klon
    126           dx(i,k)=(zmfd(i,k)+zmfu(i,k)
    127      .      +zmfa(i,k)+zmfp(i,k))*pdtime
    128      .      *RG/(paprs(i,k)-paprs(i,k+1))
    129 c          print*,'dx',k,dx(i,k)
    130         enddo
    131       enddo
    132 c
    133 c test de conservation du traceur
    134 c      conserv=0.
    135 c      do k=1, klev
    136 c        do i=1, klon
    137 c         conserv=conserv+dx(i,k)*
    138 c     .     (paprs(i,k)-paprs(i,k+1))/RG
    139 C
    140 c        enddo
    141 c      enddo
    142 c      print *,'conserv',conserv
     12  include "YOMCST.h"
     13  include "YOECUMF.h"
     14
     15! Entree
     16  REAL,INTENT(IN)                           :: pdtime
     17  REAL,DIMENSION(klon,klev),INTENT(IN)      :: da
     18  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
     19  REAL,DIMENSION(klon,klev),INTENT(IN)      :: mp
     20  REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs ! pression aux 1/2 couches (bas en haut)
     21  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pplay ! pression pour le milieu de chaque couche
     22  REAL,DIMENSION(klon,klev),INTENT(IN)      :: x     ! q de traceur (bas en haut)
     23  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
     24  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
     25
     26! Sortie
     27  REAL,DIMENSION(klon,klev),INTENT(OUT) :: dx ! tendance de traceur  (bas en haut)
     28
     29! Variables locales     
     30  REAL,DIMENSION(klon,klev)       :: zed
     31  REAL,DIMENSION(klon,klev,klev)  :: zmd
     32  REAL,DIMENSION(klon,klev,klev)  :: za
     33  REAL,DIMENSION(klon,klev)       :: zmfd,zmfa
     34  REAL,DIMENSION(klon,klev)       :: zmfp,zmfu
     35  INTEGER                         :: i,k,j
     36
     37! =========================================
     38! calcul des tendances liees au downdraft
     39! =========================================
     40  zed(:,:)=0.
     41  zmfd(:,:)=0.
     42  zmfa(:,:)=0.
     43  zmfu(:,:)=0.
     44  zmfp(:,:)=0.
     45  zmd(:,:,:)=0.
     46  za(:,:,:)=0.
     47! entrainement
     48  DO k=1,klev-1
     49     DO i=1,klon
     50        zed(i,k)=max(0.,mp(i,k)-mp(i,k+1))
     51     END DO
     52  END DO
     53
     54! calcul de la matrice d echange
     55! matrice de distribution de la masse entrainee en k
     56
     57  DO k=1,klev
     58     DO i=1,klon
     59        zmd(i,k,k)=zed(i,k)
     60     END DO
     61  END DO
     62  DO k=2,klev
     63     DO j=k-1,1,-1
     64        DO i=1,klon
     65           if(mp(i,j+1).ne.0) then
     66              zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1))
     67           ENDif
     68        END DO
     69     END DO
     70  END DO
     71  DO k=1,klev
     72     DO j=1,klev-1
     73        DO i=1,klon
     74           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
     75        END DO
     76     END DO
     77  END DO
     78!
     79! rajout du terme lie a l ascendance induite
     80!
     81  DO j=2,klev
     82     DO i=1,klon
     83        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
     84     END DO
     85  END DO
     86!
     87! tendances
     88!           
     89  DO k=1,klev
     90     DO j=1,klev
     91        DO i=1,klon
     92           zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j))
     93        END DO
     94     END DO
     95  END DO
     96!
     97! =========================================
     98! calcul des tendances liees aux flux satures
     99! =========================================
     100  DO j=1,klev
     101     DO i=1,klon
     102        zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j))
     103     END DO
     104  END DO
     105  DO k=1,klev
     106     DO j=1,klev
     107        DO i=1,klon
     108           zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j))
     109        END DO
     110     END DO
     111  END DO
     112  DO j=1,klev-1
     113     DO i=1,klon
     114        zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j))
     115     END DO
     116  END DO
     117  DO j=2,klev
     118     DO i=1,klon
     119        zmfu(i,j)=zmfu(i,j)+min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1))
     120     END DO
     121  END DO
     122
     123! =========================================
     124! calcul final des tendances
     125! =========================================
     126  DO k=1, klev
     127     DO i=1, klon
     128        dx(i,k)=(zmfd(i,k)+zmfu(i,k)       &
     129             +zmfa(i,k)+zmfp(i,k))*pdtime  &
     130             *RG/(paprs(i,k)-paprs(i,k+1))
     131        !          print*,'dx',k,dx(i,k)
     132     ENDDO
     133  ENDDO
     134
     135! test de conservation du traceur
     136!      conserv=0.
     137!      DO k=1, klev
     138!        DO i=1, klon
     139!         conserv=conserv+dx(i,k)*   &
     140!        (paprs(i,k)-paprs(i,k+1))/RG
     141!        ENDDO
     142!      ENDDO
     143!      print *,'conserv',conserv
    143144     
    144       return
    145       end
     145END SUBROUTINE cvltr

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h

    r1115 r1191  
    11!
    2 ! $Header$
     2! $Id $
    33!
    4       IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
    5 c$OMP MASTER
    6          CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
    7 c
    8          CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,
    9      .                 nhori, nid_tra)
    10          CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
    11      .                 klev, presnivs, nvert)
     4#ifdef CPP_IOIPSL
    125
     6  IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
     7!$OMP MASTER
     8     CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
     9     CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,nhori, nid_tra)
     10     CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",klev, presnivs, nvert)
    1311
     12     zsto = pdtphys
     13     zout = ecrit_tra
     14     CALL histdef(nid_tra, "phis", "Surface geop. height", "-",   &
     15          iim,jj_nb,nhori, 1,1,1, -99, 32,"once",  zsto,zout)
     16     CALL histdef(nid_tra, "aire", "Grid area", "-",              &
     17          iim,jj_nb,nhori, 1,1,1, -99, 32,"once",  zsto,zout)
    1418
    15          zsto = pdtphys
    16          zout = ecrit_tra
    17 c
    18          CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
    19      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    20      .                "once",  zsto,zout)
    21 c
    22          CALL histdef(nid_tra, "aire", "Grid area", "-",
    23      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    24      .                "once",  zsto,zout)
    25          DO it=1,nbtr
    26 C champ 2D
    27          iq=it+2
    28          iiq=niadv(iq)
    29          CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga",
    30      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    31      .                "ave(X)", zsto,zout)
    32          if (lessivage) THEN
    33          CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq),
    34      .              "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    35      .              "ave(X)", zsto,zout)
    36          endif
     19!TRACEURS
     20!----------------
     21     DO it = 1,nbtr
     22        iiq = niadv(it+2)
    3723
    38 c---Ajout Olivia
    39          CALL histdef(nid_tra, "d_tr_th_"//tname(iiq),
    40      .                "tendance thermique"// ttext(iiq), "?",
    41      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    42      .                "ave(X)", zsto,zout)
    43 c
    44          if(iflag_con.GE.2) then
    45          CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq),
    46      .                "tendance convection"// ttext(iiq), "?",
    47      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    48      .                "ave(X)", zsto,zout)
    49              endif !(iflag_con.GE.2) then
    50          CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq),
    51      .                "tendance couche limite"// ttext(iiq), "?",
    52      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    53      .                "ave(X)", zsto,zout)
    54 c---fin Olivia   
     24! CONCENTRATIONS
     25        CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga",    &
     26             iim,jj_nb,nhori, klev,1,klev,nvert, 32,"ave(X)", zsto,zout)
    5527
    56          ENDDO
     28! TD LESSIVAGE
     29        IF (lessivage .AND. aerosol(it)) THEN
     30           CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq), &
     31                "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     32                "ave(X)", zsto,zout)
     33        END IF
    5734
    58          CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
    59      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    60      .                "inst(X)",  zout,zout)
     35! TD THERMIQUES
     36        IF (iflag_thermals.gt.0) THEN
     37           CALL histdef(nid_tra, "d_tr_th_"//tname(iiq),      &
     38                "tendance thermique"// ttext(iiq), "?",       &
     39                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     40                "ave(X)", zsto,zout)
     41        ENDIF
    6142
    62          CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
    63      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    64      .                "inst(X)",  zout,zout)
    65          CALL histdef(nid_tra, "psrf1", "nature sol", "-",
    66      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    67      .                "inst(X)",  zout,zout)
    68          CALL histdef(nid_tra, "psrf2", "nature sol", "-",
    69      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    70      .                "inst(X)",  zout,zout)
    71          CALL histdef(nid_tra, "psrf3", "nature sol", "-",
    72      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    73      .                "inst(X)",  zout,zout)
    74          CALL histdef(nid_tra, "psrf4", "nature sol", "-",
    75      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    76      .                "inst(X)",  zout,zout)
    77          CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
    78      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    79      .                "inst(X)",  zout,zout)
    80          CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
    81      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    82      .                "inst(X)",  zout,zout)
    83          CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
    84      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    85      .                "inst",  zout,zout)
    86          CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
    87      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    88      .                "inst(X)",  zout,zout)
    89          CALL histdef(nid_tra, "pplay", "flux u mont","-",
    90      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    91      .                "inst(X)", zout,zout)
    92          CALL histdef(nid_tra, "t", "flux u mont","-",
    93      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    94      .                "inst(X)", zout,zout)
    95          CALL histdef(nid_tra, "mfu", "flux u mont","-",
    96      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    97      .                "ave(X)", zsto,zout)
    98          CALL histdef(nid_tra, "mfd", "flux u decen","-",
    99      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    100      .                "ave(X)", zsto,zout)
    101          CALL histdef(nid_tra, "en_u", "flux u mont","-",
    102      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    103      .                "ave(X)", zsto,zout)
    104          CALL histdef(nid_tra, "en_d", "flux u mont","-",
    105      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    106      .                "ave(X)", zsto,zout)
    107          CALL histdef(nid_tra, "de_d", "flux u mont","-",
    108      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    109      .                "ave(X)", zsto,zout)
    110          CALL histdef(nid_tra, "de_u", "flux u decen","-",
    111      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    112      .                "ave(X)", zsto,zout)
    113          CALL histdef(nid_tra, "coefh", "turbulent coef","-",
    114      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    115      .                "ave(X)", zsto,zout)
     43! TD CONVECTION
     44        IF (iflag_con.GE.2) THEN
     45           CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq),   &
     46                "tendance convection"// ttext(iiq), "?",   &
     47                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     48                "ave(X)", zsto,zout)
     49        ENDIF
    11650
    117 c
    118          CALL histend(nid_tra)
    119          ndex2d = 0
    120          ndex3d = 0
    121          ndex = 0
    122 c$OMP END MASTER
    123       END IF ! ecrit_tra>0. .AND. config_inca == 'none'
     51! TD COUCHE-LIMITE
     52        CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq),      &
     53             "tendance couche limite"// ttext(iiq), "?",   &
     54             iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     55             "ave(X)", zsto,zout)
     56     ENDDO
     57!---------------   
     58!
     59! VENT (niveau 1)
     60     CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",      &
     61          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     62          "inst(X)",  zout,zout)     
     63     CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",      &
     64          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     65          "inst(X)",  zout,zout)
     66
     67! TEMPERATURE DU SOL
     68     CALL histdef(nid_tra, "ftsol1", "temper sol", "-",    &
     69          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     70          "inst(X)",  zout,zout)
     71     CALL histdef(nid_tra, "ftsol2", "temper sol", "-",    &
     72          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     73          "inst(X)",  zout,zout)
     74     CALL histdef(nid_tra, "ftsol3", "temper sol", "-",    &
     75          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     76          "inst",  zout,zout)
     77     CALL histdef(nid_tra, "ftsol4", "temper sol", "-",    &
     78          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     79          "inst(X)",  zout,zout)
     80
     81! NATURE DU SOL
     82     CALL histdef(nid_tra, "psrf1", "nature sol", "-",     &
     83          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     84          "inst(X)",  zout,zout)
     85     CALL histdef(nid_tra, "psrf2", "nature sol", "-",     &
     86          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     87          "inst(X)",  zout,zout)
     88     CALL histdef(nid_tra, "psrf3", "nature sol", "-",     &
     89          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     90          "inst(X)",  zout,zout)
     91     CALL histdef(nid_tra, "psrf4", "nature sol", "-",     &
     92          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     93          "inst(X)",  zout,zout)
     94! DIVERS
     95     CALL histdef(nid_tra, "pplay", "flux u mont","-",     &
     96          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     97          "inst(X)", zout,zout)
     98     CALL histdef(nid_tra, "t", "flux u mont","-",         &
     99          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     100          "inst(X)", zout,zout)
     101     CALL histdef(nid_tra, "mfu", "flux u mont","-",       &
     102          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     103          "ave(X)", zsto,zout)
     104     CALL histdef(nid_tra, "mfd", "flux u decen","-",      &
     105          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     106          "ave(X)", zsto,zout)
     107     CALL histdef(nid_tra, "en_u", "flux u mont","-",      &
     108          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     109          "ave(X)", zsto,zout)
     110     CALL histdef(nid_tra, "en_d", "flux u mont","-",      &
     111          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     112          "ave(X)", zsto,zout)
     113     CALL histdef(nid_tra, "de_d", "flux u mont","-",      &
     114          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     115          "ave(X)", zsto,zout)
     116     CALL histdef(nid_tra, "de_u", "flux u decen","-",     &
     117          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     118          "ave(X)", zsto,zout)
     119     CALL histdef(nid_tra, "coefh", "turbulent coef","-",  &
     120          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     121          "ave(X)", zsto,zout)   
     122     
     123     CALL histend(nid_tra)
     124!$OMP END MASTER
     125  END IF ! ecrit_tra>0. .AND. config_inca == 'none'
     126
     127#endif
     128 

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90

    r1188 r1191  
    11!
    2 ! $Header$
     2! $Id $
    33!
    4       SUBROUTINE  initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
    5      .                   ,vdeptr,scavtr)
    6       USE dimphy
    7       USE infotrac, ONLY : nbtr
    8       IMPLICIT none
    9 c======================================================================
    10 c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
    11 c Objet: initialisation des constantes des traceurs
    12 CAA Revison pour le controle avec la temperature du sol
    13 cAA
    14 CAA   it = 1 radon ss controle de ts
    15 cAA   it = 2 plomb ss controle de ts 
    16 c======================================================================
    17 c Arguments:
    18 c nbtr------input-I- nombre de vrais traceurs (sans l'eau)
    19 c ftsol-------input-R- Temperature du sol (Kelvin)
    20 c pctsrf-----input-R-  Nature de sol (pourcentage de sol)
    21 c masktr---output-R- Masque reservoir de sol traceur (1 = reservoir)
    22 c fshtr----output-R- Flux surfacique de production dans le sol
    23 c hsoltr---output-R- Epaisseur du reservoir de sol
    24 c tautr----output-R- Constante de decroissance du traceur
    25 c vdeptr---output-R- Vitesse de depot sec dans la couche Brownienne
    26 c scavtr---output-R- Coefficient de lessivage
    27 c======================================================================
    28 cym#include "dimensions.h"
    29 cym#include "dimphy.h"
    30 #include "indicesol.h"
    31 c======================================================================
    32 C
    33       INTEGER i, it
    34       REAL pctsrf(klon,nbsrf) !Pourcentage de sol (f(nature du sol))
    35       REAL ftsol(klon,nbsrf)  ! Temperature du sol pour le controle Rn
    36 c                             ! le cas echeant
    37       REAL masktr(klon,nbtr)  ! Masque de l'echange avec la surface
    38 c                                 (possible => 1 )
    39       REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
    40       REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
    41       REAL tautr(nbtr)       ! Constante de decroissance radioactive
    42       REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
    43       REAL scavtr(nbtr)      ! Coefficient de lessivage
    44       REAL s
    45 C
    46       WRITE(*,'(''PASSAGE initrrnpb ...'',$)')
    47       print*,'nbtr= ',nbtr
    48       print*,'nbsrf= ',nbsrf
    49       print*,'klon= ',klon
    50 C
    51 C Puis les initialisation specifiques a chaque traceur (pour le moment, Rn222)
    52 C
    53 C
    54 C Radon it = 1
    55 c
    56       IF ( nbtr .LE. 0 ) STOP 'initrrnpb pas glop pas glop'
    57       it = 1
    58       s = 1.E4  !  Source: atome par m2
    59       hsoltr(it) = 0.1      ! Hauteur equivalente du reservoir :
    60 c                              1 m * porosite 0.1
    61       tautr(it) = 4.765E5  ! Decroissance du radon, secondes
    62 cAA
    63 c      tautr(it) = 4.765E55  ! Decroissance du radon,infinie
    64 cAA
    65       vdeptr(it) = 0. ! Pas de depot sec pour le radon
    66       scavtr(it) = 0. ! Pas de lessivage pour le radon
     4SUBROUTINE  initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
     5  USE dimphy
     6  USE infotrac, ONLY : nbtr
     7  IMPLICIT NONE
     8!======================================================================
     9! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
     10! Objet: initialisation des constantes des traceurs
     11!AA Revison pour le controle avec la temperature du sol
     12!AA
     13!AA   it = 1 radon ss controle de ts
     14!AA   it = 2 plomb ss controle de ts 
     15!======================================================================
     16! Arguments:
     17! nbtr.............. nombre de vrais traceurs (sans l'eau)
     18! ftsol....input-R-  Temperature du sol (Kelvin)
     19! pctsrf...input-R-  Nature de sol (pourcentage de sol)
     20! masktr...output-R- Masque reservoir de sol traceur (1 = reservoir)
     21! fshtr....output-R- Flux surfacique de production dans le reservoir de sol
     22! hsoltr...output-R- Epaisseur equivalente du reservoir de sol
     23! tautr....output-R- Constante de decroissance radioactive du traceur
     24! vdeptr...output-R- Vitesse de depot sec dans la couche Brownienne
     25! scavtr...output-R- Coefficient de lessivage
     26!======================================================================
     27  INCLUDE "indicesol.h"
     28!======================================================================
    6729
    68       print*, '-------------- SOURCE DU RADON ------------------------ '
    69       print*,'it = ',it
    70       print*,'Source : ', s
    71       print*,'Hauteur equivalente du reservoir de sol: ',hsoltr(it)
    72       print*,'Decroissance (s): ', tautr(it)
    73       print*,'Vitesse de depot sec: ',vdeptr(it)
    74       print*,'Facteur de lessivage: ',scavtr(it)
     30  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf
     31  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol
     32  REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: masktr
     33  REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: fshtr
     34  REAL,DIMENSION(nbtr),INTENT(OUT)      :: hsoltr
     35  REAL,DIMENSION(nbtr),INTENT(OUT)      :: tautr
     36  REAL,DIMENSION(nbtr),INTENT(OUT)      :: vdeptr
     37  REAL,DIMENSION(nbtr),INTENT(OUT)      :: scavtr
     38  INTEGER                               :: i, it
     39  REAL                                  :: s
    7540
    76       DO i = 1,klon
    77         masktr(i,it) = 0.
    78         IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1.
    79         fshtr(i,it) = s * masktr(i,it)
     41  WRITE(*,*)'PASSAGE initrrnpb ...'
     42!
     43! Radon it = 1
     44!----------------
     45  IF ( nbtr .LE. 0 ) STOP '**PHYTRAC:initrrnpb:** nbtr < 0; verifier RN dans traceur.def'
     46  it = 1
     47  s = 1.E4             ! Source: atome par m2
     48  hsoltr(it) = 0.1     ! Hauteur equivalente du reservoir :
     49                       ! 1 m * porosite 0.1
     50  tautr(it) = 4.765E5  ! Decroissance du radon, secondes
     51  vdeptr(it) = 0.      ! Pas de depot sec pour le radon
     52  scavtr(it) = 0.      ! Pas de lessivage pour le radon
     53 
     54  WRITE(*,*)'-------------- SOURCE DU RADON ------------------------ '
     55  WRITE(*,*)'it = ',it
     56  WRITE(*,*)'Source : ', s
     57  WRITE(*,*)'Hauteur equivalente du reservoir de sol: ',hsoltr(it)
     58  WRITE(*,*)'Decroissance (s): ', tautr(it)
     59  WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it)
     60  WRITE(*,*)'Facteur de lessivage: ',scavtr(it)
    8061
    81 cAA
    82 cAA quelques tests
    83 cAA POur l'instant le pctsrf(i,3) = 1.0
    84 cAA lorsqu'il ya de la terre mias ne prend aucune autre valeur
    85 cAA il n'est donc pas necessaire de multiplier fshtr par pctsrf
    86 cAA
    87 c       print*, '------------------------------------------ '
    88 c        print*, 'masktr(',i,it,')= ',masktr(i,it)
    89 c        print*, 'fshtr(',i,it,')= ',fshtr(i,it)
    90 c        print*, 'pctsrf(',i,',1)= ',pctsrf(i,1)
    91 c        print*, 'pctsrf(',i,',2)= ',pctsrf(i,2)
    92 c        print*, 'pctsrf(',i,',3)= ',pctsrf(i,3)
    93 c        print*, 'pctsrf(',i,',4)= ',pctsrf(i,4)
    94 c        print*, 's = ',s
    95 c        print*, '------------------------------------------ '
     62  DO i = 1,klon
     63     masktr(i,it) = 0.
     64     IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1.
     65     fshtr(i,it) = s * masktr(i,it)
     66  END DO
     67!
     68! 210Pb it = 2
     69!----------------
     70  IF ( nbtr .LE. 1 ) STOP '**PHYTRAC**:initrrnpb:** nbtr <= 1; verifier PB dans traceur.def'
     71  it = 2
     72  s = 0.                ! Pas de source
     73  hsoltr(it) = 10.      ! Hauteur equivalente du reservoir
     74                        ! a partir duquel le depot Brownien a lieu
     75  tautr(it) = 1.028E9   ! Decroissance du Pb210, secondes
     76  vdeptr(it) = 1.E-3    ! 1 mm/s pour le 210Pb
     77  scavtr(it) =  .5      ! Lessivage du Pb210
     78  DO i = 1,klon
     79     masktr(i,it) = 1.  ! Le depot sec peut avoir lieu partout
     80     fshtr(i,it) = s * masktr(i,it)
     81  END DO
     82  WRITE(*,*)'-------------- SOURCE DU PLOMB ------------------------ '
     83  WRITE(*,*)'it = ',it
     84  WRITE(*,*)'Source : ', s
     85  WRITE(*,*)'Hauteur equivalente du reservoir : ',hsoltr(it)
     86  WRITE(*,*)'Decroissance (s): ', tautr(it)
     87  WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it)
     88  WRITE(*,*)'Facteur de lessivage: ',scavtr(it)
    9689
    97       END DO
    98 C
    99 C 210Pb it = 2
    100 C
    101       IF ( nbtr .LE. 1 ) STOP 'initrrnpb pas glop pas glop'
    102       it = 2
    103       s = 0. !  Pas de source !!!
    104       hsoltr(it) = 10.     ! Hauteur equivalente du reservoir
    105 c                              a partir duquel le
    106 c                              depot Brownien a lieu
    107       tautr(it) = 1.028E9 ! Decroissance du Pb210, secondes
    108       vdeptr(it) = 1.E-3 ! 1 mm/s pour le 210Pb
    109       scavtr(it) =  .5   ! Lessivage du Pb210
    110       DO i = 1,klon
    111         masktr(i,it) = 1. ! Le depot sec peut avoir lieu partout
    112         fshtr(i,it) = s * masktr(i,it)
    113       END DO
    114       print*, '-------------- SOURCE DU PLOMB ------------------------ '
    115       print*,'it = ',it
    116       print*,'Source : ', s
    117       print*,'Hauteur equivalente du reservoir : ',hsoltr(it)
    118       print*,'Decroissance (s): ', tautr(it)
    119       print*,'Vitesse de depot sec: ',vdeptr(it)
    120       print*,'Facteur de lessivage: ',scavtr(it)
    121 c
    122       WRITE(*,*) 'initialisation rnpb ok'
    123 c
    124       RETURN
    125       END
     90  WRITE(*,*) 'Initialisation RN et PB ok'
     91
     92END SUBROUTINE initrrnpb

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90

    r1188 r1191  
    11!
    2 ! $Header$
     2! $Id$
    33!
    4       SUBROUTINE minmaxqfi(zq,qmin,qmax,comment)
    5       USE dimphy
    6       IMPLICIT none
     4SUBROUTINE minmaxqfi(zq,qmin,qmax,comment)
     5  USE dimphy
     6  IMPLICIT NONE
    77
    8 cym#include "dimensions.h"
    9 cym#include "dimphy.h"
     8! Entrees
     9  REAL,DIMENSION(klon,klev), INTENT(IN)   :: zq
     10  REAL,INTENT(IN)                         :: qmin,qmax
     11  CHARACTER(LEN=*),INTENT(IN)             :: comment
    1012
    11       CHARACTER*(*) comment
    12       real qmin,qmax
    13       real zq(klon,klev)
    14 
    15       INTEGER jadrs(klon), jbad, k, i
    16 
    17       DO k = 1, klev
    18          jbad = 0
    19          DO i = 1, klon
    20          IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
    21             jbad = jbad + 1
    22             jadrs(jbad) = i
    23          ENDIF
    24          ENDDO
    25          IF (jbad.GT.0) THEN
    26          PRINT*, comment
    27          DO i = 1, jbad
    28             PRINT*, "i,k,q=", jadrs(i),k,zq(jadrs(i),k)
    29          ENDDO
    30          ENDIF
    31       ENDDO
    32 
    33       return
    34       end
     13! Local 
     14  INTEGER,DIMENSION(klon)     :: jadrs
     15  INTEGER                     :: i, jbad, k
     16 
     17  DO k = 1, klev
     18     jbad = 0
     19     DO i = 1, klon
     20        IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
     21           jbad = jbad + 1
     22           jadrs(jbad) = i
     23        ENDIF
     24     ENDDO
     25     IF (jbad.GT.0) THEN
     26        WRITE(*,*)comment
     27        DO i = 1, jbad
     28           WRITE(*,*) "i,k,q=", jadrs(i),k,zq(jadrs(i),k)
     29        ENDDO
     30     ENDIF
     31  ENDDO
     32 
     33END SUBROUTINE minmaxqfi

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90

    r1188 r1191  
    11!
    2 ! $Header$
     2! $Id $
    33!
    4       SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d,
    5      .                 pplay,paprs,x,dx)
    6       USE dimphy
    7       IMPLICIT NONE
    8 c=====================================================================
    9 c Objet : Melange convectif de traceurs a partir des flux de masse
    10 c Date : 13/12/1996 -- 13/01/97
    11 c Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD),
    12 c         Brinkop et Sausen (1996) et Boucher et al. (1996).
    13 c ATTENTION : meme si cette routine se veut la plus generale possible,
    14 c             elle a herite de certaines notations et conventions du
    15 c             schema de Tiedtke (1993).
    16 c --En particulier, les couches sont numerotees de haut en bas !!!
    17 c   Ceci est valable pour les flux
    18 c   mais pas pour les entrees x, pplay, paprs !!!!
    19 c --pmfu est positif, pmfd est negatif
    20 c --Tous les flux d'entrainements et de detrainements sont positifs
    21 c   contrairement au schema de Tiedtke d'ou les changements de signe!!!!
    22 c=====================================================================
    23 c
    24 cym#include "dimensions.h"
    25 cym#include "dimphy.h"
    26 #include "YOMCST.h"
    27 #include "YOECUMF.h"
    28 c
    29       REAL pdtime
    30 c--les flux sont definis au 1/2 niveaux
    31 c--pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls
    32       REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
    33       REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
    34       REAL pen_u(klon,klev) ! flux entraine dans le panache montant
    35       REAL pde_u(klon,klev) ! flux detraine dans le panache montant
    36       REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
    37       REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
     4SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d,pplay,paprs,x,dx)
     5  USE dimphy
     6  IMPLICIT NONE
     7!=====================================================================
     8! Objet : Melange convectif de traceurs a partir des flux de masse
     9! Date : 13/12/1996 -- 13/01/97
     10! Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD),
     11!         Brinkop et Sausen (1996) et Boucher et al. (1996).
     12! ATTENTION : meme si cette routine se veut la plus generale possible,
     13!             elle a herite de certaines notations et conventions du
     14!             schema de Tiedtke (1993).
     15! 1. En particulier, les couches sont numerotees de haut en bas !!!
     16!    Ceci est valable pour les flux
     17!    mais pas pour les entrees x, pplay, paprs !!!!
     18! 2. pmfu est positif, pmfd est negatif
     19! 3. Tous les flux d'entrainements et de detrainements sont positifs
     20!    contrairement au schema de Tiedtke d'ou les changements de signe!!!!
     21!=====================================================================
     22!
     23  include "YOMCST.h"
     24  include "YOECUMF.h"
    3825
    39       REAL pplay(klon,klev)    ! pression aux couches (bas en haut)
    40       REAL paprs(klon,klev+1)  ! pression aux 1/2 couches (bas en haut)
    41       REAL x(klon,klev)        ! q de traceur (bas en haut)
    42       REAL dx(klon,klev)     ! tendance de traceur  (bas en haut)
    43 c
    44 c--flux convectifs mais en variables locales
    45       REAL zmfu(klon,klev+1)
    46       REAL zmfd(klon,klev+1)
    47       REAL zen_u(klon,klev)
    48       REAL zde_u(klon,klev)
    49       REAL zen_d(klon,klev)
    50       REAL zde_d(klon,klev)
    51       real zmfe
    52 c
    53 c--variables locales     
    54 c--les flux de x sont definis aux 1/2 niveaux
    55 c--xu et xd sont definis aux niveaux complets
    56       REAL xu(klon,klev)        ! q de traceurs dans le panache montant
    57       REAL xd(klon,klev)        ! q de traceurs dans le panache descendant
    58       REAL zmfux(klon,klev+1)   ! flux de x dans le panache montant
    59       REAL zmfdx(klon,klev+1)   ! flux de x dans le panache descendant
    60       REAL zmfex(klon,klev+1)   ! flux de x dans l'environnement
    61       INTEGER i, k
    62       REAL zmfmin
    63       PARAMETER (zmfmin=1.E-10)
     26  REAL,INTENT(IN) :: pdtime  ! pdtphys
     27!
     28! les flux sont definis au 1/2 niveaux
     29! => pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls
     30!
     31  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
     32  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
     33  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
     34  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
     35  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
     36  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
    6437
    65 c =========================================
    66 c
    67 c
    68 c   Extension des flux UP et DN sur klev+1 niveaux
    69 c =========================================
    70       do k=1,klev
    71          do i=1,klon
    72             zmfu(i,k)=pmfu(i,k)
    73             zmfd(i,k)=pmfd(i,k)
    74          enddo
    75       enddo
    76       do i=1,klon
    77          zmfu(i,klev+1)=0.
    78          zmfd(i,klev+1)=0.
    79       enddo
     38  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay ! pression aux couches (bas en haut)
     39  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
     40  REAL,DIMENSION(klon,klev),INTENT(IN)   :: x     ! q de traceur (bas en haut)
     41  REAL,DIMENSION(klon,klev),INTENT(INOUT) :: dx   ! tendance de traceur  (bas en haut)
    8042
    81 c--modif pour diagnostiquer les detrainements
    82 c =========================================
    83 c   on privilegie l'ajustement de l'entrainement dans l'ascendance.
     43! flux convectifs mais en variables locales
     44  REAL,DIMENSION(klon,klev+1) :: zmfu  ! copie de pmfu avec klev+1 = 0
     45  REAL,DIMENSION(klon,klev+1) :: zmfd  ! copie de pmfd avec klev+1 = 0
     46  REAL,DIMENSION(klon,klev)   :: zen_u
     47  REAL,DIMENSION(klon,klev)   :: zde_u
     48  REAL,DIMENSION(klon,klev)   :: zen_d
     49  REAL,DIMENSION(klon,klev)   :: zde_d
     50  REAL                        :: zmfe
    8451
    85       do k=1, klev
    86          do i=1, klon
    87             zen_d(i,k)=pen_d(i,k)
    88             zde_u(i,k)=pde_u(i,k)
    89             zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k)
    90             zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k)
    91          enddo
    92       enddo
    93 c
    94 c--calcul des flux dans le panache montant
    95 c =========================================
    96 c
    97 c Dans la premiere couche, on prend q comme valeur de qu
    98 c
    99       do i=1, klon
    100          zmfux(i,1)=0.0
    101       enddo
    102 c
    103 c Autres couches
    104       do k=1,klev
    105          do i=1, klon
    106             if ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN
    107                xu(i,k)=x(i,k)
    108             else
    109                xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k))
    110      s               /(zmfu(i,k+1)+zde_u(i,k))
    111             endif
    112             zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k)
    113          enddo
    114       enddo
    115 c
    116 c--calcul des flux dans le panache descendant
    117 c =========================================
    118 c   
    119       do i=1, klon
    120          zmfdx(i,klev+1)=0.0
    121       enddo
    122 c
    123       do k=klev,1,-1
    124          do i=1, klon
    125             if ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN
    126                xd(i,k)=x(i,k)
    127             else
    128                xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k)) /
    129      .               (zmfd(i,k)-zde_d(i,k))
    130             endif
    131             zmfdx(i,k)=zmfd(i,k)*xd(i,k)
    132          enddo
    133       enddo
    134 c
    135 c--introduction du flux de retour dans l'environnement
    136 c =========================================
    137 c
    138       do k=2, klev
    139          do i=1, klon
    140             zmfe=-zmfu(i,k)-zmfd(i,k)
    141             if (zmfe.le.0.) then
    142                zmfex(i,k)= zmfe*x(i,k)
    143             else
    144                zmfex(i,k)= zmfe*x(i,k-1)
    145             endif
    146          enddo
    147       enddo
     52! variables locales     
     53! les flux de x sont definis aux 1/2 niveaux
     54! xu et xd sont definis aux niveaux complets
     55  REAL,DIMENSION(klon,klev)   :: xu      ! q de traceurs dans le panache montant
     56  REAL,DIMENSION(klon,klev)   :: xd      ! q de traceurs dans le panache descendant
     57  REAL,DIMENSION(klon,klev+1) :: zmfux   ! flux de x dans le panache montant
     58  REAL,DIMENSION(klon,klev+1) :: zmfdx   ! flux de x dans le panache descendant
     59  REAL,DIMENSION(klon,klev+1) :: zmfex   ! flux de x dans l'environnement
     60  INTEGER                     :: i, k
     61  REAL,PARAMETER              :: zmfmin=1.E-10
    14862
    149       do i=1, klon
    150          zmfex(i,1)=0.
    151          zmfex(i,klev+1)=0.
    152       enddo
    153 c
    154 c--calcul final des tendances
    155 c
    156       do k=1, klev
    157          do i=1, klon
    158             dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime*
    159      .                      ( zmfux(i,k) - zmfux(i,k+1) +
    160      .                        zmfdx(i,k) - zmfdx(i,k+1) +
    161      .                        zmfex(i,k) - zmfex(i,k+1) )
    162          enddo
    163       enddo
    164 c
    165       return
    166       end
     63! ==============================================
     64! Extension des flux UP et DN sur klev+1 niveaux
     65! ==============================================
     66  DO k=1,klev
     67     DO i=1,klon
     68        zmfu(i,k)=pmfu(i,k)
     69        zmfd(i,k)=pmfd(i,k)
     70     ENDDO
     71  ENDDO
     72  DO i=1,klon
     73     zmfu(i,klev+1)=0.
     74     zmfd(i,klev+1)=0.
     75  ENDDO
     76! ==========================================
     77! modif pour diagnostiquer les detrainements
     78! ==========================================
     79!   on privilegie l'ajustement de l'entrainement dans l'ascendance.
     80
     81  DO k=1, klev
     82     DO i=1, klon
     83        zen_d(i,k)=pen_d(i,k)
     84        zde_u(i,k)=pde_u(i,k)
     85        zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k)
     86        zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k)
     87     ENDDO
     88  ENDDO
     89! =========================================
     90! calcul des flux dans le panache montant
     91! =========================================
     92!
     93! Dans la premiere couche, on prend q comme valeur de qu
     94
     95  DO i=1, klon
     96     zmfux(i,1)=0.0
     97  ENDDO
     98
     99! Autres couches
     100  DO k=1,klev
     101     DO i=1, klon
     102        IF ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN
     103           xu(i,k)=x(i,k)
     104        ELSE
     105           xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k))/(zmfu(i,k+1)+zde_u(i,k))
     106        ENDIF
     107        zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k)
     108     ENDDO
     109  ENDDO
     110! ==========================================
     111! calcul des flux dans le panache descendant
     112! ==========================================
     113   
     114  DO i=1, klon
     115     zmfdx(i,klev+1)=0.0
     116  ENDDO
     117
     118  DO k=klev,1,-1
     119     DO i=1, klon
     120        IF ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN
     121           xd(i,k)=x(i,k)
     122        ELSE
     123           xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k))/(zmfd(i,k)-zde_d(i,k))
     124        ENDIF
     125        zmfdx(i,k)=zmfd(i,k)*xd(i,k)
     126     ENDDO
     127  ENDDO
     128! ===================================================
     129! introduction du flux de retour dans l'environnement
     130! ===================================================
     131
     132  DO k=2, klev
     133     DO i=1, klon
     134        zmfe=-zmfu(i,k)-zmfd(i,k)
     135        IF (zmfe.le.0.) then
     136           zmfex(i,k)= zmfe*x(i,k)
     137        ELSE
     138           zmfex(i,k)= zmfe*x(i,k-1)
     139        ENDIF
     140     ENDDO
     141  ENDDO
     142
     143  DO i=1, klon
     144     zmfex(i,1)=0.
     145     zmfex(i,klev+1)=0.
     146  ENDDO
     147! ==========================
     148! calcul final des tendances
     149! ==========================
     150  DO k=1, klev
     151     DO i=1, klon
     152        dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime*  &
     153             ( zmfux(i,k) - zmfux(i,k+1) +            &
     154             zmfdx(i,k) - zmfdx(i,k+1) +              &
     155             zmfex(i,k) - zmfex(i,k+1) )
     156     ENDDO
     157  ENDDO
     158 
     159END SUBROUTINE nflxtr

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F

    r1054 r1191  
    1919      USE iostart
    2020      USE write_field_phy
     21      USE infotrac
     22      USE traclmdz_mod,    ONLY : traclmdz_from_restart
     23
    2124      IMPLICIT none
    2225c======================================================================
    … …  
    4851      REAL run_off_lic_0(klon)
    4952      REAL fractint(klon)
     53      REAL trs(klon,nbtr)
    5054
    5155      CHARACTER*6 ocean_in
    … …  
    6266      INTEGER length
    6367      PARAMETER (length=100)
     68      INTEGER it, iiq
    6469      REAL tab_cntrl(length), tabcntr0(length)
    6570      CHARACTER*7 str7
    … …  
    983988      PRINT*,'(ecart-type) wake_fip:', xmin, xmax
    984989c
     990c Read and send field trs to traclmdz
     991c
     992      IF (type_trac == 'lmdz') THEN
     993         DO it=1,nbtr
     994            iiq=niadv(it+2)
     995            CALL get_field("trs_"//tname(iiq),trs(:,it),found)
     996            IF (.NOT. found) THEN
     997               PRINT*,
     998     $           "phyetat0: Le champ <trs_"//tname(iiq)//"> est absent"
     999               PRINT*, "Depart legerement fausse. Mais je continue"
     1000               trs(:,it) = 0.
     1001            ENDIF
     1002            xmin = 1.0E+20
     1003            xmax = -1.0E+20
     1004            xmin = MINval(trs(:,it))
     1005            xmax = MAXval(trs(:,it))
     1006            PRINT*,"(ecart-type) trs_"//tname(iiq)//" :", xmin, xmax
     1007
     1008         END DO
     1009         
     1010         CALL traclmdz_from_restart(trs)
     1011      END IF
     1012
    9851013
    9861014c on ferme le fichier
    … …  
    10051033      CALL fonte_neige_init(run_off_lic_0)
    10061034
     1035
    10071036      RETURN
    10081037      END

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F

    r1001 r1191  
    1212      USE phys_state_var_mod
    1313      USE iostart
     14      USE traclmdz_mod, ONLY : traclmdz_to_restart
     15      USE infotrac
    1416
    1517      IMPLICIT none
    … …  
    4244      REAL agesno(klon,nbsrf)
    4345      REAL run_off_lic_0(klon)
     46      REAL trs(klon,nbtr)
    4447c
    4548      INTEGER nid, nvarid, idim1, idim2, idim3
    … …  
    5255      CHARACTER (len=7) :: str7
    5356      CHARACTER (len=2) :: str2
    54 
     57      INTEGER           :: it, iiq
     58     
    5559c======================================================================
    5660c
    … …  
    311315      CALL put_field("WAKE_FIP","",wake_fip)
    312316
     317
     318! trs from traclmdz_mod
     319      IF (type_trac == 'lmdz') THEN
     320         CALL traclmdz_to_restart(trs)
     321         DO it=1,nbtr
     322            iiq=niadv(it+2)
     323            CALL put_field("trs_"//tname(iiq),"",trs(:,it))
     324         END DO
     325      END IF
     326
    313327      CALL close_restartphy
    314328!$OMP BARRIER

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

    r1188 r1191  
    111111      PARAMETER (ok_stratus=.FALSE.)
    112112c======================================================================
    113       LOGICAL, SAVE :: rnpb=.TRUE.
    114 c$OMP THREADPRIVATE(rnpb)
    115113      REAL amn, amx
    116114      INTEGER igout
    … …  
    31333131c====================================================================
    31343132C
    3135       IF (config_inca /= 'none') rnpb=.FALSE.
    3136 
    3137       call phytrac (     rnpb,
    3138      I                   itap,
    3139      I                   julien,
    3140      I                   gmtime,
    3141      I                   debut,
    3142      I                   lafin,
    3143      I                   nlon,
    3144      I                   klev,
    3145      I                   dtime,
    3146      I                   u,
    3147      I                   v,
    3148      I                   t,
    3149      I                   paprs,
    3150      I                   pplay,
    3151      I                   pmfu,
    3152      I                   pmfd,
    3153      I                   pen_u,
    3154      I                   pde_u,
    3155      I                   pen_d,
    3156      I                   pde_d,
    3157      I                   cdragh,
    3158      I                   coefh,
    3159      I                   fm_therm,
    3160      I                   entr_therm,
    3161      I                   u1,
    3162      I                   v1,
    3163      I                   ftsol,
    3164      I                   pctsrf,
    3165      I                   rlat,
    3166      I                   frac_impa,
    3167      I                   frac_nucl,
    3168      I                   rlon,
    3169      I                   presnivs,
    3170      I                   pphis,
    3171      I                   pphi,
    3172      I                   albsol1,
    3173      I                   qx(1,1,1),
    3174      I                   rhcl,
    3175      I                   cldfra,
    3176      I                   rneb,
    3177      I                   diafra,
    3178      I                   cldliq,
    3179      I                   itop_con,
    3180      I                   ibas_con,
    3181      I                   pmflxr,
    3182      I                   pmflxs,
    3183      I                   prfl,
    3184      I                   psfl,
    3185      I                   da,
    3186      I                   phi,
    3187      I                   mp,
    3188      I                   upwd,
    3189      I                   dnwd,
    3190      I                   aerosol_couple,
    3191      I                   flxmass_w,
    3192      I                   tau_aero,
    3193      I                   piz_aero,
    3194      I                   cg_aero,
    3195      I                   ccm,
    3196      I                   rfname,
    3197      O                   tr_seri)
     3133
     3134      call phytrac (
     3135     I     itap,     julien,    gmtime,   debut,
     3136     I     lafin,    dtime,     u, v,     t,
     3137     I     paprs,    pplay,     pmfu,     pmfd,
     3138     I     pen_u,    pde_u,     pen_d,    pde_d,
     3139     I     cdragh,   coefh,     fm_therm, entr_therm,
     3140     I     u1,       v1,        ftsol,    pctsrf,
     3141     I     rlat,     frac_impa, frac_nucl,rlon,
     3142     I     presnivs, pphis,     pphi,     albsol1,
     3143     I     qx(:,:,ivap),rhcl,   cldfra,   rneb,
     3144     I     diafra,   cldliq,    itop_con, ibas_con,
     3145     I     pmflxr,   pmflxs,    prfl,     psfl,
     3146     I     da,       phi,       mp,       upwd,     
     3147     I     dnwd,     aerosol_couple,      flxmass_w,
     3148     I     tau_aero, piz_aero,  cg_aero,  ccm,
     3149     I     rfname,
     3150     O     tr_seri)
    31983151
    31993152      IF (offline) THEN

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90

    r1188 r1191  
    1 !
    2 c
    3 c
    4       SUBROUTINE phytrac (rnpb,
    5      I                    nstep,
    6      I                    julien,
    7      I                    gmtime,
    8      I                    debutphy,
    9      I                    lafin,
    10      I                    nlon,
    11      I                    nlev,
    12      I                    pdtphys,
    13      I                    u,
    14      I                    v,
    15      I                    t_seri,
    16      I                    paprs,
    17      I                    pplay,
    18      I                    pmfu,
    19      I                    pmfd,
    20      I                    pen_u,
    21      I                    pde_u,
    22      I                    pen_d,
    23      I                    pde_d,
    24      I                    cdragh,
    25      I                    coefh,
    26      I                    fm_therm,
    27      I                    entr_therm,
    28      I                    yu1,
    29      I                    yv1,
    30      I                    ftsol,
    31      I                    pctsrf,
    32      I                    xlat,
    33      I                    frac_impa,
    34      I                    frac_nucl,
    35      I                    xlon,
    36      I                    presnivs,
    37      I                    pphis,
    38      I                    pphi,
    39      I                    albsol,
    40      I                    sh,
    41      I                    rh,
    42      I                    cldfra,
    43      I                    rneb,
    44      I                    diafra,
    45      I                    cldliq,
    46      I                    itop_con,
    47      I                    ibas_con,
    48      I                    pmflxr,
    49      I                    pmflxs,
    50      I                    prfl,
    51      I                    psfl,
    52      I                    da,
    53      I                    phi,
    54      I                    mp,
    55      I                    upwd,
    56      I                    dnwd,
    57      I                    aerosol_couple,
    58      I                    flxmass_w,
    59      I                    tau_aero,
    60      I                    piz_aero,
    61      I                    cg_aero,
    62      I                    ccm,
    63      I                    rfname,
    64      O                    tr_seri)
    65 
    66       USE ioipsl
    67       USE dimphy
    68       USE infotrac
    69       USE mod_grid_phy_lmdz
    70       USE mod_phys_lmdz_para
    71       USE comgeomphy
    72       USE iophy
    73       USE vampir
    74 
    75       IMPLICIT none
    76 c======================================================================
    77 c Auteur(s) FH
    78 c Objet: Moniteur general des tendances traceurs
    79 c
    80 cAA Remarques en vrac:
    81 cAA--------------------
    82 cAA 2/ Le choix du radon et du pb se fait juste avec un data
    83 cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
    84 cAA    une variable "type de traceur"
    85 c======================================================================
    86 #include "YOMCST.h"
    87 #include "dimensions.h"
    88 #include "indicesol.h"
    89 #include "clesphys.h"
    90 #include "temps.h"
    91 #include "paramet.h"
    92 #include "control.h"
    93 #include "thermcell.h"
    94 c======================================================================
    95 
    96 c Arguments:
    97 c
    98 c   EN ENTREE:
    99 c   ==========
    100 c
    101 c   divers:
    102 c   -------
    103 c
    104       integer nlon  ! nombre de points horizontaux
    105       integer nlev  ! nombre de couches verticales
    106       integer nstep  ! appel physique
    107       integer julien !jour julien
    108       integer itop_con(nlon)
    109       integer ibas_con(nlon)
    110       real gmtime
    111       real pdtphys  ! pas d'integration pour la physique (seconde)
    112       real t_seri(nlon,nlev) ! temperature
    113       real tr_seri(nlon,nlev,nbtr) ! traceur 
    114       real u(klon,klev)
    115       real v(klon,klev)
    116       real sh(nlon,nlev)     ! humidite specifique
    117       real rh(nlon,nlev)     ! humidite relative
    118       real cldliq(nlon,nlev) ! eau liquide nuageuse
    119       real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
    120       real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
    121       real rneb(nlon,nlev)   ! fraction nuageuse (grande echelle)
    122       real albsol(nlon)  ! albedo surface
    123       real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
    124       real ps(nlon)  ! pression surface
    125       real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
    126       real pphi(nlon,klev) ! geopotentiel
    127       real pphis(klon)
    128       REAL presnivs(klev)
    129       logical debutphy       ! le flag de l'initialisation de la physique
    130       logical lafin          ! le flag de la fin de la physique
    131 c Olivia     
    132       integer nsplit
    133       REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
    134       REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale
    135       LOGICAL aerosol_couple
    136 
    137       REAL flxmass_w(klon,klev)
    138       CHARACTER(len=8) :: solsym(nbtr)
    139       integer la
    140       REAL              ::    tau_aero(klon,klev,9,2)
    141       REAL              ::    piz_aero(klon,klev,9,2)
    142       REAL              ::    cg_aero(klon,klev,9,2)
    143       character*4       ::    rfname(9)
    144       REAL              ::    ccm(klon,klev,2)
    145 c
    146 c   convection:
    147 c   -----------
    148 c
    149       REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
    150       REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
    151       REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
    152 
    153 c
    154 c   thermiques:
    155 c   -----------
    156 c
    157       real fm_therm(klon,klev+1),entr_therm(klon,klev)
    158         real fm_therm1(klon,klev)
    159 c
    160       REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
    161       REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
    162       REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
    163 c KE
    164       real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
    165       REAL upwd(nlon,nlev)      ! saturated updraft mass flux
    166       REAL dnwd(nlon,nlev)      ! saturated downdraft mass flux
    167 
    168 c
    169 c   Couche limite:
    170 c   --------------
    171 c
    172       REAL cdragh(nlon,nlev)! coeff drag pour T et Q
    173       REAL coefh(nlon,nlev) ! coeff melange CL
    174       REAL yu1(nlon)        ! vents au premier niveau
    175       REAL yv1(nlon)        ! vents au premier niveau
    176       REAL xlat(nlon)       ! latitudes pour chaque point
    177       REAL xlon(nlon)       ! longitudes pour chaque point
    178 
    179 c
    180 c   Lessivage:
    181 c   ----------
    182 c
    183 c pour le ON-LINE
    184 c
    185       REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
    186       REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
    187 c
    188 cAA
    189 cAA Arguments necessaires pour les sources et puits de traceur:
    190 cAA ----------------
    191 cAA
    192       real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
    193       real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
    194 c abder
    195       real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
    196       real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
    197 c fin
    198 cAA ----------------------------
    199 cAA  VARIABLES LOCALES TRACEURS
    200 cAA ----------------------------
    201 cAA
    202 cAA Sources et puits des traceurs:
    203 cAA ------------------------------
    204 cAA
    205 cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
    206 
    207       REAL source(klon,nbtr)       ! a voir lorsque le flux est prescrit
    208 cAA
    209 cAA Pour la source de radon et son reservoir de sol
    210 cAA ................................................
    211  
    212       REAL,save,allocatable :: trs(:,:)    ! Conc. radon ds le sol
    213 c$OMP THREADPRIVATE(trs)
    214       REAL :: trs_tmp(klon_glo)
    215      
    216       REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur
    217 c                            Masque de l'echange avec la surface
    218 c                           (1 = reservoir) ou (possible => 1 )
    219 c$OMP THREADPRIVATE(masktr)
    220       REAL,save,allocatable :: fshtr(:,:)  ! Flux surfacique dans le reservoir de sol
    221 c$OMP THREADPRIVATE(fshtr)
    222       REAL,save,allocatable :: hsoltr(:)      ! Epaisseur equivalente du reservoir de sol
    223 c$OMP THREADPRIVATE(hsoltr)
    224       REAL,save,allocatable :: tautr(:)       ! Constante de decroissance radioactive
    225 c$OMP THREADPRIVATE(tautr)
    226       REAL,save,allocatable :: vdeptr(:)      ! Vitesse de depot sec dans la couche Brownienne
    227 c$OMP THREADPRIVATE(vdeptr)
    228       REAL,save,allocatable :: scavtr(:)      ! Coefficient de lessivage
    229 c$OMP THREADPRIVATE(scavtr)
    230 cAA
    231       CHARACTER*2 itn
    232 C maf ioipsl
    233       CHARACTER*2 str2
    234       INTEGER nhori, nvert
    235       REAL zsto, zout, zjulian
    236       INTEGER, SAVE :: nid_tra
    237 c$OMP THREADPRIVATE(nid_tra)
    238       INTEGER, SAVE :: nid_tra2,nid_tra3
    239 c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
    240       INTEGER ndex(1)
    241       INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
    242       REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
    243       REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
    244       REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
    245       REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
    246 c
    247       integer itau_w   ! pas de temps ecriture = nstep + itau_phy
    248 c
    249       logical ok_sync
    250       parameter (ok_sync = .true.)
    251 C
    252 C nature du traceur
    253 c
    254       logical,save,allocatable :: aerosol(:)  ! Nature du traceur
    255 c                            ! aerosol(it) = true  => aerosol
    256 c                            ! aerosol(it) = false => gaz
    257 c                            ! nat_trac(it) = 1. aerosol
    258       logical,save,allocatable :: clsol(:)    ! clsol(it) = true => CL sol calculee
    259       logical,save,allocatable :: radio(:)    ! radio(it)=true => decroisssance radioactive
    260 c$OMP THREADPRIVATE(aerosol,clsol,radio) 
    261 C
    262 c======================================================================
    263 c
    264 c Declaration des procedures appelees
    265 c
    266 c--modif convection tiedtke
    267       INTEGER i, k, it
    268       INTEGER iq, iiq
    269       REAL delp(klon,klev)
    270 c--end modif
    271 c
    272 c Variables liees a l'ecriture de la bande histoire physique
    273 c
    274 c Variables locales pour effectuer les appels en serie
    275 c----------------------------------------------------
    276 c
    277       REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
    278       REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs  couche limite
    279       REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs  conv pour chq traceur
    280       REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
    281       REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
    282 c                                   ! radioactive du rn - > pb
    283       REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
    284 c                                          ! par impaction
    285       REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
    286 c                                          ! par nucleation
    287       REAL fluxrn(klon,klev)
    288       REAL fluxpb(klon,klev)
    289       REAL pbimpa(klon,klev)
    290       REAL pbnucl(klon,klev)
    291       REAL rn(klon,klev)
    292       REAL pb(klon,klev)
    293       REAL flestottr(klon,klev,nbtr) ! flux de lessivage
    294 c                                    ! dans chaque couche
    295       real zmasse(klon,klev)
    296       real ztra_th(klon,klev)
    297 
    298 C
    299       character*20 modname
    300       character*80 abort_message
    301 c
    302 c   Controles
    303 c-------------
    304       logical first,couchelimite,convection,lessivage,sorties,
    305      s        rnpb,inirnpb
    306       save first,couchelimite,convection,lessivage,
    307      s        sorties,inirnpb
    308 c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
    309 c$OMP+              sorties,inirnpb)
    310 c      data first,couchelimite,convection,lessivage,sorties
    311 c     s     /.true.,.true.,.false.,.true.,.true./
    312 c Olivia
    313        data first,couchelimite,convection,lessivage,
    314      s      sorties
    315      s     /.true.,.true.,.true.,.true.,.true./
    316 
    317 ! Variables needed for configuration with INCA
    318       INTEGER           :: lastgas
    319       INTEGER           :: ncsec
    320 
    321       REAL, PARAMETER   :: dry_mass = 28.966
    322       REAL, POINTER     :: hbuf(:)
    323       REAL, ALLOCATABLE :: obuf(:)
    324       REAL              :: calday
    325       REAL              :: pdel(klon,klev)
    326 c
    327 c======================================================================
    328 
    329          modname='phytrac'
    330 
    331          ps(:)=paprs(:,1)
    332 
    333          if (debutphy) then
    334            allocate( trs(klon,nbtr) )
    335            allocate( masktr(klon,nbtr))
    336            allocate( fshtr(klon,nbtr) )
    337            allocate( hsoltr(nbtr))
    338            allocate( tautr(nbtr))
    339            allocate( vdeptr(nbtr))
    340            allocate( scavtr(nbtr))
    341            allocate( aerosol(nbtr))
    342            allocate( clsol(nbtr))
    343            allocate( radio(nbtr))
    344 
    345 
    346 ! FH 2008/05/09 correction de la frequence d'ecriture des traceurs
    347 !         ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy))
    348           print*,'dans phytrac ',pdtphys,ecrit_tra
    349 
    350          inirnpb=rnpb
    351          PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
    352 C         
    353 c=============================================================
    354 c   Initialisation des sorties
    355 c=============================================================
    356 
    357 #ifdef CPP_IOIPSL
    358 #include "ini_histrac.h"
    359 #endif
    360 
    361 c======================================================================
    362 c   Initialisation de certaines variables pour le Rn et le Pb
    363 c======================================================================
    364 
    365 c Initialisation du traceur dans le sol (couche limite radonique)
    366 c
    367 c        print*,'valeur de debut dans phytrac :',debutphy
    368          trs(:,:) = 0.
    369 c$OMP MASTER         
    370        if (is_mpi_root) then
    371          trs_tmp(:)=0.
    372          open (99,file='starttrac',status='old',
    373      .         err=999,form='formatted')
    374          read(99,*) (trs_tmp(i),i=1,klon_glo)
    375 999      close(99)
    376        endif
    377 c$OMP END MASTER
    378        call Scatter(trs_tmp,trs(:,1))
    379 
    380 c         print*, 'apres starttrac'
    381 
    382 c Initialisation de la fraction d'aerosols lessivee
    383 c
    384          d_tr_lessi_impa(:,:,:) = 0.
    385          d_tr_lessi_nucl(:,:,:) = 0.
    386 c
    387 c Initialisation de la nature des traceurs
    388 c
    389          DO it = 1, nbtr
    390             aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
    391             radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
    392             clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
    393          ENDDO
    394 c
    395       ENDIF  ! fin debutphy
    396 c Initialisation du traceur dans le sol (couche limite radonique)
    397       if(inirnpb) THEN
    398 c
    399          radio(1)= .true.
    400          radio(2)= .true.
    401          clsol(1)= .true.
    402          clsol(2)= .true.
    403          aerosol(2) = .TRUE. ! le Pb est un aerosol
    404 c
    405          call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
    406      .                   ,vdeptr,scavtr)
    407          inirnpb=.false.
    408       endif
    409 
    410 
    411       IF (config_inca == 'none') THEN
    412        DO i=1,nlon
    413           pftsol1(i) = ftsol(i,1)
    414           pftsol2(i) = ftsol(i,2)
    415           pftsol3(i) = ftsol(i,3)
    416           pftsol4(i) = ftsol(i,4)
    417 
    418           ppsrf1(i) = pctsrf(i,1)
    419           ppsrf2(i) = pctsrf(i,2)
    420           ppsrf3(i) = pctsrf(i,3)
    421           ppsrf4(i) = pctsrf(i,4)
    422 
    423        ENDDO
    424 
    425       ELSE ! config_inca /=none
    426 #ifdef INCA
    427       call VTe(VTphysiq)
    428       call VTb(VTinca)
    429 
    430 !======================================================================
    431 !     Chimie
    432 !======================================================================
    433 
    434         calday = FLOAT(julien) + gmtime
    435         ncsec  = NINT (86400.*gmtime)
    436 
    437         DO k = 1, nlev
    438         pdel(:,k) = paprs(:,k) - paprs (:,k+1)
     1!$Id $
     2
     3SUBROUTINE phytrac(                            &
     4     nstep,     julien,   gmtime,   debutphy,  &
     5     lafin,     pdtphys,  u, v,     t_seri,     &
     6     paprs,     pplay,    pmfu,     pmfd,      &
     7     pen_u,     pde_u,    pen_d,    pde_d,     &
     8     cdragh,    coefh,    fm_therm, entr_therm,&
     9     yu1,       yv1,      ftsol,    pctsrf,    &
     10     xlat,      frac_impa,frac_nucl,xlon,      &
     11     presnivs,  pphis,    pphi,     albsol,    &
     12     sh,        rh,       cldfra,   rneb,      &
     13     diafra,    cldliq,   itop_con, ibas_con,  &
     14     pmflxr,    pmflxs,   prfl,     psfl,      &
     15     da,        phi,      mp,       upwd,      &
     16     dnwd,      aerosol_couple,     flxmass_w, &
     17     tau_aero,  piz_aero,  cg_aero, ccm,       &
     18     rfname,                                   &
     19     tr_seri)         
     20!
     21!======================================================================
     22! Auteur(s) FH
     23! Objet: Moniteur general des tendances traceurs
     24!======================================================================
     25
     26  USE ioipsl
     27  USE dimphy
     28  USE infotrac
     29  USE mod_grid_phy_lmdz
     30  USE mod_phys_lmdz_para
     31  USE comgeomphy
     32  USE iophy
     33  USE traclmdz_mod
     34  USE tracinca_mod
     35
     36
     37  IMPLICIT NONE
     38
     39  INCLUDE "YOMCST.h"
     40  INCLUDE "dimensions.h"
     41  INCLUDE "indicesol.h"
     42  INCLUDE "clesphys.h"
     43  INCLUDE "temps.h"
     44  INCLUDE "paramet.h"
     45  INCLUDE "control.h"
     46  INCLUDE "thermcell.h"
     47!==========================================================================
     48!                   -- ARGUMENT DESCRIPTION --
     49!==========================================================================
     50
     51! Input arguments
     52!----------------
     53!Configuration grille,temps:
     54  INTEGER,INTENT(IN) :: nstep      ! Appel physique
     55  INTEGER,INTENT(IN) :: julien     ! Jour julien
     56  REAL,INTENT(IN)    :: gmtime
     57  REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
     58  LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
     59  LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
     60 
     61  REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
     62  REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point
     63!
     64!Physique:
     65!--------
     66  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
     67  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       !
     68  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       !
     69  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
     70  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
     71  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
     72  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
     73  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
     74  REAL,DIMENSION(klon),INTENT(IN)        :: pphis
     75  REAL,DIMENSION(klev),INTENT(IN)        :: presnivs
     76  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau liquide nuageuse
     77  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
     78  REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
     79  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
     80  INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
     81  INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
     82  REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
     83!
     84!Convection:
     85!----------
     86  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
     87  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
     88  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
     89  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
     90  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
     91  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
     92
     93!...Tiedke     
     94  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
     95  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
     96
     97  LOGICAL,INTENT(IN)                       :: aerosol_couple
     98  REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
     99  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
     100  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
     101  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
     102  CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname
     103  REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm
     104!... K.Emanuel
     105  REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
     106  REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
     107  REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
     108  REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
     109  REAL,DIMENSION(klon,klev),INTENT(IN)     :: dnwd      ! saturated downdraft mass flux
     110!
     111!Thermiques:
     112!----------
     113  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: fm_therm
     114  REAL,DIMENSION(klon,klev),INTENT(IN)     :: entr_therm
     115!
     116!Couche limite:
     117!--------------
     118!
     119  REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
     120  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
     121  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
     122  REAL,DIMENSION(klon),INTENT(IN)      :: yv1    ! vents au premier niveau
     123!
     124!Lessivage:
     125!----------
     126!
     127! pour le ON-LINE
     128!
     129  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
     130  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
     131
     132! Arguments necessaires pour les sources et puits de traceur:
     133  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
     134  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
     135
     136
     137! Output argument
     138!----------------
     139  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] 
     140
     141!=======================================================================================
     142!                        -- LOCAL VARIABLES --
     143!=======================================================================================
     144
     145  INTEGER :: i, k, it
     146  INTEGER :: nsplit
     147
     148!Sources et Reservoirs de traceurs (ex:Radon):
     149!--------------------------------------------
     150!
     151  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit
     152!$OMP THREADPRIVATE(source)
     153
     154!
     155!Entrees/Sorties: (cf ini_histrac.h et write_histrac.h) 
     156!---------------
     157  INTEGER                   :: iiq, ierr
     158  INTEGER                   :: nhori, nvert
     159  REAL                      :: zsto, zout, zjulian
     160  INTEGER,SAVE              :: nid_tra     ! pointe vers le fichier histrac.nc         
     161!$OMP THREADPRIVATE(nid_tra)
     162  REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
     163  INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
     164  LOGICAL,PARAMETER :: ok_sync=.TRUE.
     165
     166!
     167! Nature du traceur
     168!------------------
     169  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol  ! aerosol(it) = true  => aerosol => lessivage
     170!$OMP THREADPRIVATE(aerosol)                        ! aerosol(it) = false => gaz
     171  REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
     172!
     173! Tendances de traceurs (Td):
     174!------------------------
     175!
     176  REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
     177  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl  ! Td couche limite/traceur
     178  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv  ! Td convection/traceur
     179  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th  ! Td thermique
     180  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction
     181  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation
     182!
     183! Physique
     184!----------   
     185  REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
     186  REAL,DIMENSION(klon,klev)      :: zmasse    ! densité atmosphérique Kg/m2
     187  REAL,DIMENSION(klon,klev)      :: ztra_th
     188 
     189!Controles:
     190!---------
     191  LOGICAL,SAVE :: couchelimite=.TRUE.
     192  LOGICAL,SAVE :: convection=.TRUE.
     193  LOGICAL,SAVE :: lessivage
     194!$OMP THREADPRIVATE(couchelimite,convection,lessivage)
     195
     196  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
     197
     198
     199!######################################################################
     200!                    -- INITIALIZATION --
     201!######################################################################
     202  IF (debutphy) THEN
     203     WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
     204     ALLOCATE( source(klon,nbtr), stat=ierr)
     205     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
     206     
     207     ALLOCATE( aerosol(nbtr), stat=ierr)
     208     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
     209     
     210
     211     ! Initialize module for specific tracers
     212     SELECT CASE(type_trac)
     213     CASE('lmdz')
     214        CALL traclmdz_init(pctsrf, ftsol, aerosol, lessivage)
     215     CASE('inca')
     216        source(:,:)=0.
     217        CALL tracinca_init(aerosol,lessivage)
     218     END SELECT
     219!
     220! Initialize diagnostic output
     221! ----------------------------
     222     INCLUDE "ini_histrac.h"
     223       
     224  END IF
     225!############################################ END INITIALIZATION #######
     226
     227!===============================================================================
     228!    -- Do specific treatment according to chemestry model or local LMDZ tracers
     229!     
     230!===============================================================================
     231  SELECT CASE(type_trac)
     232  CASE('lmdz')
     233     !    -- Traitement des traceurs avec traclmdz
     234     
     235     CALL traclmdz(&
     236          pdtphys,  t_seri,                         &
     237          paprs,    pplay,        cdragh,  coefh,   &
     238          yu1,      yv1,          ftsol,   pctsrf,  &
     239          xlat,     couchelimite,                   &
     240          tr_seri,  source,       solsym,  d_tr_cl)
     241     
     242  CASE('inca')
     243     !    -- CHIMIE INCA  config_inca = aero or chem --
     244
     245     CALL tracinca(&
     246          nstep,    julien,   gmtime,         lafin,     &
     247          pdtphys,  t_seri,   paprs,          pplay,     &
     248          pmfu,     ftsol,    pctsrf,         pphis,     &
     249          pphi,     albsol,   sh,             rh,        &
     250          cldfra,   rneb,     diafra,         cldliq,    &
     251          itop_con, ibas_con, pmflxr,         pmflxs,    &
     252          prfl,     psfl,     aerosol_couple, flxmass_w, &
     253          tau_aero, piz_aero, cg_aero,        ccm,       &
     254          rfname,                                        &
     255          tr_seri,  source,   solsym)     
     256  END SELECT
     257
     258!======================================================================
     259!       -- Calcul de l'effet de la convection --
     260!======================================================================
     261  IF (convection) THEN
     262     DO it=1, nbtr
     263        IF ( conv_flg(it) == 0 ) CYCLE
     264       
     265        IF (iflag_con.LT.2) THEN
     266           d_tr_cv(:,:,:)=0.
     267        ELSE IF (iflag_con.EQ.2) THEN
     268!..Tiedke
     269           CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
     270                pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
     271        ELSE
     272!..K.Emanuel
     273           CALL cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(:,:,it),&
     274                upwd,dnwd,d_tr_cv(:,:,it))
     275        END IF
     276
     277        DO k = 1, klev
     278           DO i = 1, klon       
     279              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
     280           END DO
    439281        END DO
    440282
    441         IF (config_inca == 'aero') THEN
    442            CALL aerosolmain (aerosol_couple,
    443      $                 tr_seri,
    444      $                 pdtphys,
    445      $                 pplay,
    446      $                 pdel,
    447      $                 prfl,
    448      $                 pmflxr,
    449      $                 psfl,
    450      $                 pmflxs,
    451      $                 pmfu,
    452      $                 itop_con,
    453      $                 ibas_con,
    454      $                 pphi,
    455      $                 airephy, ! paire,
    456      $                 nstep,
    457      $                 rneb,         ! for chimiaq
    458      $                 t_seri,       ! for chimiaq
    459      $                 rh,
    460      $                 tau_aero,
    461      $                 piz_aero,
    462      $                 cg_aero,
    463      $                 rfname,
    464      $                 ccm,
    465      $                 lafin)
     283        CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
     284             
     285     END DO ! nbtr
     286  END IF ! convection
     287
     288!======================================================================
     289!    -- Calcul de l'effet des thermiques --
     290!======================================================================
     291
     292  DO k=1,klev
     293     DO i=1,klon
     294        zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
     295     END DO
     296  END DO
     297
     298  DO it=1,nbtr
     299     DO k=1,klev
     300        DO i=1,klon
     301           d_tr_th(i,k,it)=0.
     302           tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
     303           tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
     304        END DO
     305     END DO
     306  END DO
     307 
     308  IF (iflag_thermals.GT.0) THEN   
     309     nsplit=10
     310     DO it=1, nbtr
     311        DO isplit=1,nsplit
     312
     313           CALL dqthermcell(klon,klev,pdtphys/nsplit, &
     314                fm_therm,entr_therm,zmasse, &
     315                tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
     316
     317           DO k=1,klev
     318              DO i=1,klon
     319                 d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
     320                 d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
     321                 tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr(i,k),0.)
     322              END DO
     323           END DO
     324        END DO ! nsplit
     325     END DO ! it
     326  END IF ! Thermiques
     327
     328!======================================================================
     329!     -- Calcul de l'effet de la couche limite --
     330!======================================================================
     331
     332  IF (couchelimite) THEN
     333
     334     DO k = 1, klev
     335        DO i = 1, klon
     336           delp(i,k) = paprs(i,k)-paprs(i,k+1)
     337        END DO
     338     END DO
     339
     340     DO it=1, nbtr
     341       
     342        IF( pbl_flg(it) /= 0 ) THEN
     343       
     344           CALL cltrac(pdtphys, coefh,t_seri,       &
     345                tr_seri(:,:,it), source(:,it),      &
     346                paprs, pplay, delp,                 &
     347                d_tr_cl(:,:,it))
     348           
     349           DO k = 1, klev
     350              DO i = 1, klon
     351                 tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
     352              END DO
     353           END DO
    466354        END IF
    467355
    468         CALL chemmain (tr_seri,    !mmr
    469      $                 nstep,      !nstep
    470      $                 calday,     !calday
    471      $                 julien,     !ncdate
    472      $                 ncsec,      !ncsec
    473      $                 1,          !lat
    474      $                 pdtphys,    !delt
    475      $                 paprs(1,1), !ps
    476      $                 pplay,      !pmid
    477      $                 pdel,       !pdel
    478      $                 airephy,
    479      $                 pctsrf(1,1),!oro
    480      $                 ftsol,      !tsurf
    481      $                 albsol,     !albs
    482      $                 pphi,       !zma
    483      $                 pphis,      !phis
    484      $                 cldfra,     !cldfr
    485      $                 rneb,       !cldfr_st
    486      $                 diafra,     !cldfr_cv
    487      $                 itop_con,   !cldtop
    488      $                 ibas_con,   !cldbot
    489      $                 cldliq,     !cwat
    490      $                 prfl,       !flxrst
    491      $                 pmflxr,     !flxrcv
    492      $                 psfl,       !flxsst
    493      $                 pmflxs,     !flxscv
    494      $                 pmfu,       !flxupd
    495      $                 flxmass_w,  !flxmass_w
    496      $                 t_seri,     !tfld
    497      $                 sh,         !sh
    498      $                 rh,         !rh
    499      $                 iip1,       !nx
    500      $                 jjp1,       !ny
    501      $                 source,
    502      $                 solsym)
    503 
    504 
    505       call VTe(VTinca)
    506       call VTb(VTphysiq)
    507 #endif
    508       END IF ! config_inca
    509 c======================================================================
    510 c   Calcul de l'effet de la convection
    511 c======================================================================
    512 c     print*,'Avant convection'
    513        do it=1,nbtr
    514           WRITE(itn,'(i2)') it
    515 c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
    516        enddo
    517 
    518       if (convection) then
    519 
    520 c      print*,'Pas de temps dans phytrac : ',pdtphys
    521       DO it=1, nbtr
    522 
    523       IF ( config_inca/='none') THEN
    524          IF ( conv_flg(it) == 0 ) CYCLE
    525       END IF
    526 
    527       if (iflag_con.lt.2) then
    528        d_tr_cv=0.
    529       else if (iflag_con.eq.2) then
    530 c tiedke
    531       CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
    532      .            pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
    533       else
    534 c KE
    535       call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
    536      .           upwd,dnwd,d_tr_cv(1,1,it))
    537       endif
    538 
    539        DO k = 1, nlev
    540        DO i = 1, klon
    541          tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
    542        ENDDO
    543        ENDDO
    544 
    545        IF (config_inca == 'none') THEN
    546         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
    547        ELSE
    548         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
    549      .         //solsym(it))
    550        END IF   
    551      
    552       ENDDO
    553 
    554       endif ! convection
    555 c        print*,'Apres convection'
    556 c      do it=1,nbtr
    557 c         WRITE(itn,'(i1)') it
    558 c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
    559 c      enddo
    560 
    561 
    562 c======================================================================
    563 c   Calcul de l'effet des thermiques
    564 c======================================================================
    565 
    566       do k=1,klev
    567          do i=1,klon
    568             zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
    569          enddo
    570       enddo
    571 
    572 c      print*,'masse dans ph ',zmasse
    573       do it=1,nbtr
    574          do k=1,klev
    575             do i=1,klon
    576                d_tr_th(i,k,it)=0.
    577                tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
    578                tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
    579             enddo
    580          enddo
    581       enddo
    582 
    583       if (iflag_thermals.gt.0) then
    584 c        print*,'calcul de leffet des thermiques'
    585         nsplit=10
    586         DO it=1, nbtr
    587 c        WRITE(itn,'(i1)') it
    588 c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
    589 c            print*,'avant dqthermiquesretro'
    590 c             call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI      ')
    591 
    592          do isplit=1,nsplit
    593 c  Abderr 25 11 02
    594 C Thermiques
    595 c       print*,'Avant dans phytrac'
    596             call dqthermcell(klon,klev,pdtphys/nsplit
    597      .       ,fm_therm,entr_therm,zmasse
    598      .       ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
    599 
    600             do k=1,klev
    601                do i=1,klon
    602                   d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
    603                   d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
    604                   tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
    605                enddo
    606             enddo
    607           enddo ! nsplit
    608 c          print*,'apres thermiques'
    609 c             call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th      ')
    610 c            do k=1,klev
    611 c       print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
    612 c          enddo
    613 
    614 c        WRITE(itn,'(i1)') it
    615 c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
    616        ENDDO ! it
    617        endif ! Thermiques
    618 c       print*,'ATTENTION: sdans thermniques'
    619      
    620 c======================================================================
    621 c   Calcul de l'effet de la couche limite
    622 c======================================================================
    623 c       print *,'Avant couchelimite'
    624 c      do it=1,nbtr
    625 c         WRITE(itn,'(i1)') it
    626 c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
    627 c      enddo
    628 
    629       if (couchelimite) then
    630 
    631       DO k = 1, nlev
    632       DO i = 1, klon
    633          delp(i,k) = paprs(i,k)-paprs(i,k+1)
    634       ENDDO
    635       ENDDO
    636 
    637 C maf modif pour tenir compte du cas rnpb + traceur
    638       DO it=1, nbtr
    639 
    640          IF ( config_inca/='none' ) THEN
    641             IF( pbl_flg(it) == 0 ) CYCLE
    642          END IF
    643 
    644 c     print *,'it',it,clsol(it)
    645       if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
    646           CALL cltracrn(it, pdtphys, yu1, yv1,
    647      e                    cdragh, coefh,t_seri,ftsol,pctsrf,
    648      e                    tr_seri(1,1,it),trs(1,it),
    649      e                    paprs, pplay, delp,
    650      e                    masktr(1,it),fshtr(1,it),hsoltr(it),
    651      e                    tautr(it),vdeptr(it),
    652      e                    xlat,
    653      s                    d_tr_cl(1,1,it),d_trs)
    654           DO k = 1, nlev
    655             DO i = 1, klon
    656               tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
    657             ENDDO
    658           ENDDO
    659 c
    660 c Traceur ds sol
    661 c
    662           DO i = 1, klon
    663             trs(i,it) = trs(i,it) + d_trs(i)
    664           END DO
    665 C
    666 C maf provisoire suppression des prints
    667 C         WRITE(itn,'(i1)') it
    668 C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
    669       else ! couche limite avec flux prescrit
    670 
    671          IF (config_inca == 'none') THEN
    672 Cmaf provisoire source / traceur a creer
    673             DO i=1, klon
    674                source(i,it) = 0.0 ! pas de source, pour l'instant
    675             ENDDO
    676          END IF
    677 
    678           CALL cltrac(pdtphys, coefh,t_seri,
    679      s               tr_seri(1,1,it), source(:,it),
    680      e               paprs, pplay, delp,
    681      s               d_tr_cl(1,1,it))
    682             DO k = 1, nlev
    683                DO i = 1, klon
    684                   tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
    685                ENDDO
    686             ENDDO
    687 Cmaf          WRITE(itn,'(i1)') it
    688 cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
    689       endif
    690       ENDDO
    691 c
    692       endif ! couche limite
    693 
    694 c      print*,'Apres couchelimite'
    695 c      do it=1,nbtr
    696 c         WRITE(itn,'(i1)') it
    697 c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
    698 c      enddo
    699 
    700 c======================================================================
    701 c   Calcul de l'effet du puits radioactif
    702 c======================================================================
    703 
    704 C MAF il faudrait faire une modification pour passer dans radiornpb
    705 c si radio=true mais pour l'instant radiornpb propre au cas rnpb
    706       if(rnpb) then
    707 c       print *, 'decroissance radiactive activee'
    708         call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
    709 C
    710         DO it=1,nbtr
    711             if(radio(it)) then
    712             DO k = 1, nlev
    713                DO i = 1, klon
    714                   tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
    715                ENDDO
    716             ENDDO
    717             WRITE(itn,'(i1)') it
    718             CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
    719             endif
    720         ENDDO
    721 c
    722       endif ! rnpb decroissance  radioactive
    723 C
    724 c======================================================================
    725 c   Calcul de l'effet de la precipitation
    726 c======================================================================
    727 
    728 c      print*,'LESSIVAGE =',lessivage
    729       IF (lessivage) THEN
    730 
    731 c     print*,'avant lessivage'
    732 
    733           d_tr_lessi_nucl(:,:,:) = 0.
    734           d_tr_lessi_impa(:,:,:) = 0.
    735           flestottr(:,:,:) = 0.
    736 c
    737 c tendance des aerosols nuclees et impactes
    738 c
    739        DO it = 1, nbtr
    740          IF (aerosol(it)) THEN
    741            DO k = 1, nlev
     356     END DO
     357     
     358  END IF ! couche limite
     359
     360
     361!======================================================================
     362!   Calcul de l'effet de la precipitation
     363!======================================================================
     364
     365  IF (lessivage) THEN
     366     
     367     d_tr_lessi_nucl(:,:,:) = 0.
     368     d_tr_lessi_impa(:,:,:) = 0.
     369     flestottr(:,:,:) = 0.
     370!=========================
     371! LESSIVAGE LARGE SCALE :
     372!=========================
     373
     374! Tendance des aerosols nuclees et impactes
     375! -----------------------------------------
     376     DO it = 1, nbtr
     377        IF (aerosol(it)) THEN
     378           DO k = 1, klev
    742379              DO i = 1, klon
    743                d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
    744      s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
    745                d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
    746      s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
    747               ENDDO
    748            ENDDO
    749          ENDIF
    750        ENDDO
    751 c
    752 c Mises a jour des traceurs + calcul des flux de lessivage
    753 c Mise a jour due a l'impaction et a la nucleation
    754 c
    755 c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
    756 c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
    757 c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
    758        DO it = 1, nbtr
    759 c         print*,'IT=',it,aerosol(it)
    760          IF (aerosol(it)) THEN
    761 c           print*,'IT=',it,' On lessive'
    762            DO k = 1, nlev
    763               DO i = 1, klon
    764                tr_seri(i,k,it)=tr_seri(i,k,it)
    765      s         *frac_impa(i,k)*frac_nucl(i,k)
    766               ENDDO
    767            ENDDO
    768          ENDIF
    769        ENDDO
    770 c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
    771 c
    772 c Flux lessivage total
    773 c
    774       DO it = 1, nbtr
    775            DO k = 1, nlev
    776             DO i = 1, klon
    777                flestottr(i,k,it) = flestottr(i,k,it) -
    778      s                   ( d_tr_lessi_nucl(i,k,it)   +
    779      s                     d_tr_lessi_impa(i,k,it) ) *
    780      s                   ( paprs(i,k)-paprs(i,k+1) ) /
    781      s                   (RG * pdtphys)
    782             ENDDO
    783            ENDDO
    784 c
    785 Cmaf        WRITE(itn,'(i1)') it
    786 Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
    787       ENDDO
    788 c
    789 c     print*,'apres lessivage'
    790       ENDIF
    791 Cc
    792       DO k = 1, nlev
    793          DO i = 1, klon
    794             fluxrn(i,k) = flestottr(i,k,1)
    795             fluxpb(i,k) = flestottr(i,k,2)
    796             rn(i,k) = tr_seri(i,k,1)
    797             pb(i,k) = tr_seri(i,k,2)
    798             pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
    799             pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
    800          ENDDO
    801       ENDDO
    802 
    803 c=============================================================
    804 c   Ecriture des sorties
    805 c=============================================================
    806 
    807 #ifdef CPP_IOIPSL
    808 #include "write_histrac.h"
    809 #endif
    810 
    811 c=============================================================
    812 
    813       if (lafin) then
    814          print*, 'c est la fin de la physique'
    815          call Gather(trs(:,1),trs_tmp)
    816 c$OMP MASTER     
    817          if (is_mpi_root) then
    818          
    819            open (99,file='restarttrac',  form='formatted')
    820            do i=1,klon_glo
    821                write(99,*) trs_tmp(i)
    822            enddo
    823            PRINT*, 'Ecriture du fichier restarttrac'
    824            close(99)
    825          endif
    826 c$OMP END MASTER
    827       else
    828 c         print*, 'physique pas fini'
    829       endif
    830 
    831 
    832       RETURN
    833       END
     380                 d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
     381                      ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
     382                 d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
     383                      ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
     384
     385!
     386! Flux lessivage total
     387! ------------------------------------------------------------
     388                 flestottr(i,k,it) = flestottr(i,k,it) -   &
     389                      ( d_tr_lessi_nucl(i,k,it)   +        &
     390                      d_tr_lessi_impa(i,k,it) ) *          &
     391                      ( paprs(i,k)-paprs(i,k+1) ) /        &
     392                      (RG * pdtphys)
     393!
     394! Mise a jour des traceurs due a l'impaction,nucleation
     395! ----------------------------------------------------------------------
     396                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
     397              END DO
     398           END DO
     399        END IF
     400     END DO
     401     
     402  END IF ! lessivage
     403
     404!=============================================================
     405!   Ecriture des sorties
     406!=============================================================
     407  INCLUDE "write_histrac.h"
     408
     409
     410END SUBROUTINE phytrac

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90

    r1188 r1191  
    11!
    2 ! $Header$
     2! $Id $
    33!
    4       SUBROUTINE radiornpb(tr,dtime,tautr,d_tr)
    5       USE dimphy
    6       USE infotrac, ONLY : nbtr
    7       IMPLICIT none
    8 c======================================================================
    9 c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
    10 c Objet: Decroissance radioactive d'un traceur dans l'atmosphere
    11 CG240694 : Pour un traceur, le radon
    12 CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
    13 c======================================================================
    14 c Arguments:
    15 c======================================================================
    16 cym#include "dimensions.h"
    17 cym#include "dimphy.h"
    18 c======================================================================
    19 C
    20       INTEGER i , k , it
    21       REAL tr(klon,klev,nbtr) , d_tr(klon,klev,nbtr)
    22       REAL dtime
    23       REAL tautr(nbtr)
    24 C
    25 c      WRITE(*,'(''PASSAGE radiornpb ... '',$)')
    26 C Attention, pour un pas de temps beaucoup plus petit que la decroissance!!!
     4SUBROUTINE radio_decay(radio,rnpb,dtime,tautr,tr,d_tr)
     5!
     6! Caluclate radioactive decay for all tracers with radio(it)=true
     7!
     8  USE dimphy
     9  USE infotrac, ONLY : nbtr
     10  IMPLICIT NONE
     11!-----------------------------------------------------------------------
     12! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
     13! Objet: Calcul de la tendance radioactive des traceurs type radioelements
     14!CG240694 : Pour un traceur, le radon
     15!CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
     16!-----------------------------------------------------------------------
     17!
     18! Entrees
     19!
     20  LOGICAL,DIMENSION(nbtr),INTENT(IN)        :: radio ! .true. = traceur radioactif 
     21  LOGICAL,INTENT(IN)                        :: rnpb  ! .true. = decroissance RN = source PB
     22  REAL,INTENT(IN)                           :: dtime ! Pas de temps physique (secondes)
     23  REAL,DIMENSION(nbtr),INTENT(IN)           :: tautr ! Constante de decroissance radioactive
     24  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr    ! Concentrations traceurs U/kgA
     25!
     26! Sortie
     27!
     28  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr  ! Tendance de decroissance radioactive
     29!
     30! Locales
     31!
     32  INTEGER  :: i,k,it
    2733
    28       DO it = 1,2
    29            IF ( tautr(it) .GT. 0. ) THEN
    30                 DO k = 1,klev
    31                 DO i = 1,klon
    32                 d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it)
    33                 END DO
    34                 END DO
    35            ELSE
    36                 DO k = 1,klev
    37                 DO i = 1,klon
    38                 d_tr(i,k,it) = 0.
    39                 END DO
    40                 END DO
    41            END IF
    42       END DO
    43 C
    44 CG161294 : Cas particulier radon 1 => plomb 2
    45 c
    46       DO k = 1,klev
     34
     35  DO it = 1,nbtr
     36     IF ( radio(it) ) THEN
     37        IF (tautr(it) .GT. 0.) THEN
     38           DO k = 1,klev
     39              DO i = 1,klon
     40                 d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it)
     41              END DO
     42           END DO
     43        ELSE
     44           d_tr(:,:,it) = 0.
     45        END IF
     46     ELSE
     47        d_tr(:,:,it) = 0.
     48     END IF
     49  END DO
     50!-------------------------------------------------------
     51!CG161294 : Cas particulier radon [it=1] => plomb [it=2]
     52!-------------------------------------------------------
     53  IF ( rnpb ) THEN
     54     DO k = 1,klev
    4755        DO i = 1,klon
    48           d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1)
     56           d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1)
    4957        ENDDO
    50       ENDDO
    51 c
    52 c      WRITE(*,*) ' radiornpb OK'
    53 c
    54       RETURN
    55       END
     58     ENDDO
     59  ENDIF
     60
     61END SUBROUTINE radio_decay

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h

    r1115 r1191  
     1!$Id $
     2!***************************************
     3!  ECRITURE DU FICHIER :  histrac.nc
     4!***************************************
     5#ifdef CPP_IOIPSL
     6  IF (ecrit_tra > 0. .AND. config_inca == 'none') THEN
     7     
     8     itau_w = itau_phy + nstep
     9     
     10     CALL histwrite_phy(nid_tra,"phis",itau_w,pphis)
     11     CALL histwrite_phy(nid_tra,"aire",itau_w,airephy)
     12
     13!TRACEURS
     14!----------------
     15     DO it=1,nbtr
     16        iiq=niadv(it+2)
     17
     18! CONCENTRATIONS
     19        CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it))
     20
     21! TD LESSIVAGE       
     22        IF (lessivage .AND. aerosol(it)) THEN
     23           CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w,flestottr(:,:,it))
     24        ENDIF
     25
     26! TD THERMIQUES
     27        IF (iflag_thermals.gt.0) THEN
     28           CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w,d_tr_th(:,:,it))
     29        ENDIF
     30
     31! TD CONVECTION
     32        IF (iflag_con.GE.2) THEN
     33           CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w,d_tr_cv(:,:,it))
     34        ENDIF
     35
     36! TD COUCHE-LIMITE
     37        CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w,d_tr_cl(:,:,it))
     38     ENDDO
     39!---------------
    140!
    2 ! $Header$
    341!
     42! VENT (niveau 1)   
     43     CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1)
     44     CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1)
     45!
     46! TEMPERATURE DU SOL
     47     zx_tmp_fi2d(:)=ftsol(:,1)         
     48     CALL histwrite_phy(nid_tra,"ftsol1",itau_w,zx_tmp_fi2d)
     49     zx_tmp_fi2d(:)=ftsol(:,2)
     50     CALL histwrite_phy(nid_tra,"ftsol2",itau_w,zx_tmp_fi2d)
     51     zx_tmp_fi2d(:)=ftsol(:,3)
     52     CALL histwrite_phy(nid_tra,"ftsol3",itau_w,zx_tmp_fi2d)
     53     zx_tmp_fi2d(:)=ftsol(:,4)
     54     CALL histwrite_phy(nid_tra,"ftsol4",itau_w,zx_tmp_fi2d)
     55!     
     56! NATURE DU SOL
     57     zx_tmp_fi2d(:)=pctsrf(:,1)
     58     CALL histwrite_phy(nid_tra,"psrf1",itau_w,zx_tmp_fi2d)
     59     zx_tmp_fi2d(:)=pctsrf(:,2)
     60     CALL histwrite_phy(nid_tra,"psrf2",itau_w,zx_tmp_fi2d)
     61     zx_tmp_fi2d(:)=pctsrf(:,3)
     62     CALL histwrite_phy(nid_tra,"psrf3",itau_w,zx_tmp_fi2d)
     63     zx_tmp_fi2d(:)=pctsrf(:,4)
     64     CALL histwrite_phy(nid_tra,"psrf4",itau_w,zx_tmp_fi2d)
     65 
     66! DIVERS   
     67     CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay)     
     68     CALL histwrite_phy(nid_tra,"t",itau_w,t_seri)     
     69     CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu)
     70     CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd)
     71     CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u)
     72     CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d)
     73     CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d)
     74     CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u)
     75     CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh)
    476
    5       IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
    6       ndex = 0
    7       ndex2d = 0
    8       ndex3d = 0
    9 c
    10       itau_w = itau_phy + nstep
     77     IF (ok_sync) THEN
     78!$OMP MASTER
     79        CALL histsync(nid_tra)
     80!$OMP END MASTER
     81     ENDIF
    1182
    12       CALL histwrite_phy(nid_tra,"phis",itau_w,pphis)
    13 C
    14       CALL histwrite_phy(nid_tra,"aire",itau_w,airephy)
     83  ENDIF !ecrit_tra>0. .AND. config_inca == 'none'
    1584
    16       DO it=1,nbtr
    17        iiq=niadv(it+2)
    18 
    19        CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it))
    20        if (lessivage) THEN
    21        CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w,
    22      .                                   flestottr(:,:,it))
    23       endif
    24      
    25 c----Olivia
    26        CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w,
    27      .                                           d_tr_th(:,:,it))
    28 
    29          if(iflag_con.GE.2) then
    30        CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w,
    31      .                                           d_tr_cv(:,:,it))
    32          endif !(iflag_con.GE.2) then
    33        CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w,
    34      .                                           d_tr_cl(:,:,it))
    35 c---fin Olivia     
    36      
    37       ENDDO
    38 
    39 
    40 C abder
    41          CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1)
    42 
    43          CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1)
    44 
    45          CALL histwrite_phy(nid_tra,"ftsol1",itau_w,pftsol1)
    46 
    47          CALL histwrite_phy(nid_tra,"ftsol2",itau_w,pftsol2)
    48 
    49          CALL histwrite_phy(nid_tra,"ftsol3",itau_w,pftsol3)
    50 
    51          CALL histwrite_phy(nid_tra,"ftsol4",itau_w,pftsol4)
    52 
    53          CALL histwrite_phy(nid_tra,"psrf1",itau_w,ppsrf1)
    54 
    55          CALL histwrite_phy(nid_tra,"psrf2",itau_w,ppsrf2)
    56 
    57          CALL histwrite_phy(nid_tra,"psrf3",itau_w,ppsrf3)
    58 
    59          CALL histwrite_phy(nid_tra,"psrf4",itau_w,ppsrf4)
    60         CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay)
    61 
    62         CALL histwrite_phy(nid_tra,"t",itau_w,t_seri)
    63         CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu)
    64         CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd)
    65         CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u)
    66         CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d)
    67         CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d)
    68         CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u)
    69         CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh)
    70 
    71 
    72 c abder
    73 
    74       if (ok_sync) then
    75 c$OMP MASTER
    76         call histsync(nid_tra)
    77 c$OMP END MASTER
    78        endif
    79 
    80        END IF !ecrit_tra>0. .AND. config_inca == 'none'
    81 
    82 
    83 
     85#endif
Note: See TracChangeset for help on using the changeset viewer.