Changeset 6296 for CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/LMDZ

Timestamp:

12/22/22 16:17:21 (19 months ago)

Author:

aclsce

Message:

IPSLCM7 : updated component revisions accordingly with ICOLMDZOR_v7 for DYNAMICO, ICOSA_LMDZ, ORCHIDEE and LMDZ.

Location:

CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/LMDZ

Files:

• create_etat0_limit_unstruct.F90 (deleted)
• create_etat0_unstruct.F90 (deleted)
• oasis.F90 (modified) (10 diffs)
• physiq_mod.F90 (modified) (92 diffs)

CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/LMDZ/oasis.F90

@@ -24 +24 @@
   USE mod_prism_get_proto
   USE mod_prism_put_proto
+#endif
+#ifdef CPP_CPLOCNINCA
+  USE incaoasis, ONLY : inforcv
 #endif
 #endif
@@ -59 +62 @@
   INTEGER, PARAMETER :: ids_qraiic = 28
   INTEGER, PARAMETER :: ids_qsnoic = 29
-  INTEGER, PARAMETER :: maxsend    = 29  ! Maximum number of fields to send
+  INTEGER, PARAMETER :: ids_delta_sst = 30, ids_delta_sal = 31
+  
+  INTEGER, PARAMETER :: maxsend    = 31  ! Maximum number of fields to send
   
   ! Id for fields received from ocean
+
   INTEGER, PARAMETER :: idr_sisutw = 1
   INTEGER, PARAMETER :: idr_icecov = 2
@@ -70 +76 @@
   INTEGER, PARAMETER :: idr_curenz = 7
   INTEGER, PARAMETER :: idr_oceco2 = 8
-  INTEGER, PARAMETER :: maxrecv    = 8  ! Maximum number of fields to receive
-  
+
+  INTEGER, PARAMETER :: idr_sss = 9
+  ! bulk salinity of the surface layer of the ocean, in ppt
+
+  INTEGER, PARAMETER :: maxrecv    = 9  ! Maximum number of fields to receive
+  
+#ifdef CPP_CPLOCNINCA
+  INTEGER, PARAMETER :: idr_ocedms = 1
+  INTEGER, PARAMETER :: maxrcv = 1
+#endif
 
   TYPE, PUBLIC ::   FLD_CPL            ! Type for coupling field information
@@ -110 +124 @@
     USE geometry_mod, ONLY: ind_cell_glo                    
     USE mod_phys_lmdz_mpi_data, ONLY: klon_mpi_para_nb
-    
-
+    use config_ocean_skin_m, only: activate_ocean_skin
 
 ! Local variables
@@ -127 +140 @@
     CHARACTER (len = 80)               :: abort_message 
     LOGICAL, SAVE                      :: cpl_current_omp
-    INTEGER, DIMENSION(klon_mpi)           :: ind_cell_glo_mpi
+    INTEGER, DIMENSION(klon_mpi)       :: ind_cell_glo_mpi
 
 !*    1. Initializations
@@ -164 +177 @@
 
 !************************************************************************************
-! Gather global index to be used for oasis decomposition 
+! Gather global index to be used for oasis decomposition
 !************************************************************************************
     CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
@@ -189 +202 @@
     infosend(ids_calvin)%action = .TRUE. ; infosend(ids_calvin)%name = 'COCALVIN'
     
+    if (activate_ocean_skin == 2) then
+       infosend(ids_delta_sst)%action = .TRUE.
+       infosend(ids_delta_sst)%name = 'CODELSST'
+       infosend(ids_delta_sal)%action = .TRUE.
+       infosend(ids_delta_sal)%name = 'CODELSSS'
+    end if
+           
     IF (version_ocean=='nemo') THEN
         infosend(ids_shftot)%action = .TRUE. ; infosend(ids_shftot)%name = 'COQSRMIX'
@@ -225 +245 @@
    inforecv(idr_icealw)%action = .TRUE. ; inforecv(idr_icealw)%name = 'SIICEALW'
    inforecv(idr_icetem)%action = .TRUE. ; inforecv(idr_icetem)%name = 'SIICTEMW'
+
+   if (activate_ocean_skin >= 1) then
+      inforecv(idr_sss)%action = .TRUE.
+      inforecv(idr_sss)%name = 'SISUSALW'
+   end if
    
    IF (cpl_current ) THEN
@@ -235 +260 @@
        inforecv(idr_oceco2)%action = .TRUE. ; inforecv(idr_oceco2)%name = 'SICO2FLX'
    ENDIF
+#ifdef CPP_CPLOCNINCA
+       inforcv(idr_ocedms)%action = .TRUE. ; inforcv(idr_ocedms)%name = 'SIDMSFLX'
+#endif
 
 !************************************************************************************
@@ -319 +347 @@
        ENDIF
     END DO
-    
+
+! Now, if also coupling CPL with INCA, initialize here fields to be exchanged.
+#ifdef CPP_CPLOCNINCA
+    DO jf=1,maxrcv 
+       IF (inforcv(jf)%action) THEN
+          CALL prism_def_var_proto(inforcv(jf)%nid, inforcv(jf)%name, il_part_id, &
+               il_var_nodims, PRISM_In, il_var_actual_shape, il_var_type, &
+               ierror)
+          IF (ierror .NE. PRISM_Ok) THEN
+             WRITE(lunout,*) 'inicma : Problem with prism_def_var_proto for field : ',&
+                  inforcv(jf)%name
+             abort_message=' Problem in call to prism_def_var_proto for fields to receive'
+             CALL abort_physic(modname,abort_message,1)
+          ENDIF
+       ENDIF
+    END DO
+#endif
+ 
 !************************************************************************************
 ! Atmospheric Fields to send

CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/LMDZ/physiq_mod.F90

@@ -1 +1 @@
 !
-! $Id: physiq_mod.F90 3666 2020-04-20 10:13:34Z lfalletti $
+! $Id: physiq_mod.F90 4298 2022-10-17 08:15:06Z pcadule $
 !
 !#define IO_DEBUG
@@ -16 +16 @@
        d_u, d_v, d_t, d_qx, d_ps)
 
+! For clarity, the "USE" section is now arranged in alphabetical order,
+! with a separate section for CPP keys
+! PLEASE try to follow this rule 
+
+    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
+    USE aero_mod
+    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
+  &      fl_ebil, fl_cor_ebil
     USE assert_m, only: assert
+    USE change_srf_frac_mod
+    USE conf_phys_m, only: conf_phys
+    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
+    USE CFMIP_point_locations   ! IM stations CFMIP
+    USE cmp_seri_mod
+    USE dimphy
+    USE etat0_limit_unstruct_mod
+    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
+    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
+    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg, longitude,latitude, &
+         boundslon,boundslat, dx, dy, ind_cell_glo
     USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
          histwrite, ju2ymds, ymds2ju, getin
-    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE indice_sol_mod
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, types_trac, nqCO2
+    USE readTracFiles_mod, ONLY: addPhase
+    USE strings_mod,  ONLY: strIdx
+    USE iophy
+    USE limit_read_mod, ONLY : init_limit_read
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
+    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
+    USE mod_phys_lmdz_para
+    USE netcdf95, only: nf95_close
+    USE netcdf, only: nf90_fill_real     ! IM for NMC files
+    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
+    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
+    USE pbl_surface_mod, ONLY : pbl_surface
+    USE phyaqua_mod, only: zenang_an
+    USE phystokenc_mod, ONLY: offline, phystokenc
     USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
-         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
+         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend
+!!  USE phys_local_var_mod, ONLY : a long list of variables
+!!              ==> see below, after "CPP Keys" section
+    USE phys_state_var_mod ! Variables sauvegardees de la physique
+    USE phys_output_mod
+    USE phys_output_ctrlout_mod
+    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
+    USE radlwsw_m, only: radlwsw
+    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
+    USE regr_pr_time_av_m, only: regr_pr_time_av
+    USE surface_data,     ONLY : type_ocean, ok_veget, landice_opt
+    USE time_phylmdz_mod, only: annee_ref, current_time, day_ini, day_ref, &
+          day_step_phy, itau_phy, pdtphys, raz_date, start_time, update_time, ndays
+    USE tracinca_mod, ONLY: config_inca
+    USE tropopause_m,     ONLY: dyn_tropopause
+    USE ice_sursat_mod,  ONLY: flight_init, airplane
+    USE vampir
+    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
     USE write_field_phy
-    USE dimphy
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
-    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
-    USE mod_phys_lmdz_para
-    USE iophy
-    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
-    USE phystokenc_mod, ONLY: offline, phystokenc
-    USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
-    USE vampir
-    USE pbl_surface_mod, ONLY : pbl_surface
-    USE change_srf_frac_mod
-    USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
-    USE tropopause_m,     ONLY: dyn_tropopause
+#ifdef CPP_XIOS
+    USE wxios, ONLY: g_ctx, wxios_set_context
+#endif
+    USE lscp_mod, ONLY : lscp
+    USE wake_ini_mod, ONLY : wake_ini
+    USE thermcell_ini_mod, ONLY : thermcell_ini
+
+    !USE cmp_seri_mod
+!    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
+!  &      fl_ebil, fl_cor_ebil
+
+!!!!!!!!!!!!!!!!!! "USE" section for CPP keys !!!!!!!!!!!!!!!!!!!!!!!!
+!
+!
 #ifdef CPP_Dust
-    USE phytracr_spl_mod, ONLY: phytracr_spl
+    USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
+    USE phys_output_write_spl_mod
+#else
+    USE phytrac_mod, ONLY : phytrac_init, phytrac
+    USE phys_output_write_mod
 #endif
+
+
+#ifdef REPROBUS
+    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
+         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
+         ztrop, gravit,itroprep, Z1,Z2,fac,B
+#endif
+
+
+#ifdef CPP_RRTM
+    USE YOERAD, ONLY : NRADLP
+    USE YOESW, ONLY : RSUN
+#endif
+
+
 #ifdef CPP_StratAer
     USE strataer_mod, ONLY: strataer_init
 #endif
-    USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
+
+
+#ifdef CPP_XIOS
+    USE xios, ONLY: xios_update_calendar, xios_context_finalize
+    USE xios, ONLY: xios_get_field_attr, xios_field_is_active, xios_context
+    USE xios, ONLY: xios_set_current_context
+    USE wxios, ONLY: missing_val, missing_val_omp
+#endif
+#ifndef CPP_XIOS
+    USE paramLMDZ_phy_mod
+#endif
+!
+!
+!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
+
+USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
        ! [Variables internes non sauvegardees de la physique]
        ! Variables locales pour effectuer les appels en serie
-       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
+       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
        ! Dynamic tendencies (diagnostics)
-       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
+       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
        d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
        ! Physic tendencies
@@ -66 +156 @@
        !
        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, &
        !
@@ -119 +211 @@
        cdragm, cdragh,                   &
        zustar, zu10m, zv10m, rh2m, qsat2m, &
-       zq2m, zt2m, weak_inversion, &
-       zq2m_cor,zt2m_cor,zu10m_cor,zv10m_cor, & ! pour corriger d'un bug
-       zrh2m_cor,zqsat2m_cor, &
+       zq2m, zt2m, zn2mout, weak_inversion, &
        zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
        t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
@@ -134 +224 @@
        zxrunofflic,                            &
        zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
+       delta_qsurf,                            &
        rain_lsc, rain_num,                     &
        !
@@ -139 +230 @@
        zxfluxlat_x, zxfluxlat_w, &
        !
-       d_t_vdf_x, d_t_vdf_w, &
-       d_q_vdf_x, d_q_vdf_w, &
-       pbl_tke_input, &
+       pbl_tke_input, tke_dissip, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
        cdragh_x, cdragh_w, &
@@ -168 +257 @@
        alp_bl_stat, n2, s2,  &
        proba_notrig, random_notrig,  &
-       cv_gen,  &
+!!       cv_gen,  &  !moved to phys_state_var_mod
        !
        dnwd0,  &
@@ -195 +284 @@
        ref_liq, ref_ice, theta,  &
        ref_liq_pi, ref_ice_pi,  &
-       zphi, zx_rh,  &
+       zphi, zx_rh, zx_rhl, zx_rhi,  &
        pmfd, pmfu,  &
        !
@@ -216 +305 @@
        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
        !
-    USE phys_state_var_mod ! Variables sauvegardees de la physique
-#ifdef CPP_Dust
-    USE phys_output_write_spl_mod
-#else
-    USE phys_output_var_mod ! Variables pour les ecritures des sorties
-#endif
-
-    USE phys_output_write_mod
-    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
-    USE phys_output_mod
-    USE phys_output_ctrlout_mod
-    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
-    USE regr_pr_time_av_m, only: regr_pr_time_av
-    USE netcdf95, only: nf95_close
-    !IM for NMC files
-    USE netcdf, only: nf90_fill_real
-    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
-    USE aero_mod
-    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
-    USE conf_phys_m, only: conf_phys
-    USE radlwsw_m, only: radlwsw
-    USE phyaqua_mod, only: zenang_an
-    USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
-         start_time, pdtphys, day_ini
-    USE tracinca_mod, ONLY: config_inca
-#ifdef CPP_XIOS
-    USE wxios, ONLY: missing_val, missing_val_omp
-    USE xios, ONLY: xios_get_field_attr, xios_field_is_active
-#endif
-#ifdef REPROBUS
-    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
-         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
-         ztrop, gravit,itroprep, Z1,Z2,fac,B
-#endif
-    USE indice_sol_mod
-    USE phytrac_mod, ONLY : phytrac_init, phytrac
-    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
-
-#ifdef CPP_RRTM
-    USE YOERAD, ONLY : NRADLP
-    USE YOESW, ONLY : RSUN
-#endif
-    USE ioipsl_getin_p_mod, ONLY : getin_p
-
-#ifndef CPP_XIOS
-    USE paramLMDZ_phy_mod
-#endif
-
-    USE cmp_seri_mod
-    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
-  &      fl_ebil, fl_cor_ebil
-
-    !IM stations CFMIP
-    USE CFMIP_point_locations
-    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
-    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
-    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
-    USE etat0_limit_unstruct_mod
-#ifdef CPP_XIOS
-    USE xios, ONLY: xios_update_calendar, xios_context_finalize
-#endif
-    USE limit_read_mod, ONLY : init_limit_read
-    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
-    USE readaerosol_mod, ONLY : init_aero_fromfile
-    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
 
     IMPLICIT NONE
@@ -336 +360 @@
     include "dimsoil.h"
     include "clesphys.h"
-    include "thermcell.h"
+    include "alpale.h"
     include "dimpft.h"
     !======================================================================
     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/strato
+    !$OMP THREADPRIVATE(flag_volc_surfstrat)
     LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
     PARAMETER (ok_cvl=.TRUE.)
@@ -421 +447 @@
     !======================================================================
     !
-    INTEGER ivap          ! indice de traceurs pour vapeur d'eau
-    PARAMETER (ivap=1)
-    INTEGER iliq          ! indice de traceurs pour eau liquide
-    PARAMETER (iliq=2)
-    !CR: on ajoute la phase glace
-    INTEGER isol          ! indice de traceurs pour eau glace
-    PARAMETER (isol=3)
+    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional)
+    INTEGER,SAVE :: ivap, iliq, isol, irneb
+!$OMP THREADPRIVATE(ivap, iliq, isol, irneb)
     !
     !
@@ -469 +491 @@
     REAL dtadd(klon,klev)
 
+!#ifdef CPP_XIOS
+!    TYPE(xios_context), SAVE :: g_ctx
+!#endif
+
 #ifndef CPP_XIOS
     REAL, SAVE :: missing_val=nf90_fill_real
@@ -602 +628 @@
                                                         ! gust-front in the grid cell.
     !$OMP THREADPRIVATE(iflag_alp_wk_cond)
-
-    INTEGER,  SAVE               :: iflag_bug_t2m_ipslcm61=1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_ipslcm61)
-    INTEGER,  SAVE               :: iflag_bug_t2m_stab_ipslcm61=-1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_stab_ipslcm61)
 
     REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
@@ -842 +863 @@
     real zqsat(klon,klev)
     !
-    INTEGER i, k, iq, j, nsrf, ll, l
+    INTEGER i, k, iq, j, nsrf, ll, l, itr
     !
     REAL t_coup
@@ -949 +970 @@
     !IM cf. AM 081204 BEG
     LOGICAL ptconvth(klon,klev)
+
+    REAL picefra(klon,klev)
     !IM cf. AM 081204 END
     !
@@ -1020 +1043 @@
 !JLD    REAL zstophy, zout
 
-    CHARACTER*20 modname
+    CHARACTER (LEN=20) :: modname='physiq_mod'
     CHARACTER*80 abort_message
     LOGICAL, SAVE ::  ok_sync, ok_sync_omp
@@ -1172 +1195 @@
     integer iostat
 
+    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
     REAL zzz
     !albedo SB >>>
@@ -1186 +1210 @@
     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
@@ -1196 +1223 @@
     phys_tstep=NINT(pdtphys)
 #ifdef CPP_XIOS
-    IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
+! switch to XIOS LMDZ physics context 
+!!!!$OMP MASTER
+!!!!    WRITE(*,*)'PHYSICS XIOS Context :', g_ctx
+!!!!    CALL wxios_set_context()
+!!!!$OMP END MASTER
+    IF (.NOT. debut .AND. is_omp_master) THEN
+       CALL wxios_set_context()
+       CALL xios_update_calendar(itap+1)
+    ENDIF
 #endif
 
@@ -1230 +1265 @@
 
     IF (first) THEN 
+       ivap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
+       iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
+       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
+       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
+!       CALL init_etat0_limit_unstruct
+!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
        !CR:nvelles variables convection/poches froides
 
@@ -1245 +1286 @@
             fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
             iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
-            chemistry_couple, &
-            flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
+            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
+            chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
             flag_bc_internal_mixture, bl95_b0, bl95_b1, &
                                 ! nv flags pour la convection et les
@@ -1264 +1304 @@
 #endif
 
+       !!CALL flight_init
+
        print*, '================================================='
        !
@@ -1269 +1311 @@
        IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
           WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
-               '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
+               '(H2O_g, H2O_l, H2O_s) but nqo=', nqo, '. Might as well stop here.'
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
+          WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
+          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
+               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
+          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
+          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1303 +1370 @@
     forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
 
-    modname = 'physiq'
 
     IF (debut) THEN
@@ -1314 +1380 @@
        tau_gl=86400.*tau_gl
        WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
-
-       iflag_bug_t2m_ipslcm61 = 1
-       CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
-       iflag_bug_t2m_stab_ipslcm61 = -1
-       CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
 
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
@@ -1347 +1408 @@
        iflag_phytrac = 1 ! by default we do want to call phytrac
        CALL getin_p('iflag_phytrac',iflag_phytrac)
+#ifdef CPP_Dust
+       IF (iflag_phytrac.EQ.0) THEN 
+         WRITE(lunout,*) 'In order to run with SPLA, iflag_phytrac will be forced to 1'
+         iflag_phytrac = 1
+       ENDIF
+#endif
        nvm_lmdz = 13
        CALL getin_p('NVM',nvm_lmdz)
@@ -1404 +1471 @@
        tau_overturning_th(:)=0.
 
-       IF (type_trac == 'inca') THEN
+       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN
           ! jg : initialisation jusqu'au ces variables sont dans restart
           ccm(:,:,:) = 0.
@@ -1609 +1676 @@
 #ifdef CPP_COSP
       IF (ok_cosp) THEN
-           DO k = 1, klev
-             DO i = 1, klon
-               phicosp(i,k) = pphi(i,k) + pphis(i)
-             ENDDO
-           ENDDO
+!           DO k = 1, klev
+!             DO i = 1, klon
+!               phicosp(i,k) = pphi(i,k) + pphis(i)
+!             ENDDO
+!           ENDDO
         CALL phys_cosp(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
@@ -1631 +1698 @@
 #ifdef CPP_COSP2
         IF (ok_cosp) THEN
-           DO k = 1, klev
-             DO i = 1, klon
-               phicosp(i,k) = pphi(i,k) + pphis(i)
-             ENDDO
-           ENDDO
+!           DO k = 1, klev
+!             DO i = 1, klon
+!               phicosp(i,k) = pphi(i,k) + pphis(i)
+!             ENDDO
+!           ENDDO
           CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
@@ -1677 +1744 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Nouvelle initialisation pour le rayonnement RRTM
-       !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
-       ! Initialisation des champs dans phytrac qui sont utilisés par phys_output_write
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL wake_ini(rg,rd,rv,prt_level)
+       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
+   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       ! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+#ifdef CPP_Dust
+       ! Quand on utilise SPLA, on force iflag_phytrac=1
+       CALL phytracr_spl_out_init()
+       CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,                  &
+                                pplay, lmax_th, aerosol_couple,                 &
+                                ok_ade, ok_aie, ivap, ok_sync,                  &
+                                ptconv, read_climoz, clevSTD,                   &
+                                ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
+                                flag_aerosol, flag_aerosol_strat, ok_cdnc)
+#else
+       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1 
+       ! donc seulement dans ce cas on doit appeler phytrac_init()
        IF (iflag_phytrac == 1 ) THEN
           CALL phytrac_init()
-        ENDIF
-
+       ENDIF
        CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
                               pplay, lmax_th, aerosol_couple,                 &
@@ -1692 +1782 @@
                               ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
                               flag_aerosol, flag_aerosol_strat, ok_cdnc)
+#endif
+
 
 #ifdef CPP_XIOS
@@ -1713 +1805 @@
             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
@@ -1757 +1849 @@
           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
@@ -1806 +1903 @@
              d_s_wk(:) = 0.
              d_dens_wk(:) = 0.
-          ENDIF
+          ENDIF  !  (iflag_wake>=1)
 
           !        do i = 1,klon
@@ -1817 +1914 @@
        !   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)
@@ -1887 +1984 @@
        !$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
@@ -1904 +2038 @@
        !c         ENDDO
        !
-       IF (type_trac == 'inca') THEN
+       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN ! ModThL
 #ifdef INCA
           CALL VTe(VTphysiq)
@@ -1911 +2045 @@
           WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
-          CALL chemini(  &
-               rg, &
-               ra, &
-               cell_area, &
-               latitude_deg, &
-               longitude_deg, &
-               presnivs, &
-               calday, &
-               klon, &
-               nqtot, &
-               nqo, &
-               pdtphys, &
-               annee_ref, &
-               year_cur, &
-               day_ref,  &
-               day_ini, &
-               start_time, &
-               itau_phy, &
-               date0, &
-               io_lon, &
-               io_lat, &
-               chemistry_couple, &
-               init_source, &
-               init_tauinca, &
-               init_pizinca, &
-               init_cginca, &
-               init_ccminca)
+          call init_const_lmdz( &
+          ndays, nbsrf, is_oce,is_sic, is_ter,is_lic, calend, &
+          config_inca)
+
+          CALL init_inca_geometry( & 
+               longitude, latitude, &
+               boundslon, boundslat, &
+               cell_area, ind_cell_glo) 
+
+          if (grid_type==unstructured) THEN 
+             CALL chemini(  pplay, &
+                  nbp_lon, nbp_lat, &
+                  latitude_deg, &
+                  longitude_deg, &
+                  presnivs, &
+                  calday, &
+                  klon, &
+                  nqtot, &
+                  nqo+nqCO2, &
+                  pdtphys, &
+                  annee_ref, &
+                  year_cur, &
+                  day_ref,  &
+                  day_ini, &
+                  start_time, &
+                  itau_phy, &
+                  date0, &
+                  chemistry_couple, &
+                  init_source, &
+                  init_tauinca, &
+                  init_pizinca, &
+                  init_cginca, &
+                  init_ccminca)
+          ELSE
+             CALL chemini(  pplay, &
+                  nbp_lon, nbp_lat, &
+                  latitude_deg, &
+                  longitude_deg, &
+                  presnivs, &
+                  calday, &
+                  klon, &
+                  nqtot, &
+                  nqo+nqCO2, &
+                  pdtphys, &
+                  annee_ref, &
+                  year_cur, &
+                  day_ref,  &
+                  day_ini, &
+                  start_time, &
+                  itau_phy, &
+                  date0, &
+                  chemistry_couple, &
+                  init_source, &
+                  init_tauinca, &
+                  init_pizinca, &
+                  init_cginca, &
+                  init_ccminca, &
+                  io_lon, &
+                  io_lat)
+          ENDIF
 
 
@@ -1952 +2119 @@
 #endif
        ENDIF
-       IF (type_trac == 'repr') THEN
+       !
+       IF (ANY(types_trac == 'repr')) THEN
 #ifdef REPROBUS
           CALL chemini_rep(  &
@@ -2000 +2168 @@
           SFRWL(6)=3.02191470E-02
        END SELECT
-
-
        !albedo SB <<<
 
@@ -2040 +2206 @@
 
 
+
     ENDIF
     !
@@ -2055 +2222 @@
       ENDIF
     ENDIF
+! switch to XIOS LMDZ physics context just in case
+!$OMP MASTER
+!!!!#ifdef CPP_XIOS
+!!!!    WRITE(*,*)'PHYSICS XIOS Context :', g_ctx
+!!!!    CALL xios_set_current_context(g_ctx)
+!!!!#endif
+!$OMP END MASTER
+
     !
     !
@@ -2066 +2241 @@
 
     ! Update time and other variables in Reprobus
-    IF (type_trac == 'repr') THEN
+    IF (ANY(types_trac == 'repr')) THEN
 #ifdef REPROBUS
        CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
@@ -2124 +2299 @@
       ! RomP <<<
     ENDIF
-
     !
     ! Ne pas affecter les valeurs entrees de u, v, h, et q
@@ -2136 +2310 @@
           ql_seri(i,k) = qx(i,k,iliq)
           !CR: ATTENTION, on rajoute la variable glace
-          IF (nqo.eq.2) THEN
+          IF (nqo.EQ.2) THEN             !--vapour and liquid only
              qs_seri(i,k) = 0.
-          ELSE IF (nqo.eq.3) THEN
+             rneb_seri(i,k) = 0.
+          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
              qs_seri(i,k) = qx(i,k,isol)
+             rneb_seri(i,k) = 0.
+          ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
+             qs_seri(i,k) = qx(i,k,isol)
+             rneb_seri(i,k) = qx(i,k,irneb)
           ENDIF
        ENDDO
@@ -2155 +2334 @@
 
     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(.NOT.tracers(iq)%isInPhysics) 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,addPhase('H2O','g'))) = 0.0
     ENDIF
 !
@@ -2179 +2355 @@
     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(.NOT.tracers(iq)%isInPhysics) CYCLE
+         itr = itr+1
+         tr_ancien(:,:,itr)=tr_seri(:,:,itr)       
+       enddo
     ENDIF
     !
@@ -2212 +2391 @@
        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 <<<
+       !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
+       d_rneb_dyn(:,:)=0.0
     ELSE
        d_u_dyn(:,:)  = 0.0
@@ -2229 +2406 @@
        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 <<<
+       d_rneb_dyn(:,:)=0.0
        ancien_ok = .TRUE.
     ENDIF
@@ -2463 +2637 @@
     !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
     !   zu10m,     zv10m,   fder,
-    !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+    !   zxqsurf,   delta_qsurf,
+    !   rh2m,      zxfluxu, zxfluxv,
     !   frugs,     agesno,    fsollw,  fsolsw,
     !   d_ts,      fevap,     fluxlat, t2m,
@@ -2513 +2688 @@
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
-            zsig,      sollwdown, pphi,    cldt,      &
-            rain_fall, snow_fall, solsw,   sollw,     &
+            sollwdown,    cldt,      &
+            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &
             t_seri,    q_seri,    u_seri,  v_seri,    &
@@ -2524 +2699 @@
                                 !albedo SB <<<
             cdragh,    cdragm,  u1,    v1,            &
+            beta_aridity, &
                                 !albedo SB >>>
                                 ! albsol1,   albsol2,   sens,    evap,      &
@@ -2529 +2705 @@
                                 !albedo SB <<<
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
-            zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
+            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
             d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
                                 !nrlmd<
@@ -2550 +2726 @@
             s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
             zustar, zu10m,     zv10m,   fder, &
-            zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
+            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
             z0m, z0h,     agesno,    fsollw,  fsolsw, &
             d_ts,      fevap,     fluxlat, t2m, &
@@ -2576 +2752 @@
 !>jyg
        ENDIF
-
-!add limitation for t,q at and wind at 10m
-        if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
-          CALL borne_var_surf( klon,klev,nbsrf,                 &
-            iflag_bug_t2m_stab_ipslcm61,                        &
-            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),    &
-            ftsol,zxqsurf,pctsrf,paprs,                         &
-            t2m, q2m, u10m, v10m,                               &
-            zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor,           &
-            zrh2m_cor, zqsat2m_cor)
-        ELSE
-          zt2m_cor(:)=zt2m(:)
-          zq2m_cor(:)=zq2m(:)
-          zu10m_cor(:)=zu10m(:)
-          zv10m_cor(:)=zv10m(:)
-          zqsat2m_cor=999.999
-        ENDIF
 
        !---------------------------------------------------------------------
@@ -2797 +2956 @@
          ENDDO
        ELSE
-               t_w(:,:) = t_seri(:,:)
+                t_w(:,:) = t_seri(:,:)
                 q_w(:,:) = q_seri(:,:)
                 t_x(:,:) = t_seri(:,:)
@@ -2873 +3032 @@
           !
           !>jyg
-          IF (type_trac == 'repr') THEN
+          IF (ANY(types_trac == 'repr')) THEN
              nbtr_tmp=ntra
           ELSE
@@ -3013 +3172 @@
 
        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
@@ -3439 +3605 @@
     ! Computation of ratqs, the width (normalized) of the subrid scale 
     ! water distribution
+
+    tke_dissip_ave(:,:)=0.
+    l_mix_ave(:,:)=0.
+    wprime_ave(:,:)=0.
+
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
+          l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
+          wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
     CALL  calcratqs(klon,klev,prt_level,lunout,        &
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
          ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
-         tau_ratqs,fact_cldcon,   &
+         tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
-         paprs,pplay,q_seri,zqsat,fm_therm, &
-         ratqs,ratqsc)
-
+         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
+         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
+         ratqs,ratqsc,ratqs_inter)
 
     !
@@ -3456 +3635 @@
     ENDIF
     !
+
+    picefra(:,:)=0.
+
+    IF (ok_new_lscp) THEN
+
+    !--mise à jour de flight_m et flight_h2o dans leur module
+    IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
+      CALL airplane(debut,pphis,pplay,paprs,t_seri)
+    ENDIF
+
+    CALL lscp(phys_tstep,missing_val,paprs,pplay, &
+         t_seri, q_seri,ptconv,ratqs, &
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, & 
+         cldliq, picefra, rain_lsc, snow_lsc, &
+         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
+         frac_impa, frac_nucl, beta_prec_fisrt, &
+         prfl, psfl, rhcl,  &
+         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
+         iflag_ice_thermo, ok_ice_sursat)
+
+    ELSE
+
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
@@ -3465 +3666 @@
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
          iflag_ice_thermo)
+
+    ENDIF
     !
     WHERE (rain_lsc < 0) rain_lsc = 0.
@@ -3487 +3690 @@
        ENDDO
     ENDDO
-    IF (nqo==3) THEN
+    IF (nqo >= 3) THEN
     DO k = 1, klev
        DO i = 1, klon
@@ -3704 +3907 @@
           ENDIF
           zx_rh(i,k) = q_seri(i,k)/zx_qs
+            IF (iflag_ice_thermo .GT. 0) THEN
+          zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
+          zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
+            ENDIF
           zqsat(i,k)=zx_qs
        ENDDO
@@ -3730 +3937 @@
     ENDDO
 
-    IF (type_trac == 'inca') THEN
+    IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN ! ModThL
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -3773 +3980 @@
             nbp_lon, &
             nbp_lat-1, &
-            tr_seri, &
+            tr_seri(:,:,1+nqCO2:nbtr), &
             ftsol, &
             paprs, &
@@ -3784 +3991 @@
        CALL VTe(VTinca)
        CALL VTb(VTphysiq)
-#endif 
-    ENDIF !type_trac = inca
-    IF (type_trac == 'repr') THEN
+#endif
+    ENDIF !type_trac = inca or inco
+    IF (ANY(types_trac == 'repr')) THEN
 #ifdef REPROBUS
     !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
@@ -3956 +4163 @@
 
        IF (ok_newmicro) then
-          IF (iflag_rrtm.NE.0) THEN 
+! AI          IF (iflag_rrtm.NE.0) THEN 
+          IF (iflag_rrtm.EQ.1) THEN
 #ifdef CPP_RRTM
              IF (ok_cdnc.AND.NRADLP.NE.3) THEN
@@ -3970 +4178 @@
           ENDIF
           CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
-               paprs, pplay, t_seri, cldliq, cldfra, &
+               paprs, pplay, t_seri, cldliq, picefra, cldfra, &
                cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
                flwp, fiwp, flwc, fiwc, &
                mass_solu_aero, mass_solu_aero_pi, &
-               cldtaupi, re, fl, ref_liq, ref_ice, &
+               cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
                ref_liq_pi, ref_ice_pi)
        ELSE
           CALL nuage (paprs, pplay, &
-               t_seri, cldliq, cldfra, cldtau, cldemi, &
+               t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
                cldh, cldl, cldm, cldt, cldq, &
                ok_aie, &
@@ -4114 +4322 @@
           ! 
           !--interactive CO2 in ppm from carbon cycle
-          IF (carbon_cycle_rad.AND..NOT.debut) THEN
-            RCO2=RCO2_glo
-          ENDIF
+          IF (carbon_cycle_rad) RCO2=RCO2_glo
           !
           IF (prt_level .GE.10) THEN
@@ -4130 +4336 @@
                t_seri,q_seri,wo, &
                cldfrarad, cldemirad, cldtaurad, &
-               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
-               flag_aerosol, &
-               flag_aerosol_strat, flag_aer_feedback, &
+               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
+               flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
                tau_aero, piz_aero, cg_aero, &
                tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & 
@@ -4143 +4348 @@
                heat,heat0,cool,cool0,albpla, &
                heat_volc,cool_volc, &
-               topsw,toplw,solsw,sollw, &
+               topsw,toplw,solsw,solswfdiff,sollw, &
                sollwdown, &
                topsw0,toplw0,solsw0,sollw0, &
@@ -4217 +4422 @@
                      t_seri,q_seri,wo, &
                      cldfrarad, cldemirad, cldtaurad, &
-                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
-                     flag_aerosol, &
-                     flag_aerosol_strat, flag_aer_feedback, &
+                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
+                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
                      tau_aero, piz_aero, cg_aero, &
                      tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
@@ -4230 +4434 @@
                      heatp,heat0p,coolp,cool0p,albplap, &
                      heat_volc,cool_volc, &
-                     topswp,toplwp,solswp,sollwp, &
+                     topswp,toplwp,solswp,solswfdiffp,sollwp, &
                      sollwdownp, &
                      topsw0p,toplw0p,solsw0p,sollw0p, &
@@ -4248 +4452 @@
                      ZLWFT0_i, ZFLDN0, ZFLUP0, &
                      ZSWFT0_i, ZFSDN0, ZFSUP0)
-          endif !ok_4xCO2atm
+          ENDIF !ok_4xCO2atm
        ENDIF ! aerosol_couple
        itaprad = 0
@@ -4629 +4833 @@
 
     CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
-
+   !
+   ! Prevent pbl_tke_w from becoming negative
+    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
+   !
 
        ENDIF
@@ -4700 +4907 @@
 #ifdef CPP_COSPV2
        IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
+!        IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
 
           IF (prt_level .GE.10) THEN
              print*,'freq_cosp',freq_cosp
           ENDIF
+           DO k = 1, klev
+             DO i = 1, klon
+               phicosp(i,k) = pphi(i,k) + pphis(i)
+             ENDDO
+           ENDDO
           mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
                  print*,'Dans physiq.F avant appel '
@@ -4755 +4968 @@
     !
 
-    IF (type_trac=='repr') THEN
+    IF (ANY(types_trac=='repr')) THEN
 !MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
 !MM                               dans Reprobus
@@ -4766 +4979 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-       ch_in(:,:) = qx(:,:,iliq)
-    ENDIF
-
-    IF (iflag_phytrac == 1 ) THEN
+       IF (nqo >= 3) THEN 
+          ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
+       ELSE 
+          ch_in(:,:) = qx(:,:,iliq)
+       ENDIF
+    ENDIF
 
 #ifdef CPP_Dust
-      CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
+    !  Avec SPLA, iflag_phytrac est forcé =1 
+    CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
                       pdtphys,ftsol,                                   &  ! I
                       t,q_seri,paprs,pplay,RHcl,                  &  ! I
@@ -4788 +5004 @@
 
 #else
-
-    CALL phytrac ( &
+    IF (iflag_phytrac == 1 ) THEN
+      CALL phytrac ( &
          itap,     days_elapsed+1,    jH_cur,   debut, &
          lafin,    phys_tstep,     u, v,     t, &
@@ -4826 +5042 @@
 
 #endif
+    ENDIF    ! (iflag_phytrac=1)
 
 #endif
-    ENDIF    ! (iflag_phytrac=1)
+    !ENDIF    ! (iflag_phytrac=1)
 
     IF (offline) THEN
@@ -4849 +5066 @@
     ! Calculer le transport de l'eau et de l'energie (diagnostique)
     !
-    CALL transp (paprs,zxtsol, &
-         t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
-         ve, vq, ue, uq, vwat, uwat)
+    CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
+                 ue, ve, uq, vq, uwat, vwat)
     !
     !IM global posePB BEG
     IF(1.EQ.0) THEN
        !
-       CALL transp_lay (paprs,zxtsol, &
-            t_seri, q_seri, u_seri, v_seri, zphi, &
+       CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
             ve_lay, vq_lay, ue_lay, uq_lay)
        !
     ENDIF !(1.EQ.0) THEN
     !IM global posePB END
+    !
     ! Accumuler les variables a stocker dans les fichiers histoire:
     !
@@ -4872 +5088 @@
     d_t_ec(:,:)=0.
     forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
-    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
+    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
          u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
          zmasse,exner,d_t_ec)
@@ -4905 +5121 @@
     ENDDO
     !
-    IF (type_trac == 'inca') THEN
+    IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -4914 +5130 @@
             pplay, &
             t_seri, &
-            tr_seri, &
+            tr_seri(:,:,1+nqCO2:nbtr), &
             nbtr, &
             paprs, &
@@ -4922 +5138 @@
             pphis, &
             zx_rh, &
-            aps, bps, ap, bp)
+            aps, bps, ap, bp, lafin)
 
        CALL VTe(VTinca)
@@ -4929 +5145 @@
     ENDIF
 
+    IF (ANY(types_trac == 'repr')) THEN
+#ifdef REPROBUS
+        CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
+#endif
+    ENDIF
 
     !
@@ -4952 +5173 @@
           d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
           !CR: on ajoute le contenu en glace
-          IF (nqo.eq.3) THEN
+          IF (nqo >= 3) THEN
              d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
+          ENDIF
+          !--ice_sursat: nqo=4, on ajoute rneb
+          IF (nqo == 4) THEN
+             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
           ENDIF
        ENDDO
     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) THEN
+       itr = 0
+       DO iq = 1, nqtot
+          IF(.NOT.tracers(iq)%isInPhysics) 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
@@ -5003 +5227 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
+    rneb_ancien(:,:) = rneb_seri(:,:)
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
     CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
     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 <<<
     !==========================================================================
@@ -5128 +5348 @@
   CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,  &
        pplay, lmax_th, aerosol_couple,                 &
-       ok_ade, ok_aie, ivap, ok_sync,         &
+       ok_ade, ok_aie, ivap, ok_sync,                  &
        ptconv, read_climoz, clevSTD,                   &
        ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
@@ -5147 +5367 @@
 #endif
 
-! On remet des variables a .false. apres un premier appel
-    IF (debut) THEN
-#ifdef CPP_XIOS
-      swaero_diag=.FALSE.
-      swaerofree_diag=.FALSE.
-      dryaod_diag=.FALSE.
-      ok_4xCO2atm= .FALSE.
-!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
-
-      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)
-!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
-#endif
-    ENDIF
-
     !====================================================================
     ! Arret du modele apres hgardfou en cas de detection d'un
@@ -5214 +5384 @@
     !
 
+    ! Disabling calls to the prt_alerte function
+    alert_first_call = .FALSE.
+    
     IF (lafin) THEN
        itau_phy = itau_phy + itap
@@ -5232 +5405 @@
 #ifdef CPP_XIOS
        IF (is_omp_master) CALL xios_context_finalize
+
+#ifdef INCA
+       if (ANY(types_trac == 'inca' )) then 
+          IF (is_omp_master .and. grid_type==unstructured) THEN 
+             CALL finalize_inca
+          ENDIF
+       endif
+#endif
+
 #endif
        WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1