Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4050)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4056)
@@ -39,7 +39,7 @@
     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
@@ -61,5 +61,6 @@
     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
@@ -123,6 +124,5 @@
 #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, &
@@ -188,6 +188,6 @@
        !
        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, &
        !
@@ -262,6 +262,4 @@
        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, &
@@ -939,6 +937,5 @@
     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)
@@ -1133,5 +1130,5 @@
 !JLD    REAL zstophy, zout
 
-    CHARACTER*20 modname
+    CHARACTER (LEN=20) :: modname='physiq_mod'
     CHARACTER*80 abort_message
     LOGICAL, SAVE ::  ok_sync, ok_sync_omp
@@ -1306,4 +1303,7 @@
     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
@@ -1424,5 +1424,4 @@
     forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
 
-    modname = 'physiq'
 
     IF (debut) THEN
@@ -1853,6 +1852,10 @@
 
        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
@@ -1899,5 +1902,5 @@
             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
@@ -1949,7 +1952,12 @@
           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
@@ -2004,5 +2012,5 @@
              d_deltaxt_ajs_cv(:,:,:) = 0.
 #endif
-          ENDIF
+          ENDIF  !  (iflag_wake>=1)
 
           !        do i = 1,klon
@@ -2015,6 +2023,6 @@
        !   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)
@@ -2085,14 +2093,51 @@
        !$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
@@ -2150,4 +2195,5 @@
 #endif
        ENDIF
+       !
        IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
@@ -2198,6 +2244,4 @@
           SFRWL(6)=3.02191470E-02
        END SELECT
-
-
        !albedo SB <<<
 
@@ -2331,5 +2375,4 @@
       ! RomP <<<
     ENDIF
-
     !
     ! Ne pas affecter les valeurs entrees de u, v, h, et q
@@ -2408,17 +2451,10 @@
 
     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
@@ -2427,17 +2463,13 @@
          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
@@ -2445,9 +2477,8 @@
     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
@@ -2519,10 +2550,5 @@
        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 <<<
 
@@ -2627,12 +2653,5 @@
 
        ! !! 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.
@@ -3363,5 +3382,5 @@
             ENDDO
          ENDDO
-       ELSE !IF (iflag_wake>=1) THEN
+       ELSE
                 t_w(:,:) = t_seri(:,:)
                 q_w(:,:) = q_seri(:,:)
@@ -3752,10 +3771,9 @@
           ! 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
@@ -4484,5 +4502,5 @@
     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, &
@@ -4492,5 +4510,4 @@
          pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
          ratqs,ratqsc,ratqs_inter)
-
 
     !
@@ -5608,5 +5625,5 @@
                      ZLWFT0_i, ZFLDN0, ZFLUP0, &
                      ZSWFT0_i, ZFSDN0, ZFSUP0)
-          endif !ok_4xCO2atm
+          ENDIF !ok_4xCO2atm
        ENDIF ! aerosol_couple
        itaprad = 0
@@ -6485,20 +6502,15 @@
 #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
@@ -6558,11 +6570,5 @@
     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 <<<
     !==========================================================================
@@ -6795,4 +6801,7 @@
 ! ISO
 
+    ! Disabling calls to the prt_alerte function
+    alert_first_call = .FALSE.
+    
     IF (lafin) THEN
        itau_phy = itau_phy + itap
