Index: LMDZ6/branches/Ocean_skin/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmdiso/physiq_mod.F90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/phylmdiso/physiq_mod.F90	(revision 4369)
@@ -39,7 +39,7 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, types_trac, nqCO2
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, types_trac
     USE readTracFiles_mod, ONLY: addPhase
-    USE strings_mod,  ONLY: strIdx, strStack, int2str
+    USE strings_mod,  ONLY: strIdx
     USE iophy
     USE limit_read_mod, ONLY : init_limit_read
@@ -56,5 +56,5 @@
     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
@@ -69,14 +69,16 @@
     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
+    USE surface_data,     ONLY : type_ocean, ok_veget
+    USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time
     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
+#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
 
@@ -97,8 +99,16 @@
 
 
+#ifdef INCA
+    USE geometry_mod,      ONLY: longitude, latitude, boundslon, boundslat, ind_cell_glo
+    USE time_phylmdz_mod,  ONLY: ndays
+    USE infotrac_phy,      ONLY: nqCO2
+#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
+    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
+#if defined INCA || defined REPROBUS
+    USE time_phylmdz_mod,    ONLY: annee_ref, day_ini, day_ref, start_time
+    USE vertical_layers_mod, ONLY: aps, bps, ap, bp
 #endif
 
@@ -106,5 +116,5 @@
 #ifdef CPP_RRTM
     USE YOERAD, ONLY : NRADLP
-    USE YOESW, ONLY : RSUN
+!    USE YOESW, ONLY : RSUN
 #endif
 
@@ -116,6 +126,7 @@
 
 #ifdef CPP_XIOS
-    USE xios, ONLY: xios_update_calendar, xios_context_finalize, &
-            xios_get_field_attr, xios_field_is_active
+    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
@@ -180,6 +191,6 @@
        d_t_ajsb,d_q_ajsb, &
        d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
-       d_t_ajs_w,d_q_ajs_w, &
-       d_t_ajs_x,d_q_ajs_x, &
+!       d_t_ajs_w,d_q_ajs_w, &
+!       d_t_ajs_x,d_q_ajs_x, &
        !
        d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
@@ -194,5 +205,5 @@
        d_ts, &
        !
-       d_t_oli,d_u_oli,d_v_oli, &
+!       d_t_oli,d_u_oli,d_v_oli, &
        d_t_oro,d_u_oro,d_v_oro, &
        d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
@@ -366,5 +377,4 @@
 #endif
        !
-
 
     IMPLICIT NONE
@@ -579,5 +589,4 @@
     !
     !
-    INTEGER debug
     INTEGER n
     !ym      INTEGER npoints
@@ -636,5 +645,5 @@
     ! Upmost level reached by deep convection and related variable (jyg)
     !
-    INTEGER izero
+!    INTEGER izero
     INTEGER k_upper_cv
     !------------------------------------------------------------------
@@ -814,5 +823,4 @@
     REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
     ! RomP <<<
-    REAL          :: calday
 
     !IM cf FH pour Tiedtke 080604
@@ -934,5 +942,5 @@
     real zqsat(klon,klev)
     !
-    INTEGER i, k, iq, j, nsrf, ll, l, itr
+    INTEGER i, k, iq, nsrf, l, itr
 #ifdef ISO
     real zxt_apres(ntraciso,klon)
@@ -1056,4 +1064,5 @@
 
     REAL picefra(klon,klev)
+    REAL zrel_mount(klon)
     !IM cf. AM 081204 END
     !
@@ -1265,4 +1274,5 @@
 
 #ifdef INCA
+    REAL :: calday, zxsnow_dummy(klon)
     ! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
     REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
@@ -1313,5 +1323,9 @@
     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 
+    IF (.NOT. debut .AND. is_omp_master) THEN
+       CALL wxios_set_context()
+       CALL xios_update_calendar(itap+1)
+    ENDIF
 #endif
 
@@ -1403,4 +1417,16 @@
           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)
@@ -1847,13 +1873,15 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Nouvelle initialisation pour le rayonnement RRTM
-       !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
 
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL wake_ini(rg,rd,rv,prt_level)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
        !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2151,5 +2179,5 @@
        !c         ENDDO
        !
-       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN  ! ModThL
+       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN ! ModThL
 #ifdef INCA
           CALL VTe(VTphysiq)
@@ -2158,31 +2186,64 @@
           WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
-          CALL chemini(  &
-               rg, &
-               ra, &
-               cell_area, &
-               latitude_deg, &
-               longitude_deg, &
-               presnivs, &
-               calday, &
-               klon, &
-               nqtot, &
-               nqo+nqCO2, &
-               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
 
 
@@ -2286,4 +2347,5 @@
 
 
+
     ENDIF
     !
@@ -2390,9 +2452,11 @@
           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)
-          ELSE IF (nqo.eq.4) THEN
+             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)
@@ -2658,6 +2722,6 @@
        ! !! RomP >>>   td dyn traceur
        IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
+       ! !! RomP <<<
        d_rneb_dyn(:,:)=0.0
-       ! !! RomP <<<
        ancien_ok = .TRUE.
     ENDIF
@@ -4608,4 +4672,9 @@
     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, &
@@ -4619,4 +4688,5 @@
 
     ELSE
+
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
@@ -4690,5 +4760,5 @@
        ENDDO
     ENDDO
-    IF (nqo==3) THEN
+    IF (nqo >= 3) THEN
     DO k = 1, klev
        DO i = 1, klon
@@ -5732,6 +5802,9 @@
        DO i=1,klon
           itest(i)=0
-          !        IF ((zstd(i).gt.10.0)) THEN
-          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
+          zrel_mount(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
+          !zrel_mount: relative mountain height wrt relief explained by mean slope
+          ! -> condition on zrel_mount can deactivate the drag on tilted planar terrains
+          !    such as ice sheets (work by V. Wiener)
+          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0).AND.(zrel_mount(i).GE.zrel_mount_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -5786,5 +5859,9 @@
        DO i=1,klon
           itest(i)=0
-          IF ((zpic(i)-zmea(i)).GT.100.) THEN
+          !zrel_mount: relative mountain height wrt relief explained by mean slope
+          ! -> condition on zrel_mount can deactivate the lifting on tilted planar terrains
+          !    such as ice sheets (work by V. Wiener)
+          zrel_mount(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
+          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zrel_mount(i).GE.zrel_mount_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6273,5 +6350,5 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-       IF (nqo .EQ. 3) THEN 
+       IF (nqo >= 3) THEN 
           ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
        ELSE 
@@ -6360,5 +6437,5 @@
     ! 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, &
+    CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
                  ue, ve, uq, vq, uwat, vwat)
     !
@@ -6366,6 +6443,5 @@
     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)
        !
@@ -6445,5 +6521,5 @@
             pphis, &
             zx_rh, &
-            aps, bps, ap, bp)
+            aps, bps, ap, bp, lafin)
 
        CALL VTe(VTinca)
@@ -6452,4 +6528,9 @@
     ENDIF
 
+    IF (ANY(types_trac == 'repr')) THEN
+#ifdef REPROBUS
+        CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
+#endif
+    ENDIF
 
     !
@@ -6475,9 +6556,9 @@
           d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
           !CR: on ajoute le contenu en glace
-          IF (nqo.ge.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.eq.4) THEN
+          IF (nqo == 4) THEN
              d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
           ENDIF
@@ -6501,11 +6582,5 @@
       enddo ! DO k = 1, klev
     enddo !do ixt=1,ntraciso
-!#ifdef ISOVERIF
-!        write(*,*) 'physiq 6120: d_qx(1,1,:)=',d_qx(1,1,:)
-!        write(*,*) 'qx(1,1,:)=',qx(1,1,:)
-!        write(*,*) 'xt_seri(:,1,1)=',xt_seri(:,1,1)
-!#endif
-#endif
-! #ifdef ISO
+#endif
     ! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
     itr = 0
@@ -6567,4 +6642,5 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
+    rneb_ancien(:,:) = rneb_seri(:,:)
 #ifdef ISO
     xt_ancien(:,:,:)=xt_seri(:,:,:)
@@ -6826,4 +6902,13 @@
 #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
