Changeset 4056

LMDZ6/trunk/libf/dyn3d/advtrac.F90

-                      r4050
+                      r4056
 ! $Id$
+SUBROUTINE advtrac(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
+  !     Auteur :  F. Hourdin
+  !
+  !     Modif. P. Le Van     (20/12/97)
+  !            F. Codron     (10/99)
+  !            D. Le Croller (07/2001)
+  !            M.A Filiberti (04/2002)
+  !
+  USE infotrac, ONLY: nqtot, tracers, nqperes,ok_iso_verif
+  USE control_mod, ONLY: iapp_tracvl, day_step
+  USE comconst_mod, ONLY: dtvr
+  IMPLICIT NONE
+  !
+  include "dimensions.h"
+  include "paramet.h"
+  include "comdissip.h"
+  include "comgeom2.h"
+  include "description.h"
+  include "iniprint.h"
+  !-------------------------------------------------------------------
+  !     Arguments
+  !-------------------------------------------------------------------
+  INTEGER,INTENT(OUT) :: iapptrac
+  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
+  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
+  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
+  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
+  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
+  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
+  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
+  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
+  !-------------------------------------------------------------------
+  !     Ajout PPM
+  !--------------------------------------------------------
+  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
+  !-------------------------------------------------------------
+  !     Variables locales
+  !-------------------------------------------------------------
+  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
+  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
+  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
+  INTEGER iadvtr
+  INTEGER ij,l,iq,iiq
+  REAL zdpmin, zdpmax
+  EXTERNAL  minmax
+  SAVE iadvtr, massem, pbaruc, pbarvc
+  DATA iadvtr/0/
+  !----------------------------------------------------------
+  !     Rajouts pour PPM
+  !----------------------------------------------------------
+  INTEGER indice,n
+  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
+  REAL CFLmaxz,aaa,bbb ! CFL maximum
+  REAL psppm(iim,jjp1) ! pression  au sol
+  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
+  REAL qppm(iim*jjp1,llm,nqtot)
+  REAL fluxwppm(iim,jjp1,llm)
+  REAL apppm(llmp1), bpppm(llmp1)
+  LOGICAL dum,fill
+  DATA fill/.true./
+  DATA dum/.true./
+  integer,save :: countcfl=0
+  real cflx(ip1jmp1,llm)
+  real cfly(ip1jm,llm)
+  real cflz(ip1jmp1,llm)
+  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
+  INTEGER :: iadv
+  IF(iadvtr.EQ.0) THEN
+     pbaruc(:,:)=0
+     pbarvc(:,:)=0
+  ENDIF
+  !   accumulation des flux de masse horizontaux
+  DO l=1,llm
+     DO ij = 1,ip1jmp1
+        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+     ENDDO
+     DO ij = 1,ip1jm
+        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+     ENDDO
+  ENDDO
+  !   selection de la masse instantannee des mailles avant le transport.
+  IF(iadvtr.EQ.0) THEN
+#define DEBUG_IO
+#undef DEBUG_IO
+SUBROUTINE advtrac(pbaru, pbarv, p, masse,q,iapptrac,teta, flxw, pk)
+   !     Auteur :  F. Hourdin
+   !
+   !     Modif. P. Le Van     (20/12/97)
+   !            F. Codron     (10/99)
+   !            D. Le Croller (07/2001)
+   !            M.A Filiberti (04/2002)
+   !
+   USE infotrac,     ONLY: nqtot, tracers,ok_iso_verif
+   USE control_mod,  ONLY: iapp_tracvl, day_step
+   USE comconst_mod, ONLY: dtvr
+   IMPLICIT NONE
+   !
+   include "dimensions.h"
+   include "paramet.h"
+   include "comdissip.h"
+   include "comgeom2.h"
+   include "description.h"
+   include "iniprint.h"
+   !---------------------------------------------------------------------------
+   !     Arguments
+   !---------------------------------------------------------------------------
+   INTEGER, INTENT(OUT) :: iapptrac
+   REAL, INTENT(IN) :: pbaru(ip1jmp1,llm)
+   REAL, INTENT(IN) :: pbarv(ip1jm,  llm)
+   REAL, INTENT(INOUT) ::  q(ip1jmp1,llm,nqtot)
+   REAL, INTENT(IN) :: masse(ip1jmp1,llm)
+   REAL, INTENT(IN) ::     p(ip1jmp1,llmp1 )
+   REAL, INTENT(IN) ::  teta(ip1jmp1,llm)
+   REAL, INTENT(IN) ::    pk(ip1jmp1,llm)
+   REAL, INTENT(OUT) :: flxw(ip1jmp1,llm)
+   !---------------------------------------------------------------------------
+   !     Ajout PPM
+   !---------------------------------------------------------------------------
+   REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm)
+   !---------------------------------------------------------------------------
+   !     Variables locales
+   !---------------------------------------------------------------------------
+   INTEGER :: ij, l, iq, iiq, iadv
+   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
+   REAL :: zdp(ip1jmp1), zdpmin, zdpmax
+   INTEGER, SAVE :: iadvtr=0
+   REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg
+   REAL, DIMENSION(ip1jm,  llm) :: pbarvc, pbarvg
+   EXTERNAL  minmax
+   SAVE massem, pbaruc, pbarvc
+   !---------------------------------------------------------------------------
+   !     Rajouts pour PPM
+   !---------------------------------------------------------------------------
+   INTEGER indice, n
+   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
+   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
+   REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm
+   REAL ::    qppm(iim*jjp1,llm,nqtot)
+   REAL ::   psppm(iim,jjp1)           ! pression  au sol
+   REAL, DIMENSION(llmp1) :: apppm, bpppm
+   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
+   INTEGER, SAVE :: countcfl=0
+   REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz
+   REAL, DIMENSION(ip1jm  ,llm) :: cfly
+   REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax
+   IF(iadvtr == 0) THEN
+      pbaruc(:,:)=0
+      pbarvc(:,:)=0
+   END IF
+   !--- Accumulation des flux de masse horizontaux
+   DO l=1,llm
+      DO ij = 1,ip1jmp1
+         pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+      END DO
+      DO ij = 1,ip1jm
+         pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+      END DO
+   END DO
+   !--- Selection de la masse instantannee des mailles avant le transport.
+   IF(iadvtr == 0) THEN
      CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
+     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
+     !
+  ENDIF
+  iadvtr   = iadvtr+1
+  iapptrac = iadvtr
+  !   Test pour savoir si on advecte a ce pas de temps
+  IF ( iadvtr.EQ.iapp_tracvl ) THEN
+     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
+     !c
+     !   traitement des flux de masse avant advection.
+     !     1. calcul de w
+     !     2. groupement des mailles pres du pole.
+     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
+     ! ... Flux de masse diaganostiques traceurs
+     flxw = wg / REAL(iapp_tracvl)
+     !  test sur l'eventuelle creation de valeurs negatives de la masse
+     DO l=1,llm-1
+        DO ij = iip2+1,ip1jm
+           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
+                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                +       wg(ij,l+1)  - wg(ij,l)
+        ENDDO
+        CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+        DO ij = iip2,ip1jm
+           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
+        ENDDO
+        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
+           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
+                '   MAX:', zdpmax
+        ENDIF
+     ENDDO
+     !-------------------------------------------------------------------
+     ! Calcul des criteres CFL en X, Y et Z
+     !-------------------------------------------------------------------
+     if (countcfl == 0. ) then
+        cflxmax(:)=0.
+        cflymax(:)=0.
+        cflzmax(:)=0.
+     endif
+     countcfl=countcfl+iapp_tracvl
+     cflx(:,:)=0.
+     cfly(:,:)=0.
+     cflz(:,:)=0.
+     do l=1,llm
+        do ij=iip2,ip1jm-1
+           if (pbarug(ij,l)>=0.) then
+              cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
+           else
+              cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
+           endif
+        enddo
+     enddo
+     do l=1,llm
+        do ij=iip2,ip1jm-1,iip1
+           cflx(ij+iip1,l)=cflx(ij,l)
+        enddo
+     enddo
+     do l=1,llm
+        do ij=1,ip1jm
+           if (pbarvg(ij,l)>=0.) then
+              cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
+           else
+              cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
+           endif
+        enddo
+     enddo
+     do l=2,llm
+        do ij=1,ip1jm
+           if (wg(ij,l)>=0.) then
+              cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
+           else
+              cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
+           endif
+        enddo
+     enddo
+     do l=1,llm
+        cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
+        cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
+        cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
+     enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! Par defaut, on sort le diagnostic des CFL tous les jours.
+     ! Si on veut le sortir a chaque pas d'advection en cas de plantage
+     !     if (countcfl==iapp_tracvl) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     if (countcfl==day_step) then
+        do l=1,llm
+           write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
+                cflzmax(l)
+        enddo
+        countcfl=0
+     endif
+     !-------------------------------------------------------------------
+     !   Advection proprement dite (Modification Le Croller (07/2001)
+     !-------------------------------------------------------------------
+     !----------------------------------------------------
+     !        Calcul des moyennes basées sur la masse
+     !----------------------------------------------------
+     call massbar(massem,massebx,masseby)
+     !-----------------------------------------------------------
+     !     Appel des sous programmes d'advection
+     !-----------------------------------------------------------
+     if (ok_iso_verif) then
+           write(*,*) 'advtrac 227'
+           call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
+     endif !if (ok_iso_verif) then
+     do iq=1,nqperes
+        !        call clock(t_initial)
+       iadv = tracers(iq)%iadv
+       SELECT CASE(iadv)
+   ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
+   END IF
+   iadvtr   = iadvtr+1
+   iapptrac = iadvtr
+   !--- Test pour savoir si on advecte a ce pas de temps
+   IF(iadvtr /= iapp_tracvl) RETURN
+   !   ..  Modif P.Le Van  ( 20/12/97 )  ....
+   !
+   !   traitement des flux de masse avant advection.
+   !       1. calcul de w
+   !       2. groupement des mailles pres du pole.
+   CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg)
+   !--- Flux de masse diaganostiques traceurs
+   flxw = wg / REAL(iapp_tracvl)
+   !--- Test sur l'eventuelle creation de valeurs negatives de la masse
+   DO l=1,llm-1
+      DO ij = iip2+1,ip1jm
+         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
+                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                 +     wg(ij,l+1)    -     wg(ij,l)
+      END DO
+! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
+      CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+      DO ij = iip2,ip1jm
+         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
+      END DO
+      CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+      IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) &
+         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,' MAX:', zdpmax
+   END DO
+   !-------------------------------------------------------------------------
+   ! Calcul des criteres CFL en X, Y et Z
+   !-------------------------------------------------------------------------
+   IF(countcfl == 0. ) then
+      cflxmax(:)=0.
+      cflymax(:)=0.
+      cflzmax(:)=0.
+   END IF
+   countcfl=countcfl+iapp_tracvl
+   cflx(:,:)=0.
+   cfly(:,:)=0.
+   cflz(:,:)=0.
+   DO l=1,llm
+      DO ij=iip2,ip1jm-1
+         IF(pbarug(ij,l)>=0.) then
+            cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
+         END IF
+      END DO
+   END DO
+   DO l=1,llm
+      DO ij=iip2,ip1jm-1,iip1
+         cflx(ij+iip1,l)=cflx(ij,l)
+      END DO
+   END DO
+   DO l=1,llm
+      DO ij=1,ip1jm
+         IF(pbarvg(ij,l)>=0.) then
+            cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
+         END IF
+      END DO
+   END DO
+   DO l=2,llm
+      DO ij=1,ip1jm
+         IF(wg(ij,l) >= 0.) THEN
+            cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
+         END IF
+      END DO
+   END DO
+   DO l=1,llm
+      cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
+      cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
+      cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
+   END DO
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Par defaut, on sort le diagnostic des CFL tous les jours.
+   ! Si on veut le sortir a chaque pas d'advection en cas de plantage
+   !       IF(countcfl==iapp_tracvl) then
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   IF(countcfl==day_step) then
+      DO l=1,llm
+         WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l)
+      END DO
+      countcfl=0
+   END IF
+   !---------------------------------------------------------------------------
+   !   Advection proprement dite (Modification Le Croller (07/2001)
+   !---------------------------------------------------------------------------
+   !---------------------------------------------------------------------------
+   !   Calcul des moyennes basees sur la masse
+   !---------------------------------------------------------------------------
+   CALL massbar(massem,massebx,masseby)
+#ifdef DEBUG_IO
+   CALL WriteField_u('massem',massem)
+   CALL WriteField_u('wg',wg)
+   CALL WriteField_u('pbarug',pbarug)
+   CALL WriteField_v('pbarvg',pbarvg)
+   CALL WriteField_u('p_tmp',p)
+   CALL WriteField_u('pk_tmp',pk)
+   CALL WriteField_u('teta_tmp',teta)
+   DO iq=1,nqtot
+      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
+   END DO
+#endif
+   IF(ok_iso_verif) THEN
+      WRITE(*,*) 'advtrac 227'
+      CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162')
+   END IF
+   !-------------------------------------------------------------------------
+   !       Appel des sous programmes d'advection
+   !-------------------------------------------------------------------------
+   DO iq = 1, nqtot
+!     CALL clock(t_initial)
+      IF(tracers(iq)%parent /= 'air') CYCLE
+      iadv = tracers(iq)%iadv
+      !-----------------------------------------------------------------------
+      SELECT CASE(iadv)
+      !-----------------------------------------------------------------------
          CASE(0); CYCLE
+         CASE(10)
+        !   ----------------------------------------------------------------
+        !   Schema de Van Leer I MUSCL
+        !   ----------------------------------------------------------------
+           ! CRisi: on fait passer tout q pour avoir acces aux fils
+           !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
+           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
+         CASE(14)
+           !   ----------------------------------------------------------------
+           !   Schema "pseudo amont" + test sur humidite specifique
+           !    pour la vapeur d'eau. F. Codron
+           !   ----------------------------------------------------------------
+           !
+           !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
+           CALL vlspltqs( q, 2., massem, wg , &
+                pbarug,pbarvg,dtvr,p,pk,teta,iq)
+         CASE(12)
+           !   ----------------------------------------------------------------
+           !   Schema de Frederic Hourdin
+           !   ----------------------------------------------------------------
+           !            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+           end do
+         CASE(13)
+           !            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+           end do
+         CASE(20)
+           !   ----------------------------------------------------------------
+           !   Schema de pente SLOPES
+           !   ----------------------------------------------------------------
+           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+         CASE(30)
+           !   ----------------------------------------------------------------
+           !   Schema de Prather
+           !   ----------------------------------------------------------------
+           !            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
+                n,dtbon)
+         CASE(11,16,17,18)
+           !   ----------------------------------------------------------------
+           !   Schemas PPM Lin et Rood
+           !   ----------------------------------------------------------------
+           !        Test sur le flux horizontal
+           !        Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           !        Test sur le flux vertical
+           CFLmaxz=0.
+           do l=2,llm
+              do ij=iip2,ip1jm
+                 aaa=wg(ij,l)*dtvr/massem(ij,l)
+                 CFLmaxz=max(CFLmaxz,aaa)
+                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                 CFLmaxz=max(CFLmaxz,bbb)
+              enddo
+           enddo
+           if (CFLmaxz.GE.1) then
+              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+           endif
+           !-----------------------------------------------------------
+           !        Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
+           !-----------------------------------------------------------
+           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                unatppm,vnatppm,psppm)
+           do indice=1,n
+              !----------------------------------------------------------------
+              !                         VL (version PPM) horiz. et PPM vert.
+              !----------------------------------------------------------------
+             SELECT CASE(iadv)
+               CASE(11)
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+               CASE(16)
+                 !-------------------------------------------------------------
+                 !                           Monotonic PPM
+                 !-------------------------------------------------------------
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+               CASE(17)
+                 !-------------------------------------------------------------
+                 !                           Semi Monotonic PPM
+                 !-------------------------------------------------------------
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+               CASE(18)
+                 !-------------------------------------------------------------
+                 !                         Positive Definite PPM
+                 !-------------------------------------------------------------
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+             END SELECT
+           enddo
+           !-----------------------------------------------------------------
+           !               Ss-prg interface PPM3d-LMDZ.4
+           !-----------------------------------------------------------------
+           call interpost(q(1,1,iq),qppm(1,1,iq))
+         END SELECT
+        !----------------------------------------------------------------------
+        !-----------------------------------------------------------------
+        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
+        ! et Nord j=1
+        !-----------------------------------------------------------------
+        !                  call traceurpole(q(1,1,iq),massem)
+        ! calcul du temps cpu pour un schema donne
+        !                  call clock(t_final)
+        !ym                  tps_cpu=t_final-t_initial
+        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
+     end DO
+     if (ok_iso_verif) then
+           write(*,*) 'advtrac 402'
+           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
+     endif !if (ok_iso_verif) then
+     !------------------------------------------------------------------
+     !   on reinitialise a zero les flux de masse cumules
+     !---------------------------------------------------
+     iadvtr=0
+  ENDIF ! if iadvtr.EQ.iapp_tracvl
+         !--------------------------------------------------------------------
+         CASE(10)  !--- Schema de Van Leer I MUSCL
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
+            CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
+         !--------------------------------------------------------------------
+         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
+                   !--- pour la vapeur d'eau. F. Codron
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
+            CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq)
+         !--------------------------------------------------------------------
+         CASE(12)  !--- Schema de Frederic Hourdin
+         !--------------------------------------------------------------------
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+            END DO
+         !--------------------------------------------------------------------
+         CASE(13)  !--- Pas de temps adaptatif
+         !--------------------------------------------------------------------
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+            END DO
+         !--------------------------------------------------------------------
+         CASE(20)  !--- Schema de pente SLOPES
+         !--------------------------------------------------------------------
+            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+         !--------------------------------------------------------------------
+         CASE(30)  !--- Schema de Prather
+         !--------------------------------------------------------------------
+            ! Pas de temps adaptatif
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
+         !--------------------------------------------------------------------
+         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
+         !--------------------------------------------------------------------
+            ! Test sur le flux horizontal
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            ! Test sur le flux vertical
+            CFLmaxz=0.
+            DO l=2,llm
+               DO ij=iip2,ip1jm
+                  aaa=wg(ij,l)*dtvr/massem(ij,l)
+                  CFLmaxz=max(CFLmaxz,aaa)
+                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                  CFLmaxz=max(CFLmaxz,bbb)
+               END DO
+            END DO
+            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+            !----------------------------------------------------------------
+            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
+            !----------------------------------------------------------------
+            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                 unatppm,vnatppm,psppm)
+            !----------------------------------------------------------------
+            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
+            !----------------------------------------------------------------
+               SELECT CASE(iadv)
+                  !----------------------------------------------------------
+                  CASE(11)
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(16) !--- Monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(17) !--- Semi monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(18) !--- Positive Definite PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+               END SELECT
+            !----------------------------------------------------------------
+            END DO
+            !----------------------------------------------------------------
+            !     Ss-prg interface PPM3d-LMDZ.4
+            !----------------------------------------------------------------
+            CALL interpost(q(1,1,iq),qppm(1,1,iq))
+      !----------------------------------------------------------------------
+      END SELECT
+      !----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
+      !----------------------------------------------------------------------
+      !  CALL traceurpole(q(1,1,iq),massem)
+      !--- Calcul du temps cpu pour un schema donne
+      !  CALL clock(t_final)
+      !ym  tps_cpu=t_final-t_initial
+      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
+   END DO
+   IF(ok_iso_verif) then
+      WRITE(*,*) 'advtrac 402'
+      CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397')
+   END IF
+   !-------------------------------------------------------------------------
+   !   on reinitialise a zero les flux de masse cumules
+   !-------------------------------------------------------------------------
+   iadvtr=0
 END SUBROUTINE advtrac

LMDZ6/trunk/libf/dyn3d/iniacademic.F90

-                      r4052
+                      r4056
   USE filtreg_mod, ONLY: inifilr
   USE infotrac,    ONLY: nqtot,niso_possibles,ok_isotopes,ok_iso_verif,tnat,alpha_ideal, &
         & iqiso,iso_indnum, tracers
+  USE infotrac,    ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, &
+                         iqiso, tracers, iso_indnum
   USE control_mod, ONLY: day_step,planet_type
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
 #ifdef CPP_IOIPSL
   USE IOIPSL, ONLY: getin
 …
 #endif
   USE Write_Field
-  use exner_hyb_m, only: exner_hyb
-  use exner_milieu_m, only: exner_milieu
   USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
   USE logic_mod, ONLY: iflag_phys, read_start
 …
   real tetastrat ! potential temperature in the stratosphere, in K
   real tetajl(jjp1,llm)
   INTEGER i,j,l,lsup,ij, iq
+  INTEGER i,j,l,lsup,ij, iq, iName, iZone, iPhase, iqParent
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
 …
   real zz,ran1
   integer idum, iqParent, iName, iZone, iPhase
+  integer idum
   REAL zdtvr
 …
   character(len=*),parameter :: modname="iniacademic"
   character(len=80) :: abort_message
   ! Sanity check: verify that options selected by user are not incompatible
 …
         ps(:)=preff
      endif
      ! ground geopotential
      phis(:)=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
      if (pressure_exner) then
        CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
 …
      CALL getin('tetanoise',tetanoise)
      ! 2. Initialize fields towards which to relax
      ! Friction
 …
      IF (.NOT. read_start) THEN
         ! bulk initialization of temperature
         IF (iflag_phys>10000) THEN
         ! Particular case to impose a constant temperature T0=0.01*iflag_physx
 …
            teta(:,:)=tetarappel(:,:)
         ENDIF
         ! geopotential
         CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
 …
               iName = tracers(iq)%iso_iName
               if (.NOT.ok_isotopes .OR. iName <= 0) CYCLE
               iZone = tracers(iq)%iso_iZone
               iPhase= tracers(iq)%iso_iPhase
+              iZone    = tracers(iq)%iso_iZone
+              iPhase   = tracers(iq)%iso_iPhase
               iqParent = tracers(iq)%iqParent
               if (iZone == 0) q(:,:,iq) = q(:,:,iqParent)*tnat(iName) &
 …
         endif ! of if (planet_type=="earth")
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
+        endif !if (ok_iso_verif) then
+        if (ok_iso_verif) call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
         ! add random perturbation to temperature

LMDZ6/trunk/libf/dyn3d_common/infotrac.F90

-                      r4052
+                      r4056
   INTEGER, SAVE :: nbtr
+! CRisi: on retranche les isotopes des traceurs habituels
+! On fait un tableaux d'indices des traceurs qui passeront dans phytrac
+! Nombre de traceurs passes a phytrac
   INTEGER, SAVE :: nqtottr
-  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: itr_indice
-! CRisi: nb traceurs peres= directement advectes par l'air
-  INTEGER, SAVE :: nqperes
-! ThL: nb traceurs INCA
-  INTEGER, SAVE :: nqINCA
 ! ThL: nb traceurs CO2
 …
   TYPE(trac_type), TARGET,  SAVE, ALLOCATABLE ::  tracers(:)    !=== TRACERS DESCRIPTORS VECTOR
   TYPE(isot_type), TARGET,  SAVE, ALLOCATABLE :: isotopes(:)    !=== ISOTOPES PARAMETERS VECTOR
-! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the
-!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
-  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
   REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi
 …
   INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
-#ifdef CPP_StratAer
-!--CK/OB for stratospheric aerosols
-  INTEGER, SAVE :: nbtr_bin
-  INTEGER, SAVE :: nbtr_sulgas
-  INTEGER, SAVE :: id_OCS_strat
-  INTEGER, SAVE :: id_SO2_strat
-  INTEGER, SAVE :: id_H2SO4_strat
-  INTEGER, SAVE :: id_BIN01_strat
-  INTEGER, SAVE :: id_TEST_strat
-#endif
 CONTAINS
 …
     LOGICAL :: continu,nouveau_traceurdef
     INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
     CHARACTER(len=maxlen) :: tchaine
+    CHARACTER(len=maxlen) :: tchaine, msg1
     INTEGER, ALLOCATABLE  :: iqfils(:,:)
+    INTEGER :: nqINCA
     character(len=*),parameter :: modname="infotrac_init"
 …
     descrq(30)='PRA'
+    ! Coherence test between parameter type_trac, config_inca and preprocessing keys
+    IF (type_trac=='inca') THEN
+       WRITE(lunout,*) 'You have chosen to couple with INCA chemistry model : type_trac=', &
+            type_trac,' config_inca=',config_inca
+       IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
+          WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
+          CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
+       ENDIF
+    !--- MESSAGE ABOUT THE CHOSEN CONFIGURATION
+    WRITE(lunout,*)'type_trac='//TRIM(type_trac)
+    msg1 = 'For type_trac = "'//TRIM(type_trac)//'":'
+    SELECT CASE(type_trac)
+      CASE('inca'); WRITE(lunout,*)TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca
+      CASE('inco'); WRITE(lunout,*)TRIM(msg1)//' coupling jointly with INCA and CO2 cycle'
+      CASE('repr'); WRITE(lunout,*)TRIM(msg1)//' coupling with REPROBUS chemistry model'
+      CASE('co2i'); WRITE(lunout,*)TRIM(msg1)//' you have chosen to run with CO2 cycle'
+      CASE('coag'); WRITE(lunout,*)TRIM(msg1)//' tracers are treated for COAGULATION tests'
+      CASE('lmdz'); WRITE(lunout,*)TRIM(msg1)//' tracers are treated in LMDZ only'
+      CASE DEFAULT
+        CALL abort_gcm(modname,'type_trac='//TRIM(type_trac)//' not possible yet.',1)
+    END SELECT
+    !--- COHERENCE TEST BETWEEN "type_trac", "config_inca" AND PREPROCESSING KEYS
+    SELECT CASE(type_trac)
+      CASE('inca','inco'); IF(ALL(['aero', 'aeNP', 'chem']/=config_inca)) &
+        CALL abort_gcm(modname, 'Incoherence between type_trac and config_inca. Please modify "run.def"',1)
 #ifndef INCA
+       WRITE(lunout,*) 'To run this option you must add cpp key INCA and compile with INCA code'
+       CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
+        CALL abort_gcm(modname, 'You must add cpp key INCA and compile with INCA code',1)
 #endif
+    ELSE IF (type_trac=='repr') THEN
+       WRITE(lunout,*) 'You have chosen to couple with REPROBUS chemestry model : type_trac=', type_trac
+      CASE('repr')
 #ifndef REPROBUS
+       WRITE(lunout,*) 'To run this option you must add cpp key REPROBUS and compile with REPRPBUS code'
+       CALL abort_gcm('infotrac_init','You must compile with cpp key REPROBUS',1)
+        CALL abort_gcm(modname, 'You must add cpp key REPROBUS and compile with REPROBUS code',1)
 #endif
+    ELSE IF (type_trac == 'co2i') THEN
+       WRITE(lunout,*) 'You have chosen to run with CO2 cycle: type_trac=', type_trac
+    ELSE IF (type_trac == 'coag') THEN
+       WRITE(lunout,*) 'Tracers are treated for COAGULATION tests : type_trac=', type_trac
+      CASE('coag')
 #ifndef CPP_StratAer
+       WRITE(lunout,*) 'To run this option you must add cpp key StratAer and compile with StratAer code'
+       CALL abort_gcm('infotrac_init','You must compile with cpp key StratAer',1)
+        CALL abort_gcm(modname, 'You must add cpp key StratAer and compile with StratAer code',1)
 #endif
+    ELSE IF (type_trac == 'lmdz') THEN
+       WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac
+    ELSE IF (type_trac == 'inco') THEN ! ThL
+       WRITE(lunout,*) 'Using jointly INCA and CO2 cycle: type_trac =', type_trac
+       IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
+          WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
+          CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
+       ENDIF
+#ifndef INCA
+       WRITE(lunout,*) 'To run this option you must add cpp key INCA and compilewith INCA code'
+       CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
+#endif
+    ELSE
+       WRITE(lunout,*) 'type_trac=',type_trac,' not possible. Model stops'
+       CALL abort_gcm('infotrac_init','bad parameter',1)
+    ENDIF
+    ! Test if config_inca is other then none for run without INCA
+    IF (type_trac/='inca' .AND. type_trac/='inco' .AND. config_inca/='none') THEN
+       WRITE(lunout,*) 'config_inca will now be changed to none as you do not couple with INCA model'
+       config_inca='none'
+    ENDIF
+    END SELECT
+    !--- Disable "config_inca" option for a run without INCA if it differs from "none"
+    IF (ALL(['inca', 'inco', 'none'] /= config_inca)) THEN
+      WRITE(lunout,*)'setting config_inca="none" as you do not couple with INCA model'
+      config_inca = 'none'
+    END IF
 !-----------------------------------------------------------------------
 …
+!
 !-----------------------------------------------------------------------
+    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
+       IF (type_trac=='co2i') THEN
+          nqCO2 = 1
+       ELSE
+          nqCO2 = 0
+       ENDIF
+       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
+       IF(ierr.EQ.0) THEN
+          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
+          READ(90,*) nqtrue
+          write(lunout,*) 'nqtrue=',nqtrue
+       ELSE
+          WRITE(lunout,*) trim(modname),': Failed opening traceur.def'
+          CALL abort_gcm(modname,"file traceur.def not found!",1)
+       ENDIF
+!jyg<
+!!       if ( planet_type=='earth') then
+!!         ! For Earth, water vapour & liquid tracers are not in the physics
+!!         nbtr=nqtrue-2
+!!       else
+!!         ! Other planets (for now); we have the same number of tracers
+!!         ! in the dynamics than in the physics
+!!         nbtr=nqtrue
+!!       endif
+!>jyg
+    ELSE ! type_trac=inca or inco
+       IF (type_trac=='inco') THEN
+          nqCO2 = 1
+       ELSE
+          nqCO2 = 0
+       ENDIF
+!jyg<
+       ! The traceur.def file is used to define the number "nqo" of water phases
+       ! present in the simulation. Default : nqo = 2.
+       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
+       IF(ierr.EQ.0) THEN
+          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
+          READ(90,*) nqo
+       ELSE
+          WRITE(lunout,*) trim(modname),': Failed opening traceur.def'
+          CALL abort_gcm(modname,"file traceur.def not found!",1)
+       ENDIF
+       IF (nqo /= 2 .AND. nqo /= 3 ) THEN
+          IF (nqo == 4 .AND. type_trac=='inco') THEN ! ThL
+             WRITE(lunout,*) trim(modname),': you are coupling with INCA, and also using CO2i.'
+             nqo = 3    ! A ameliorier... je force 3 traceurs eau...  ThL
+             WRITE(lunout,*) trim(modname),': nqo = ',nqo
+          ELSE
+          WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowed. Only 2 or 3 water phases allowed'
+          CALL abort_gcm('infotrac_init','Bad number of water phases',1)
+          ENDIF
+       ENDIF
+       ! nbtr has been read from INCA by init_const_lmdz() in gcm.F
+    OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
+    IF(ierr.EQ.0) THEN
+       WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
+    ELSE
+       WRITE(lunout,*) trim(modname),': Failed opening traceur.def'
+       CALL abort_gcm(modname,"file traceur.def not found!",1)
+    ENDIF
+    nqCO2 = 0; IF(type_trac == 'inco') nqCO2 = 1
+    SELECT CASE(type_trac)
+       CASE('lmdz','repr','coag','co2i'); READ(90,*) nqtrue
+       CASE('inca','inco');               READ(90,*) nqo
+          ! The traceur.def file is used to define the number "nqo" of water phases
+          ! present in the simulation. Default : nqo = 2.
+          IF (nqo == 4 .AND. type_trac=='inco') nqo = 3
+          IF(ALL([2,3] /= nqo)) THEN
+             WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowed. Only 2 or 3 water phases allowed'
+             CALL abort_gcm('infotrac_init','Bad number of water phases',1)
+          END IF
+          ! nbtr has been read from INCA by init_const_lmdz() in gcm.F
 #ifdef INCA
        CALL Init_chem_inca_trac(nqINCA)
+          CALL Init_chem_inca_trac(nqINCA)
 #else
        nqINCA=0
+          nqINCA=0
 #endif
        nbtr=nqINCA+nqCO2
        nqtrue=nbtr+nqo
        WRITE(lunout,*) trim(modname),': nqo = ',nqo
        WRITE(lunout,*) trim(modname),': nbtr = ',nbtr
        WRITE(lunout,*) trim(modname),': nqtrue = ',nqtrue
        WRITE(lunout,*) trim(modname),': nqCO2 = ',nqCO2
        WRITE(lunout,*) trim(modname),': nqINCA = ',nqINCA
        ALLOCATE(hadv_inca(nqINCA), vadv_inca(nqINCA), conv_flg_inca(nqINCA), pbl_flg_inca(nqINCA), solsym_inca(nqINCA))
     ENDIF   ! type_trac 'inca' ou 'inco'
+          nbtr=nqINCA+nqCO2
+          nqtrue=nbtr+nqo
+          WRITE(lunout,*) trim(modname),': nqo    = ',nqo
+          WRITE(lunout,*) trim(modname),': nbtr   = ',nbtr
+          WRITE(lunout,*) trim(modname),': nqtrue = ',nqtrue
+          WRITE(lunout,*) trim(modname),': nqCO2 = ',nqCO2
+          WRITE(lunout,*) trim(modname),': nqINCA = ',nqINCA
+          ALLOCATE(hadv_inca(nqINCA), vadv_inca(nqINCA), conv_flg_inca(nqINCA), pbl_flg_inca(nqINCA), solsym_inca(nqINCA))
+    END SELECT
 !>jyg
 …
 !    Get choice of advection schema from file tracer.def or from INCA
 !---------------------------------------------------------------------
+    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
+          ! Continue to read tracer.def
+          DO iq=1,nqtrue
+             write(*,*) 'infotrac 237: iq=',iq
+             ! CRisi: ajout du nom du fluide transporteur
+             ! mais rester retro compatible
+             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+             write(lunout,*) 'tchaine=',trim(tchaine)
+             write(*,*) 'infotrac 238: IOstatus=',IOstatus
+             if (IOstatus.ne.0) then
+                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+             endif
+             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
+             ! espace ou pas au milieu de la chaine.
+             continu=.true.
+             nouveau_traceurdef=.false.
+             iiq=1
+             do while (continu)
+                if (tchaine(iiq:iiq).eq.' ') then
+                  nouveau_traceurdef=.true.
+                  continu=.false.
+                else if (iiq.lt.LEN_TRIM(tchaine)) then
+                  iiq=iiq+1
+                else
+                  continu=.false.
+                endif
+             enddo
+             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+             if (nouveau_traceurdef) then
+                write(lunout,*) 'C''est la nouvelle version de traceur.def'
+                tnom_0(iq)=TRIM(tchaine(1:iiq-1))
+                tnom_transp(iq)=TRIM(tchaine(iiq+1:))
+             else
+                write(lunout,*) 'C''est l''ancienne version de traceur.def'
+                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
+                tnom_0(iq)=tchaine
+                tnom_transp(iq) = 'air'
+             endif
+             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
+          ENDDO!DO iq=1,nqtrue
+          CLOSE(90)
+    IF (ANY(['lmdz', 'repr', 'coag', 'co2i'] == type_trac)) THEN
+       ! Continue to read tracer.def
+       DO iq=1,nqtrue
+          write(*,*) 'infotrac 237: iq=',iq
+          ! CRisi: ajout du nom du fluide transporteur
+          ! mais rester retro compatible
+          READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+          write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+          write(lunout,*) 'tchaine=',trim(tchaine)
+          write(*,*) 'infotrac 238: IOstatus=',IOstatus
+          if (IOstatus.ne.0) CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+          ! Y-a-t-il 1 ou 2 noms de traceurs separes par un espace ?
+          iiq = INDEX(tchaine(1:LEN_TRIM(tchaine)),' ')
+          nouveau_traceurdef=iiq/=0
+          write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+          if (nouveau_traceurdef) then
+             IF(iq==1) write(lunout,*) "Nouvelle version de traceur.def"
+             tnom_0     (iq) = tchaine(1:iiq-1)
+             tnom_transp(iq) = tchaine(iiq+1:)
+          else
+             IF(iq==1) write(lunout,*) "Nouvelle version de traceur.def: traceurs d'air uniquement"
+             tnom_0     (iq) = tchaine
+             tnom_transp(iq) = 'air'
+          endif
+          write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+          write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
+       ENDDO!DO iq=1,nqtrue
+       CLOSE(90)
        WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
 …
        IF ( planet_type=='earth') THEN
          !CR: nombre de traceurs de l eau
+         IF (tnom_0(3) == 'H2Oi') THEN
+            nqo=3
+         ELSE
+            nqo=2
+         ENDIF
+         nqo=2; IF (tnom_0(3) == 'H2Oi') nqo=3
          ! For Earth, water vapour & liquid tracers are not in the physics
          nbtr=nqtrue-nqo
 …
        ENDIF
+#ifdef CPP_StratAer
+       IF (type_trac == 'coag') THEN
+         nbtr_bin=0
+         nbtr_sulgas=0
+         DO iq=1,nqtrue
+           IF (tnom_0(iq)(1:3)=='BIN') THEN !check if tracer name contains 'BIN'
+             nbtr_bin=nbtr_bin+1
+           ENDIF
+           IF (tnom_0(iq)(1:3)=='GAS') THEN !check if tracer name contains 'GAS'
+             nbtr_sulgas=nbtr_sulgas+1
+           ENDIF
+         ENDDO
+         print*,'nbtr_bin=',nbtr_bin
+         print*,'nbtr_sulgas=',nbtr_sulgas
+         DO iq=1,nqtrue
+           IF (tnom_0(iq)=='GASOCS') THEN
+             id_OCS_strat=iq-nqo
+           ENDIF
+           IF (tnom_0(iq)=='GASSO2') THEN
+             id_SO2_strat=iq-nqo
+           ENDIF
+           IF (tnom_0(iq)=='GASH2SO4') THEN
+             id_H2SO4_strat=iq-nqo
+           ENDIF
+           IF (tnom_0(iq)=='BIN01') THEN
+             id_BIN01_strat=iq-nqo
+           ENDIF
+           IF (tnom_0(iq)=='GASTEST') THEN
+             id_TEST_strat=iq-nqo
+           ENDIF
+         ENDDO
+         print*,'id_OCS_strat  =',id_OCS_strat
+         print*,'id_SO2_strat  =',id_SO2_strat
+         print*,'id_H2SO4_strat=',id_H2SO4_strat
+         print*,'id_BIN01_strat=',id_BIN01_strat
+       ENDIF
+    ENDIF
+!jyg<
+!
+! Transfert number of tracers to Reprobus
+#ifdef REPROBUS
+    IF (type_trac == 'repr') CALL Init_chem_rep_trac(nbtr,nqo,tnom_0)
 #endif
-    ENDIF  ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag' .OR. type_trac = 'co2i')
-!jyg<
+!
-! Transfert number of tracers to Reprobus
-    IF (type_trac == 'repr') THEN
-#ifdef REPROBUS
-       CALL Init_chem_rep_trac(nbtr,nqo,tnom_0)
-#endif
-    ENDIF
+!
 ! Allocate variables depending on nbtr
 …
     IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN   ! config_inca='aero' ou 'chem'
 !>jyg
+! le module de chimie fournit les noms des traceurs
+! et les schemas d'advection associes. excepte pour ceux lus
+! dans traceur.def
+          DO iq=1,nqo+nqCO2
+             write(*,*) 'infotrac 237: iq=',iq
+             ! CRisi: ajout du nom du fluide transporteur
+             ! mais rester retro compatible
+             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+             write(lunout,*) 'tchaine=',trim(tchaine)
+             write(*,*) 'infotrac 238: IOstatus=',IOstatus
+             if (IOstatus.ne.0) then
+                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+             endif
+             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
+             ! espace ou pas au milieu de la chaine.
+             continu=.true.
+             nouveau_traceurdef=.false.
+             iiq=1
+             do while (continu)
+                if (tchaine(iiq:iiq).eq.' ') then
+                  nouveau_traceurdef=.true.
+                  continu=.false.
+                else if (iiq.lt.LEN_TRIM(tchaine)) then
+                  iiq=iiq+1
+                else
+                  continu=.false.
+                endif
+             enddo
+             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+             if (nouveau_traceurdef) then
+                write(lunout,*) 'C''est la nouvelle version de traceur.def'
+                tnom_0(iq)=tchaine(1:iiq-1)
+                tnom_transp(iq)=tchaine(iiq+1:)
+             else
+                write(lunout,*) 'C''est l''ancienne version de traceur.def'
+                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
+                tnom_0(iq)=tchaine
+                tnom_transp(iq) = 'air'
+             endif
+             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
+          ENDDO  !DO iq=1,nqo
+          CLOSE(90)
+! le module de chimie fournit les noms des traceurs et les schemas d'advection associes.
+! excepte pour ceux lus dans traceur.def
 #ifdef INCA
+       CALL init_transport( &
+            hadv_inca, &
+            vadv_inca, &
+            conv_flg_inca, &
+            pbl_flg_inca,  &
+            solsym_inca)
+       conv_flg(1+nqCO2:nbtr) = conv_flg_inca
+       pbl_flg(1+nqCO2:nbtr) = pbl_flg_inca
+       solsym(1+nqCO2:nbtr) = solsym_inca
+       IF (type_trac == 'inco') THEN
+          conv_flg(1:nqCO2) = 1
+          pbl_flg(1:nqCO2) = 1
+          solsym(1:nqCO2) = 'CO2'
+       ENDIF
+       CALL init_transport(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
+       ! DC passive CO2 tracer is at index 1: H2O was removed ; nqCO2/=0 in "inco" case only
+       conv_flg(1:nqCO2) = 1;   conv_flg(1+nqCO2:nbtr) = conv_flg_inca
+        pbl_flg(1:nqCO2) = 1;    pbl_flg(1+nqCO2:nbtr) =  pbl_flg_inca
+         solsym(1:nqCO2) = 'CO2'; solsym(1+nqCO2:nbtr) =   solsym_inca
 #endif
+!jyg<
+       DO iq = nqo+nqCO2+1, nqtrue
+          hadv(iq) = hadv_inca(iq-nqo-nqCO2)
+          vadv(iq) = vadv_inca(iq-nqo-nqCO2)
+          tnom_0(iq)=solsym_inca(iq-nqo-nqCO2)
+          tnom_transp(iq) = 'air'
+       itr = 0
+       DO iq = 1, nqtot
+          IF(iq > nqo+nqCO2) THEN
+             itr = itr+1
+             hadv  (iq) = hadv_inca  (itr)
+             vadv  (iq) = vadv_inca  (itr)
+             tnom_0(iq) = solsym_inca(itr)
+             tnom_transp(iq) = 'air'
+             CYCLE
+          END IF
+          write(*,*) 'infotrac 237: iq=',iq
+          ! CRisi: ajout du nom du fluide transporteur en restant retro-compatible
+          READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+          write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+          write(lunout,*) 'tchaine=',trim(tchaine)
+          write(*,*) 'infotrac 238: IOstatus=',IOstatus
+          if (IOstatus.ne.0) CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+          ! Y-a-t-il 1 ou 2 noms de traceurs separes par un espace ?
+          iiq = INDEX(tchaine(1:LEN_TRIM(tchaine)),' ')
+          nouveau_traceurdef=iiq/=0
+          write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+          if (nouveau_traceurdef) then
+             IF(iq==1) write(lunout,*) "Nouvelle version de traceur.def"
+             tnom_0     (iq) = tchaine(1:iiq-1)
+             tnom_transp(iq) = tchaine(iiq+1:)
+          else
+             IF(iq==1) write(lunout,*) "Nouvelle version de traceur.def: traceurs d'air uniquement"
+             tnom_0     (iq) = tchaine
+             tnom_transp(iq) = 'air'
+          endif
+          write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+          write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
        END DO
+       CLOSE(90)
     ENDIF ! (type_trac == 'inca' or 'inco')
 …
+!
     ALLOCATE(tracers(nqtot))
-    ALLOCATE(niadv(nqtot))
 !-----------------------------------------------------------------------
 …
     END DO
+!
+! Find vector keeping the correspodence between true and total tracers
+!
+    niadv(:)=0
+    iiq=0
+    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
+    WRITE(lunout,*) trim(modname),': iadv  name  long_name :'
     DO iq=1,nqtot
+       IF(tracers(iq)%iadv.GE.0) THEN
+          ! True tracer
+          iiq=iiq+1
+          niadv(iiq)=iq
+       ENDIF
+    END DO
+    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
+    WRITE(lunout,*) trim(modname),': iadv  niadv  name  long_name :'
+    DO iq=1,nqtot
+       WRITE(lunout,*) tracers(iq)%iadv,niadv(iq), ' ',trim(tracers(iq)%name),' ',trim(tracers(iq)%longName)
+       WRITE(lunout,*) tracers(iq)%iadv,' ',trim(tracers(iq)%name),' ',trim(tracers(iq)%longName)
     END DO
 …
 ! + verifier que tous les peres sont ecrits en premieres positions
     ALLOCATE(iqfils(nqtot,nqtot))
-    nqperes=0
     iqfils(:,:)=0
     tracers(:)%iqParent=0
 …
         ! ceci est un traceur pere
         WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
-        nqperes=nqperes+1
         tracers(iq)%iqParent=0
       else !if (tnom_transp(iq) == 'air') then
 …
       endif !if (tnom_transp(iq) == 'air') then
     enddo !DO iq=1,nqtot
     WRITE(lunout,*) 'infotrac: nqperes=',nqperes
+    WRITE(lunout,*) 'infotrac: nqGen0=',COUNT(tracers(:)%parent == 'air')
     WRITE(lunout,*) 'nqChilds=',tracers(:)%nqChilds
     WRITE(lunout,*) 'iqParent=',tracers(:)%iqParent
 …
     tracers(:)%isH2Ofamily = [(tracers(iq)%gen0Name(1:3) == 'H2O', iq=1, nqtot)]
 ! detecter quels sont les traceurs isotopiques parmi des traceurs
     call infotrac_isoinit(tnom_0,nqtrue)
 …
         write(lunout,*) 'infotrac 704: nqtottr,nqtot,nqo=',nqtottr,nqtot,nqo
         ! Rq: nqtottr n'est pas forcement egal a nbtr dans le cas ou new_iq /= nqtrue
+        ALLOCATE (itr_indice(nqtottr))
+        itr_indice(:)=0
+        itr=0
+        do iq=nqo+1, nqtot
+          if (tracers(iq)%iso_iName.eq.0) then
+            itr=itr+1
+            write(*,*) 'itr=',itr
+            itr_indice(itr)=iq
+          endif !if (tracers(iq)%iso_iName.eq.0) then
+        enddo
+        if (itr.ne.nqtottr) then
+        if (COUNT(tracers(:)%iso_iName == 0) /= nqtottr) &
             CALL abort_gcm('infotrac_init','pb dans le calcul de nqtottr',1)
-        endif
-        write(lunout,*) 'itr_indice=',itr_indice
 !    endif !if (ntraciso.gt.0) then

LMDZ6/trunk/libf/dyn3dmem/advtrac_loc.F90

-                      r4055
+                      r4056
+!
+! $Id$
+!
+c
+c
 #define DEBUG_IO
 #undef DEBUG_IO
+      SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
+     *                   p,  massem,q,teta,
+     *                   pk   )
+c     Auteur :  F. Hourdin
+c
+c     Modif. P. Le Van     (20/12/97)
+c            F. Codron     (10/99)
+c            D. Le Croller (07/2001)
+c            M.A Filiberti (04/2002)
+c
+      USE parallel_lmdz
+      USE Write_Field_loc
+      USE Write_Field
+      USE Bands
+      USE mod_hallo
+      USE Vampir
+      USE times
+      USE infotrac, ONLY: nqtot, tracers, ok_iso_verif
+      USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
+      USE advtrac_mod, ONLY: finmasse
+      USE comconst_mod, ONLY: dtvr
+      IMPLICIT NONE
+c
+      include "dimensions.h"
+      include "paramet.h"
+      include "comdissip.h"
+      include "comgeom2.h"
+      include "description.h"
+c-------------------------------------------------------------------
+c     Arguments
+c-------------------------------------------------------------------
+c     Ajout PPM
+c--------------------------------------------------------
+      REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
+c--------------------------------------------------------
+      INTEGER iapptrac
+      REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
+      REAL wg(ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
+      REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
+      REAL pk(ijb_u:ije_u,llm)
+c-------------------------------------------------------------
+c     Variables locales
+c-------------------------------------------------------------
+      REAL zdp(ijb_u:ije_u)
+      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
+      INTEGER,SAVE :: iadvtr=0
+c$OMP THREADPRIVATE(iadvtr)
+      INTEGER ij,l,iq,iiq
+      REAL zdpmin, zdpmax
+c----------------------------------------------------------
+c     Rajouts pour PPM
+c----------------------------------------------------------
+      INTEGER indice,n
+      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
+      REAL CFLmaxz,aaa,bbb ! CFL maximum
+      REAL psppm(iim,jjb_u:jje_u) ! pression  au sol
+      REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
+      REAL qppm(iim*jjnb_u,llm,nqtot)
+      REAL fluxwppm(iim,jjb_u:jje_u,llm)
+      REAL apppm(llmp1), bpppm(llmp1)
+      LOGICAL dum,fill
+      DATA fill/.true./
+      DATA dum/.true./
+      integer ijb,ije,ijbu,ijbv,ijeu,ijev,j, iadv
+      type(Request),SAVE :: testRequest
+SUBROUTINE advtrac_loc(pbarug, pbarvg, wg, p, massem, q, teta, pk)
+   !     Auteur :  F. Hourdin
+   !
+   !     Modif. P. Le Van     (20/12/97)
+   !            F. Codron     (10/99)
+   !            D. Le Croller (07/2001)
+   !            M.A Filiberti (04/2002)
+   !
+   USE infotrac,     ONLY: nqtot, tracers,ok_iso_verif
+   USE control_mod,  ONLY: iapp_tracvl, day_step, planet_type
+   USE comconst_mod, ONLY: dtvr
+   USE parallel_lmdz
+   USE Write_Field_loc
+   USE Write_Field
+   USE Bands
+   USE mod_hallo
+   USE Vampir
+   USE times
+   USE advtrac_mod, ONLY: finmasse
+   IMPLICIT NONE
+   !
+   include "dimensions.h"
+   include "paramet.h"
+   include "comdissip.h"
+   include "comgeom2.h"
+   include "description.h"
+!   include "iniprint.h"
+   !---------------------------------------------------------------------------
+   !     Arguments
+   !---------------------------------------------------------------------------
+   REAL, INTENT(IN) ::  pbarug(ijb_u:ije_u,llm)
+   REAL, INTENT(IN) ::  pbarvg(ijb_v:ije_v,llm)
+   REAL, INTENT(IN) ::      wg(ijb_u:ije_u,llm)
+   REAL, INTENT(IN) ::       p(ijb_u:ije_u,llmp1)
+   REAL, INTENT(IN) ::  massem(ijb_u:ije_u,llm)
+   REAL, INTENT(INOUT) ::    q(ijb_u:ije_u,llm,nqtot)
+   REAL, INTENT(IN) ::    teta(ijb_u:ije_u,llm)
+   REAL, INTENT(IN) ::      pk(ijb_u:ije_u,llm)
+   !---------------------------------------------------------------------------
+   !     Ajout PPM
+   !---------------------------------------------------------------------------
+   REAL :: massebx(ijb_u:ije_u,llm), masseby(ijb_v:ije_v,llm)
+   !---------------------------------------------------------------------------
+   !     Variables locales
+   !---------------------------------------------------------------------------
+   INTEGER :: ij, l, iq, iiq, iadv
+   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
+   REAL :: zdp(ijb_u:ije_u), zdpmin, zdpmax
+   INTEGER, SAVE :: iadvtr=0
+!$OMP THREADPRIVATE(iadvtr)
+   EXTERNAL  minmax
+   !---------------------------------------------------------------------------
+   !     Rajouts pour PPM
+   !---------------------------------------------------------------------------
+   INTEGER :: indice, n
+   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
+   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
+   REAL, DIMENSION(iim,jjb_u:jje_u,llm) :: unatppm, vnatppm, fluxwppm
+   REAL ::    qppm(iim*jjnb_u,llm,nqtot)
+   REAL ::   psppm(iim,jjb_u:jje_u)    ! pression  au sol
+   REAL, DIMENSION(llmp1) :: apppm, bpppm
+   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
+   INTEGER :: ijb, ije, ijbu, ijbv, ijeu, ijev, j
+   TYPE(Request),SAVE :: testRequest
 !$OMP THREADPRIVATE(testRequest)
+c  test sur l''eventuelle creation de valeurs negatives de la masse
+         ijb=ij_begin
+         ije=ij_end
+         if (pole_nord) ijb=ij_begin+iip1
+         if (pole_sud) ije=ij_end-iip1
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm-1
+            DO ij = ijb+1,ije
+              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
+     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
+     s                  +       wg(ij,l+1)  - wg(ij,l)
+            ENDDO
+c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
+            do ij=ijb,ije-iip1+1,iip1
+              zdp(ij)=zdp(ij+iip1-1)
+            enddo
+            DO ij = ijb,ije
+               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
+            ENDDO
+c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+c  ym ---> eventuellement a revoir
+            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
+            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
+            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
+     s        '   MAX:', zdpmax
+            ENDIF
+         ENDDO
+c$OMP END DO NOWAIT
+c-------------------------------------------------------------------
+c   Advection proprement dite (Modification Le Croller (07/2001)
+c-------------------------------------------------------------------
+c----------------------------------------------------
+c        Calcul des moyennes basées sur la masse
+c----------------------------------------------------
+cym      ----> Normalement, inutile pour les schémas classiques
+cym      ----> Revérifier lors de la parallélisation des autres schemas
+cym          call massbar_p(massem,massebx,masseby)
+! Test sur l'eventuelle creation de valeurs negatives de la masse
+   ijb = ij_begin; IF(pole_nord) ijb = ij_begin+iip1
+   ije = ij_end;   IF(pole_sud)  ije = ij_end-iip1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+   DO l=1,llm-1
+      DO ij = ijb+1,ije
+         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
+                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                 +     wg(ij,l+1)    -     wg(ij,l)
+      END DO
+! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
+!     CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+      DO ij = ijb,ije-iip1+1,iip1
+         zdp(ij)=zdp(ij+iip1-1)
+      END DO
+      DO ij = ijb,ije
+         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
+      END DO
+!     CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+! ym ---> eventuellement a revoir
+      CALL minmax( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
+      IF(MAX(ABS(zdpmin),ABS(zdpmax)) >0.5) &
+         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,'   MAX:', zdpmax
+   END DO
+!$OMP END DO NOWAIT
+   !---------------------------------------------------------------------------
+   !   Advection proprement dite (Modification Le Croller (07/2001)
+   !---------------------------------------------------------------------------
+   !---------------------------------------------------------------------------
+   !   Calcul des moyennes basees sur la masse
+   !---------------------------------------------------------------------------
+!ym   CALL massbar_p(massem,massebx,masseby)
+!ym   ----> Normalement, inutile pour les schémas classiques
+!ym   ----> Revérifier lors de la parallélisation des autres schemas
 #ifdef DEBUG_IO
+          CALL WriteField_u('massem',massem)
+          CALL WriteField_u('wg',wg)
+          CALL WriteField_u('pbarug',pbarug)
+          CALL WriteField_v('pbarvg',pbarvg)
+          CALL WriteField_u('p_tmp',p)
+          CALL WriteField_u('pk_tmp',pk)
+          CALL WriteField_u('teta_tmp',teta)
+          do j=1,nqtot
+            call WriteField_u('q_adv'//trim(int2str(j)),
+     .                q(:,:,j))
+          enddo
+   CALL WriteField_u('massem',massem)
+   CALL WriteField_u('wg',wg)
+   CALL WriteField_u('pbarug',pbarug)
+   CALL WriteField_v('pbarvg',pbarvg)
+   CALL WriteField_u('p_tmp',p)
+   CALL WriteField_u('pk_tmp',pk)
+   CALL WriteField_u('teta_tmp',teta)
+   DO iq=1,nqtot
+      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
+   END DO
 #endif
+!
+!       call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
+!
+!       call SendRequest(TestRequest)
+!c$OMP BARRIER
+!       call WaitRequest(TestRequest)
+c$OMP BARRIER
+!  CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
+!  CALL SendRequest(TestRequest)
+!!$OMP BARRIER
+!  CALL WaitRequest(TestRequest)
+!$OMP BARRIER
+          !write(*,*) 'advtrac 157: appel de vlspltgen_loc'
+          call vlspltgen_loc( q, 2., massem, wg,pbarug,pbarvg,dtvr,p,
+     *                        pk,teta )
+!  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
+   CALL vlspltgen_loc(q, tracers(:)%iadv, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
 #ifdef DEBUG_IO
+          do j=1,nqtot
+            call WriteField_u('q_adv'//trim(int2str(j)),
+     .                q(:,:,j))
+          enddo
+   DO iq = 1, nqtot
+      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
+   END DO
 #endif
+          GOTO 1234
+c-----------------------------------------------------------
+c     Appel des sous programmes d'advection
+c-----------------------------------------------------------
+      do iq=1,nqtot
+c        call clock(t_initial)
+        iadv = tracers(iq)%iadv
+        SELECT CASE(iadv)
+          CASE(0); CYCLE
+          CASE(10)
+c   ----------------------------------------------------------------
+c   Schema de Van Leer I MUSCL
+c   ----------------------------------------------------------------
+!LF            call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+          CASE(14)
+c   ----------------------------------------------------------------
+c   Schema "pseudo amont" + test sur humidite specifique
+C    pour la vapeur d'eau. F. Codron
+c   ----------------------------------------------------------------
+c
+cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
+!LF           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
+!LF     *                 pbarug,pbarvg,dtvr,p,pk,teta )
+          CASE(12)
+c   ----------------------------------------------------------------
+c   Schema de Frederic Hourdin
+c   ----------------------------------------------------------------
+          stop 'advtrac : schema non parallelise'
+c            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
+     s             dtvr,'n=',n
+           endif
+           do indice=1,n
+            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+           end do
+          CASE(13)
+          stop 'advtrac : schema non parallelise'
+c            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
+     s             dtvr,'n=',n
+           endif
+          do indice=1,n
+            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+          end do
+          CASE(20)
+c   ----------------------------------------------------------------
+c   Schema de pente SLOPES
+c   ----------------------------------------------------------------
+          stop 'advtrac : schema non parallelise'
+            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+          CASE(30)
+c   ----------------------------------------------------------------
+c   Schema de Prather
+c   ----------------------------------------------------------------
+          stop 'advtrac : schema non parallelise'
+c            Pas de temps adaptatif
+           call adaptdt(iadv,dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
+     s             dtvr,'n=',n
+           endif
+           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
+     s                     n,dtbon)
+          CASE(11,16,17,18)
+c   ----------------------------------------------------------------
+c   Schemas PPM Lin et Rood
+c   ----------------------------------------------------------------
+           stop 'advtrac : schema non parallelise'
+c        Test sur le flux horizontal
+c        Pas de temps adaptatif
+         call adaptdt(iadv,dtbon,n,pbarug,massem)
+         if (n.GT.1) then
+           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
+     s             dtvr,'n=',n
+         endif
+c        Test sur le flux vertical
+         CFLmaxz=0.
+         do l=2,llm
+           do ij=iip2,ip1jm
+            aaa=wg(ij,l)*dtvr/massem(ij,l)
+            CFLmaxz=max(CFLmaxz,aaa)
+            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+            CFLmaxz=max(CFLmaxz,bbb)
+           enddo
+         enddo
+         if (CFLmaxz.GE.1) then
+            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+         endif
+c-----------------------------------------------------------
+c        Ss-prg interface LMDZ.4->PPM3d
+c-----------------------------------------------------------
+          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
+     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
+     s                 unatppm,vnatppm,psppm)
+          do indice=1,n
+c---------------------------------------------------------------------
+c                         VL (version PPM) horiz. et PPM vert.
+c---------------------------------------------------------------------
+            SELECT CASE(iadv)
+              CASE(11)
+c                  Ss-prg PPM3d de Lin
+                  call ppm3d(1,qppm(1,1,iq),
+     s                       psppm,psppm,
+     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
+     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
+     s                       fill,dum,220.)
+              CASE(16)
+c----------------------------------------------------------------------
+c                           Monotonic PPM
+c----------------------------------------------------------------------
+c                  Ss-prg PPM3d de Lin
+                  call ppm3d(1,qppm(1,1,iq),
+     s                       psppm,psppm,
+     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
+     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
+     s                       fill,dum,220.)
+c---------------------------------------------------------------------
+              CASE(17)
+c---------------------------------------------------------------------
+c                           Semi Monotonic PPM
+c---------------------------------------------------------------------
+c                  Ss-prg PPM3d de Lin
+                  call ppm3d(1,qppm(1,1,iq),
+     s                       psppm,psppm,
+     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
+     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
+     s                       fill,dum,220.)
+c---------------------------------------------------------------------
+              CASE(18)
+c---------------------------------------------------------------------
+c                         Positive Definite PPM
+c---------------------------------------------------------------------
+c                  Ss-prg PPM3d de Lin
+                  call ppm3d(1,qppm(1,1,iq),
+     s                       psppm,psppm,
+     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
+     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
+     s                       fill,dum,220.)
+c---------------------------------------------------------------------
+            END SELECT
+            enddo
+c-----------------------------------------------------------------
+c               Ss-prg interface PPM3d-LMDZ.4
+c-----------------------------------------------------------------
+                  call interpost(q(1,1,iq),qppm(1,1,iq))
+        END SELECT
+c----------------------------------------------------------------------
+c-----------------------------------------------------------------
+c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
+c et Nord j=1
+c-----------------------------------------------------------------
+c                  call traceurpole(q(1,1,iq),massem)
+c calcul du temps cpu pour un schema donne
+       end DO
+  CONTINUE
+c$OMP BARRIER
+      if (planet_type=="earth") then
+   GOTO 1234
+   !-------------------------------------------------------------------------
+   !       Appel des sous programmes d'advection
+   !-------------------------------------------------------------------------
+   DO iq = 1, nqtot
+!     CALL clock(t_initial)
+      IF(tracers(iq)%parent /= 'air') CYCLE
+      iadv = tracers(iq)%iadv
+      !-----------------------------------------------------------------------
+      SELECT CASE(iadv)
+      !-----------------------------------------------------------------------
+         CASE(0); CYCLE
+         !--------------------------------------------------------------------
+         CASE(10)  !--- Schema de Van Leer I MUSCL
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
+!LF         CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+         !--------------------------------------------------------------------
+         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
+                   !--- pour la vapeur d'eau. F. Codron
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
+            STOP 'advtrac : appel a vlspltqs :schema non parallelise'
+!LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
+         !--------------------------------------------------------------------
+         CASE(12)  !--- Schema de Frederic Hourdin
+         !--------------------------------------------------------------------
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+            END DO
+         !--------------------------------------------------------------------
+         CASE(13)  !--- Pas de temps adaptatif
+         !--------------------------------------------------------------------
+            STOP 'advtrac : schema non parallelise'
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+            END DO
+         !--------------------------------------------------------------------
+         CASE(20)  !--- Schema de pente SLOPES
+         !--------------------------------------------------------------------
+            STOP 'advtrac : schema non parallelise'
+            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+         !--------------------------------------------------------------------
+         CASE(30)  !--- Schema de Prather
+         !--------------------------------------------------------------------
+            STOP 'advtrac : schema non parallelise'
+            ! Pas de temps adaptatif
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
+         !--------------------------------------------------------------------
+         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
+         !--------------------------------------------------------------------
+            STOP 'advtrac : schema non parallelise'
+            ! Test sur le flux horizontal
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            ! Test sur le flux vertical
+            CFLmaxz=0.
+            DO l=2,llm
+               DO ij=iip2,ip1jm
+                  aaa=wg(ij,l)*dtvr/massem(ij,l)
+                  CFLmaxz=max(CFLmaxz,aaa)
+                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                  CFLmaxz=max(CFLmaxz,bbb)
+               END DO
+            END DO
+            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+            !----------------------------------------------------------------
+            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
+            !----------------------------------------------------------------
+            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                 unatppm,vnatppm,psppm)
+            !----------------------------------------------------------------
+            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
+            !----------------------------------------------------------------
+               SELECT CASE(iadv)
+                  !----------------------------------------------------------
+                  CASE(11)
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(16) !--- Monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(17) !--- Semi monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(18) !--- Positive Definite PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+               END SELECT
+            !----------------------------------------------------------------
+            END DO
+            !----------------------------------------------------------------
+            !     Ss-prg interface PPM3d-LMDZ.4
+            !----------------------------------------------------------------
+            CALL interpost(q(1,1,iq),qppm(1,1,iq))
+      !----------------------------------------------------------------------
+      END SELECT
+      !----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 et Nord j=1
+      !----------------------------------------------------------------------
+      !  CALL traceurpole(q(1,1,iq),massem)
+      !--- Calcul du temps cpu pour un schema donne
+      !  CALL clock(t_final)
+      !ym  tps_cpu=t_final-t_initial
+      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
+   END DO
+CONTINUE
+!$OMP BARRIER
+   IF(planet_type=="earth") THEN
       ijb=ij_begin
       ije=ij_end
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+       DO l = 1, llm
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+      DO l = 1, llm
          DO ij = ijb, ije
+           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
+         ENDDO
+       ENDDO
+c$OMP END DO
+        ! CRisi: on passe nqtot et non nq
+       CALL qminimum_loc( q, nqtot, finmasse )
+      endif ! of if (planet_type=="earth")
+       RETURN
+       END
+            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
+         END DO
+      END DO
+!$OMP END DO
+      CALL qminimum_loc( q, nqtot, finmasse )
+   END IF ! of if (planet_type=="earth")
+END SUBROUTINE advtrac_loc

LMDZ6/trunk/libf/dyn3dmem/call_calfis_mod.F90

-                      r4050
+                      r4056
     TYPE(Request),SAVE :: Request_physic
 !$OMP THREADPRIVATE(Request_physic )
     INTEGER :: ijb,ije,l,j
+    INTEGER :: ijb,ije,l,iq
 …
     CALL WriteField_u('pfi',p)
     CALL WriteField_u('pkfi',pk)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q(:,:,iq))
     ENDDO
 #endif
 …
     CALL WriteField_u('pfi',p)
     CALL WriteField_u('pkfi',pk)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q(:,:,iq))
     ENDDO
 #endif
 …
     CALL Register_Hallo_u(dpfi,1,1,0,0,1,Request_physic)
     DO j=1,nqtot
       CALL Register_Hallo_u(dqfi(:,:,j),llm,1,0,0,1,Request_physic)
+    DO iq=1,nqtot
+      CALL Register_Hallo_u(dqfi(:,:,iq),llm,1,0,0,1,Request_physic)
     ENDDO
 …
     CALL WriteField_u('dtetafi',dtetafi)
     CALL WriteField_u('dpfi',dpfi)
     DO j=1,nqtot
       CALL WriteField_u('dqfi'//trim(int2str(j)),dqfi(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('dqfi'//trim(int2str(iq)),dqfi(:,:,iq))
     ENDDO
 #endif
 …
     CALL WriteField_u('tetafi',teta)
     CALL WriteField_u('psfi',ps)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q(:,:,iq))
     ENDDO
 #endif
 …
     CALL WriteField_u('tetafi',teta)
     CALL WriteField_u('psfi',ps)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q(:,:,iq))
     ENDDO
 #endif
 …
     CALL WriteField_u('tetafi',teta)
     CALL WriteField_u('psfi',ps)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q(:,:,iq))
     ENDDO
 #endif
 …
     CALL WriteField_u('tetafi',teta_dyn)
     CALL WriteField_u('psfi',ps_dyn)
     DO j=1,nqtot
       CALL WriteField_u('qfi'//trim(int2str(j)),q_dyn(:,:,j))
+    DO iq=1,nqtot
+      CALL WriteField_u('qfi'//trim(int2str(iq)),q_dyn(:,:,iq))
     ENDDO
 #endif

LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90

-                      r4052
+                      r4056
   USE filtreg_mod, ONLY: inifilr
+  USE infotrac,    ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, &
+                         iqiso, tracers, iso_indnum
+  USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
   use exner_milieu_m, only: exner_milieu
-  USE infotrac, ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, &
-                      iqiso, tracers, iso_indnum
-  USE control_mod, ONLY: day_step,planet_type
   USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v
 #ifdef CPP_IOIPSL
 …
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
   !   Author:    Frederic Hourdin      original: 15/01/93
 …
   ztot0      = 0.
   stot0      = 0.
   ang0       = 0.
+  ang0       = 0.
   if (llm == 1) then
 …
         ps_glo(:)=preff
      endif
      ! ground geopotential
      phis_glo(:)=0.
      CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
      if (pressure_exner) then
 …
      CALL massdair(p,masse_glo)
   ENDIF
   if (llm == 1) then
 …
         if (planet_type=="earth") then
            ! Earth: first two tracers will be water
            do iq=1,nqtot
               q(ijb_u:ije_u,:,iq)=0.
 !              IF(tracers(iq)%name == 'H2O'//Phases_sep//'g') q(ijb_u:ije_u,:,iq)=1.e-10
+!              IF(tracers(iq)%name == 'H2O'//phases_sep//'g') q(ijb_u:ije_u,:,iq)=1.e-10
 !              IF(tracers(iq)%name == 'H2O'//phases_sep//'l') q(ijb_u:ije_u,:,iq)=1.e-15
               if (tracers(iq)%name == 'H2Ov') q(ijb_u:ije_u,:,iq)=1.e-10
               if (tracers(iq)%name == 'H2Ol') q(ijb_u:ije_u,:,iq)=1.e-15
+              IF(tracers(iq)%name == 'H2Ov') q(ijb_u:ije_u,:,iq)=1.e-10
+              IF(tracers(iq)%name == 'H2Ol') q(ijb_u:ije_u,:,iq)=1.e-15
               ! CRisi: init des isotopes
 …
               iName = tracers(iq)%iso_iName
               if (.NOT.ok_isotopes .OR. iName <= 0) CYCLE
               iZone = tracers(iq)%iso_iZone
               iPhase= tracers(iq)%iso_iPhase
+              iZone    = tracers(iq)%iso_iZone
+              iPhase   = tracers(iq)%iso_iPhase
               iqParent = tracers(iq)%iqParent
               if (iZone == 0) q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) &
 …
         endif ! of if (planet_type=="earth")
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
+        endif !if (ok_iso_verif) then
+        if (ok_iso_verif) call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
         ! add random perturbation to temperature

LMDZ6/trunk/libf/dyn3dmem/leapfrog_loc.F

-                      r3947
+                      r4056
         call WriteField_u('pkf',pkf)
         call WriteField_u('phis',phis)
         do j=1,nqtot
           call WriteField_u('q'//trim(int2str(j)),
      .                q(:,:,j))
+        do iq=1,nqtot
+          call WriteField_u('q'//trim(int2str(iq)),
+     .                q(:,:,iq))
         enddo
       endif

LMDZ6/trunk/libf/dyn3dmem/vlspltgen_loc.F

-                      r4052
+                      r4056
       USE VAMPIR
       ! CRisi: on rajoute variables utiles d'infotrac
       USE infotrac, ONLY : nqtot,nqperes, tracers,ok_iso_verif
+      USE infotrac, ONLY : nqtot, tracers,ok_iso_verif
       USE vlspltgen_mod
       USE comconst_mod, ONLY: cpp
 …
 c$OMP BARRIER
+!      DO iq=1,nqtot
+      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
+       !write(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
+      DO iq=1,nqtot
+        ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
+        IF(tracers(iq)%parent /= 'air') CYCLE
+        !write(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
 #ifdef DEBUG_IO
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
+        CALL WriteField_u('zq',zq(:,:,iq))
+        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
         SELECT CASE(tracers(iq)%iadv)
 …
         END SELECT
       enddo !DO iq=1,nqperes
+      enddo !DO iq=1,nqtot
 …
       endif !if (ok_iso_verif) then
+      do iq=1,nqperes
+      do iq=1,nqtot
+        IF(tracers(iq)%parent /= 'air') CYCLE
         !write(*,*) 'vlspltgen 279: iq=',iq
 …
       if (ok_iso_verif) then
            call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
-      endif !if (ok_iso_verif) then
-      if (ok_iso_verif) then
            ijb=ij_begin-2*iip1
            ije=ij_end+2*iip1
 …
       endif !if (ok_iso_verif) then
+      do iq=1,nqperes
+      do iq = 1, nqtot
+       IF(tracers(iq)%parent /= 'air') CYCLE
        !write(*,*) 'vlspltgen 321: iq=',iq
 #ifdef DEBUG_IO
 …
       endif !if (ok_iso_verif) then
+      do iq=1,nqperes
+      do iq = 1, nqtot
+       IF(tracers(iq)%parent /= 'air') CYCLE
       !write(*,*) 'vlspltgen 349: iq=',iq
 #ifdef DEBUG_IO
 …
 c$OMP BARRIER
+      do iq=1,nqperes
+      do iq=1,nqtot
+        IF(tracers(iq)%parent /= 'air') CYCLE
       !write(*,*) 'vlspltgen 409: iq=',iq
 …
       endif !if (ok_iso_verif) then
+      do iq=1,nqperes
+      do iq=1,nqtot
+        IF(tracers(iq)%parent /= 'air') CYCLE
       !write(*,*) 'vlspltgen 449: iq=',iq
 #ifdef DEBUG_IO
 …
         END SELECT
        enddo !do iq=1,nqperes
+       enddo !do iq=1,nqtot
       if (ok_iso_verif) then
 …
       endif !if (ok_iso_verif) then
+      do iq=1,nqperes
+      do iq=1,nqtot
+        IF(tracers(iq)%parent /= 'air') CYCLE
       !write(*,*) 'vlspltgen 477: iq=',iq
 #ifdef DEBUG_IO
 …
         END SELECT
        enddo !do iq=1,nqperes
+       enddo !do iq=1,nqtot
       !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'

LMDZ6/trunk/libf/dynphy_lonlat/calfis.F

-                      r4046
+                      r4056
 c    Auteur :  P. Le Van, F. Hourdin
 c   .........
       USE infotrac, ONLY: nqtot, niadv, tracers
+      USE infotrac, ONLY: nqtot, tracers
       USE control_mod, ONLY: planet_type, nsplit_phys
 #ifdef CPP_PHYS
 …
 c    -----------------
       INTEGER i,j,l,ig0,ig,iq,iiq
+      INTEGER i,j,l,ig0,ig,iq,itr
       REAL zpsrf(ngridmx)
       REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
 …
 c   ---------------
+c
+      itr=0
       DO iq=1,nqtot
+          iiq=niadv(iq)
+         IF(.NOT.tracers(iq)%isAdvected) CYCLE
+         itr = itr + 1
          DO l=1,llm
             zqfi(1,l,iq) = pq(1,1,l,iiq)
             ig0          = 2
+            zqfi(1,l,itr) = pq(1,1,l,iq)
+            ig0           = 2
             DO j=2,jjm
                DO i = 1, iim
                   zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
+                  zqfi(ig0,l,itr) = pq(i,j,l,iq)
                   ig0             = ig0 + 1
                ENDDO
             ENDDO
             zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
+            zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
          ENDDO
       ENDDO
 …
          lafin_split=lafin.and.isplit==nsplit_phys
+!      if (planet_type=="earth") then
         CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
      &                   debut_split,lafin_split,
 …
      &                   flxwfi,pducov,
      &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
-!      if (planet_type=="earth") then
+!
-!         CALL physiq (ngridmx,
-!     .             llm,
-!     .             debut_split,
-!     .             lafin_split,
-!     .             jD_cur,
-!     .             jH_cur_split,
-!     .             zdt_split,
-!     .             zplev,
-!     .             zplay,
-!     .             zphi,
-!     .             zphis,
-!     .             presnivs,
-!     .             zufi,
-!     .             zvfi, zrfi,
-!     .             ztfi,
-!     .             zqfi,
-!     .             flxwfi,
-!     .             zdufi,
-!     .             zdvfi,
-!     .             zdtfi,
-!     .             zdqfi,
-!     .             zdpsrf,
-!     .             pducov)
+!
 !      else if ( planet_type=="generic" ) then
 …
       pdqfi(:,:,:,:)=0.
+C
+      itr = 0
       DO iq=1,nqtot
+         iiq=niadv(iq)
+         IF(.NOT.tracers(iq)%isAdvected) CYCLE
+         itr = itr + 1
          DO l=1,llm
             DO i=1,iip1
                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
                pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
+               pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
+               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
             ENDDO
             DO j=2,jjm
                ig0=1+(j-2)*iim
                DO i=1,iim
                   pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
+                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
                ENDDO
                pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
+               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
             ENDDO
          ENDDO

LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F

-                      r4046
+                      r4056
       USE Times
 #endif
       USE infotrac, ONLY: nqtot, niadv, tracers
+      USE infotrac, ONLY: nqtot, tracers
       USE control_mod, ONLY: planet_type, nsplit_phys
 #ifdef CPP_PHYS
 …
 c    -----------------
       INTEGER i,j,l,ig0,ig,iq,iiq
+      INTEGER i,j,l,ig0,ig,iq,itr
       REAL,ALLOCATABLE,SAVE :: zpsrf(:)
       REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
 …
+c
+      itr = 0
       DO iq=1,nqtot
+         iiq=niadv(iq)
+         IF(.NOT.tracers(iq)%isAdvected) CYCLE
+         itr = itr + 1
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
 …
              i=index_i(ig0)
              j=index_j(ig0)
              zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
+             zqfi(ig0,l,itr)  = pq(i,j,l,iq)
            enddo
          ENDDO
 …
+C
 !cdir NODEP
+      itr = 0
       DO iq=1,nqtot
+         iiq=niadv(iq)
+         IF(.NOT.tracers(iq)%isAdvected) CYCLE
+         itr = itr + 1
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
 …
               i=index_i(ig0)
               j=index_j(ig0)
               pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
               if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
+              pdqfi(i,j,l,iq) = zdqfi(ig0,l,itr)
+              if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,itr)
             ENDDO
             IF (is_north_pole_dyn) then
               DO i=1,iip1
                 pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
+                pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
               ENDDO
             ENDIF
 …
             IF (is_south_pole_dyn) then
               DO i=1,iip1
                 pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
+                pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,itr)
               ENDDO
             ENDIF

LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

-                      r4050
+                      r4056
   USE vertical_layers_mod, ONLY : init_vertical_layers
   USE infotrac, ONLY: nqtot,nqo,nbtr,nqCO2,tracers,type_trac,&
                       niadv,conv_flg,pbl_flg,solsym,&
+                      conv_flg,pbl_flg,solsym,&
                       ok_isotopes,ok_iso_verif,ok_isotrac,&
                       ok_init_iso,niso_possibles,tnat,&
                       alpha_ideal,use_iso,iqiso,iso_indnum,&
                       indnum_fn_num,index_trac,&
                       niso,ntraceurs_zone,ntraciso,nqtottr,itr_indice
+                      niso,ntraceurs_zone,ntraciso,nqtottr
 #ifdef CPP_StratAer
   USE infotrac, ONLY: nbtr_bin, nbtr_sulgas, id_OCS_strat, &
+  USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, id_OCS_strat, &
                       id_SO2_strat, id_H2SO4_strat, id_BIN01_strat
 …
   ! Initialize tracer names, numbers, etc. for physics
   CALL init_infotrac_phy(nqtot,nqo,nbtr,nqtottr,nqCO2,tracers,type_trac,&
                          niadv,conv_flg,pbl_flg,solsym,&
+                         conv_flg,pbl_flg,solsym,&
                          ok_isotopes,ok_iso_verif,ok_isotrac,&
                          ok_init_iso,niso_possibles,tnat,&
                          alpha_ideal,use_iso,iqiso,iso_indnum,&
                          indnum_fn_num,index_trac,&
+                         niso,ntraceurs_zone,ntraciso,itr_indice &
+#ifdef CPP_StratAer
+                         ,nbtr_bin,nbtr_sulgas&
+                         ,id_OCS_strat,id_SO2_strat,id_H2SO4_strat,id_BIN01_strat&
+#endif
+                         )
+                         niso,ntraceurs_zone,ntraciso)
   ! Initializations for Reprobus

LMDZ6/trunk/libf/phylmd/Dust/phys_output_write_spl_mod.F90

-                      r3805
+                      r4056
     USE pbl_surface_mod, ONLY: snow
     USE indice_sol_mod, ONLY: nbsrf
     USE infotrac, ONLY: nqtot, nqo, nbtr, type_trac
+    USE infotrac, ONLY: nqtot, nbtr, type_trac
     USE geometry_mod, ONLY: cell_area
     USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, landice_opt
 …
     INTEGER, PARAMETER :: jjmp1=jjm+1-1/jjm
     INTEGER :: itau_w
     INTEGER :: i, iinit, iinitend=1, iff, iq, nsrf, k, ll, naero
+    INTEGER :: i, iinit, iinitend=1, iff, iq, itr, nsrf, k, ll, naero
     REAL, DIMENSION (klon) :: zx_tmp_fi2d
     REAL, DIMENSION (klon,klev) :: zx_tmp_fi3d, zpt_conv
 …
 #endif
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       IF (nqtot.GE.nqo+1) THEN
+         !AS: type_trac = 'lmdz' par defaut dans libf/dyn3d/conf_gcm.F90
+         !Changé par inca, repr(obus), coag(ulation), co2i(nteractif), PAS par SPLA
+         !Cet "if" est donc inutile : IF (type_trac == 'lmdz' .OR. type_trac == 'repr') THEN
+           DO iq=nqo+1,nqtot
+             CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_lessi_impa(iq-nqo),d_tr_lessi_impa(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_lessi_nucl(iq-nqo),d_tr_lessi_nucl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_insc(iq-nqo),d_tr_insc(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_bcscav(iq-nqo),d_tr_bcscav(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_evapls(iq-nqo),d_tr_evapls(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_ls(iq-nqo),d_tr_ls(:,:,iq-nqo))
+!            CALL histwrite_phy(o_dtr_dyn(iq-nqo),d_tr_dyn(:,:,iq-nqo))
+!            CALL histwrite_phy(o_dtr_cl(iq-nqo),d_tr_cl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_trsp(iq-nqo),d_tr_trsp(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_sscav(iq-nqo),d_tr_sscav(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_sat(iq-nqo),d_tr_sat(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo))
+           itr = 0
+           DO iq = 1, nqtot
+             IF(tracers(iq)%isH2Ofamily) CYCLE
+             itr = itr+1
+             CALL histwrite_phy(o_trac(itr), tr_seri(:,:,itr))
+             CALL histwrite_phy(o_dtr_vdf(itr),d_tr_cl(:,:,itr))
+             CALL histwrite_phy(o_dtr_the(itr),d_tr_th(:,:,itr))
+             CALL histwrite_phy(o_dtr_con(itr),d_tr_cv(:,:,itr))
+             CALL histwrite_phy(o_dtr_lessi_impa(itr),d_tr_lessi_impa(:,:,itr))
+             CALL histwrite_phy(o_dtr_lessi_nucl(itr),d_tr_lessi_nucl(:,:,itr))
+             CALL histwrite_phy(o_dtr_insc(itr),d_tr_insc(:,:,itr))
+             CALL histwrite_phy(o_dtr_bcscav(itr),d_tr_bcscav(:,:,itr))
+             CALL histwrite_phy(o_dtr_evapls(itr),d_tr_evapls(:,:,itr))
+             CALL histwrite_phy(o_dtr_ls(itr),d_tr_ls(:,:,itr))
+!            CALL histwrite_phy(o_dtr_dyn(itr),d_tr_dyn(:,:,itr))
+!            CALL histwrite_phy(o_dtr_cl(itr),d_tr_cl(:,:,itr))
+             CALL histwrite_phy(o_dtr_trsp(itr),d_tr_trsp(:,:,itr))
+             CALL histwrite_phy(o_dtr_sscav(itr),d_tr_sscav(:,:,itr))
+             CALL histwrite_phy(o_dtr_sat(itr),d_tr_sat(:,:,itr))
+             CALL histwrite_phy(o_dtr_uscav(itr),d_tr_uscav(:,:,itr))
              zx_tmp_fi2d=0.
              IF (vars_defined) THEN
                 DO k=1,klev
                    zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo)
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,itr)
                 ENDDO
              ENDIF
              CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
+             CALL histwrite_phy(o_trac_cum(itr), zx_tmp_fi2d)
            ENDDO
-         !ENDIF
-       ENDIF
        IF (.NOT.vars_defined) THEN

LMDZ6/trunk/libf/phylmd/Dust/phytracr_spl_mod.F90

-                      r4046
+                      r4056
       REAL, intent(in) ::  rlon(klon)       ! longitudes pour chaque point
+!
       INTEGER i, k, it, j, ig
+      INTEGER i, k, iq, itr, j, ig
+!
 ! DEFINITION OF DIAGNOSTIC VARIABLES
 …
 #ifdef IOPHYS_DUST
+      do it=1,nbtr
+         write(str2,'(i2.2)') it
+         call iophys_ecrit('TRA'//str2,klev,'SOURCE','',tr_seri(:,:,it))
+      itr = 0
+      DO iq = 1, nqtot
+         IF(tracers(iq)%isH2Ofamily) CYCLE
+         itr = itr+1
+         write(str2,'(i2.2)') itr
+         call iophys_ecrit('TRA'//str2,klev,'SOURCE','',tr_seri(:,:,itr))
       enddo
 #endif
 …
         id_codu=-1
         id_scdu=-1
+       !print *,nbtr
+       do it=1,nbtr
+        print *, it, tracers(it+nqo)%name
+        SELECT CASE(tracers(it+nqo)%name)
+          CASE('PREC'); id_prec=it
+          CASE('FINE'); id_fine=it
+          CASE('COSS'); id_coss=it
+          CASE('CODU'); id_codu=it
+          CASE('SCDU'); id_scdu=it
+        END SELECT
+       enddo
+       ! check consistency with dust emission scheme:
+       if (ok_chimeredust) then
+        itr = 0
+        do iq=1,nqtot
+          IF(tracers(iq)%isH2Ofamily) CYCLE
+          itr = itr+1
+          print *, itr, TRIM(tracers(iq)%name)
+          SELECT CASE(tracers(iq)%name)
+            CASE('PREC'); id_prec=itr
+            CASE('FINE'); id_fine=itr
+            CASE('COSS'); id_coss=itr
+            CASE('CODU'); id_codu=itr
+            CASE('SCDU'); id_scdu=itr
+          END SELECT
+        enddo
+        ! check consistency with dust emission scheme:
+        if (ok_chimeredust) then
           if (.not.( id_scdu>0 .and. id_codu>0 .and. id_fine>0)) then
              call abort_gcm('phytracr_mod', 'pb in ok_chimdust 0',1)
           endif
        else
+        else
           if (id_scdu>0) then
        call abort_gcm('phytracr_mod', 'pb in ok_chimdust 1 SCDU',1)
 …
 ! JE before put in zero
       IF (lminmax) THEN
         DO it=1,nbtr
         CALL checknanqfi(tr_seri(:,:,it),qmin,qmax,'nan init phytracr')
         ENDDO
         DO it=1,nbtr
         CALL minmaxqfi2(tr_seri(:,:,it),qmin,qmax,'minmax init phytracr')
+        DO itr=1,nbtr
+        CALL checknanqfi(tr_seri(:,:,itr),qmin,qmax,'nan init phytracr')
+        ENDDO
+        DO itr=1,nbtr
+        CALL minmaxqfi2(tr_seri(:,:,itr),qmin,qmax,'minmax init phytracr')
         ENDDO
         CALL minmaxsource(source_tr,qmin,qmax,'maxsource init phytracr')

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_5wv_rrtm.F90

-                      r4046
+                      r4056
   USE DIMPHY
   USE aero_mod
   USE infotrac_phy
+  USE infotrac_phy, ONLY: nqtot, nbtr, tracers
   USE phys_local_var_mod, ONLY: od550aer,od865aer,ec550aer,od550lt1aer
+  !
 …
   LOGICAL :: soluble
   INTEGER :: i, k, m, itr, irh, aerindex
+  INTEGER :: i, k, m, iq, itr, irh, aerindex
   INTEGER :: spsol, spinsol, la
   INTEGER :: RH_num(klon,klev)
 …
   ENDDO
+  DO itr=1,nbtr    !--loop over tracers
+    SELECT CASE(tracers(itr+nqo)%name)
+  itr = 0
+  DO iq = 1, nqtot
+    IF(tracers(iq)%isH2Ofamily) CYCLE
+    itr = itr+1
+    SELECT CASE(tracers(iq)%name)
       CASE('PREC'); CYCLE                                  !--precursor
       CASE('FINE'); soluble=.TRUE.;  spsol=1; aerindex=1   !--fine mode accumulation mode
 …
       CASE('CODU'); soluble=.FALSE.; spsol=1; aerindex=3   !--coarse mode dust
       CASE('SCDU'); soluble=.FALSE.; spsol=2; aerindex=4   !--super coarse mode dust
       CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(itr+nqo)%name,1)
+      CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(iq)%name,1)
     END SELECT

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_6bands_rrtm.F90

-                      r4046
+                      r4056
   USE dimphy
   USE aero_mod
   USE infotrac_phy
+  USE infotrac_phy, ONLY: nqtot, nbtr, tracers
   USE phys_local_var_mod, ONLY: abs550aer
 …
+  !
   LOGICAL :: soluble
   INTEGER :: i, k, irh, itr, inu
+  INTEGER :: i, k, irh, iq, itr, inu
   INTEGER :: aerindex, spsol, spinsol
   INTEGER :: RH_num(klon,klev)
 …
   cg_ae(:,:,:,:)=0.
+  DO itr=1,nbtr    !--loop over tracers
+    SELECT CASE(tracers(itr+nqo)%name)
+  itr = 0
+  DO iq = 1, nqtot
+    IF(tracers(iq)%isH2Ofamily) CYCLE
+    itr = itr+1
+    SELECT CASE(tracers(iq)%name)
       CASE('PREC'); CYCLE                                  !--precursor
       CASE('FINE'); soluble=.TRUE.;  spsol=1; aerindex=1   !--fine mode accumulation mode
 …
       CASE('CODU'); soluble=.FALSE.; spsol=1; aerindex=3   !--coarse mode dust
       CASE('SCDU'); soluble=.FALSE.; spsol=2; aerindex=4   !--super coarse mode dust
       CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(itr+nqo)%name,1)
+      CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(iq)%name,1)
     END SELECT

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_lw_rrtm.F90

-                      r4046
+                      r4056
   USE dimphy
   USE aero_mod
   USE infotrac_phy
+  USE infotrac_phy, ONLY: nqtot, nbtr, tracers
   USE phys_state_var_mod, ONLY : tau_aero_lw_rrtm
   USE YOERAD, ONLY : NLW
 …
   INTEGER, PARAMETER :: naero=naero_soluble+naero_insoluble
+  !
   INTEGER inu, itr, spinsol
+  INTEGER inu, itr, iq, spinsol
   CHARACTER*20 modname
+  !
 …
     tau_aero_lw_rrtm = 0.0
+    !
+    DO itr=1,nbtr
+      SELECT CASE(tracers(itr+nqo)%name)
+        CASE('PREC','FINE''COSS'); CYCLE                   !--precursor or fine/coarde accumulation mode
+    itr = 0
+    DO iq = 1, nqtot
+      IF(tracers(iq)%isH2Ofamily) CYCLE
+      itr = itr+1
+      SELECT CASE(tracers(iq)%name)
+        CASE('PREC','FINE','COSS'); CYCLE                  !--precursor or fine/coarde accumulation mode
         CASE('CODU'); spinsol=1                            !--coarse mode dust
         CASE('SCDU'); spinsol=2                            !--super coarse mode dust
         CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(itr+nqo)%name,1)
+        CASE DEFAULT; CALL abort_physic(modname,'I cannot do aerosol optics for '//tracers(iq)%name,1)
       END SELECT
+      !

LMDZ6/trunk/libf/phylmd/Dust/splaerosol_optic_rrtm.F90

-                      r4046
+                      r4056
   USE dimphy
   USE aero_mod
   USE infotrac_phy
+  USE infotrac_phy, ONLY: nbtr, nqtot, tracers
   USE YOMCST, ONLY: RD, RG
 …
   REAL, DIMENSION(klon,klev,nwave,naero_tot), INTENT(OUT)  :: tau3d_aero
   INTEGER i, k, itr
+  INTEGER i, k, iq, itr
   REAL, DIMENSION(klon,klev) :: zdm, zdh
   REAL zrho, pdel
 …
   mass_solu_aero_pi(:,:) = 0.0
+  !
+  DO itr=1,nbtr
+    IF (tracers(itr+nqo)%name=='FINE') THEN
+  itr = 0
+  DO iq = 1, nqtot
+    IF(tracers(iq)%isH2Ofamily) CYCLE
+    itr = itr+1
+    IF(tracers(iq)%name/='FINE') THEN
       mass_solu_aero(:,:)    = tr_seri(:,:,itr)
       mass_solu_aero_pi(:,:) = tr_seri(:,:,itr)

LMDZ6/trunk/libf/phylmd/infotrac_phy.F90

-                      r4052
+                      r4056
 !$OMP THREADPRIVATE(nqtottr)
 ! ThL : number of CO2 tracers                   ModThL
+! ThL : number of CO2 tracers   ModThL
   INTEGER, SAVE :: nqCO2
 !$OMP THREADPRIVATE(nqCO2)
 #ifdef CPP_StratAer
+! nbtr_bin: number of aerosol bins for StratAer model
+! nbtr_sulgas: number of sulfur gases for StratAer model
+  INTEGER, SAVE :: nbtr_bin, nbtr_sulgas
+!$OMP THREADPRIVATE(nbtr_bin,nbtr_sulgas)
+  INTEGER, SAVE :: id_OCS_strat, id_SO2_strat, id_H2SO4_strat, id_BIN01_strat
+!$OMP THREADPRIVATE(id_OCS_strat, id_SO2_strat, id_H2SO4_strat, id_BIN01_strat)
+  !=== SPECIFIC TO STRATOSPHERIC AEROSOLS (CK/OB)
+  INTEGER, SAVE ::  nbtr_bin, nbtr_sulgas         !--- number of aerosols bins and sulfur gases for StratAer model
+!$OMP THREADPRIVATE(nbtr_bin, nbtr_sulgas)
+  INTEGER, SAVE ::  id_OCS_strat, id_SO2_strat, id_H2SO4_strat, id_BIN01_strat, id_TEST_strat
+!$OMP THREADPRIVATE(id_OCS_strat, id_SO2_strat, id_H2SO4_strat, id_BIN01_strat, id_TEST_strat)
 #endif
-! CRisi: nb traceurs pères= directement advectés par l'air
-  INTEGER, SAVE :: nqperes
-!$OMP THREADPRIVATE(nqperes)
 ! Tracers parameters
   TYPE(trac_type), TARGET, ALLOCATABLE, SAVE :: tracers(:)
 !$OMP THREADPRIVATE(tracers)
-! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the
-!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
-  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
-!$OMP THREADPRIVATE(niadv)
 ! conv_flg(it)=0 : convection desactivated for tracer number it
 …
 !$OMP THREADPRIVATE(niso,ntraceurs_zone,ntraciso)
-    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  itr_indice ! numéro iq entre 1 et nqtot qui correspond au traceur itr entre 1 et nqtottr
-!$OMP THREADPRIVATE(itr_indice)
 CONTAINS
   SUBROUTINE init_infotrac_phy(nqtot_,nqo_,nbtr_,nqtottr_,nqCO2_,tracers_,type_trac_,&
                                niadv_,conv_flg_,pbl_flg_,solsym_,&
+                               conv_flg_,pbl_flg_,solsym_,&
                                ok_isotopes_,ok_iso_verif_,ok_isotrac_,&
                                ok_init_iso_,niso_possibles_,tnat_,&
                                alpha_ideal_,use_iso_,iqiso_,iso_indnum_,&
                                indnum_fn_num_,index_trac_,&
+                               niso_,ntraceurs_zone_,ntraciso_,itr_indice_&
+#ifdef CPP_StratAer
+                               ,nbtr_bin_,nbtr_sulgas_&
+                               ,id_OCS_strat_,id_SO2_strat_,id_H2SO4_strat_,id_BIN01_strat_&
+#endif
+                               )
+                               niso_,ntraceurs_zone_,ntraciso_)
     ! transfer information on tracers from dynamics to physics
 …
     INTEGER,INTENT(IN) :: nqtottr_
     INTEGER,INTENT(IN) :: nqCO2_
-#ifdef CPP_StratAer
-    INTEGER,INTENT(IN) :: nbtr_bin_
-    INTEGER,INTENT(IN) :: nbtr_sulgas_
-    INTEGER,INTENT(IN) :: id_OCS_strat_
-    INTEGER,INTENT(IN) :: id_SO2_strat_
-    INTEGER,INTENT(IN) :: id_H2SO4_strat_
-    INTEGER,INTENT(IN) :: id_BIN01_strat_
-#endif
     TYPE(trac_type), INTENT(IN) :: tracers_(nqtot_) ! tracers descriptors
     CHARACTER(len=*),INTENT(IN) :: type_trac_
-    INTEGER,INTENT(IN) :: niadv_ (nqtot_) ! equivalent dyn / physique
     INTEGER,INTENT(IN) :: conv_flg_(nbtr_)
     INTEGER,INTENT(IN) :: pbl_flg_(nbtr_)
 …
     INTEGER,INTENT(IN) :: ntraceurs_zone_
     INTEGER,INTENT(IN) :: ntraciso_
+    INTEGER,INTENT(IN) :: itr_indice_(nqtottr_)
+    CHARACTER(LEN=30) :: modname="init_infotrac_phy"
+    INTEGER :: iq, itr
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: tnames(:)
+    CHARACTER(LEN=maxlen) :: modname="init_infotrac_phy"
     nqtot=nqtot_
 …
     ALLOCATE(tracers(nqtot)); tracers(:) = tracers_(:)
 #ifdef CPP_StratAer
+    nbtr_bin=nbtr_bin_
+    nbtr_sulgas=nbtr_sulgas_
+    id_OCS_strat=id_OCS_strat_
+    id_SO2_strat=id_SO2_strat_
+    id_H2SO4_strat=id_H2SO4_strat_
+    id_BIN01_strat=id_BIN01_strat_
+    IF (type_trac == 'coag') THEN
+      nbtr_bin    = COUNT([(tracers(iq)%name(1:3)=='BIN', iq=1, nqtot)])
+      nbtr_sulgas = COUNT([(tracers(iq)%name(1:3)=='GAS', iq=1, nqtot)])
+      tnames = PACK(tracers(:)%name, MASK=.NOT.tracers(:)%isH2Ofamily)
+      id_BIN01_strat = strIdx(tnames, 'BIN01'   )
+      id_OCS_strat   = strIdx(tnames, 'GASOSC'  )
+      id_SO2_strat   = strIdx(tnames, 'GASSO2'  )
+      id_H2SO4_strat = strIdx(tnames, 'GASH2SO4')
+      id_TEST_strat  = strIdx(tnames, 'GASTEST' )
+      WRITE(lunout,*)'nbtr_bin       =', nbtr_bin
+      WRITE(lunout,*)'nbtr_sulgas    =', nbtr_sulgas
+      WRITE(lunout,*)'id_BIN01_strat =', id_BIN01_strat
+      WRITE(lunout,*)'id_OCS_strat   =',   id_OCS_strat
+      WRITE(lunout,*)'id_SO2_strat   =',   id_SO2_strat
+      WRITE(lunout,*)'id_H2SO4_strat =', id_H2SO4_strat
+      WRITE(lunout,*)'id_TEST_strat  =',  id_TEST_strat
+    END IF
 #endif
     type_trac = type_trac_
-    ALLOCATE(niadv(nqtot))
-    niadv(:)=niadv_(:)
     ALLOCATE(conv_flg(nbtr))
     conv_flg(:)=conv_flg_(:)
 …
     ENDIF ! of IF(ok_isotopes)
+    ALLOCATE(itr_indice(nqtottr))
+    itr_indice(:)=itr_indice_(:)
+    WRITE(*,*) 'infotrac_phy 207: nqtottr=',nqtottr
+    WRITE(*,*) 'ntraciso,niso=',ntraciso,niso
+#ifdef ISOVERIF
+    ! DC: the "1" will be replaced by iH2O (H2O isotopes group index)
+    WRITE(*,*) 'iso_iName=',PACK(tracers(:)%iso_iName, MASK=tracers(:)%iso_iGroup==1)
+#endif
   END SUBROUTINE init_infotrac_phy

LMDZ6/trunk/libf/phylmd/init_be.F90

r2351	r4056
5	5
6	6	USE dimphy
7		~~USE infotrac_phy, ONLY : nbtr~~
8	7	USE indice_sol_mod
9	8	USE geometry_mod, ONLY : longitude, latitude

LMDZ6/trunk/libf/phylmd/phyetat0.F90

-                      r4046
+                      r4056
        ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, ratqs_inter
 !FC
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
   USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
   USE infotrac_phy, only: nbtr, nqo, type_trac, tracers, niadv
   USE traclmdz_mod,    ONLY : traclmdz_from_restart
   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
   USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
   USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
+  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
+  USE iostart,          ONLY: close_startphy, get_field, get_var, open_startphy
+  USE infotrac_phy,     ONLY: nqtot, nbtr, type_trac, tracers
+  USE traclmdz_mod,     ONLY: traclmdz_from_restart
+  USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_cpl, co2_send
+  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
+  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
 #ifdef CPP_XIOS
 …
   INTEGER length
   PARAMETER (length=100)
   INTEGER it, iiq, isw
+  INTEGER it, iq, isw
   REAL tab_cntrl(length), tabcntr0(length)
   CHARACTER*7 str7
 …
   IF (type_trac == 'lmdz') THEN
+     DO it=1, nbtr
+!!        iiq=niadv(it+2)                                                           ! jyg
+        iiq=niadv(it+nqo)                                                           ! jyg
+        found=phyetat0_get(1,trs(:,it),"trs_"//TRIM(tracers(iiq)%name), &
+                                  "Surf trac"//TRIM(tracers(iiq)%name),0.)
+     ENDDO
+     it = 0
+     DO iq = 1, nqtot
+        IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+        it = it+1
+        found=phyetat0_get(1,trs(:,it),"trs_"//TRIM(tracers(iq)%name), &
+                                  "Surf trac"//TRIM(tracers(iq)%name),0.)
+     END DO
      CALL traclmdz_from_restart(trs)
   ENDIF
+!--OB now this is for co2i - ThL: and therefore also for inco
   IF (type_trac == 'co2i' .OR. type_trac == 'inco') THEN
      IF (carbon_cycle_cpl) THEN
 …
    CALL get_field(name, field, found)
    IF (.NOT. found) THEN
      WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
+     WRITE(lunout,*) "phyetat0: Le champ <",TRIM(name),"> est absent"
      WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
      field(:,:)=default

LMDZ6/trunk/libf/phylmd/phyredem.F90

-                      r4046
+                      r4056
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
   USE traclmdz_mod, ONLY : traclmdz_to_restart
   USE infotrac_phy, ONLY: type_trac, niadv, tracers, nbtr, nqo
+  USE infotrac_phy, ONLY: type_trac, nqtot, tracers, nbtr
   USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
 …
   CHARACTER (len=2) :: str2
   CHARACTER (len=256) :: nam, lnam
   INTEGER           :: it, iiq, pass
+  INTEGER           :: it, iq, pass
   !======================================================================
 …
     IF (type_trac == 'lmdz') THEN
        CALL traclmdz_to_restart(trs)
+       DO it=1, nbtr
+!!        iiq=niadv(it+2)                                                           ! jyg
+          iiq=niadv(it+nqo)                                                           ! jyg
+          CALL put_field(pass,"trs_"//tracers(iiq)%name, "", trs(:, it))
+       it = 0
+       DO iq = 1, nqtot
+          IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+          it = it+1
+          CALL put_field(pass,"trs_"//tracers(iq)%name, "", trs(:, it))
        END DO
     END IF
 …
   IMPLICIT NONE
   INTEGER, INTENT(IN)            :: pass
+  INTEGER, INTENT(IN)           :: pass
   CHARACTER(LEN=*),  INTENT(IN) :: nam, lnam
   REAL,              INTENT(IN) :: field(:,:)

LMDZ6/trunk/libf/phylmd/phys_output_mod.F90

-                      r4046
+                      r4056
     USE iophy
     USE dimphy
+    USE infotrac_phy, ONLY: nqtot, nqo, niadv, tracers, type_trac
+    USE strings_mod,  ONLY: maxlen
+    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso, maxlen
     USE ioipsl
     USE phys_cal_mod, only : hour, calend
 …
     CHARACTER(LEN=4), DIMENSION(nlevSTD)  :: clevSTD
     REAL, DIMENSION(nlevSTD)              :: rlevSTD
     INTEGER                               :: nsrf, k, iq, iiq, iff, i, j, ilev, jq
+    INTEGER                               :: nsrf, k, iq, iff, i, j, ilev, itr, ixt, iiso, izone
     INTEGER                               :: naero
     LOGICAL                               :: ok_veget
 …
     LOGICAL, DIMENSION(nfiles)            :: phys_out_filestations
-    CHARACTER(LEN=50) :: outiso
-    CHARACTER(LEN=20) :: unit
     CHARACTER(LEN=maxlen) :: tnam, lnam, dn
     INTEGER :: flag(nfiles)
 …
 !!!!!!!!!! stockage dans une region limitee pour chaque fichier !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !                 entre [phys_out_lonmin,phys_out_lonmax] et [phys_out_latmin,phys_out_latmax]
+    LOGICAL, DIMENSION(nfiles), SAVE  :: phys_out_regfkey       = (/ .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., &
+                                                                     .FALSE., .FALSE., .FALSE., .FALSE., .FALSE. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_lonmin        = (/ -180., -180., -180., -180., -180., &
+                                                                     -180., -180., -180., -180., -180. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_lonmax        = (/  180.,  180.,  180.,  180.,  180., &
+.,  180.,  180.,  180.,  180. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_latmin        = (/  -90.,  -90.,  -90.,  -90.,  -90., &
+                                                                      -90.,  -90.,  -90.,  -90.,  -90. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_latmax        = (/   90.,   90.,   90.,   90.,   90., &
+.,   90.,   90.,   90.,   90. /)
+    LOGICAL, DIMENSION(nfiles), SAVE :: &
+      phys_out_regfkey = [.FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE.]
+    REAL,    DIMENSION(nfiles), SAVE :: &
+      phys_out_lonmin  = [  -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.], &
+      phys_out_lonmax  = [   180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.], &
+      phys_out_latmin  = [   -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.], &
+      phys_out_latmax  = [    90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.]
     REAL, DIMENSION(klev,2) :: Ahyb_bounds, Bhyb_bounds
     REAL, DIMENSION(klev+1)   :: lev_index
 …
     ALLOCATE(o_dtr_dry(nqtot),o_dtr_vdf(nqtot))
     levmax = (/ klev, klev, klev, klev, klev, klev, nlevSTD, nlevSTD, nlevSTD, klev /)
+    levmax = [klev, klev, klev, klev, klev, klev, nlevSTD, nlevSTD, nlevSTD, klev]
     phys_out_filenames(1) = 'histmth'
 …
     CALL wxios_add_vaxis("bnds", 2, (/1.,2./))
      CALL wxios_add_vaxis("Alt", &
+    CALL wxios_add_vaxis("Alt", &
             levmax(iff) - levmin(iff) + 1, pseudoalt)
+    IF (NSW.EQ.6) THEN
+!
+!wl1_sun: minimum bound of wavelength (in um)
+!
+      wl1_sun(1)=0.180
+      wl1_sun(2)=0.250
+      wl1_sun(3)=0.440
+      wl1_sun(4)=0.690
+      wl1_sun(5)=1.190
+      wl1_sun(6)=2.380
+!
+!wl2_sun: maximum bound of wavelength (in um)
+!
+      wl2_sun(1)=0.250
+      wl2_sun(2)=0.440
+      wl2_sun(3)=0.690
+      wl2_sun(4)=1.190
+      wl2_sun(5)=2.380
+      wl2_sun(6)=4.000
+!
+    ELSE IF(NSW.EQ.2) THEN
+!
+!wl1_sun: minimum bound of wavelength (in um)
+!
+      wl1_sun(1)=0.250
+      wl1_sun(2)=0.690
+!
+!wl2_sun: maximum bound of wavelength (in um)
+!
+      wl2_sun(1)=0.690
+      wl2_sun(2)=4.000
+    ENDIF
+    ! wl1_sun/wl2_sun: minimum/maximum bound of wavelength (in um)
+    SELECT CASE(NSW)
+      CASE(6)
+        wl1_sun(1:6) = [0.180, 0.250, 0.440, 0.690, 1.190, 2.380]
+        wl2_sun(1:6) = [0.250, 0.440, 0.690, 1.190, 2.380, 4.000]
+      CASE(2)
+        wl1_sun(1:2) = [0.250, 0.690]
+        wl2_sun(1:2) = [0.690, 4.000]
+    END SELECT
     DO ISW=1, NSW
 …
      ENDIF ! clef_files
        IF (nqtot>=nqo+1) THEN
+!
           DO iq=nqo+1,nqtot
             iiq=niadv(iq); jq = iq-nqo
             dn = 'd'//TRIM(tracers(iiq)%name)//'_'
+          itr = 0
+          DO iq = 1, nqtot
+            IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+            itr = itr + 1
+            dn = 'd'//TRIM(tracers(iq)%name)//'_'
             flag = [1, 5, 5, 5, 10, 10, 11, 11, 11, 11]
+            lnam = 'Tracer '//TRIM(tracers(iiq)%longName)
+            tnam = TRIM(tracers(iiq)%name); o_trac          (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'Tracer '//TRIM(tracers(iq)%longName)
+            tnam = TRIM(tracers(iq)%name);  o_trac          (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [4, 7, 7, 7, 10, 10, 11, 11, 11, 11]
             lnam = 'Tendance tracer '//TRIM(tracers(iiq)%longName)
             tnam = TRIM(dn)//'vdf';         o_dtr_vdf       (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'Tendance tracer '//TRIM(tracers(iq)%longName)
+            tnam = TRIM(dn)//'vdf';         o_dtr_vdf       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [5, 7, 7, 7, 10, 10, 11, 11, 11, 11]
             tnam = TRIM(dn)//'the';         o_dtr_the       (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'con';         o_dtr_con       (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'the';         o_dtr_the       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'con';         o_dtr_con       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [7, 7, 7, 7, 10, 10, 11, 11, 11, 11]
             tnam = TRIM(dn)//'lessi_impa';  o_dtr_lessi_impa(jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'lessi_nucl';  o_dtr_lessi_nucl(jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'insc';        o_dtr_insc      (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'bcscav';      o_dtr_bcscav    (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'evapls';      o_dtr_evapls    (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'ls';          o_dtr_ls        (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'trsp';        o_dtr_trsp      (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'sscav';       o_dtr_sscav     (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'sat';         o_dtr_sat       (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             tnam = TRIM(dn)//'uscav';       o_dtr_uscav     (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             lnam = 'tracer tendency dry deposition'//TRIM(tracers(iiq)%longName)
             tnam = 'cum'//TRIM(dn)//'dry';  o_dtr_dry       (jq) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'lessi_impa';  o_dtr_lessi_impa(itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'lessi_nucl';  o_dtr_lessi_nucl(itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'insc';        o_dtr_insc      (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'bcscav';      o_dtr_bcscav    (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'evapls';      o_dtr_evapls    (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'ls';          o_dtr_ls        (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'trsp';        o_dtr_trsp      (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'sscav';       o_dtr_sscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'sat';         o_dtr_sat       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            tnam = TRIM(dn)//'uscav';       o_dtr_uscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'tracer tendency dry deposition'//TRIM(tracers(iq)%longName)
+            tnam = 'cum'//TRIM(dn)//'dry';  o_dtr_dry       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [1, 4, 10, 10, 10, 10, 11, 11, 11, 11]
             lnam = 'Cumulated tracer '//TRIM(tracers(iiq)%longName)
             tnam = 'cum'//TRIM(tracers(iiq)%name); o_trac_cum(jq)= ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'Cumulated tracer '//TRIM(tracers(iq)%longName)
+            tnam = 'cum'//TRIM(tracers(iq)%name); o_trac_cum(itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
           ENDDO
-       ENDIF
    ENDDO !  iff
 …
     ENDIF
+!  DO iq=nqo+1,nqtot
+!    iiq=niadv(iq)
+!    dn = 'd'//TRIM(tracers(iiq)%name)//'_'
+!    WRITE(*,'(a,i1,a,10i3)')'trac(',iiq,')%flag = ',o_trac(iiq)%flag
+!    WRITE(*,'(a,i1,a)')'trac(',iiq,')%tnam = '//TRIM(o_trac(iiq)%name)
+!    WRITE(*,'(a,i1,a)')'trac(',iiq,')%lnam = '//TRIM(o_trac(iiq)%description)
+!  DO iq=1,nqtot
+!    IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+!    WRITE(*,'(a,i1,a,10i3)')'trac(',iq,')%flag = ',o_trac(iq)%flag
+!    WRITE(*,'(a,i1,a)')'trac(',iq,')%name = '//TRIM(o_trac(iq)%name)
+!    WRITE(*,'(a,i1,a)')'trac(',iq,')%description = '//TRIM(o_trac(iq)%description)
 !  END DO

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

-                      r4046
+                      r4056
     USE dimphy, ONLY: klon, klev, klevp1
     USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niadv
+    USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niso, ntraciso, maxlen
     USE strings_mod,  ONLY: maxlen
     USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy
 …
          o_map_emis_Anv, o_map_pcld_Anv, o_map_tcld_Anv, &
          o_map_ntot, o_map_hc,o_map_hist,o_map_Cb,o_map_ThCi,o_map_Anv, &
+#ifdef ISO
+! Isotopes
+         o_xtprecip,o_xtplul,o_xtpluc,o_xtovap,o_xtoliq,o_xtcond, &
+         o_xtevap,o_dxtdyn,o_dxtldyn,o_dxtcon,o_dxtlsc,o_dxteva, &
+         o_dxtajs,o_dxtvdf,o_dxtthe, o_dxtch4, &
+         o_dxtprod_nucl,o_dxtcosmo,o_dxtdecroiss, &
+#endif
+! Tropopause
          o_alt_tropo, &
-! Tropopause
          o_p_tropopause, o_z_tropopause, o_t_tropopause,  &
          o_col_O3_strato, o_col_O3_tropo,                 &
 …
          rhlevSTD, O3STD, O3daySTD, uvSTD, vqSTD, vTSTD, wqSTD, vphiSTD, &
          wTSTD, u2STD, v2STD, T2STD, missing_val_nf90, delta_sal, ds_ns, &
+#ifdef ISO
+        xtrain_con, xtsnow_con, xtrain_fall, xtsnow_fall, &
+#endif
          dt_ns, delta_sst
 …
          east_gwstress, west_gwstress, &
          d_q_ch4, pmfd, pmfu, ref_liq, ref_ice, rhwriteSTD, &
+#ifdef ISO
+        xtrain_lsc, xtsnow_lsc, xt_seri, xtl_seri,xts_seri,xtevap, &
+        d_xt_dyn,d_xtl_dyn,d_xt_con,d_xt_vdf,d_xt_ajsb, &
+        d_xt_lsc,d_xt_eva,d_xt_ch4, &
+        d_xt_ajs, d_xt_ajsb, &
+        d_xt_prod_nucl,d_xt_cosmo,d_xt_decroiss, &
+#endif
          ep, epmax_diag, &  ! epmax_cape
          p_tropopause, t_tropopause, z_tropopause
 …
     USE pbl_surface_mod, ONLY: snow
     USE indice_sol_mod, ONLY: nbsrf
+#ifdef ISO
+    USE isotopes_mod, ONLY: iso_HTO
+#endif
     USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
     USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, landice_opt
 …
     ! Local
     INTEGER :: itau_w
     INTEGER :: i, iinit, iinitend=1, iff, iq, iiq, nsrf, k, ll, naero
+    INTEGER :: i, iinit, iinitend=1, iff, iq, nsrf, k, ll, naero
     REAL, DIMENSION (klon) :: zx_tmp_fi2d, zpt_conv2d, wind100m
     REAL, DIMENSION (klon,klev) :: zx_tmp_fi3d, zpt_conv
 …
 #endif
     REAL, PARAMETER :: un_jour=86400.
     INTEGER ISW
+    INTEGER :: ISW, itr, ixt, it
     CHARACTER*1 ch1
     CHARACTER(LEN=maxlen) :: varname, dn
 …
     REAL,DIMENSION(klon,klev) :: z, dz
     REAL,DIMENSION(klon)      :: zrho, zt
-    INTEGER :: nqup
     ! On calcul le nouveau tau:
 …
           CALL xios_get_handle("fields_strataer_trac_3D", group_handle)
           ! On boucle sur les traceurs pour les ajouter au groupe puis fixer les attributs
           DO iq=nqo+1, nqtot
             iiq=niadv(iq)
             dn = 'd'//TRIM(tracers(iiq)%name)//'_'
             WRITE (lunout,*) 'XIOS var=', nqo, iq, nqtot, tracers(iiq)%name
+          DO iq = 1, nqtot
+            IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+            dn = 'd'//TRIM(tracers(iq)%name)//'_'
+            WRITE (lunout,*) 'XIOS var=', nqo, iq, nqtot, tracers(iq)%name
             unt = "kg kg-1"
             varname=trim(tracers(iiq)%name)
+            varname=trim(tracers(iq)%name)
             CALL xios_add_child(group_handle, child, varname)
             CALL xios_set_attr(child, name=varname, unit=unt)
 …
             CALL xios_add_child(group_handle, child, varname)
             CALL xios_set_attr(child, name=varname, unit=unt)
           ENDDO
+          END DO
           !On ajoute les variables 2D traceurs par l interface fortran
           CALL xios_get_handle("fields_strataer_trac_2D", group_handle)
           ! On boucle sur les traceurs pour les ajouter au groupe puis fixer les attributs
           DO iq=nqo+1, nqtot
             iiq=niadv(iq)
+          DO iq = 1, nqtot
+            IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
             unt = "kg m-2"
             varname='cum'//trim(tracers(iiq)%name)
+            varname='cum'//trim(tracers(iq)%name)
             WRITE (lunout,*) 'XIOS var=', iq, nqtot, varname
             CALL xios_add_child(group_handle, child, varname)
 …
             unt = "kg m-2 s-1"
             varname='cumd'//trim(tracers(iiq)%name)//'_dry'
+            varname='cumd'//trim(tracers(iq)%name)//'_dry'
             CALL xios_add_child(group_handle, child, varname)
             CALL xios_set_attr(child, name=varname, unit=unt)
 …
        ENDIf
        CALL histwrite_phy(o_aire, zx_tmp_fi2d)
        IF (vars_defined) THEN
           DO i=1, klon
 …
        IF (iflag_phytrac == 1 ) then
          IF (type_trac == 'lmdz' .OR. type_trac == 'coag') THEN
+           DO iq=nqo+1, nqtot
+           itr = 0
+           DO iq = 1, nqtot
+             IF(tracers(iq)%isH2Ofamily) CYCLE
+             itr = itr + 1
+!            write(*,*) 'phys_output_write_mod 2337: itr=',itr
              !--3D fields
              CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_lessi_impa(iq-nqo),d_tr_lessi_impa(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_lessi_nucl(iq-nqo),d_tr_lessi_nucl(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_insc(iq-nqo),d_tr_insc(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_bcscav(iq-nqo),d_tr_bcscav(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_evapls(iq-nqo),d_tr_evapls(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_ls(iq-nqo),d_tr_ls(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_trsp(iq-nqo),d_tr_trsp(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_sscav(iq-nqo),d_tr_sscav(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_sat(iq-nqo),d_tr_sat(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo))
+             CALL histwrite_phy(o_trac(itr), tr_seri(:,:,itr))
+             CALL histwrite_phy(o_dtr_vdf(itr),d_tr_cl(:,:,itr))
+             CALL histwrite_phy(o_dtr_the(itr),d_tr_th(:,:,itr))
+             CALL histwrite_phy(o_dtr_con(itr),d_tr_cv(:,:,itr))
+             CALL histwrite_phy(o_dtr_lessi_impa(itr),d_tr_lessi_impa(:,:,itr))
+             CALL histwrite_phy(o_dtr_lessi_nucl(itr),d_tr_lessi_nucl(:,:,itr))
+             CALL histwrite_phy(o_dtr_insc(itr),d_tr_insc(:,:,itr))
+             CALL histwrite_phy(o_dtr_bcscav(itr),d_tr_bcscav(:,:,itr))
+             CALL histwrite_phy(o_dtr_evapls(itr),d_tr_evapls(:,:,itr))
+             CALL histwrite_phy(o_dtr_ls(itr),d_tr_ls(:,:,itr))
+             CALL histwrite_phy(o_dtr_trsp(itr),d_tr_trsp(:,:,itr))
+             CALL histwrite_phy(o_dtr_sscav(itr),d_tr_sscav(:,:,itr))
+             CALL histwrite_phy(o_dtr_sat(itr),d_tr_sat(:,:,itr))
+             CALL histwrite_phy(o_dtr_uscav(itr),d_tr_uscav(:,:,itr))
             !--2D fields
              CALL histwrite_phy(o_dtr_dry(iq-nqo), flux_tr_dry(:,iq-nqo))
+             CALL histwrite_phy(o_dtr_dry(itr), flux_tr_dry(:,itr))
              zx_tmp_fi2d=0.
              IF (vars_defined) THEN
                 DO k=1,klev
                    zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo)
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,itr)
                 ENDDO
              ENDIF
              CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
+             CALL histwrite_phy(o_trac_cum(itr), zx_tmp_fi2d)
            ENDDO !--iq
          ENDIF   !--type_trac
+!
          IF (type_trac == 'co2i') THEN
+           DO iq=nqo+1, nqtot
+           itr = 0
+           DO iq = 1, nqtot
+             IF(tracers(iq)%isH2Ofamily) CYCLE
+             itr = itr + 1
+!            write(*,*) 'phys_output_write_mod 2370: itr=',itr
              !--3D fields
              CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
+             CALL histwrite_phy(o_trac(itr), tr_seri(:,:,itr))
+             CALL histwrite_phy(o_dtr_vdf(itr),d_tr_cl(:,:,itr))
+             CALL histwrite_phy(o_dtr_the(itr),d_tr_th(:,:,itr))
+             CALL histwrite_phy(o_dtr_con(itr),d_tr_cv(:,:,itr))
              !--2D fields
              !--CO2 burden
 …
              IF (vars_defined) THEN
                 DO k=1,klev
                    zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo)
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,itr)
                 ENDDO
              ENDIF
              CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
+             CALL histwrite_phy(o_trac_cum(itr), zx_tmp_fi2d)
            ENDDO !--iq
            !--CO2 net fluxes
 …
          IF (type_trac == 'inco') THEN
+           nqup = nqo+1
+           DO iq=nqo+1, nqup
+           DO iq = 1, nqtot
+             IF(tracers(iq)%isH2Ofamily .OR. .NOT.tracers(iq)%isAdvected) CYCLE
+             itr = itr+1
+             IF(tracers(iq)%component /= 'co2i') CYCLE
              !--3D fields
              CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
+             CALL histwrite_phy(o_trac   (itr),tr_seri(:,:,itr))
+             CALL histwrite_phy(o_dtr_vdf(itr),d_tr_cl(:,:,itr))
+             CALL histwrite_phy(o_dtr_the(itr),d_tr_th(:,:,itr))
+             CALL histwrite_phy(o_dtr_con(itr),d_tr_cv(:,:,itr))
              !--2D fields
              !--CO2 burden
 …
              IF (vars_defined) THEN
                 DO k=1,klev
                    zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo)
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,itr)
                 ENDDO
              ENDIF
              CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
+             CALL histwrite_phy(o_trac_cum(itr), zx_tmp_fi2d)
            ENDDO !--iq
            !--CO2 net fluxes
 …
        end if
+#ifdef ISO
+    do ixt=1,ntraciso
+!        write(*,*) 'ixt'
+        IF (vars_defined) zx_tmp_fi2d(:) = xtrain_fall(ixt,:) + xtsnow_fall(ixt,:)
+        CALL histwrite_phy(o_xtprecip(ixt), zx_tmp_fi2d)
+        IF (vars_defined) zx_tmp_fi2d(:) = xtrain_lsc(ixt,:) + xtsnow_lsc(ixt,:)
+        CALL histwrite_phy(o_xtplul(ixt), zx_tmp_fi2d)
+        IF (vars_defined) zx_tmp_fi2d(:) = xtrain_con(ixt,:) + xtsnow_con(ixt,:)
+        CALL histwrite_phy(o_xtpluc(ixt), zx_tmp_fi2d)
+        CALL histwrite_phy(o_xtevap(ixt),   xtevap(ixt,:))
+        CALL histwrite_phy(o_xtovap(ixt),  xt_seri(ixt,:,:))
+        CALL histwrite_phy(o_xtoliq(ixt), xtl_seri(ixt,:,:))
+        IF (vars_defined) zx_tmp_fi3d(:,:)=xtl_seri(ixt,:,:)+xts_seri(ixt,:,:)
+        CALL histwrite_phy(o_xtcond(ixt), zx_tmp_fi3d)
+        CALL histwrite_phy(o_dxtdyn(ixt),   d_xt_dyn(ixt,:,:))
+        CALL histwrite_phy(o_dxtldyn(ixt), d_xtl_dyn(ixt,:,:))
+        IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_con(ixt,:,:)/pdtphys
+        CALL histwrite_phy(o_dxtcon(ixt), zx_tmp_fi3d)
+        IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_lsc(ixt,:,:)/pdtphys
+        CALL histwrite_phy(o_dxtlsc(ixt), zx_tmp_fi3d)
+        IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_eva(ixt,:,:)/pdtphys
+        CALL histwrite_phy(o_dxteva(ixt), zx_tmp_fi3d)
+        IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_vdf(ixt,:,:)/pdtphys
+        CALL histwrite_phy(o_dxtvdf(ixt), zx_tmp_fi3d)
+        IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_ajsb(ixt,:,:)/pdtphys
+        CALL histwrite_phy(o_dxtajs(ixt), zx_tmp_fi3d)
+        IF (vars_defined) zx_tmp_fi3d(:,:)=(d_xt_ajs(ixt,:,:)-d_xt_ajsb(ixt,:,:))/pdtphys
+        CALL histwrite_phy(o_dxtthe(ixt), zx_tmp_fi3d)
+        IF (ok_qch4) THEN
+          IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_ch4(ixt,:,:)/pdtphys
+          CALL histwrite_phy(o_dxtch4(ixt), zx_tmp_fi3d)
+        END IF
+        IF (ixt == iso_HTO) THEN
+          IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_prod_nucl(ixt,:,:)/pdtphys
+          CALL histwrite_phy(o_dxtprod_nucl(ixt), zx_tmp_fi3d)
+          IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_cosmo(ixt,:,:)/pdtphys
+          CALL histwrite_phy(o_dxtcosmo(ixt), zx_tmp_fi3d)
+          IF (vars_defined) zx_tmp_fi3d(:,:)=d_xt_decroiss(ixt,:,:)/pdtphys
+          CALL histwrite_phy(o_dxtdecroiss(ixt), zx_tmp_fi3d)
+        END IF
+    !write(*,*) 'phys_output_write_mod 2531'
+    enddo !do ixt=1,ntraciso
+#endif
        IF (.NOT.vars_defined) THEN
           !$OMP MASTER

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

-                      r4009
+                      r4056
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac, nqCO2
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, nqCO2
+    USE readTracFiles_mod, ONLY: phases_sep
+    USE strings_mod,  ONLY: strIdx
     USE iophy
     USE limit_read_mod, ONLY : init_limit_read
 …
+       !
        d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
+       d_t_vdf_x, d_t_vdf_w, &
+       d_q_vdf_x, d_q_vdf_w, &
        d_ts, &
+       !
 …
        zxfluxlat_x, zxfluxlat_w, &
+       !
-       d_t_vdf_x, d_t_vdf_w, &
-       d_q_vdf_x, d_q_vdf_w, &
        pbl_tke_input, tke_dissip, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
 …
     LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
     !$OMP THREADPRIVATE(ok_volcan)
     INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf ou dans la strato
+    INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato
     !$OMP THREADPRIVATE(flag_volc_surfstrat)
     LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
 …
     real zqsat(klon,klev)
+    !
     INTEGER i, k, iq, j, nsrf, ll, l
+    INTEGER i, k, iq, j, nsrf, ll, l, itr
+    !
     REAL t_coup
 …
     CHARACTER (LEN=20) :: modname='physiq_mod'
     CHARACTER*80 message, abort_message
+    CHARACTER*80 abort_message
     LOGICAL, SAVE ::  ok_sync, ok_sync_omp
     !$OMP THREADPRIVATE(ok_sync)
 …
          iflag_phytrac = 1
        ENDIF
 #endif
+#endif
        nvm_lmdz = 13
        CALL getin_p('NVM',nvm_lmdz)
 …
     tke0(:,:)=pbl_tke(:,:,is_ave)
+    !CR:Nombre de traceurs de l'eau: nqo
+    !  IF (nqtot.GE.3) THEN
+    IF (nqtot.GE.(nqo+1)) THEN
+       !     DO iq = 3, nqtot
+       DO iq = nqo+1, nqtot
+    IF (nqtot > nqo) THEN
+       ! water isotopes are not included in tr_seri
+       itr = 0
+       DO iq = 1, nqtot
+         IF(tracers(iq)%isH2Ofamily) CYCLE
+         itr = itr+1
           DO  k = 1, klev
              DO  i = 1, klon
+                !              tr_seri(i,k,iq-2) = qx(i,k,iq)
+                tr_seri(i,k,iq-nqo) = qx(i,k,iq)
+                tr_seri(i,k,itr) = qx(i,k,iq)
              ENDDO
           ENDDO
        ENDDO
     ELSE
+       DO k = 1, klev
+          DO i = 1, klon
+             tr_seri(i,k,1) = 0.0
+          ENDDO
+       ENDDO
+! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
+!       tr_seri(:,:,strIdx(tracers(:)%name,'H2O'//phases_sep//'g')) = 0.0
+       tr_seri(:,:,strIdx(tracers(:)%name,'H2Ov')) = 0.0
     ENDIF
+!
 …
     IF (debut) THEN
       WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
+      DO iq = nqo+1, nqtot
+           tr_ancien(:,:,iq-nqo)=tr_seri(:,:,iq-nqo)
+      ENDDO
+       itr = 0
+       do iq = 1, nqtot
+         IF(tracers(iq)%isH2Ofamily) CYCLE
+         itr = itr+1
+         tr_ancien(:,:,itr)=tr_seri(:,:,itr)
+       enddo
     ENDIF
+    !
 …
        d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
        ! !! RomP >>>   td dyn traceur
+       IF (nqtot.GT.nqo) THEN     ! jyg
+          DO iq = nqo+1, nqtot      ! jyg
+              d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
+          ENDDO
+       ENDIF
+       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
        ! !! RomP <<<
     ELSE
 …
        d_qs_dyn2d(:) = 0.0
        ! !! RomP >>>   td dyn traceur
+       IF (nqtot.GT.nqo) THEN                                       ! jyg
+          DO iq = nqo+1, nqtot                                      ! jyg
+              d_tr_dyn(:,:,iq-nqo)= 0.0                             ! jyg
+          ENDDO
+       ENDIF
+       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
        ! !! RomP <<<
        ancien_ok = .TRUE.
 …
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
              sollwdown,    cldt,      &
+            sollwdown,    cldt,      &
             rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &
 …
          ENDDO
        ELSE
                t_w(:,:) = t_seri(:,:)
+                t_w(:,:) = t_seri(:,:)
                 q_w(:,:) = q_seri(:,:)
                 t_x(:,:) = t_seri(:,:)
 …
        DO i = 1, klon
+          ema_pcb(i)  = paprs(i,ibas_con(i))
+          ! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
+          ! où i n'est pas un point convectif et donc ibas_con(i)=0
+          ! c'est un pb indépendant des isotopes
+          if (ibas_con(i) > 0) then
+             ema_pcb(i)  = paprs(i,ibas_con(i))
+          else
+             ema_pcb(i)  = 0.0
+          endif
        ENDDO
        DO i = 1, klon
 …
     wprime_ave(:,:)=0.
     DO nsrf = 1, nbsrf
        DO i = 1, klon
 …
        ENDDO
     ENDDO
     CALL  calcratqs(klon,klev,prt_level,lunout,        &
 …
        print *,'itap, ->fisrtilp ',itap
     ENDIF
+    !
     picefra(:,:)=0.
 …
          iflag_ice_thermo)
     ENDIF
+    !
     WHERE (rain_lsc < 0) rain_lsc = 0.
     WHERE (snow_lsc < 0) snow_lsc = 0.
 …
 #ifndef CPP_XIOS
+          !--OB 30/05/2016 modified 21/10/2016
+          !--here we return swaero_diag and dryaod_diag to FALSE
+          !--and histdef will switch it back to TRUE if necessary
+          !--this is necessary to get the right swaero at first step
+          !--but only in the case of no XIOS as XIOS is covered elsewhere
+          IF (debut) swaerofree_diag = .FALSE.
+          IF (debut) swaero_diag = .FALSE.
+          IF (debut) dryaod_diag = .FALSE.
+          !--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
+          !--as for swaero_diag, see above
+          IF (debut) ok_4xCO2atm = .FALSE.
+          !
           !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
           !IM des taux doit etre different du taux actuel
 …
     ENDDO
+    !
     !CR: nb de traceurs eau: nqo
     !  IF (nqtot.GE.3) THEN
     IF (nqtot.GE.(nqo+1)) THEN
        !     DO iq = 3, nqtot
        DO iq = nqo+1, nqtot
+    IF (nqtot > nqo+1) THEN
+       itr = 0
+       DO iq = 1, nqtot
+          IF(tracers(iq)%isH2Ofamily) CYCLE
+          itr = itr+1
           DO  k = 1, klev
              DO  i = 1, klon
+                ! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
+                d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
+                d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
              ENDDO
           ENDDO
 …
     CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
     ! !! RomP >>>
+    !CR: nb de traceurs eau: nqo
+    IF (nqtot.GT.nqo) THEN
+       DO iq = nqo+1, nqtot
+          tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
+       ENDDO
+    ENDIF
+    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
     ! !! RomP <<<
     !==========================================================================

LMDZ6/trunk/libf/phylmd/phytrac_mod.F90

-                      r4012
+                      r4056
   SUBROUTINE phytrac_init()
     USE dimphy
     USE infotrac_phy, ONLY: nbtr, nqCO2, type_trac
+    USE infotrac_phy, ONLY: nbtr, type_trac
     USE tracco2i_mod, ONLY: tracco2i_init
     IMPLICIT NONE
 …
     USE phys_local_var_mod, ONLY: budg_dep_dry_h2so4, budg_dep_wet_h2so4
     USE phys_local_var_mod, ONLY: budg_dep_dry_part,  budg_dep_wet_part
     USE infotrac, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
+    USE infotrac_phy, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
     USE aerophys
 #endif
 …
           iflag_con_trac= 1
        CASE('inco')
           source(:,1:nqCO2) = 0.                          ! from CO2i   ModThL
           source(:,nqCO2+1:nbtr)=init_source(:,:)         ! from INCA   ModThL
           aerosol(1:nqCO2) = .FALSE.                      ! from CO2i   ModThL
           CALL tracinca_init(aerosol(nqCO2+1:nbtr),lessivage)     ! from INCA   ModThL
           pbl_flg(1:nqCO2) = 1              ! From CO2i         ModThL
           iflag_the_trac= 1           ! From CO2i
           iflag_vdf_trac= 1           ! From CO2i
           iflag_con_trac= 1           ! From CO2i
+          source(:,1:nqCO2) = 0.                          ! from CO2i ModThL
+          source(:,nqCO2+1:nbtr)=init_source(:,:)         ! from INCA ModThL
+          aerosol(1:nqCO2) = .FALSE.                      ! from CO2i ModThL
+          CALL tracinca_init(aerosol(nqCO2+1:nbtr),lessivage)     ! from INCA ModThL
+          pbl_flg(1:nqCO2) = 1                            ! From CO2i ModThL
+          iflag_the_trac = 1                              ! From CO2i
+          iflag_vdf_trac = 1                              ! From CO2i
+          iflag_con_trac = 1                              ! From CO2i
 #ifdef CPP_StratAer
        CASE('coag')

LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90

-                      r4050
+                      r4056
     USE dimphy
     USE infotrac_phy, ONLY: nbtr, tracers, niadv, solsym
+    USE infotrac_phy, ONLY: nbtr
     ! Input argument
 …
     ! Initialization of the tracers should be done here only for those not found in the restart file.
     USE dimphy
     USE infotrac_phy, ONLY: nbtr, nqo, tracers, pbl_flg, conv_flg, niadv
+    USE infotrac_phy, ONLY: nbtr, nqtot, tracers, pbl_flg, conv_flg
     USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
     USE press_coefoz_m, ONLY: press_coefoz
 …
 ! Local variables
     INTEGER :: ierr, it, iiq, i, k
+    INTEGER :: ierr, it, iq, i, k
     REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global
     REAL, DIMENSION(klev)          :: mintmp, maxtmp
 …
     id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
     id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
+    DO it=1,nbtr
+!!       iiq=niadv(it+2)                                                            ! jyg
+       iiq=niadv(it+nqo)                                                            ! jyg
+       SELECT CASE(strLower(tracers(iiq)%name))
+    it = 0
+    DO iq = 1, nqtot
+       IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+       it = it+1
+       SELECT CASE(strLower(tracers(iq)%name))
          CASE("rn");      id_rn     = it ! radon
          CASE("pb");      id_pb     = it ! plomb
 …
          CASE("pcq0");    id_pcq0   = it
          CASE DEFAULT
            WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iiq)%name)
+           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
        END SELECT
        SELECT CASE(strLower(tracers(iiq)%name))
+       SELECT CASE(strLower(tracers(iq)%name))
          CASE("pb")                        !--- RomP >>> profil initial de PB210
            OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
 …
 ! Check if all tracers have restart values
 ! ----------------------------------------------
+    DO it=1,nbtr
+!!       iiq=niadv(it+2)                                                            ! jyg
+       iiq=niadv(it+nqo)                                                            ! jyg
+    it = 0
+    DO iq = 1, nqtot
+       IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+       it = it+1
        ! Test if tracer is zero everywhere.
        ! Done by master process MPI and master thread OpenMP
 …
        IF (zero) THEN
           ! The tracer was not found in restart file or it was equal zero everywhere.
           WRITE(lunout,*) "The tracer ",trim(tracers(iiq)%name)," will be initialized"
+          WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
           IF (it==id_pcsat .OR. it==id_pcq .OR. &
                it==id_pcs0 .OR. it==id_pcq0) THEN

LMDZ6/trunk/libf/phylmdiso/phyetat0.F90

-                      r4050
+                      r4056
        ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, ratqs_inter
 !FC
+  USE geometry_mod, ONLY : longitude_deg, latitude_deg
+  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
+  USE infotrac_phy, only: nbtr, nqo, type_trac, tracers, niadv, &
+        itr_indice ! C Risi
+  USE traclmdz_mod,    ONLY : traclmdz_from_restart
+  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
+  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
+  USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
+  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
+  USE iostart,          ONLY: close_startphy, get_field, get_var, open_startphy
+  USE infotrac_phy,     ONLY: nqtot, nbtr, type_trac, tracers
+  USE traclmdz_mod,     ONLY: traclmdz_from_restart
+  USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_cpl, co2_send
+  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
+  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
 #ifdef CPP_XIOS
 …
   INTEGER length
   PARAMETER (length=100)
   INTEGER it, iiq, isw
+  INTEGER it, iq, isw
   REAL tab_cntrl(length), tabcntr0(length)
   CHARACTER*7 str7
 …
   REAL Rland_ice(niso,klon)
 #endif
-  INTEGER iq
   ! FH1D
   !     real iolat(jjm+1)
 …
   IF (type_trac == 'lmdz') THEN
+     DO it=1, nbtr
+!!        iiq=niadv(it+2)                                                           ! jyg
+!        iiq=niadv(it+nqo)   C Risi: on efface pour remplacer
+        iq=itr_indice(it)
+        iiq=niadv(iq)                                                        ! jyg
+        found=phyetat0_get(1,trs(:,it),"trs_"//TRIM(tracers(iiq)%name), &
+                                  "Surf trac"//TRIM(tracers(iiq)%name),0.)
+     ENDDO
+     it = 0
+     DO iq = 1, nqtot
+        IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+        it = it+1
+        found=phyetat0_get(1,trs(:,it),"trs_"//TRIM(tracers(iq)%name), &
+                                  "Surf trac"//TRIM(tracers(iq)%name),0.)
+     END DO
      CALL traclmdz_from_restart(trs)
   ENDIF
 …
    CALL get_field(name, field, found)
    IF (.NOT. found) THEN
      WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
+     WRITE(lunout,*) "phyetat0: Le champ <",TRIM(name),"> est absent"
      WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
      field(:,:)=default

LMDZ6/trunk/libf/phylmdiso/phyredem.F90

-                      r4046
+                      r4056
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
                                 wake_deltat, wake_deltaq, wake_s, wake_dens, &
+                                awake_dens, cv_gen,                          &
                                 wake_cstar,                                  &
                                 wake_pe, wake_fip, fm_therm, entr_therm,     &
 …
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
   USE traclmdz_mod, ONLY : traclmdz_to_restart
   USE infotrac_phy, ONLY: type_trac, niadv, tracers, nbtr, nqo,itr_indice
+  USE infotrac_phy, ONLY: type_trac, nqtot, tracers, nbtr, niso, ntraciso
 #ifdef ISO
-  USE infotrac_phy, ONLY: itr_indice,niso,ntraciso
 #ifdef ISOVERIF
   USE isotopes_verif_mod
 …
   REAL Rland_ice(niso,klon)
 #endif
-  INTEGER iq ! C Risi
   INTEGER nid, nvarid, idim1, idim2, idim3
 …
   CHARACTER (len=2) :: str2
   CHARACTER (len=256) :: nam, lnam
   INTEGER           :: it, iiq, pass
+  INTEGER           :: it, iq, pass
   !======================================================================
 …
     CALL put_field_srf1(pass,"TS","Temperature",ftsol(:,:))
+!!    CALL put_field_srf1(pass,"DELTA_TS","w-x surface temperature difference", delta_tsurf(:,:))
     CALL put_field_srf1(pass,"DELTATS","w-x surface temperature difference", delta_tsurf(:,:))
 !    CALL put_field_srf1(pass,"BETA_S","Aridity factor", beta_aridity(:,:))
     CALL put_field_srf1(pass,"BETAS","Aridity factor", beta_aridity(:,:))
+    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then
+       CALL put_field_srf1(pass, "DELTATS", &
+                      "w-x surface temperature difference",  delta_tsurf(:,:))
+       CALL put_field_srf1(pass, "BETAS", "Aridity factor", beta_aridity(:,:))
+    end IF
 !    End surface variables
 …
     CALL put_field(pass,"WAKE_DENS", "Wake num. /unit area", wake_dens)
+    CALL put_field(pass,"AWAKE_DENS", "Active Wake num. /unit area", awake_dens)
+    CALL put_field(pass,"CV_GEN", "CB birth rate", cv_gen)
     CALL put_field(pass,"WAKE_CSTAR", "WAKE_CSTAR", wake_cstar)
 …
     IF (type_trac == 'lmdz') THEN
        CALL traclmdz_to_restart(trs)
+       DO it=1, nbtr
+!!        iiq=niadv(it+2)                                                           ! jyg
+          !iiq=niadv(it+nqo) ! C Risi: on efface pour remplacer:
+          iq=itr_indice(it)                                                           ! jyg
+          iiq=niadv(iq)                                                           ! jyg
+          CALL put_field(pass,"trs_"//tracers(iiq)%name, "", trs(:, it))
+       it = 0
+       DO iq = 1, nqtot
+          IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+          it = it+1
+          CALL put_field(pass,"trs_"//tracers(iq)%name, "", trs(:, it))
        END DO
+    END IF
+    IF (type_trac == 'co2i' .OR. type_trac == 'inco') THEN
        IF (carbon_cycle_cpl) THEN
           IF (.NOT. ALLOCATED(co2_send)) THEN
 …
   IMPLICIT NONE
   INTEGER, INTENT(IN)            :: pass
+  INTEGER, INTENT(IN)           :: pass
   CHARACTER(LEN=*),  INTENT(IN) :: nam, lnam
   REAL,              INTENT(IN) :: field(:,:)
 …
       CHARACTER*7 str7
       CHARACTER*2 str2
       CHARACTER*50 striso_sortie
+      CHARACTER*50 outiso
       integer lnblnk
 #ifdef ISOTRAC
 …
      if (ixt.le.niso) then
         striso_sortie=striso(ixt)
+        outiso=striso(ixt)
      else
 #ifdef ISOTRAC
         iiso=index_iso(ixt)
         izone=index_zone(ixt)
         striso_sortie=striso(iiso)//strtrac(izone)
+        outiso=striso(iiso)//strtrac(izone)
 #else
         write(*,*) 'phyredem 546: ixt,ntraciso=', ixt,ntraciso
 …
 #endif
      endif !if (ixt.le.niso) then
       write(*,*) 'phyredem 550: ixt,striso_sortie=',ixt,striso_sortie(1:lnblnk(striso_sortie))
+      write(*,*) 'phyredem 550: ixt,outiso=',ixt,TRIM(outiso)
       iso_tmp_lonsrf(:,:)=fxtevap(ixt,:,:)
+      CALL put_field_srf1(pass,"XTEVAP"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &           "Evaporation de surface",iso_tmp_lonsrf)
+      CALL put_field_srf1(pass, "XTEVAP"//TRIM(outiso), "Evaporation de surface",iso_tmp_lonsrf)
       iso_tmp_lonsrf(:,:)=xtsnow(ixt,:,:)
+      CALL put_field_srf1(pass,"XTSNOW"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &           "NEIGE",iso_tmp_lonsrf)
+      CALL put_field_srf1(pass, "XTSNOW"//TRIM(outiso), "NEIGE",       iso_tmp_lonsrf)
       iso_tmp(:)=xtrain_fall(ixt,:)
+      CALL put_field(pass,"xtrain_f"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "precipitation liquide",iso_tmp)
+      CALL put_field(pass,    "xtrain_f"//TRIM(outiso), "precipitation liquide",iso_tmp)
       iso_tmp(:)=xtsnow_fall(ixt,:)
+      CALL put_field(pass,"xtsnow_f"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &    "precipitation solide",iso_tmp)
+      CALL put_field(pass,    "xtsnow_f"//TRIM(outiso), "precipitation solide",iso_tmp)
       iso_tmp_lonlev(:,:)=xt_ancien(ixt,:,:)
+      CALL put_field(pass,"XTANCIEN"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "QANCIEN",iso_tmp_lonlev)
+      CALL put_field(pass,    "XTANCIEN"//TRIM(outiso), "QANCIEN",     iso_tmp_lonlev)
       iso_tmp_lonlev(:,:)=xtl_ancien(ixt,:,:)
+      CALL put_field(pass,"XTLANCIEN"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "QLANCIEN",iso_tmp_lonlev)
+      CALL put_field(pass,   "XTLANCIEN"//TRIM(outiso), "QLANCIEN",    iso_tmp_lonlev)
       iso_tmp_lonlev(:,:)=xts_ancien(ixt,:,:)
+      CALL put_field(pass,"XTSANCIEN"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "QSANCIEN",iso_tmp_lonlev)
+      CALL put_field(pass,   "XTSANCIEN"//TRIM(outiso), "QSANCIEN",    iso_tmp_lonlev)
       iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
+      CALL put_field(pass,"WAKE_DELTAXT"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "WAKE_DELTAQ",iso_tmp_lonlev)
+      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAQ", iso_tmp_lonlev)
       iso_tmp(:)=xtrun_off_lic_0(ixt,:)
+      CALL put_field(pass,"XTRUNOFFLIC0"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "Runofflic0",iso_tmp)
+      CALL put_field(pass,"XTRUNOFFLIC0"//TRIM(outiso), "Runofflic0",  iso_tmp)
       iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
+      CALL put_field(pass,"WAKE_DELTAXT"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "WAKE_DELTAXT",iso_tmp_lonlev)
+      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAXT",iso_tmp_lonlev)
       ! variables seulement pour niso:
 …
       iso_tmp(:)=xtsol(ixt,:)
+      CALL put_field(pass,"XTSOL"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "Eau dans le sol (mm)",iso_tmp)
+      CALL put_field(pass,      "XTSOL"//TRIM(outiso), "Eau dans le sol (mm)",iso_tmp)
       iso_tmp(:)=Rland_ice(ixt,:)
+      CALL put_field(pass,"Rland_ice"//striso_sortie(1:lnblnk(striso_sortie)), &
+     &   "ratio land ice",iso_tmp)
+      CALL put_field(pass,  "Rland_ice"//TRIM(outiso), "ratio land ice",      iso_tmp)
       endif ! if (ixt.le.niso) then

LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90

-                      r4050
+                      r4056
     USE iophy
     USE dimphy
+    USE infotrac_phy, ONLY: nqtot, nqo, niadv, tracers, type_trac, maxlen, &
+        nqtottr,itr_indice ! C Risi
+    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso, maxlen
     USE ioipsl
     USE phys_cal_mod, only : hour, calend
 …
 #endif
 #ifdef ISO
-    USE infotrac_phy,ONLY: niso, ntraciso
     USE isotopes_mod, ONLY: striso,iso_HTO
 #ifdef ISOTRAC
 …
     CHARACTER(LEN=4), DIMENSION(nlevSTD)  :: clevSTD
     REAL, DIMENSION(nlevSTD)              :: rlevSTD
+    INTEGER                               :: nsrf, k, iq, iiq, iff, i, j, ilev
+    INTEGER                               :: itr ! C Risi
+    INTEGER                               :: nsrf, k, iq, iff, i, j, ilev, itr, ixt, iiso, izone
     INTEGER                               :: naero
     LOGICAL                               :: ok_veget
 …
 #ifdef ISO
-    INTEGER  :: ixt,iiso,izone
     CHARACTER(LEN=LEN(striso)) :: outiso
     CHARACTER(LEN=20) :: unit
 …
 !!!!!!!!!! stockage dans une region limitee pour chaque fichier !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !                 entre [phys_out_lonmin,phys_out_lonmax] et [phys_out_latmin,phys_out_latmax]
+    LOGICAL, DIMENSION(nfiles), SAVE  :: phys_out_regfkey       = (/ .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., &
+                                                                     .FALSE., .FALSE., .FALSE., .FALSE., .FALSE. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_lonmin        = (/ -180., -180., -180., -180., -180., &
+                                                                     -180., -180., -180., -180., -180. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_lonmax        = (/  180.,  180.,  180.,  180.,  180., &
+.,  180.,  180.,  180.,  180. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_latmin        = (/  -90.,  -90.,  -90.,  -90.,  -90., &
+                                                                      -90.,  -90.,  -90.,  -90.,  -90. /)
+    REAL, DIMENSION(nfiles), SAVE     :: phys_out_latmax        = (/   90.,   90.,   90.,   90.,   90., &
+.,   90.,   90.,   90.,   90. /)
+    LOGICAL, DIMENSION(nfiles), SAVE :: &
+      phys_out_regfkey = [.FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE.]
+    REAL,    DIMENSION(nfiles), SAVE :: &
+      phys_out_lonmin  = [  -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.,   -180.], &
+      phys_out_lonmax  = [   180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.,    180.], &
+      phys_out_latmin  = [   -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.,    -90.], &
+      phys_out_latmax  = [    90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.,     90.]
     REAL, DIMENSION(klev,2) :: Ahyb_bounds, Bhyb_bounds
     REAL, DIMENSION(klev+1)   :: lev_index
 …
     CALL wxios_add_vaxis("bnds", 2, (/1.,2./))
      CALL wxios_add_vaxis("Alt", &
+    CALL wxios_add_vaxis("Alt", &
             levmax(iff) - levmin(iff) + 1, pseudoalt)
+    IF (NSW.EQ.6) THEN
+!
+!wl1_sun: minimum bound of wavelength (in um)
+!
+      wl1_sun(1)=0.180
+      wl1_sun(2)=0.250
+      wl1_sun(3)=0.440
+      wl1_sun(4)=0.690
+      wl1_sun(5)=1.190
+      wl1_sun(6)=2.380
+!
+!wl2_sun: maximum bound of wavelength (in um)
+!
+      wl2_sun(1)=0.250
+      wl2_sun(2)=0.440
+      wl2_sun(3)=0.690
+      wl2_sun(4)=1.190
+      wl2_sun(5)=2.380
+      wl2_sun(6)=4.000
+!
+    ELSE IF(NSW.EQ.2) THEN
+!
+!wl1_sun: minimum bound of wavelength (in um)
+!
+      wl1_sun(1)=0.250
+      wl1_sun(2)=0.690
+!
+!wl2_sun: maximum bound of wavelength (in um)
+!
+      wl2_sun(1)=0.690
+      wl2_sun(2)=4.000
+    ENDIF
+    ! wl1_sun/wl2_sun: minimum/maximum bound of wavelength (in um)
+    SELECT CASE(NSW)
+      CASE(6)
+        wl1_sun(1:6) = [0.180, 0.250, 0.440, 0.690, 1.190, 2.380]
+        wl2_sun(1:6) = [0.250, 0.440, 0.690, 1.190, 2.380, 4.000]
+      CASE(2)
+        wl1_sun(1:2) = [0.250, 0.690]
+        wl2_sun(1:2) = [0.690, 4.000]
+    END SELECT
     DO ISW=1, NSW
 …
      ENDIF ! clef_files
+        write(lunout,*) 'phys_output_mid 496: nqtottr=',nqtottr
+        write(lunout,*) 'itr_indice=',itr_indice
+!       IF (nqtot>=nqo+1) THEN
+!
+            !DO iq=nqo+1,nqtot
+            ! C Risi: on modifie la boucle
+          DO itr=1,nqtottr ! C Risi
+            iq=itr_indice(itr)  ! C Risi
+            write(*,*) 'phys_output_mid 503: itr=',itr
+            iiq=niadv(iq)
+            dn = 'd'//TRIM(tracers(iiq)%name)//'_'
+          itr = 0
+          DO iq = 1, nqtot
+            IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
+            itr = itr + 1
+            dn = 'd'//TRIM(tracers(iq)%name)//'_'
             flag = [1, 5, 5, 5, 10, 10, 11, 11, 11, 11]
             lnam = 'Tracer '//TRIM(tracers(iiq)%longName)
             tnam = TRIM(tracers(iiq)%name); o_trac          (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'Tracer '//TRIM(tracers(iq)%longName)
+            tnam = TRIM(tracers(iq)%name);  o_trac          (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [4, 7, 7, 7, 10, 10, 11, 11, 11, 11]
             lnam = 'Tendance tracer '//TRIM(tracers(iiq)%longName)
+            lnam = 'Tendance tracer '//TRIM(tracers(iq)%longName)
             tnam = TRIM(dn)//'vdf';         o_dtr_vdf       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
 …
             tnam = TRIM(dn)//'uscav';       o_dtr_uscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             lnam = 'tracer tendency dry deposition'//TRIM(tracers(iiq)%longName)
+            lnam = 'tracer tendency dry deposition'//TRIM(tracers(iq)%longName)
             tnam = 'cum'//TRIM(dn)//'dry';  o_dtr_dry       (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
             flag = [1, 4, 10, 10, 10, 10, 11, 11, 11, 11]
             lnam = 'Cumulated tracer '//TRIM(tracers(iiq)%longName)
             tnam = 'cum'//TRIM(tracers(iiq)%name); o_trac_cum(itr)= ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
+            lnam = 'Cumulated tracer '//TRIM(tracers(iq)%longName)
+            tnam = 'cum'//TRIM(tracers(iq)%name); o_trac_cum(itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
           ENDDO
    ENDDO !  iff
-        write(*,*) 'phys_output_mid 589'
 #ifdef ISO
+    write(*,*) 'phys_output_mid 589'
     do ixt=1,ntraciso
       if (ixt.le.niso) then
+      if (ixt <= niso) then
         outiso=striso(ixt)
       else
 …
       END IF
     enddo !do ixt=1,niso
+#endif
+        write(*,*) 'phys_output_mid 596'
+    write(*,*) 'phys_output_mid 596'
+#endif

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

-                      r4050
+                      r4056
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac,ok_isotopes, &
         nqtottr,itr_indice ! C Risi
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, nqCO2, ok_isotopes
+    USE readTracFiles_mod, ONLY: phases_sep
+    USE strings_mod,  ONLY: strIdx
     USE iophy
     USE limit_read_mod, ONLY : init_limit_read
 …
     USE phys_output_mod
     USE phys_output_ctrlout_mod
+    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
+    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level, &
+         alert_first_call, call_alert, prt_alerte
     USE readaerosol_mod, ONLY : init_aero_fromfile
     USE readaerosolstrato_m, ONLY : init_readaerosolstrato
 …
 #ifdef ISO
     USE infotrac_phy, ONLY:  &
+        iqiso,iso_indnum,tracers,ok_isotrac, &
+        niso,ntraciso,nqtottr,itr_indice ! ajout C Risi pour isos
+        iqiso,iso_indnum,ok_isotrac,niso, ntraciso
      USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO, &
         & bidouille_anti_divergence,ok_bidouille_wake, &
 …
+       !
        d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
        d_t_vdf_w,d_q_vdf_w, &
        d_t_vdf_x,d_q_vdf_x, &
+       d_t_vdf_x, d_t_vdf_w, &
+       d_q_vdf_x, d_q_vdf_w, &
        d_ts, &
+       !
 …
        zxfluxlat_x, zxfluxlat_w, &
+       !
-       d_t_vdf_x, d_t_vdf_w, &
-       d_q_vdf_x, d_q_vdf_w, &
        pbl_tke_input, tke_dissip, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
 …
     real zqsat(klon,klev)
+    !
+    INTEGER i, k, iq, j, nsrf, ll, l
+    INTEGER itr ! C Risi
+    INTEGER i, k, iq, j, nsrf, ll, l, itr
 #ifdef ISO
     real zxt_apres(ntraciso,klon)
 …
 !JLD    REAL zstophy, zout
     CHARACTER*20 modname
+    CHARACTER (LEN=20) :: modname='physiq_mod'
     CHARACTER*80 abort_message
     LOGICAL, SAVE ::  ok_sync, ok_sync_omp
 …
     pi = 4. * ATAN(1.)
+    ! set-up call to alerte function
+    call_alert = (alert_first_call .AND. is_master)
     ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
     jjmp1=nbp_lat
 …
     forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
-    modname = 'physiq'
     IF (debut) THEN
 …
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
+       ! Initialisation des champs dans phytrac* qui sont utilisés par phys_output_write*
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       ! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 #ifdef CPP_Dust
        ! Quand on utilise SPLA, on force iflag_phytrac=1
 …
             ENDDO
           ENDDO
         ELSE
+       ELSE
           pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
 !>jyg
 …
           CALL abort_physic(modname,abort_message,1)
        ENDIF
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       ! Initialisation pour la convection de K.E. et pour les poches froides
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
+       WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
+            ok_cvl
+       WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", ok_cvl
+       !
        !KE43
 …
              d_deltaxt_ajs_cv(:,:,:) = 0.
 #endif
           ENDIF
+          ENDIF  !  (iflag_wake>=1)
           !        do i = 1,klon
 …
        !   ALLOCATE(lonGCM(0), latGCM(0))
        !   ALLOCATE(iGCM(0), jGCM(0))
        ENDIF
+       ENDIF  !  (iflag_con.GE.3)
+       !
        DO i=1,klon
           rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
 …
        !$OMP BARRIER
        missing_val=missing_val_omp
+       !
+       ! Now we activate some double radiation call flags only if some
+       ! diagnostics are requested, otherwise there is no point in doing this
+       IF (is_master) THEN
+         !--setting up swaero_diag to TRUE in XIOS case
+         IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
+            xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
+            xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  &
+              (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
+                                  xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  &
+            !!!--for now these fields are not in the XML files so they are omitted
+            !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
+            swaero_diag=.TRUE.
+         !--setting up swaerofree_diag to TRUE in XIOS case
+         IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
+            xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
+            xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
+            xios_field_is_active("LWupTOAcleanclr")) &
+            swaerofree_diag=.TRUE.
+         !--setting up dryaod_diag to TRUE in XIOS case
+         DO naero = 1, naero_tot-1
+          IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
+         ENDDO
+         !
+         !--setting up ok_4xCO2atm to TRUE in XIOS case
+         IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
+            xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
+            xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
+            xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
+            xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
+            xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
+            ok_4xCO2atm=.TRUE.
+       ENDIF
+       !$OMP BARRIER
+       CALL bcast(swaero_diag)
+       CALL bcast(swaerofree_diag)
+       CALL bcast(dryaod_diag)
+       CALL bcast(ok_4xCO2atm)
 #endif
+       !
        CALL printflag( tabcntr0,radpas,ok_journe, &
             ok_instan, ok_region )
+       !
+       !
+       !
        ! Prescrire l'ozone dans l'atmosphere
+       !
+       !
        !c         DO i = 1, klon
 …
 #endif
        ENDIF
+       !
        IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
 …
           SFRWL(6)=3.02191470E-02
        END SELECT
        !albedo SB <<<
 …
       ! RomP <<<
     ENDIF
+    !
     ! Ne pas affecter les valeurs entrees de u, v, h, et q
 …
     tke0(:,:)=pbl_tke(:,:,is_ave)
+    !C Risi:Nombre de traceurs de l'eau: nqo
+    !  IF (nqtot.GE.3) THEN
+    IF (nqtot.GE.(nqo+1)) THEN
+       !     DO iq = 3, nqtot
+!       DO iq = nqo+1, nqtot
+       ! CR: on modifie directement le code ici.
+       ! les isotopes ne sont pas dispatchés dans tr_seri, il faut les enlever.
+       ! on a prévu pour ça un tableau d'indice dans infotrac
+#ifdef ISOVERIF
+       write(*,*) 'physiq 1971: nqtottr=',nqtottr
+#endif
+       do itr=1,nqtottr
+         iq=itr_indice(itr)
+    IF (nqtot > nqo) THEN
+       ! water isotopes are not included in tr_seri
+       itr = 0
+       DO iq = 1, nqtot
+         IF(tracers(iq)%isH2Ofamily) CYCLE
+         itr = itr+1
 !#ifdef ISOVERIF
 !         write(*,*) 'physiq 1973: itr,iq=',itr,iq
 …
          DO  k = 1, klev
              DO  i = 1, klon
                 tr_seri(i,k,itr) = qx(i,k,iq) ! modif C Risi
+                tr_seri(i,k,itr) = qx(i,k,iq)
              ENDDO
-          ENDDO !DO  k = 1, klev
-          !write(*,*) 'physiq 1980'
-       enddo !do itr=1,nqtottr
-    ELSE !IF (nqtot.GE.(nqo+1)) THEN
-       DO k = 1, klev
-          DO i = 1, klon
-             tr_seri(i,k,1) = 0.0
           ENDDO
        ENDDO
+    ENDIF !IF (nqtot.GE.(nqo+1)) THEN
+    ELSE
+! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
+!       tr_seri(:,:,strIdx(tracers(:)%name,'H2O'//phases_sep//'g')) = 0.0
+       tr_seri(:,:,strIdx(tracers(:)%name,'H2Ov')) = 0.0
+    ENDIF
+!
 ! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
 …
     IF (debut) THEN
       WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
+!      DO iq = nqo+1, nqtot
+!           tr_ancien(:,:,iq-nqo)=tr_seri(:,:,iq-nqo)
+!      ENDDO
+       ! modif CRisi:
+       do itr=1,nqtottr
+       itr = 0
+       do iq = 1, nqtot
+         IF(tracers(iq)%isH2Ofamily) CYCLE
+         itr = itr+1
          tr_ancien(:,:,itr)=tr_seri(:,:,itr)
        enddo
 …
        d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
        ! !! RomP >>>   td dyn traceur
+       IF (nqtot.GT.nqo) THEN     ! jyg
+!          DO iq = nqo+1, nqtot      ! jyg
+          DO itr=1,nqtottr      ! C Risi modif directe
+              d_tr_dyn(:,:,itr)=(tr_seri(:,:,itr)-tr_ancien(:,:,itr))/phys_tstep ! jyg
+          ENDDO
+       ENDIF
+       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
        ! !! RomP <<<
 …
        ! !! RomP >>>   td dyn traceur
+       IF (nqtot.GT.nqo) THEN                                       ! jyg
+!          DO iq = nqo+1, nqtot                                      ! jyg
+!              d_tr_dyn(:,:,iq-nqo)= 0.0                             ! jyg
+! Modif C Risi:
+          DO itr=1,nqtottr
+                d_tr_dyn(:,:,itr)= 0.0
+          ENDDO
+       ENDIF
+       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
        ! !! RomP <<<
        ancien_ok = .TRUE.
 …
             ENDDO
          ENDDO
        ELSE !IF (iflag_wake>=1) THEN
+       ELSE
                 t_w(:,:) = t_seri(:,:)
                 q_w(:,:) = q_seri(:,:)
 …
           ! où i n'est pas un point convectif et donc ibas_con(i)=0
           ! c'est un pb indépendant des isotopes
+          if (ibas_con(i).gt.0) then
+            ema_pcb(i)  = paprs(i,ibas_con(i))
+          else ! if (ibas_con(i).gt.0) then
+              ema_pcb(i)  = 0.0
+          endif ! if (ibas_con(i).gt.0) then
+          if (ibas_con(i) > 0) then
+             ema_pcb(i)  = paprs(i,ibas_con(i))
+          else
+             ema_pcb(i)  = 0.0
+          endif
        ENDDO
        DO i = 1, klon
 …
     ENDDO
    CALL  calcratqs(klon,klev,prt_level,lunout,        &
+    CALL  calcratqs(klon,klev,prt_level,lunout,        &
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
          ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
 …
          pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
          ratqs,ratqsc,ratqs_inter)
+    !
 …
                      ZLWFT0_i, ZFLDN0, ZFLUP0, &
                      ZSWFT0_i, ZFSDN0, ZFSUP0)
           endif !ok_4xCO2atm
+          ENDIF !ok_4xCO2atm
        ENDIF ! aerosol_couple
        itaprad = 0
 …
 #endif
 ! #ifdef ISO
+    !
+    !CR: nb de traceurs eau: nqo
+    !  IF (nqtot.GE.3) THEN
+    IF (nqtot.GE.(nqo+1)) THEN
+       !     DO iq = 3, nqtot
+!       DO iq = nqo+1, nqtot ! modif C Risi
+        do itr=1,nqtottr
+         iq=itr_indice(itr)
+          DO  k = 1, klev
+             DO  i = 1, klon
+                ! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
+                d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
+             ENDDO
+    ! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
+    itr = 0
+    DO iq = 1, nqtot
+       IF(tracers(iq)%isH2Ofamily) CYCLE
+       itr = itr+1
+       DO  k = 1, klev
+          DO  i = 1, klon
+             d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
           ENDDO
        ENDDO ! !do itr=1,nqtottr
     ENDIF
+       ENDDO
+    ENDDO
+    !
     !IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
 …
     CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
     ! !! RomP >>>
+    !CR: nb de traceurs eau: nqo
+    IF (nqtot.GT.nqo) THEN
+       ! DO iq = nqo+1, nqtot ! modif C Risi
+       do itr=1,nqtottr
+          tr_ancien(:,:,itr) = tr_seri(:,:,itr)
+       ENDDO
+    ENDIF
+    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
     ! !! RomP <<<
     !==========================================================================
 …
 ! ISO
+    ! Disabling calls to the prt_alerte function
+    alert_first_call = .FALSE.
     IF (lafin) THEN
        itau_phy = itau_phy + itap

Context Navigation

Legend:

LMDZ6/trunk/libf/dyn3d/advtrac.F90

LMDZ6/trunk/libf/dyn3d/iniacademic.F90

LMDZ6/trunk/libf/dyn3d_common/infotrac.F90

LMDZ6/trunk/libf/dyn3dmem/advtrac_loc.F90

LMDZ6/trunk/libf/dyn3dmem/call_calfis_mod.F90

LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90

LMDZ6/trunk/libf/dyn3dmem/leapfrog_loc.F

LMDZ6/trunk/libf/dyn3dmem/vlspltgen_loc.F

LMDZ6/trunk/libf/dynphy_lonlat/calfis.F

LMDZ6/trunk/libf/dynphy_lonlat/calfis_loc.F

LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

LMDZ6/trunk/libf/phylmd/Dust/phys_output_write_spl_mod.F90

LMDZ6/trunk/libf/phylmd/Dust/phytracr_spl_mod.F90

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_5wv_rrtm.F90

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_6bands_rrtm.F90

LMDZ6/trunk/libf/phylmd/Dust/splaeropt_lw_rrtm.F90

LMDZ6/trunk/libf/phylmd/Dust/splaerosol_optic_rrtm.F90

LMDZ6/trunk/libf/phylmd/infotrac_phy.F90

LMDZ6/trunk/libf/phylmd/init_be.F90

LMDZ6/trunk/libf/phylmd/phyetat0.F90

LMDZ6/trunk/libf/phylmd/phyredem.F90

LMDZ6/trunk/libf/phylmd/phys_output_mod.F90

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

LMDZ6/trunk/libf/phylmd/phytrac_mod.F90

LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90

LMDZ6/trunk/libf/phylmdiso/phyetat0.F90

LMDZ6/trunk/libf/phylmdiso/phyredem.F90

LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

Download in other formats: