Index: /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90	(revision 5197)
+++ /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90	(revision 5198)
@@ -1,53 +1,11 @@
+! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
 MODULE physiqex_mod
+
 IMPLICIT NONE
+
 CONTAINS
 
-! ------------------------------------------------------------------------------------------------
-!
-! Main authors : F. Hourdin, S. Riette, Q. Libois, A. Idelkadi
-! --------------
-!
-! Object : A version of the LMDZ physics that include PHYEX 
-! -------  Tu be used in particular to be coupled to non hydrostatic version of the dynamical
-!          core dynamico or its limited area version.
-!
-!
-! Development steps :
-! -------------------
-!
-! physiqex_mod.F90 has been introducued in phylmd as an empty driver for the physics (with only
-!    a newtonian relaxation for temperature and drag at the surface).
-! This driver has been used to introduce the PHYEX package, an externailzed physics package coming
-!    from MesoNH, representing turbulence, convection, clouds and precipitations.
-! This work has been done for a large part in Dephy sessions (Oleron, 2022, Fréjus 2023, Lège 2024).
-!
-! Fréjus 2023 : first version runing in 1D for Arm cumulus
-! Lège 2024   : first versions with a generic soil model (developped at this occasion) and ECrad
-!
-!
-! To do list
-! ---------
-! To do list pour un branchement plus propre
-! * PHYEX considère toutes les espèces microphysiques et la tke comme pronostiques, des termes de tendances 
-!   sont calculés. L'avance temporelle est faite en fin de pas de temps mais pourrait être déplacée au dessus si
-!   les tendances de ces variables passaient par l'interface
-! * PHYEX a besoin de variables avec mémoire d'un pas de temps à l'autre (en plus des variables pronostiques
-!   décrites juste au dessus). Ce sont les tableaux ALLOCATABLE. Ces variables pourraient devenir
-!   des traceurs.
-! * La variable ZDZMIN est ici calculée avec un MINVAL qui conduira à des résultats différents
-!   lorsque la distribution (noeuds/procs) sera différente. Cette variable est utile pour déterminer
-!   un critère CFL pour la sédimentation des précipitations. Il suffit donc d'avoir une valeur
-!   approchante.
-! * Certains allocatable sont en klev, d'autres en klev+2. Il est possible de changer ceci mais il faut gérer
-!   les recopies pour que les params voient des tableaux en klev+2 (en effectuant des recopies des
-!   niveaux extrêmes comme fait pour le vent par exemple)
-! * L'eau tombée en surface (précipitations, sédimentation du nuage, terme de dépôt) se trouve dans
-!   les variables ZINPRC, ZINPRR, ZINPRS, ZINPRG
-!
-!-------------------------------------------------------------------------------------------------
-
-
       SUBROUTINE physiqex (nlon,nlev, &
-     &            debut,lafin,pdtphys_, &
+     &            debut,lafin,pdtphys, &
      &            paprs,pplay,pphi,pphis,presnivs, &
      &            u,v,rot,t,qx, &
@@ -57,40 +15,10 @@
       USE dimphy, only : klon,klev
       USE infotrac_phy, only : nqtot
-      USE geometry_mod, only : latitude, cell_area, latitude_deg, longitude_deg
+      USE geometry_mod, only : latitude
 !      USE comcstphy, only : rg
       USE ioipsl, only : ymds2ju
-      !USE phys_state_var_mod, only : phys_state_var_init
-      USE phys_state_var_mod
+      USE phys_state_var_mod, only : phys_state_var_init
       USE phyetat0_mod, only: phyetat0
       USE output_physiqex_mod, ONLY: output_physiqex
-!!!!!! AI 04-2024
-      USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time
-      USE limit_read_mod, ONLY : init_limit_read
-      USE conf_phys_m, only: conf_phys
-      USE ozonecm_m, only: ozonecm
-      USE aero_mod
-      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, calend
-      USE radlwsw_m, only: radlwsw
-      USE MODD_CST
-      USE phys_output_var_mod, ONLY : phys_output_var_init, cloud_cover_sw
-      USE write_field_phy
-      USE phyaqua_mod, only: zenang_an
-!!!!!!      
-      ! PHYEX internal modules
-      USE MODE_INIT_PHYEX, ONLY: INIT_PHYEX, FILL_DIMPHYEX
-      USE MODD_DIMPHYEX,   ONLY: DIMPHYEX_t
-      USE MODD_PHYEX,      ONLY: PHYEX_t
-      USE MODI_INI_PHYEX,  ONLY: INI_PHYEX
-      USE MODI_ICE_ADJUST, ONLY: ICE_ADJUST
-      USE MODI_RAIN_ICE, ONLY: RAIN_ICE
-      USE MODI_TURB
-      USE MODI_SHALLOW_MF
-
-      USE MODD_CST
-      USE ioipsl_getin_p_mod, ONLY : getin_p
-      USE generic_soil_ini, ONLY : soil_ini
-      USE generic_soil_conduct, ONLY : soil_conduct
-      USE generic_soil_fluxes, ONLY : soil_fluxes
 
       IMPLICIT none
@@ -98,9 +26,10 @@
 ! Routine argument:
 !
+
       integer,intent(in) :: nlon ! number of atmospheric colums
       integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
       logical,intent(in) :: debut ! signals first call to physics
       logical,intent(in) :: lafin ! signals last call to physics
-      real,intent(in) :: pdtphys_ ! physics time step (s)
+      real,intent(in) :: pdtphys ! physics time step (s)
       real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
       real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
@@ -119,157 +48,6 @@
       real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
       real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
-      real :: d_qr(klon, klev), d_qs(klon, klev), d_qg(klon, klev) ! tendency for rain, snow, graupel
-      real :: d_tke(klon, klev)
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! AI 04-2024 a nettoyer
-  !
-      CHARACTER (LEN=20) :: modname='physiq_mod'
-      CHARACTER*80 abort_message
-  ! Declarations pas de temps
-      INTEGER, SAVE :: itap=0         ! compteur pour la physique
-      !$OMP THREADPRIVATE(itap)
-      INTEGER lmt_pas
-      SAVE lmt_pas                ! frequence de mise a jour
-      !$OMP THREADPRIVATE(lmt_pas)
-      REAL zdtime, zdtime1, zdtime2, zlongi
-   ! Declarations pour appel de conf_phys 
-      LOGICAL ok_hf
-      SAVE ok_hf
-      !$OMP THREADPRIVATE(ok_hf)
-      LOGICAL ok_journe ! sortir le fichier journalier
-      SAVE ok_journe
-      !$OMP THREADPRIVATE(ok_journe)
-      LOGICAL ok_mensuel ! sortir le fichier mensuel
-      SAVE ok_mensuel
-      !$OMP THREADPRIVATE(ok_mensuel)
-      LOGICAL ok_instan ! sortir le fichier instantane
-      SAVE ok_instan
-      !$OMP THREADPRIVATE(ok_instan)
-      LOGICAL ok_LES ! sortir le fichier LES
-      SAVE ok_LES
-      !$OMP THREADPRIVATE(ok_LES)
-      LOGICAL callstats ! sortir le fichier stats
-      SAVE callstats
-      !$OMP THREADPRIVATE(callstats)
-      LOGICAL ok_region ! sortir le fichier regional
-      PARAMETER (ok_region=.FALSE.)
-      REAL seuil_inversion
-      SAVE seuil_inversion
-      !$OMP THREADPRIVATE(seuil_inversion)
-      ! Parametres lies au nouveau schema de nuages (SB, PDF)
-      REAL, SAVE :: fact_cldcon
-      REAL, SAVE :: facttemps
-      !$OMP THREADPRIVATE(fact_cldcon,facttemps)
-      LOGICAL, SAVE :: ok_newmicro
-      !$OMP THREADPRIVATE(ok_newmicro)
-      INTEGER, SAVE :: iflag_cld_th
-      !$OMP THREADPRIVATE(iflag_cld_th)
-      REAL ratqsbas,ratqshaut,tau_ratqs
-      SAVE ratqsbas,ratqshaut,tau_ratqs
-      !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
-      ! ozone 
-      INTEGER,SAVE :: read_climoz                ! Read ozone climatology
-      !     (let it keep the default OpenMP shared attribute)
-      !     Allowed values are 0, 1 and 2
-      !     0: do not read an ozone climatology
-      !      1: read a single ozone climatology that will be used day and night
-      !     2: read two ozone climatologies, the average day and night
-      !     climatology and the daylight climatology
-      INTEGER,SAVE :: ncid_climoz
-      REAL zzzo
-      REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
-                 ! the ozone fields, old method.
-      !           
-      ! aerosols
-      ! params et flags aerosols
-      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
-      LOGICAL ok_alw            ! Apply aerosol LW effect or not
-      LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
-      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
-      SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
-      !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
-      LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
-      ! false : lecture des aerosol dans un fichier
-      INTEGER, SAVE :: flag_aerosol
-      !$OMP THREADPRIVATE(flag_aerosol)
-      LOGICAL, SAVE :: flag_bc_internal_mixture
-      !$OMP THREADPRIVATE(flag_bc_internal_mixture)
-      !
-      !--STRAT AEROSOL
-      INTEGER, SAVE :: flag_aerosol_strat
-      !$OMP THREADPRIVATE(flag_aerosol_strat)
-      !
-      !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
-      LOGICAL, SAVE :: flag_aer_feedback
-      !$OMP THREADPRIVATE(flag_aer_feedback)
-      !
-      LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3
-      ! false : use offline chemistry O3
-      !$OMP THREADPRIVATE(chemistry_couple)
-      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)
-      !
-      INTEGER, SAVE :: iflag_radia=0     ! active ou non le rayonnement (MPL)
-      !$OMP THREADPRIVATE(iflag_radia)
-      !
-      REAL, SAVE :: alp_offset
-      !$OMP THREADPRIVATE(alp_offset)
-! Declarations pour le rayonnement
-! Le rayonnement n'est pas calcule tous les pas, il faut donc
-! sauvegarder les sorties du rayonnement
-     ! Surface
-     REAL zxtsol(KLON)
-     INTEGER itaprad
-     SAVE itaprad
-     !$OMP THREADPRIVATE(itaprad)
-     ! Date equinoxe de printemps
-     !LOGICAL, parameter ::  ok_rad=.true. ! a nettoyer
-     REAL,SAVE ::  solarlong0
-     !$OMP THREADPRIVATE(solarlong0)
-     ! 6 bandes SW : argument non utilise dans radlwsw a nettoyer ?
-     REAL,DIMENSION(6), parameter :: SFRWL = (/ 1.28432794E-03, 0.12304168, 0.33106142, &
-                                         & 0.32870591, 0.18568763, 3.02191470E-02 /)
-     LOGICAL, parameter :: new_orbit = .TRUE.
-     INTEGER, parameter :: mth_eq=3, day_eq=21
-     REAL :: day_since_equinox, jD_eq
-     !REAL zdtime, zdtime1, zdtime2, zlongi
-     REAL dist, rmu0(klon), fract(klon), zrmu0(klon), zfract(klon)
-     CHARACTER(len=512) :: namelist_ecrad_file
-     ! 
-     REAL, dimension(klon, klev) :: zqsat
-     REAL, dimension(klon, klev) :: ref_liq, ref_ice, ref_liq_pi, ref_ice_pi ! rayons effectifs des gout
-     REAL, dimension(klon, klev) :: cldtaurad   ! epaisseur optique
-     REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
-     REAL, dimension(klon, klev) :: cldemirad   ! emissivite pour
-     REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
-     ! aerosols AI attention a les declarer ailleurs dans phylmd : phys_output ou phy_stats et non phys_local
-     REAL, dimension(klon,klev+1)  :: ZLWFT0_i, ZSWFT0_i,ZFLDN0, ZFLUP0, ZFSDN0, ZFSUP0
-     REAL, dimension(klon)       :: oplwad_aero, sollwad_aero,&
-                                    toplwai_aero, sollwai_aero, &
-                                    toplwad0_aero, sollwad0_aero
-     REAL, dimension(klon,9)       :: solsw_aero, solsw0_aero,topsw_aero, topsw0_aero
-
-     REAL, dimension(klon)       :: &
-                                    topswcf_aero, solswcf_aero, &
-                                    solswad_aero, &
-                                    solswad0_aero, &
-                                    topswai_aero, solswai_aero, &
-                                    toplwad_aero, topswad0_aero, &
-                                    topswad_aero
-     REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer
-
-     REAL,dimension(klon,klev)   :: d_t_swr, d_t_sw0, d_t_lwr, d_t_lw0
-     !
-     EXTERNAL suphel ! AI attention 
-     !
-     include "YOMCST.h"
-     include "clesphys.h"
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Fin AI
-
-    include "flux_arp.h"
-
+!    include "clesphys.h"
     INTEGER        length
     PARAMETER    ( length = 100 )
@@ -279,143 +57,16 @@
     !$OMP THREADPRIVATE(clesphy0)
 
-! Saved variables
-REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PSIGS !variance of s
-REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PCF_MF, PRC_MF, PRI_MF !shallow convection cloud
-REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ZQR, ZQS, ZQG !rain, snow, graupel specifiq contents
-REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  PTKEM       ! TKE
-TYPE(DIMPHYEX_t),SAVE    :: D
-TYPE(PHYEX_t), SAVE      :: PHYEX
-!
-INTEGER, PARAMETER       :: KRR=6, KRRL=2, KRRI=3, KSV=0
-INTEGER                  :: JRR
-! Time-State variables and Sources variables
-REAL, DIMENSION(klon,klev+2)     ::  ZUT, ZVT ! wind component on klev+2
-REAL, DIMENSION(klon,klev+2)     ::  ZPABST   ! absolute pressure
-REAL, DIMENSION(klon,klev+2,KRR) ::  ZRX,ZRXS,ZRXS0 ! qx and source of q from LMDZ to rt
-REAL, DIMENSION(klon,klev+2)     ::  ZRHOD    ! rho_dry
-REAL, DIMENSION(klon,klev+2,0)   ::  ZSVT     ! passive scal. var.
-REAL, DIMENSION(klon,klev+2)     :: PWT ! vertical wind velocity (used only for diagnostic)
-REAL, DIMENSION(klon,klev+2) ::  ZRUS,ZRVS,ZRWS,ZRTHS,ZRTKES ! sources of momentum, conservative potential temperature, Turb. Kin. Energy
-REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS !tendency
-REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS0
-REAL, DIMENSION(KLON,KLEV+2) :: ZTKES
-REAL, DIMENSION(KLON,KLEV+2) :: ZTKES0
-! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative ! mixing ratio
-REAL, DIMENSION(klon,klev+2,KRR) ::  ZRRS
-REAL, DIMENSION(klon,klev+2)   ::  PTHVREF  ! Virtual Potential Temperature of the reference state
-! Adjustment variables
-REAL, DIMENSION(KLON,KLEV+2) :: zqdm !1-qt=1/(1+rt)
-REAL, DIMENSION(KLON,KLEV+2) :: zqt !mixing ratio and total specifiq content
-REAL, DIMENSION(KLON,KLEV+2) :: ZTHETA, ZEXN !theta and exner function
-REAL, DIMENSION(KLON,KLEV+2) :: zz_mass !altitude above ground of mass points
-REAL, DIMENSION(KLON,KLEV+2) :: zz_flux !altitude above ground of flux points
-REAL, DIMENSION(KLON,KLEV+2) :: zdzm !distance between two consecutive mass points (expressed on flux points)
-REAL, DIMENSION(KLON,KLEV+2) :: zdzf !distance between two consecutive flux points (expressed on mass points)
-REAL, DIMENSION(KLON) :: zs !surface orography
-REAL, DIMENSION(KLON) :: zsigqsat !coeff for the extra term for s variance
-REAL, DIMENSION(0) :: ZMFCONV
-REAL, DIMENSION(KLON,KLEV+2) :: ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR, ZICE_CLD_WGT !used only in HIRLAM config
-REAL, DIMENSION(KLON,KLEV+2) :: ZSRC ! bar(s'rc')
-REAL, DIMENSION(KLON,KLEV+2) :: ZCLDFR !cloud fraction
-REAL, DIMENSION(KLON,KLEV+2) :: ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF !subgrid autoconversion
-! Rain_ice variables
-real :: ZDZMIN
-real, dimension(klon, klev+2) :: ZCIT !ice concentration
-real, dimension(klon) :: ZINPRC, ZINPRR, ZINPRS, ZINPRG !precipitation flux at ground
-real, dimension(klon, klev+2) :: ZEVAP3D !evaporation (diag)
-real, dimension(klon, klev+2) :: ZRAINFR
-real, dimension(klon) :: ZINDEP !deposition flux, already contained in ZINPRC
-real, dimension(klon)         :: ZSEA, ZTOWN !sea and town fractions in the frid cell
-! Turbulence variables
-REAL, DIMENSION(klon,klev+2)  ::  PDXX,PDYY,PDZX,PDZY ! metric coefficients
-REAL, DIMENSION(klon,klev+2)  ::  PZZ       !  physical distance between 2 succesive grid points along the K direction
-REAL, DIMENSION(klon)         ::  PDIRCOSXW, PDIRCOSYW, PDIRCOSZW ! Director Cosinus along x, y and z directions at surface w-point
-REAL, DIMENSION(klon)         ::  PCOSSLOPE       ! cosinus of the angle between i and the slope vector
-REAL, DIMENSION(klon)         ::  PSINSLOPE       ! sinus of the angle   between i and the slope vector
-REAL, DIMENSION(klon,klev+2)  ::  PRHODJ   ! dry density * Grid size
-REAL, DIMENSION(0,0)          ::  MFMOIST  ! moist mass flux dual scheme
-REAL, DIMENSION(0,0,0)        ::  PHGRAD      ! horizontal gradients
-REAL, DIMENSION(klon)         ::  PSFTH,PSFRV,PSFU,PSFV ! normal surface fluxes of theta, Rv, (u,v) parallel to the orography
-REAL, DIMENSION(klon,0)       ::  PSFSV ! normal surface fluxes of Scalar var. KSV=0
-REAL, DIMENSION(klon)         ::  ZBL_DEPTH  ! BL height for TOMS
-REAL, DIMENSION(klon)         ::  ZSBL_DEPTH ! SBL depth for RMC01
-REAL, DIMENSION(klon,klev+2)  ::  ZCEI ! Cloud Entrainment instability index to emphasize localy turbulent fluxes
-REAL, DIMENSION(klon,klev+2,0)::  ZRSVS ! Source terms for all passive scalar variables
-REAL, DIMENSION(klon,klev+2)  ::  ZFLXZTHVMF ! MF contribution for vert. turb. transport used in the buoy. prod. of TKE
-REAL, DIMENSION(klon,klev+2)  ::  ZWTH       ! heat flux
-REAL, DIMENSION(klon,klev+2)  ::  ZWRC       ! cloud water flux
-REAL, DIMENSION(klon,klev+2,0)::  ZWSV       ! scalar flux
-REAL, DIMENSION(klon,klev+2)  ::  ZTP        ! Thermal TKE production (MassFlux + turb)
-REAL, DIMENSION(klon,klev+2)  ::  ZDP        ! Dynamic TKE production
-REAL, DIMENSION(klon,klev+2)  ::  ZTDIFF     ! Diffusion TKE term
-REAL, DIMENSION(klon,klev+2)  ::  ZTDISS     ! Dissipation TKE term
-REAL, DIMENSION(0,0)          ::  PLENGTHM, PLENGTHH ! length scale from vdfexcu (HARMONIE-AROME)
-REAL :: ZTHVREFZIKB ! for electricity scheme interface
-REAL, DIMENSION(klon,klev+2)  ::  ZDXX,ZDYY,ZDZX,ZDZY,ZZZ
-REAL, DIMENSION(klon)         ::  ZDIRCOSXW,ZDIRCOSYW,ZDIRCOSZW,ZCOSSLOPE,ZSINSLOPE
-! Shallow variables 
-REAL, DIMENSION(klon,klev+2) ::  PDUDT_MF     ! tendency of U   by massflux scheme
-REAL, DIMENSION(klon,klev+2) ::  PDVDT_MF     ! tendency of V   by massflux scheme
-REAL, DIMENSION(klon,klev+2) ::  PDTHLDT_MF   ! tendency of thl by massflux scheme
-REAL, DIMENSION(klon,klev+2) ::  PDRTDT_MF    ! tendency of rt  by massflux scheme
-REAL, DIMENSION(klon,klev+2,KSV)::  PDSVDT_MF    ! tendency of Sv  by massflux scheme
-REAL, DIMENSION(klon,klev+2) ::  PSIGMF
-REAL, DIMENSION(klon,klev+2) ::  ZFLXZTHMF
-REAL, DIMENSION(klon,klev+2) ::  ZFLXZRMF
-REAL, DIMENSION(klon,klev+2) ::  ZFLXZUMF
-REAL, DIMENSION(klon,klev+2) ::  ZFLXZVMF
-REAL, DIMENSION(klon,klev+2) ::  PTHL_UP   ! Thl updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PRT_UP    ! Rt  updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PRV_UP    ! Vapor updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PU_UP     ! U wind updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PV_UP     ! V wind updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PRC_UP    ! cloud content updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PRI_UP    ! ice content   updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PTHV_UP   ! Thv   updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PW_UP     ! vertical speed updraft characteristics
-REAL, DIMENSION(klon,klev+2) ::  PFRAC_UP  ! updraft fraction
-REAL, DIMENSION(klon,klev+2) ::  PEMF      ! updraft mass flux
-REAL, DIMENSION(klon,klev+2) ::  PDETR     ! updraft detrainment
-REAL, DIMENSION(klon,klev+2) ::  PENTR     ! updraft entrainment
-INTEGER,DIMENSION(klon) ::IKLCL,IKETL,IKCTL ! level of LCL,ETL and CTL
-! Values after saturation adjustement
-REAL, DIMENSION(klon,klev) :: t_adj ! Adjusted temperature
-REAL, DIMENSION(klon,klev) :: qv_adj, ql_adj, qi_adj, qr_adj, qs_adj, qg_adj !specific contents after adjustement
-!
+
 real :: temp_newton(klon,klev)
 integer :: k
 logical, save :: first=.true.
 !$OMP THREADPRIVATE(first)
-!real,save :: rg=9.81
-!!$OMP THREADPRIVATE(rg)
+
+real,save :: rg=9.81
+!$OMP THREADPRIVATE(rg)
 
 ! For I/Os
 integer :: itau0
-
-real, dimension(nlon) :: swdnsrf,lwdnsrf,lwupsrf,sensible,latent,qsats,capcal,fluxgrd, fluxu,fluxv,fluxt,fluxq
-
-integer, parameter :: nsoil=11
-real, dimension(nsoil) :: zout
-real, dimension(nlon,nlev) :: tsoil_out
-
-integer :: isoil
-real,dimension(:), allocatable, save :: tsrf
-real, dimension(:,:), allocatable, save :: tsoil
-!$OMP THREADPRIVATE(tsrf,tsoil)
-
-real :: thermal_inertia0,z0m0,z0h0,z0q0,betasoil
-
-!AI
-!real, save :: rpi, stefan, karman, fsw0, timesec
-!!$OMP THREADPRIVATE(rpi, stefan, karman, fsw0, timesec)
-real, save :: stefan, fsw0, timesec
-!$OMP THREADPRIVATE(stefan, fsw0, timesec)
-integer, save :: iecri=1,iflag_surface=0
-!$OMP THREADPRIVATE(iecri,iflag_surface)
-real,save :: zjulian
-!$OMP THREADPRIVATE(zjulian)
-! AI
-!integer, save :: itap=0
-!!$OMP THREADPRIVATE(iecri,iflag_surface,rpi, stefan, karman, fsw0, timesec,tsrf, tsoil,itap)
+real :: zjulian
 
 
@@ -425,124 +76,28 @@
 
 print*,'Debut physiqex',debut
-! AI
-pdtphys=pdtphys_
-CALL update_time(pdtphys)
-phys_tstep=NINT(pdtphys)
-
 ! initializations
 if (debut) then ! Things to do only for the first call to physics 
-  print*,'Debut physiqex IN'
-  print*,'NOUVEAU PHYSIQEX 222222222222222222222'
-  CALL suphel ! initialiser constantes et parametres phys.
-  ! Read parameters and flags in files .def  
-  CALL conf_phys(ok_journe, ok_mensuel, &
-            ok_instan, ok_hf, &
-            ok_LES, &
-            callstats, &
-            solarlong0,seuil_inversion, &
-            fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
-            iflag_cld_th,ratqsbas,ratqshaut,tau_ratqs, &
-            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, &
-            read_climoz, alp_offset)
+print*,'Debut physiqex IN'
 
-    ! IF (.NOT. create_etat0_limit) 
-  CALL init_limit_read(days_elapsed)
+! load initial conditions for physics (including the grid)
+  call phys_state_var_init(1) ! some initializations, required before calling phyetat0
+  call phyetat0("startphy.nc", clesphy0, tabcntr0)
 
-  ! load initial conditions for physics (including the grid)
-    ! AI
-    !call phys_state_var_init(1) ! some initializations, required before calling phyetat0
-    call phys_state_var_init(read_climoz) ! some initializations, required before calling phyetat0
-    call phyetat0("startphy.nc", clesphy0, tabcntr0)
-    CALL phys_output_var_init
-    
-  ! AI
-  ! Init increments et pas
-    !itap    = 0
-    itaprad = 0
-    lmt_pas = NINT(86400./phys_tstep * 1.0)   ! tous les jours
-    WRITE(*,*)'La frequence de lecture surface est de ',  &
-            lmt_pas
-    IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
-          radpas = NINT( 86400./phys_tstep)/nbapp_rad
-    ELSE
-       WRITE(*,*) 'le nombre de pas de temps physique doit etre un ', &
-            'multiple de nbapp_rad'
-       WRITE(*,*) 'changer nbapp_rad ou alors commenter ce test ', &
-            'mais 1+1<>2'
-       abort_message='nbre de pas de temps physique n est pas multiple ' &
-            // 'de nbapp_rad'
-       CALL abort_physic(modname,abort_message,1)
-    ENDIF
+! Initialize outputs:
+  itau0=0
+  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
+  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
+  call ymds2ju(1979, 1, 1, 0.0, zjulian)
 
-  ! Initialize outputs:
-    itau0=0
-    ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
-    !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
-    !call ymds2ju(1979, 1, 1, 0.0, zjulian)
-    !print*,'zjulian=',zjulian
-  
-  ZDZMIN=MINVAL((pphi(:,2:) - pphi(:,1:klev-1))/9.81)
-  CALL INIT_PHYEX(pdtphys, ZDZMIN, PHYEX)
-  CALL FILL_DIMPHYEX(KLON, KLEV, D)
-
-  !  
-  ! Variables saved
-  ALLOCATE(PTKEM(klon,klev+2))
-  ALLOCATE(PSIGS(klon,klev+2))
-  ALLOCATE(PCF_MF(klon,klev+2))
-  ALLOCATE(PRC_MF(klon,klev+2))
-  ALLOCATE(PRI_MF(klon,klev+2))
-  ALLOCATE(ZQR(klon, klev))
-  ALLOCATE(ZQS(klon, klev))
-  ALLOCATE(ZQG(klon, klev))
-  PSIGS=0.
-  PCF_MF=0.
-  PRC_MF=0.
-  PRI_MF=0.
-  ZQR=0.
-  ZQS=0.
-  ZQG=0.
-  PTKEM(:,:) = PHYEX%TURBN%XTKEMIN ! TODO: init from TKE at stationnary state
-
-  ! ----------------------------------------------------------------
-  ! Initialisation du sol generique
-  ! ----------------------------------------------------------------
-  allocate(tsrf(nlon))
-  allocate(tsoil(nlon,nsoil))
-  thermal_inertia0=2000.
-  z0m0=0.05
-  z0h0=0.05
-  z0q0=0.05
-  betasoil=0.2
-  call soil_ini(klon,nsoil,PHYEX%CST%XCPD,PHYEX%CST%XLVTT,thermal_inertia0,z0m0,z0h0,z0q0,betasoil)
-!
 #ifndef CPP_IOIPSL_NO_OUTPUT
   ! Initialize IOIPSL output file
 #endif
-! AI Partie Surface depplacee
-! compute tendencies to return to the dynamics:
-! "friction" on the first layer
-!   tsrf(1:nlon)=t(1:nlon,1)
-!   do isoil=1,nsoil
-!      tsoil(1:nlon,isoil)=tsrf(1:nlon)
-!      !tsoil(1:nlon,isoil)=302.
-!   enddo
-!   rpi=2.*asin(1.)
-!   timesec=5*3600
-!   stefan=5.67e-8
-!   fsw0=900.
-!   call getin_p('iecri',iecri)
-!   call getin_p('iflag_surface',iflag_surface)
 
 endif ! of if (debut)
 
-! AI
-! Increment
-!itap   = itap + 1
 !------------------------------------------------------------
 ! Initialisations a chaque pas de temps
 !------------------------------------------------------------
+
 
 ! set all tendencies to zero
@@ -551,672 +106,28 @@
 d_t(1:klon,1:klev)=0.
 d_qx(1:klon,1:klev,1:nqtot)=0.
-d_qr(1:klon,1:klev)=0.
-d_qs(1:klon,1:klev)=0.
-d_qg(1:klon,1:klev)=0.
 d_ps(1:klon)=0.
-d_tke(1:klon,1:klev)=0.
-!
-! AI attention
-d_t_swr = 0.
-d_t_sw0 = 0.
-d_t_lwr = 0.
-d_t_lw0 = 0.
-!
-ZDXX(:,:) = 0.
-ZDYY(:,:) = 0.
-ZDZX(:,:) = 0.
-ZDZY(:,:) = 0.
-ZDIRCOSXW(:) = 1.
-ZDIRCOSYW(:) = 1.
-ZDIRCOSZW(:) = 1.
-ZCOSSLOPE(:) = 0.
-ZSINSLOPE(:) = 1.
-PHGRAD(:,:,:) = 0.
-ZBL_DEPTH(:) = 0. ! needed only with LRMC01 key (correction in the surface boundary layer)
-ZSBL_DEPTH(:) = 0.
-ZCEI(:,:) = 0.  ! needed only if HTURBLEN_CL /= 'NONE' modification of mixing lengh inside clouds
-ZSVT(:,:,:) = 0.
-PWT(:,:) = 0.
-ZUT(:,2:klev+1) = u(:,:)
-ZVT(:,2:klev+1) = v(:,:)
-!
-!------------------------------------------------------------
-! Conversions and extra levels
-!------------------------------------------------------------
-!TODO check in Meso-NH how values are extrapolated outside of the physical domain
-zqt(:,2:klev+1) = qx(:,:,1) + qx(:,:,2) + qx(:,:,3)
-zqt(:,1)=0.
-zqt(:,klev+2)=0.
-zqdm(:,:)=1.-zqt(:,:) !equal to 1/(1+rt)
-ZRX(:,2:klev+1,1) = qx(:,:,1) / zqdm(:,2:klev+1)
-ZRX(:,2:klev+1,2) = qx(:,:,2) / zqdm(:,2:klev+1)
-ZRX(:,2:klev+1,4) = qx(:,:,3) / zqdm(:,2:klev+1)
-ZRX(:,2:klev+1,3) = ZQR(:,:)
-ZRX(:,2:klev+1,5) = ZQS(:,:)
-ZRX(:,2:klev+1,6) = ZQG(:,:)
-DO JRR=1,KRR
-  CALL VERTICAL_EXTEND(ZRX(:,:,JRR),klev)
-END DO
-!
-ZEXN(:,2:klev+1) = (pplay(:,:) / PHYEX%CST%XP00) ** (PHYEX%CST%XRD/PHYEX%CST%XCPD)
-ZTHETA(:,2:klev+1) = t(:,:) / ZEXN(:,2:klev+1)
-CALL VERTICAL_EXTEND(ZEXN,klev)
-CALL VERTICAL_EXTEND(ZTHETA,klev)
-
-!TODO check in Meso-NH how zz_mass and zz_flux are initialized outside of the physical domain
-zs(:) = pphis(:)/PHYEX%CST%XG
-zz_mass(:,2:klev+1) = pphi(:,:) / PHYEX%CST%XG
-zz_mass(:,1) = 2*zs-zz_mass(:,2)
-zz_mass(:,klev+2)=2.*zz_mass(:,klev+1)-zz_mass(:,klev)
-
-do k=2, klev+2
-  zz_flux(:,k)=(zz_mass(:,k-1)+zz_mass(:,k))/2.
-enddo
-zz_flux(:,1)=2*zz_mass(:,1)-zz_flux(:,2)
-
-!zdzf is the distance between two consecutive flux points (expressed on mass points)
-do k=1,klev+1
-  zdzf(:,k)=zz_flux(:,k+1)-zz_flux(:,k)
-enddo
-zdzf(:,klev+2)=(zz_mass(:,klev+2)-zz_flux(:,klev+2))*2.
-
-!zdzm distance between two consecutive mass points (expressed on flux points)
-do k=2,klev+2
-  zdzm(:,k)=zz_mass(:,k)-zz_mass(:,k-1)
-enddo
-zdzm(:,1)=(zz_mass(:,1)-zz_flux(:,1))*2.
-
-ZPABST(:,2:klev+1) = pplay(:,:)
-ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV))
-DO k=2,klev+1
-  PRHODJ(:,k) = ZRHOD(:,k) * (zdzf(:,k)*cell_area(:))
-END DO
-PTHVREF(:,:) = ZTHETA(:,:) * (1. + PHYEX%CST%XRV/PHYEX%CST%XRD * ZRX(:,:,1)) * ZQDM(:,:)
-
-CALL VERTICAL_EXTEND(ZPABST,klev)
-CALL VERTICAL_EXTEND(PRHODJ,klev)
-CALL VERTICAL_EXTEND(ZRHOD,klev)
-CALL VERTICAL_EXTEND(ZUT,klev)
-CALL VERTICAL_EXTEND(ZVT,klev)
- 
 
 !------------------------------------------------------------
-! Tendencies
+! Calculs
 !------------------------------------------------------------
-!For Meso-NH, initialia values for the tendencies are filled with
-!a pseudo-tendecy computed by dividing the state variable by the time step
-!This mechanism enables the possibility for the different parametrisations
-!to guess the value at the end of the time step.
-!For the wind components, we could do the same way but it is not needed
-!as the parametrisations don't use the S varaible to guess the futur value of the wind.
-ZRXS(:,:,:) = ZRX(:,:,:)/pdtphys
-ZTHETAS(:,:)=ZTHETA(:,:)/pdtphys
-ZTKES(:,:)=PTKEM(:,:)/pdtphys
-!To compute the actual tendency, we save the initial values of these variables
-ZRXS0(:,:,:) = ZRXS(:,:,:)
-ZTHETAS0(:,:)=ZTHETAS(:,:)
-ZTKES0(:,:)=ZTKES(:,:)
-!------------------------------------------------------------
-! Adjustment
-!------------------------------------------------------------
-!
-ZSRC(:,:) = 0.
-ZSIGQSAT=PHYEX%NEBN%VSIGQSAT
-CALL ICE_ADJUST (D, PHYEX%CST, PHYEX%RAIN_ICE_PARAMN, PHYEX%NEBN, PHYEX%TURBN, PHYEX%PARAM_ICEN,    &
-                &PHYEX%MISC%TBUCONF, KRR,                                                           &
-                &'ADJU',                                                                            &
-                &pdtphys, ZSIGQSAT,                                                                 &
-                &PRHODJ, ZEXN, ZRHOD, PSIGS, .FALSE., zmfconv,                                      &
-                &ZPABST, ZZ_MASS,                                                                   &
-                &ZEXN, PCF_MF, PRC_MF, PRI_MF,                                                      &
-                &ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR,                                              &
-                &ZRX(:,:,1), ZRX(:,:,2), ZRXS(:,:,1), ZRXS(:,:,2), ZTHETA, ZTHETAS,                 &
-                &PHYEX%MISC%COMPUTE_SRC, ZSRC, ZCLDFR,                                              &
-                &ZRX(:,:,3), ZRX(:,:,4), ZRXS(:,:,4), ZRX(:,:,5), ZRX(:,:,6),                       &
-                &PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET,                                           &
-                &ZICE_CLD_WGT,                                                                      &
-                &PHLC_HRC=ZHLC_HRC, PHLC_HCF=ZHLC_HCF, PHLI_HRI=ZHLI_HRI, PHLI_HCF=ZHLI_HCF         )
-!
-!Variables are updated with their adjusted values (to be used by the other parametrisations)
-ZTHETA(:,:)=ZTHETAS(:,:)*pdtphys
-ZRX(:,:,:)=ZRXS(:,:,:)*pdtphys
-t_adj=ZTHETA(:,2:klev+1)*ZEXN(:,2:klev+1)
-qv_adj=ZRX(:,2:klev+1,1)*zqdm(:,2:klev+1)
-ql_adj=ZRX(:,2:klev+1,2)*zqdm(:,2:klev+1)
-qi_adj=ZRX(:,2:klev+1,4)*zqdm(:,2:klev+1)
-qr_adj=ZRX(:,2:klev+1,3)*zqdm(:,2:klev+1)
-qs_adj=ZRX(:,2:klev+1,5)*zqdm(:,2:klev+1)
-qg_adj=ZRX(:,2:klev+1,6)*zqdm(:,2:klev+1)
-ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV))
-CALL VERTICAL_EXTEND(ZRHOD,klev)
-!
-!------------------------------------------------------------
-! Radiation
-!------------------------------------------------------------
-!
-!AI a nettoyer apres
-! calcul de qsat Fred mai 2024
-zqsat(:,:) = EXP( XALPW - XBETAW / t(:,:) - XGAMW * ALOG( t(:,:) ) )
-zqsat(:,:) = XRD / XRV * zqsat(:,:) / ( pplay(:,:) - zqsat(:,:))
-!
-if (iflag_radia.ge.1) then
-!
-!ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors
-!Cloud fraction is available in the ZCLDFR variable
-!
-! AI attention : dans un 1er temps read_climoz=0
-! Update ozone if day change
-     !print*,'solarlong0, itap, lmt_pas, days_elapsed : ', solarlong0, itap, lmt_pas, days_elapsed
-     !stop
-    IF (MOD(itap-1,lmt_pas) == 0) THEN
-!       IF (read_climoz <= 0) THEN
-          ! Once per day, update ozone from Royer:
-          IF (solarlong0<-999.) then
-             ! Generic case with evolvoing season
-             zzzo=real(days_elapsed+1)
-          ELSE IF (abs(solarlong0-1000.)<1.e-4) then
-             ! Particular case with annual mean insolation
-             zzzo=real(90) ! could be revisited
-             IF (read_climoz/=-1) THEN
-                abort_message ='read_climoz=-1 is recommended when ' &
-                               // 'solarlong0=1000.'
-                CALL abort_physic (modname,abort_message,1)
-             ENDIF
-          ELSE
-             ! Case where the season is imposed with solarlong0
-             zzzo=real(90) ! could be revisited
-          ENDIF
-          zzzo=real(days_elapsed+1)
-          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzzo)
-          !print*,'wo = ', wo
-!       ELSE
-!          !--- ro3i = elapsed days number since current year 1st january, 0h
-!          ro3i=days_elapsed+jh_cur-jh_1jan
-!          !--- scaling for old style files (360 records)
-!          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
-!          IF(adjust_tropopause) THEN
-!             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
-!                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
-!                      time_climoz ,  longitude_deg,   latitude_deg,          &
-!                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
-!          ELSE
-!             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
-!                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
-!                      time_climoz )
-!          ENDIF
-!          ! Convert from mole fraction of ozone to column density of ozone in a
-!          ! cell, in kDU:
-!          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
-!               * zmasse / dobson_u / 1e3
-!          ! (By regridding ozone values for LMDZ only once a day, we
-!          ! have already neglected the variation of pressure in one
-!          ! day. So do not recompute "wo" at each time step even if
-!          ! "zmasse" changes a little.)
-!       ENDIF
-    ENDIF
-!=========================================================================
-! Calculs de l'orbite.
-! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
-! doit donc etre plac\'e avant radlwsw et pbl_surface
-  !CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
-  CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
-  print*,'mth_eq, day_eq, jD_eq, jH_cur =', mth_eq, day_eq, jD_eq, jH_cur
-  day_since_equinox = (jD_cur + jH_cur) - jD_eq
-  !choix entre calcul de la longitude solaire vraie ou valeur fixee = solarlong0
-  IF (solarlong0<-999.) THEN
-     IF (new_orbit) THEN
-        ! calcul selon la routine utilisee pour les planetes
-        CALL solarlong(day_since_equinox, zlongi, dist)
-        print*,'day_since_equinox, zlongi, dist : ', day_since_equinox, zlongi, dist
-     ELSE
-        ! calcul selon la routine utilisee pour l'AR4
-        CALL orbite(REAL(days_elapsed+1),zlongi,dist)
-     ENDIF
-  ELSE
-     zlongi=solarlong0  ! longitude solaire vraie
-     dist=1.            ! distance au soleil / moyenne
-  ENDIF
-    !IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
 
-! AI a refaire en definissant les itap, itaprad    
-    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        ! ============================
-    ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
-    ! l'annee a partir d'une formule analytique.
-    ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
-    ! non nul aux poles.
-    IF (abs(solarlong0-1000.)<1.e-4) THEN
-       CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
-            latitude_deg,longitude_deg,rmu0,fract)
-       swradcorr(:) = 1.0
-       ! JrNt(:) = 1.0
-       zrmu0(:) = rmu0(:)
-    ELSE
-       ! recode par Olivier Boucher en sept 2015
-       SELECT CASE (iflag_cycle_diurne)
-       CASE(0)
-          !  Sans cycle diurne
-          CALL angle(zlongi, latitude_deg, fract, rmu0)
-          swradcorr = 1.0
-          ! JrNt = 1.0
-          zrmu0 = rmu0
-       CASE(1)
-          !  Avec cycle diurne sans application des poids
-          !  bit comparable a l ancienne formulation cycle_diurne=true
-          !  on integre entre gmtime et gmtime+radpas
-          zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
-          CALL zenang(zlongi,jH_cur,0.0,zdtime, &
-               latitude_deg,longitude_deg,rmu0,fract)
-          zrmu0 = rmu0
-          swradcorr = 1.0
-          ! Calcul du flag jour-nuit
-          ! JrNt = 0.0
-          !WHERE (fract.GT.0.0) JrNt = 1.0
-       CASE(2)
-          !  Avec cycle diurne sans application des poids
-          !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
-          !  Comme cette routine est appele a tous les pas de temps de
-          !  la physique meme si le rayonnement n'est pas appele je
-          !  remonte en arriere les radpas-1 pas de temps
-          !  suivant. Petite ruse avec MOD pour prendre en compte le
-          !  premier pas de temps de la physique pendant lequel
-          !  itaprad=0
-          zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
-          zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
-          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
-               latitude_deg,longitude_deg,rmu0,fract)
-          !
-          ! Calcul des poids
-          !
-          zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
-          zdtime2=0.0    !--pas de temps de la physique qui se termine
-          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
-               latitude_deg,longitude_deg,zrmu0,zfract)
-          swradcorr = 0.0
-          WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
-               swradcorr=zfract/fract*zrmu0/rmu0
-          ! Calcul du flag jour-nuit
-          !JrNt = 0.0
-          !WHERE (zfract.GT.0.0) JrNt = 1.0
-       END SELECT
-    ENDIF
-    !print*,'fract = ',fract
-    !sza_o = ACOS (rmu0) *180./pi
-
-    ! Calcul de l'ensoleillement :
- IF (MOD(itaprad,radpas).EQ.0) THEN
-   namelist_ecrad_file='namelist_ecrad'
-   !
- ! Q0 gestion des pas de temps :
-   ! itap, itaprad ?????
-   ! 
- ! Q1 champs lies a la surface ? sorties appel au modele de surface
-   ! zxtsol temp au sol : dvrait etre recuperee de la partie surface
-   ! albsol_dir, albsol_dif aussi
-   albsol_dir = 1. ! AI attention
-   albsol_dif = 1. ! AI attention
-   zxtsol(:) = t(:,1) ! AI attention
-   !
- ! Q2 ozone ? lmdz ok a completer
-   ! wo : pour l'instant appel ozonecm
-   !
- ! Q3 orbite ? lmdz ok a comleter
-   ! dist, rmu0 et fract
-   ! pour l'instant cas solarlong0<-999.
-   ! comment introduire les pas de temps rayonnement
-   !  
- ! Q3 proprietes optiques des nuages ? call cloud_optics_prop.F90 ? non necessaires pour ecrad
-   ! epaisseur optique, emissivite : cldfrarad, cldemirad, cldtaurad , non necessaire pour ecrad
-   ! rayons effectifs : ref_liq, ref_ice ?????
-   ! cloud liq/Ice water mixing contenent (kg/kg) ????
-   ! les rayons effectifs seront calcules en utilisant la routine de l'IFS
-   ! en fonct des contenus en eau, .....
-   cldemirad = 1.
-   cldtaurad = 5.
-   cldtaupirad = 5.
-   ref_liq = 1.
-   ref_ice = 1.
-   ref_liq_pi = 1.
-   ref_ice_pi = 1.
-   !
- ! Q4 aerosols ?
-   ! simplifier l'ecriture : mettre tous les flags et parametres dans aero_mod.F90 ?
-   ! 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, &
-   ! tau_aero_lw_rrtm, &
-   !
- ! Q5 thermo : zqsat
-   ! zqsat specific humidity at liq saturation => 
-   ! calcul de Fred ou a  calculer dans ecrad enf fct de t et p
-   ! flwc, fiwc : qx( 2) et qx(  3) 
-
-   print*,'Entree radlwsw a itap : ', itap
-   CALL radlwsw &
-              (debut, dist, rmu0, fract,  &
-              & paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif,  &
-              & t,qx(:,:,1),wo, &
-              & ZCLDFR, cldemirad, cldtaurad, &
-              & 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, &
-               ! Rajoute par OB pour RRTM
-              & tau_aero_lw_rrtm, &
-              & cldtaupirad, m_allaer, &
-               !
-              ! & zqsat, flwc, fiwc, &
-              & zqsat, qx(:,:,2), qx(:,:,2), &
-              & ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
-              & namelist_ecrad_file, &
-              & heat,heat0,cool,cool0,albpla, &
-              & heat_volc,cool_volc, &
-              & topsw,toplw,solsw,solswfdiff,sollw, &
-              & sollwdown, &
-              & topsw0,toplw0,solsw0,sollw0, &
-              & lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,  &
-              & swdnc0, swdn0, swdn, swupc0, swup0, swup, &
-              & topswad_aero, solswad_aero, &
-              & topswai_aero, solswai_aero, &
-              & topswad0_aero, solswad0_aero, &
-              & topsw_aero, topsw0_aero, &
-              & solsw_aero, solsw0_aero, &
-              & topswcf_aero, solswcf_aero, &
-               !-C. Kleinschmitt for LW diagnostics
-              & toplwad_aero, sollwad_aero,&
-              & toplwai_aero, sollwai_aero, &
-              & toplwad0_aero, sollwad0_aero,&
-               !-end
-              & ZLWFT0_i, ZFLDN0, ZFLUP0, &
-              & ZSWFT0_i, ZFSDN0, ZFSUP0, &
-              & cloud_cover_sw)
- itaprad = 0
- ENDIF ! IF (MOD(itaprad,radpas).EQ.0)    
- itaprad = itaprad + 1
-
- IF (iflag_radia.eq.0) THEN
-    heat=0.
-    cool=0.
-    sollw=0.   ! MPL 01032011
-    solsw=0.
-    radsol=0.
-    swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
-    swup0=0.
-    lwup=0.
-    lwup0=0.
-    lwdn=0.
-    lwdn0=0.
-  ENDIF
- ! Ajouter la tendance des rayonnements (tous les pas)
- ! avec une correction pour le cycle diurne dans le SW
- !
- DO k=1, klev
- ! AI attention 
-    !d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY
-    !d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY
-    !d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
-    !d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
-    d_t_swr(:,k)=swradcorr(:)*heat(:,k)/RDAY
-    d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)/RDAY
-    d_t_lwr(:,k)=-cool(:,k)/RDAY
-    d_t_lw0(:,k)=-cool0(:,k)/RDAY
- ENDDO
-
- !t_seri(:,:) = t_seri(:,:) + d_t_swr(:,:) + d_t_lwr(:,:)
- d_t = d_t + d_t_swr + d_t_lwr
-
- !CALL writefield_phy('d_t_swr',d_t_swr,klev)
- !CALL writefield_phy('d_t_lwr',d_t_lwr,klev)
- !CALL writefield_phy('temp',t,klev)
-
-endif   ! iflag_radia a nettoyer
-!
-! AI passer peut-etre la surface avant radiation
-!    albedo, temp surf, ... pour le rayonnement
-!
-!ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors
-!Cloud fraction is available in the ZCLDFR variable
-
-!
-!------------------------------------------------------------
-! Surface
-!------------------------------------------------------------
-!
-! AI 
 ! compute tendencies to return to the dynamics:
 ! "friction" on the first layer
-if (debut) then
-   tsrf(1:nlon)=t(1:nlon,1)
-   do isoil=1,nsoil
-      tsoil(1:nlon,isoil)=tsrf(1:nlon)
-      !tsoil(1:nlon,isoil)=302.
-   enddo
-   rpi=2.*asin(1.)
-   timesec=5*3600
-   stefan=5.67e-8
-   fsw0=900.
-   call getin_p('iecri',iecri)
-   call getin_p('iflag_surface',iflag_surface)
-endif
+d_u(1:klon,1)=-u(1:klon,1)/86400.
+d_v(1:klon,1)=-v(1:klon,1)/86400.
+! newtonian relaxation towards temp_newton()
+do k=1,klev
+  temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
+  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
+enddo
 
-timesec=timesec+pdtphys
 
-! Computation of qsats. Done even it is not used for ouptuts
-! AI
-!qsats(:) = EXP( XALPW - XBETAW / t(:,1) - XGAMW * ALOG( t(:,1) ) )
-!qsats(:) = XRD / XRV * qsats(:) / ( pplay(:,1) - qsats(:) )
-qsats(:) = zqsat(:,1)
-
-if (iflag_surface == 0 ) then
-      ! Latent and sensible heat fluxes imposed in 1D (exemple Arm cumulus)
-      ! radiative fluxes set to 0.
-      ! capcal and fluxgrd=0. replace generic_soil and allow to compute a surface temperature
-      ! A case close to Rico can be obtained by imposing fsens=-5E-3 et flat=-6E-5
-      capcal(:)=1.e6
-      fluxgrd(:)=0.
-      sensible(:)=-fsens
-      latent(:)=-flat
-      swdnsrf(:)=0.
-      lwdnsrf(:)=0.
-      fluxt(1:klon)=sensible/PHYEX%CST%XCPD
-      fluxq(1:klon)=latent(:)/PHYEX%CST%XLVTT
-      fluxu(1:klon)=-0.01*u(1:klon,1)
-      fluxv(1:klon)=-0.01*v(1:klon,1)
-else
-    ! AI
-    if (iflag_radia.eq.0) then    
-      ! Using the generic soil model
-      if ( nlon == 1 ) then
-          ! Imposed radiative fluxes to simulate a 1D case close to Arm Cumulus
-          swdnsrf(:)=fsw0*max(cos(2*rpi*(timesec-43200.)/(1.15*86400.)),0.)
-          lwdnsrf(:)=400.
-      else
-          ! Model forced in latitude by an idelaized radiative flux that depends on latitude
-          ! only. Used for 3D tests without radiation.
-          swdnsrf(:)=600*cos(latitude(:))+100
-          lwdnsrf(:)=0.
-      endif
-    else
-         swdnsrf(:)=solsw(:)
-         lwdnsrf(:)=sollwdown(:)
-    endif     
-      ! updating the soil temperature (below the surface) from the surface temperature
-      ! en returning the surface conductive flux and a heat capacity (accounting for the
-      ! sentivity of the conductive flux to the surface temperature.
-      call soil_conduct(klon,nsoil,debut,pdtphys,tsrf,tsoil,capcal,fluxgrd,zout)
-
-      ! tsoil_out = tsoil with the vertical dimension of atmosphere outputs
-
-      ! Computing the turbulent fluxes Fqx = rho w'x' = rho (kappa/ln(z/z0x))^2 [ x_surface - x ]. Upward
-      call soil_fluxes(klon,tsrf,qsats,pphi(:,1)/rg,ZRHOD(:,2),u(:,1),v(:,1),t(:,1),qx(:,1,1),fluxu,fluxv,fluxt,fluxq)
-      sensible(:)=fluxt(:)*PHYEX%CST%XCPD
-      latent(:)=fluxq(:)*PHYEX%CST%XLVTT
-endif
-tsoil_out=0. ; tsoil_out(:,1:min(nsoil,nlev))=tsoil(:,1:min(nsoil,nlev))
-
-! Computing the surface temperature
-lwupsrf=stefan*tsrf(:)**4
-tsrf(:)=tsrf(:)+pdtphys*(-latent-sensible+lwdnsrf+swdnsrf-lwupsrf+fluxgrd)/capcal(:)
-
-PSFTH(:) = fluxt(:)
-PSFRV(:) = fluxq(:)
-PSFU(:) = fluxu(:)
-PSFV(:) = fluxv(:)
-PSFSV(:,:) = 0.
-!
-!------------------------------------------------------------
-! Shallow convection
-!------------------------------------------------------------
-!
-CALL SHALLOW_MF(D, PHYEX%CST, PHYEX%NEBN, PHYEX%PARAM_MFSHALLN, PHYEX%TURBN, PHYEX%CSTURB,           &
-     &KRR=KRR, KRRL=KRRL, KRRI=KRRI, KSV=KSV,                                                        &
-     &ONOMIXLG=PHYEX%MISC%ONOMIXLG,KSV_LGBEG=PHYEX%MISC%KSV_LGBEG,KSV_LGEND=PHYEX%MISC%KSV_LGEND,    &
-     &PTSTEP=pdtphys,                                                                                &
-     &PDZZ=zdzm(:,:),PZZ=zz_mass(:,:),                                                               &
-     &PRHODJ=PRHODJ(:,:),PRHODREF=ZRHOD(:,:),                                                        &
-     &PPABSM=ZPABST(:,:),PEXNM=ZEXN(:,:),                                                            &
-     &PSFTH=PSFTH(:),PSFRV=PSFRV(:),                                                                 &
-     &PTHM=ZTHETA(:,:),PRM=ZRX(:,:,:),PUM=ZUT(:,:),PVM=ZVT(:,:),                                     &
-     &PTKEM=PTKEM(:,:),PSVM=ZSVT(:,:,:),                                                             &
-     &PDUDT_MF=PDUDT_MF(:,:),PDVDT_MF=PDVDT_MF(:,:),                                                 &
-     &PDTHLDT_MF=PDTHLDT_MF(:,:),PDRTDT_MF=PDRTDT_MF(:,:),PDSVDT_MF=PDSVDT_MF(:,:,:),                &
-     &PSIGMF=PSIGMF(:,:),PRC_MF=PRC_MF(:,:),PRI_MF=PRI_MF(:,:),PCF_MF=PCF_MF(:,:),                   &
-     &PFLXZTHVMF=ZFLXZTHVMF(:,:),                                                                    &
-     &PFLXZTHMF=ZFLXZTHMF(:,:),PFLXZRMF=ZFLXZRMF(:,:),PFLXZUMF=ZFLXZUMF(:,:),PFLXZVMF=ZFLXZVMF(:,:), &
-     &PTHL_UP=PTHL_UP(:,:),PRT_UP=PRT_UP(:,:),PRV_UP=PRV_UP(:,:),                                    &
-     &PRC_UP=PRC_UP(:,:),PRI_UP=PRI_UP(:,:),                                                         &
-     &PU_UP=PU_UP(:,:), PV_UP=PV_UP(:,:), PTHV_UP=PTHV_UP(:,:), PW_UP=PW_UP(:,:),                    &
-     &PFRAC_UP=PFRAC_UP(:,:),PEMF=PEMF(:,:),PDETR=PDETR(:,:),PENTR=PENTR(:,:),                       &
-     &KKLCL=IKLCL(:),KKETL=IKETL(:),KKCTL=IKCTL(:),PDX=1000.0,PDY=1000.0,KBUDGETS=PHYEX%MISC%NBUDGET )
-
-! Add tendencies of shallow to total physics tendency
-d_u(:,1:klev) = d_u(:,1:klev) + PDUDT_MF(:,2:klev+1)
-d_v(:,1:klev) = d_v(:,1:klev) + PDVDT_MF(:,2:klev+1) 
-ZRXS(:,:,1)=ZRXS(:,:,1)+PDRTDT_MF(:,:)
-ZTHETAS(:,:)=ZTHETAS(:,:)+PDTHLDT_MF(:,:)
-! TODO add SV tendencies
-!
-!------------------------------------------------------------
-! Turbulence
-!------------------------------------------------------------
-! out tendencies
-ZRUS(:,:) = 0.
-ZRVS(:,:) = 0.
-ZRWS(:,:) = 0.
-ZRSVS(:,:,:) = 0.
-ZRTKES(:,:) =  ZTKES(:,:) * PRHODJ(:,:)
-DO JRR=1, KRR
-  ZRRS(:,:,JRR) = ZRXS(:,:,JRR) * PRHODJ(:,:)
-ENDDO
-ZRTHS(:,:) = ZTHETAS(:,:) * PRHODJ(:,:)
-CALL TURB(PHYEX%CST, PHYEX%CSTURB, PHYEX%MISC%TBUCONF, PHYEX%TURBN, PHYEX%NEBN, D, PHYEX%MISC%TLES,               &
-   & KRR, KRRL, KRRI, PHYEX%MISC%HLBCX, PHYEX%MISC%HLBCY, PHYEX%MISC%KGRADIENTS, PHYEX%MISC%KHALO,                &
-   & PHYEX%TURBN%NTURBSPLIT, PHYEX%TURBN%LCLOUDMODIFLM, KSV, PHYEX%MISC%KSV_LGBEG, PHYEX%MISC%KSV_LGEND,          &
-   & PHYEX%MISC%KSV_LIMA_NR, PHYEX%MISC%KSV_LIMA_NS, PHYEX%MISC%KSV_LIMA_NG, PHYEX%MISC%KSV_LIMA_NH,              &
-   & PHYEX%MISC%O2D, PHYEX%MISC%ONOMIXLG, PHYEX%MISC%OFLAT, PHYEX%MISC%OCOUPLES,                                  &
-   & PHYEX%MISC%OBLOWSNOW,PHYEX%MISC%OIBM,                                                                        &
-   & PHYEX%MISC%OFLYER, PHYEX%MISC%COMPUTE_SRC, PHYEX%MISC%PRSNOW,                                                &
-   & PHYEX%MISC%OOCEAN, PHYEX%MISC%ODEEPOC, PHYEX%MISC%ODIAG_IN_RUN,                                              &
-   & PHYEX%TURBN%CTURBLEN_CLOUD, PHYEX%MISC%CMICRO, PHYEX%MISC%CELEC,                                             &
-   & pdtphys,PHYEX%MISC%ZTFILE,                                                                 &
-   & ZDXX(:,:),ZDYY(:,:),zdzm(:,:),                                                             &
-   & ZDZX(:,:),ZDZY(:,:),zz_flux(:,:),                                                          &
-   & ZDIRCOSXW(:),ZDIRCOSYW(:),ZDIRCOSZW(:),ZCOSSLOPE(:),ZSINSLOPE(:),                          &
-   & PRHODJ(:,:),PTHVREF(:,:), PHGRAD(:,:,:), zs(:),                                            &
-   & PSFTH(:),PSFRV(:),PSFSV(:,:),PSFU(:),PSFV(:),                                              &
-   & ZPABST(:,:),ZUT(:,:),ZVT(:,:),PWT(:,:),PTKEM(:,:),ZSVT(:,:,:),ZSRC(:,:),                   &
-   & PLENGTHM(:,:),PLENGTHH(:,:),MFMOIST(:,:),                                                  &
-   & ZBL_DEPTH(:),ZSBL_DEPTH(:),                                                                &
-   & ZCEI(:,:), PHYEX%TURBN%XCEI_MIN, PHYEX%TURBN%XCEI_MAX, PHYEX%TURBN%XCOEF_AMPL_SAT,         &
-   & ZTHETA(:,:),ZRX(:,:,:),                                                                    &
-   & ZRUS(:,:),ZRVS(:,:),ZRWS(:,:),ZRTHS(:,:),ZRRS(:,:,:),ZRSVS(:,:,:),ZRTKES(:,:),             &
-   & PSIGS(:,:),                                                                                &
-   & ZFLXZTHVMF(:,:),ZWTH(:,:),ZWRC(:,:),ZWSV(:,:,:),ZDP(:,:),ZTP(:,:),ZTDIFF(:,:),ZTDISS(:,:), &
-   & PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET                                                    )
-DO JRR=1, KRR
-  ZRXS(:,:,JRR) = ZRRS(:,:,JRR) / PRHODJ(:,:)
-ENDDO
-ZTHETAS(:,:) = ZRTHS(:,:) / PRHODJ(:,:)
-ZTKES(:,:) = ZRTKES(:,:) / PRHODJ(:,:)
-! Add tendencies of turb to total physics tendency
-d_u(:,1:klev) = d_u(:,1:klev) + ZRUS(:,2:klev+1)/PRHODJ(:,2:klev+1)
-d_v(:,1:klev) = d_v(:,1:klev) + ZRVS(:,2:klev+1)/PRHODJ(:,2:klev+1)
-IF(PHYEX%PARAM_MFSHALLN%CMF_CLOUD=='STAT') THEN
-  PSIGS(:,:)=SQRT(PSIGS(:,:)**2 + PSIGMF(:,:)**2)
-ENDIF
-!------------------------------------------------------------
-! Microphysics
-!------------------------------------------------------------
-ZSEA=1.
-ZTOWN=0.
-ZCIT=0.
-ZTHVREFZIKB=0
-CALL RAIN_ICE (D, PHYEX%CST, PHYEX%PARAM_ICEN, PHYEX%RAIN_ICE_PARAMN, PHYEX%RAIN_ICE_DESCRN,                      &
-               PHYEX%ELEC_PARAM, PHYEX%ELEC_DESCR, PHYEX%MISC%TBUCONF,                                            &
-               PHYEX%MISC%OELEC, PHYEX%MISC%OSEDIM_BEARD,                                                         &
-               ZTHVREFZIKB,                                                                                       &
-               pdtphys, KRR, ZEXN,                                                                                &
-               zdzf, PRHODJ, ZRHOD, ZEXN, ZPABST, ZCIT, ZCLDFR,                                                   &
-               ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF,                                                            &
-               ztheta, ZRX(:,:,1), ZRX(:,:,2), ZRX(:,:,3), ZRX(:,:,4), ZRX(:,:,5),                                &
-               ZRX(:,:,6), zthetas, ZRXS(:,:,1), ZRXS(:,:,2), ZRXS(:,:,3), ZRXS(:,:,4), ZRXS(:,:,5), ZRXS(:,:,6), &
-               ZINPRC, ZINPRR, ZEVAP3D,                                                                           &
-               ZINPRS, ZINPRG, ZINDEP, ZRAINFR, PSIGS,                                                            &
-               PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET,                                                           &
-               ZSEA, ZTOWN                                                                                        )
-
-!------------------------------------------------------------
-! Tendencies and time evolution (values for next time step)
-!------------------------------------------------------------
-! Tendencies, mixing ratio -> specific
-d_qx(:,1:klev,1)=d_qx(:,1:klev,1) + (ZRXS(:,2:klev+1,1)-ZRXS0(:,2:klev+1,1))*ZQDM(:,2:klev+1)
-d_qx(:,1:klev,2)=d_qx(:,1:klev,2) + (ZRXS(:,2:klev+1,2)-ZRXS0(:,2:klev+1,2))*ZQDM(:,2:klev+1)
-d_qx(:,1:klev,3)=d_qx(:,1:klev,3) + (ZRXS(:,2:klev+1,4)-ZRXS0(:,2:klev+1,4))*ZQDM(:,2:klev+1)
-d_qr(:,1:klev)=d_qr(:,1:klev) + (ZRXS(:,2:klev+1,3)-ZRXS0(:,2:klev+1,3))*ZQDM(:,2:klev+1)
-d_qs(:,1:klev)=d_qs(:,1:klev) + (ZRXS(:,2:klev+1,5)-ZRXS0(:,2:klev+1,5))*ZQDM(:,2:klev+1)
-d_qg(:,1:klev)=d_qg(:,1:klev) + (ZRXS(:,2:klev+1,6)-ZRXS0(:,2:klev+1,6))*ZQDM(:,2:klev+1)
-! Tendency, theta -> T
-d_t(:,1:klev)=d_t(:,1:klev) + (zthetas(:,2:klev+1)-zthetas0(:,2:klev+1))*ZEXN(:,2:klev+1)
-! TKE
-d_tke(:,1:klev)=d_tke(:,1:klev) + (ZTKES(:,2:klev+1) - ZTKES0(:,2:klev+1))
-
-!Time evolution
-ZQR(:,:)=ZQR(:,:)+d_qr(:,:)*pdtphys
-ZQS(:,:)=ZQS(:,:)+d_qs(:,:)*pdtphys
-ZQG(:,:)=ZQG(:,:)+d_qg(:,:)*pdtphys
-PTKEM(:,2:klev+1)=PTKEM(:,2:klev+1)+d_tke(:,:)*pdtphys
-!
 !------------------------------------------------------------
 ! Entrees sorties
 !------------------------------------------------------------
 
-itap=itap+1
-!zjulian=zjulian+pdtphys/86400.
-!print*,'avant mod(itap,iecri) , zjulian itap',zjulian,itap,iecri
-if ( mod(itap,iecri) == 1 .or. iecri == 1 ) then
-print*,'ECRITURE ITAP=',itap
-    call output_physiqex(debut,zjulian,pdtphys_,presnivs,paprs,u,v,t,qx,ZCLDFR,ZQR,ZQS,ZQG,PTKEM,ZTHETA)
-    call iophys_ecrit('tsrf',       1,'tsrf',       ' ',tsrf)
-    call iophys_ecrit('capcal',     1,'capcal',     ' ',capcal)
-    call iophys_ecrit('swdnsrf',    1,'swdnsrf',    ' ',swdnsrf)
-    call iophys_ecrit('lwdnsrf',    1,'lwdnsrf',    ' ',lwdnsrf)
-    call iophys_ecrit('lwupsrf',    1,'lwupsrf',    ' ',lwupsrf)
-    call iophys_ecrit('sensible',   1,'sensible',   ' ',sensible)
-    call iophys_ecrit('latent',     1,'latent',     ' ',latent)
-    call iophys_ecrit('qsats',      1,'qsats',      ' ',qsats)
-    call iophys_ecrit('tsoil',   nlev,'tsoi',      ' ',tsoil_out)
-!AI
-if (iflag_radia.ge.1) then
-    call iophys_ecrit('lwdn0',       1,'lwdn0',       ' ',lwdn0)
-    call iophys_ecrit('lwdn',       1,'lwdn',       ' ',lwdn)
-    call iophys_ecrit('lwup0',       1,'lwup0',       ' ',lwup0)
-    call iophys_ecrit('lwup',       1,'lwup',       ' ',lwup)
-    call iophys_ecrit('swdn0',       1,'swdn0',       ' ',swdn0)
-    call iophys_ecrit('swdn',       1,'swdn',       ' ',swdn)
-    call iophys_ecrit('swup0',       1,'swup0',       ' ',swup0)
-    call iophys_ecrit('swup',       1,'swup',       ' ',swup)
-  endif
-endif
+
+call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)
+
 
 ! if lastcall, then it is time to write "restartphy.nc" file
@@ -1227,14 +138,4 @@
 
 end subroutine physiqex
-!
-SUBROUTINE VERTICAL_EXTEND(PX,KLEV)
-
- ! fill extra vetical levels to fit MNH interface
-
-REAL, DIMENSION(:,:),   INTENT(INOUT)   :: PX
-INTEGER, INTENT(IN) :: KLEV
-PX(:,1     )= PX(:,2)
-PX(:,KLEV+2)= PX(:,KLEV+1)
-END SUBROUTINE VERTICAL_EXTEND
 
 END MODULE physiqex_mod
