Context Navigation

physiq.F90 @ 2212

Last change on this file since 2212 was 2104, checked in by fhourdin, 10 years ago
Corrections de phymar par Gilles Delaygue.
File size: 67.6 KB

Rev	Line
[2089]	1	! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
	2	!#define IO_DEBUG
	3
	4	!========================================================================================================================
	5	! Adaptation de physiq.F90 de phydev pour mettre en oeuvre une interface dynamique de LMDZ/physique de MAR.
	6	! Martin Ménégoz (Février 2014) martin.menegoz@lgge.obs.ujf-grenoble.fr
	7	! Contact LGGE: gallee@lgge.obs.ujf-grenoble.fr (modèle MAR) ; contact LMD: Frederic.hourdin@lmd.jussieu.fr (modèle LMDZ)
	8	!========================================================================================================================
	9
	10	SUBROUTINE physiq (nlon,nlev, &
	11	& debut,lafin,jD_cur, jH_cur,pdtphys, &
	12	& paprs,pplay,pphi,pphis,ppresnivs,clesphy0, &
	13	& u,v,t,qx, &
	14	& flxmass_w, &
	15	& d_u, d_v, d_t, d_qx, d_ps &
	16	& , dudyn )
	17
	18	USE dimphy, only : klon,klev,klevp1
	19	USE infotrac, only : nqtot
	20	USE comgeomphy, only : rlatd,rlond,airephy
	21	!USE comcstphy, only : rg
	22	USE iophy, only : histbeg_phy,histwrite_phy
	23	USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju,ju2ymds,histsync
	24	USE mod_phys_lmdz_para, only : jj_nb
	25	USE phys_state_var_mod, only : phys_state_var_init
	26	! Modules LMDZ utilisés pour LMDZ-MAR:
	27	USE control_mod
	28
	29	! Modules utilisé pour MAR dans DYN_to_PHY_MAR.f90
	30	use Mod_Real
	31	use Mod_PHY____dat
	32	use Mod_PHY____grd
	33	use Mod_PHY_S0_grd
	34	use Mod_SISVAT_grd
[2104]	35	use Mod_SISVAT_dat
[2089]	36	use Mod_SISVAT_gpt ! surface heat flxes in outputs
	37	use Mod_PHY_DY_kkl
	38	use Mod_PHY_CM_kkl ! pour écrire les snowfall et rainfall dans les outputs
	39	use Mod_PHY_CP_kkl ! pour écrire les snowfall et rainfall dans les outputs
	40
	41	#ifdef CPP_XIOS
	42	USE wxios
	43	#endif
	44
	45	IMPLICIT none
	46
	47	#include "dimensions.h"
	48	! Include LMDZ utilisés ici pour MAR-LMDZ:
	49	#include "comvert.h"
	50	#include "YOMCST.h"
	51
	52	!======================================================================
	53	! Arguments du moniteur general de la physique dans le modele LMDZ
	54	!======================================================================
	55	!
	56	! Arguments:
	57	!
	58	! nlon----input-I-nombre de points horizontaux
	59	! nlev----input-I-nombre de couches verticales, doit etre egale a klev
	60	! debut---input-L-variable logique indiquant le premier passage
	61	! lafin---input-L-variable logique indiquant le dernier passage
	62	! jD_cur -R-jour courant a l'appel de la physique (jour julien)
	63	! jH_cur -R-heure courante a l'appel de la physique (jour julien)
	64	! pdtphys-input-R-pas d'integration pour la physique (seconde)
	65	! paprs---input-R-pression pour chaque inter-couche (en Pa)
	66	! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	67	! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	68	! pphis---input-R-geopotentiel du sol
	69	! ppresnivs-input_R_pressions approximat. des milieux couches ( en PA)
	70	! u-------input-R-vitesse dans la direction X (de O a E) en m/s
	71	! v-------input-R-vitesse Y (de S a N) en m/s
	72	! t-------input-R-temperature (K)
	73	! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	74	! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	75	! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
	76	! flxmass_w -input-R- flux de masse verticale
	77	! d_u-----output-R-tendance physique de "u" (m/s/s)
	78	! d_v-----output-R-tendance physique de "v" (m/s/s)
	79	! d_t-----output-R-tendance physique de "t" (K/s)
	80	! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	81	! d_ps----output-R-tendance physique de la pression au sol
	82	!IM
	83
	84	!======================================================================
	85	! Variables ajoutés pour l'utilisation de la physique de MAR
	86	!======================================================================
	87
	88	integer,parameter :: jjmp1=jjm+1-1/jjm
	89	integer,parameter :: iip1=iim+1
	90
	91	! Indices des traceurs (qx)
	92	! NB: Tous les traceurs de MAR (vapeur d'eau et hydrométéores, solides et liquides) sont transportés via la variables qx (tendances: d_qx)
	93	! Il doivent être déclarés dans un fichier traceur.def avant de lancer la simulation
	94
	95	INTEGER, PARAMETER :: ivap = 1 ! indice de traceurs pour vapeur d'eau
	96	INTEGER, PARAMETER :: iliq = 2 ! indice de traceurs pour eau liquide
	97	INTEGER, PARAMETER :: iice = 3 ! indice de traceurs pour eau solide
	98	INTEGER, PARAMETER :: icin = 4 ! indice de traceurs pour CIN
	99	INTEGER, PARAMETER :: isnow = 5 ! indice de traceurs pour flocons de neige
	100	INTEGER, PARAMETER :: irain = 6 ! indice de traceurs pour gouttes de pluie
	101	INTEGER, PARAMETER :: itke = 7 ! indice de traceurs pour tke
	102	INTEGER, PARAMETER :: itked = 8 ! indice de traceurs pour tke dissipation
	103	! #cw INTEGER, PARAMETER :: iccn = 9 ! indice de traceurs pour CCN
	104
	105	!======================================================================
	106	! Arguments de la routine physiq.F90:
	107	!======================================================================
	108	integer,intent(in) :: nlon ! number of atmospheric colums
	109	integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
	110	real,intent(in) :: jD_cur ! current day number (Julian day)
	111	real,intent(in) :: jH_cur ! current time of day (as fraction of day)
	112	logical,intent(in) :: debut ! signals first call to physics
	113	logical,intent(in) :: lafin ! signals last call to physics
	114	real,intent(in) :: pdtphys ! physics time step (s)
	115	real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
	116	real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
	117	real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
	118	real,intent(in) :: pphis(klon) ! surface geopotential
	119	real,intent(in) :: ppresnivs(klev) ! pseudo-pressure (Pa) of mid-layers
	120	integer,parameter :: longcles=20
	121	real,intent(in) :: clesphy0(longcles) ! Not used
	122	real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
	123	real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
	124	real,intent(in) :: t(klon,klev) ! temperature (K)
	125	real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
	126	real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
	127	real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
	128	real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
	129	real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
	130	real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
	131	real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
	132	real,intent(in) :: dudyn(iim+1,jjmp1,klev) ! Not used
	133
	134	!======================================================================
	135	! Variables locales temporelles
	136	!======================================================================
	137	integer,save :: itau=0 ! counter to count number of calls to physics
	138	!$OMP THREADPRIVATE(itau)
	139	!real :: temp_newton(klon,klev)
	140	integer :: k
	141	logical, save :: first=.true.
	142	!$OMP THREADPRIVATE(first)
	143
	144	!======================================================================
	145	! Inputs/Outputs
	146	!======================================================================
	147	integer :: itau0
	148	real :: zjulian
	149	real :: dtime
	150	integer :: nhori ! horizontal coordinate ID
	151	integer,save :: nid_hist ! output file ID
	152	!$OMP THREADPRIVATE(nid_hist)
	153	integer :: zvertid ! vertical coordinate ID
	154	integer,save :: iwrite_phys ! output every iwrite_phys physics step
	155	!$OMP THREADPRIVATE(iwrite_phys)
	156	integer :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
	157	real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
	158	real :: t_wrt ! frequency of the IOIPSL outputs
	159
	160
	161	!======================================================================
	162	! Arguments d'entrée pour la physique de MAR
	163	!======================================================================
	164
	165	! NB : Pour plus de visibilité pour les développeurs de LMDZ et MAR, on conserve les noms de variables propres à MAR et LMDZ.
	166	! La correspondance entre les variables se fait dans cette interface.
	167
	168	! Flags pour MAR (A ranger ultérieurement dans un fichier .def, selon le protocole LMDZ).
	169	! Attention, ces Flags pilotent l'initialisation de nombreuses variables. Toutes les configurations n'ont pas été testées.
	170	! Ici, on trouve les déclarations utilisées pour une experience aquaplanète.
	171	! Contacter Hubert Gallée pour connaître la signification des Flags.
	172
	173
	174	logical :: FlagSV = .TRUE. ! Surface
	175	logical :: FlagSV_Veg = .TRUE. ! Végétation
	176	logical :: FlagSV_SNo = .TRUE. ! Neige
	177	logical :: FlagSV_BSn = .FALSE. !
	178	logical :: FlagSV_KzT = .TRUE. !
	179	logical :: FlagSV_SWD = .TRUE. ! (T,F) == SW INPUT is (downward,absorbed)
	180	! FlagSV_SBC est utilisé pour lire des options de conditions aux limites. Ici, on le met à FALSE pour l'aquaplanète
	181	logical :: FlagSV_SBC = .FALSE. ! Flag pour utiliser le PHY_SISVAT_INP propre à notre experience
	182	logical :: FlagSV_UBC = .FALSE. !
	183	logical :: FlagAT = .TRUE. !
	184	character(len=1) :: TypeAT = 'e' ! Turbulence Model (e, K, L or H)
	185	logical :: FlagAT_TKE = .TRUE. !
	186	logical :: FlagCM = .TRUE. !
	187	!Gilles : T > F car LMDz utilise les tendances des traceurs (sauf frac.nuage)
	188	logical :: FlagCM_UpD = .FALSE. ! Update de la microphysique immédiat
	189	logical :: FlagCP = .TRUE. ! Avec ou sans convection
	190	logical :: FlagRT = .TRUE. !
	191	logical :: FlagS0_SLO = .FALSE. !
	192	logical :: FlagS0_MtM = .FALSE. !
	193	logical :: Flag_O = .TRUE. !
	194	!Gilles: FlagVR T > F pas de sortie interessante
	195	logical :: FlagVR = .FALSE. !
	196
	197	!======================================================================
	198	! Arguments pour la physique de MAR :
	199	!======================================================================
	200
	201	! Pas de temps par défaut. Ils sont modifiés plus bas en fonction de pdtphys défini dans les fichiers .def.
	202
	203	real :: dt0DYn = 60. ! Time Step, Dynamics [s] recalculé en fonction de day_step
	204	real :: dt0_SV = 60. ! Time Step, SISVAT [s] f(pdtphys)
	205	real :: dt0_AT = 60. ! Time Step, Atm_AT [s] f(pdtphys)
	206	real :: dt0_CM = 60. ! Time Step, CMiPhy [s] f(pdtphys)
	207	real :: dt0_CP = 600. ! Time Step, CVAmnh [s] f(pdtphys)
	208	real :: dt0_RT = 900. ! Time Step, radCEP [s] f(pdtphys)
	209	real :: dx = 4.0e+4 ! Grid Spacing, MAX authorised for CMiPhy [m] Limite max tolérée pour la microphysique MAR (40)
	210	real(kind=real8), dimension(klon) :: dx___HOST ! Grid Spacing max along x axis.
	211	real(kind=real8), dimension(klon) :: dy___HOST ! Grid Size [m]
	212	real :: DD_AxX = 90.00 ! x-Axis Direction [degree] Standard (90dg = Eastward) [dg] utilise juste pour les conditions aux limite et la position du soleil.
	213	real, dimension(klevp1) :: s_HOST ! MAR ne fonctionne que avec des coordonnées sigma: sigma=(p-Psommet)/(psurface -Psommet)
	214	real, dimension(klon) :: sh___HOST ! Surface Height [m] =0 en aquaplanète. A calculer ultérieurement en fonction de l'environnement LMDZ.
	215	real(kind=real8), dimension(klon) :: sh_a_HOST ! Topography Anomaly [m]
	216	real(kind=real8), dimension(klon) :: slopxHOST ! Slope, x-direction [-]
	217	real(kind=real8), dimension(klon) :: slopyHOST ! Slope, y-direction [-]
	218	real(kind=real8), dimension(klon) :: slopeHOST ! Slope [-]
	219	real(kind=real8), allocatable :: MMaskHOST(:,:) ! Mountain Mask [-]
	220	real, dimension(klon) :: lonh_HOST ! Longitude [hour]
	221	real, dimension(klon) :: latr_HOST ! Latitude [radian]
	222	real :: ptop_HOST ! Pressure Model Top [kPa] 0 in LMDZ
	223	real, dimension(klon,klevp1) :: pkta_HOST ! Reduced Potential Temperature [KX/s] ! Temperature * p ** (R/Cp) [KX]
	224	real, dimension(klon) :: psa__HOST ! Pressure Thickness [kPa] a recalculer en fonction du geopotentiel de surface: P*=Psurface-Ptop
	225	real, dimension(klon,klevp1) :: gZa__HOST ! Geopotential Height [m2/s2] pphi dans LMDZ
	226	real, dimension(klon,klevp1) :: gZam_HOST ! Geopotential Height, mid-level [m2/s2] On fait la moyenne entre couche du dessus et couche du dessous.
	227	real, dimension(klon,klev) :: Ua___HOST ! Wind Speed, x-direction [m/s]
	228	real, dimension(klon,klev) :: Va___HOST ! Wind Speed, y-direction [m/s]
	229	real, dimension(klon,klev) :: Wa___HOST ! Wind Speed, z-direction [m/s] Non echangé entre phys et dyn dans LMDZ. Necessaire pour MAR. On le recalcule ci-aprés !
	230	real, dimension(klon,klevp1) :: qv___HOST ! Specific Humidity [kg/kg] qx(:,:,ivap)
	231	real, dimension(klon,klev) :: qw___HOST ! Cloud Droplets Concentration [kg/kg] qx(:,:,iliq)
	232	! #cw real, dimension(klon,klev) :: CCN__HOST ! CCN Concentration [-/kg] qx(:,:,iccn)
	233	real, dimension(klon,klev) :: qi___HOST ! Cloud Crystals Concentration [kg/kg] qx(:,:,iice)
	234	real, dimension(klon,klev) :: CIN__HOST ! CIN Concentration [-/kg] qx(:,:,icin)
	235	real, dimension(klon,klev) :: CF___HOST ! Cloud Fraction [-/kg] Dans le futur, la CF sera transportée dans la dynamique, mais pour l'instant, on n'en a pas besoin.
	236
	237	real, dimension(klon,klev) :: qs___HOST ! Snow Particles Concentration [kg/kg] qx(:,:,isnow)
	238	real, dimension(klon,klev) :: qr___HOST ! Rain Drops Concentration [kg/kg] qx(:,:,irain)
	239	real, dimension(klon,klev) :: TKE__HOST ! Turbulent Kinetic Energy [m2/s2] qx(:,:,itke)
	240	real, dimension(klon,klev) :: eps__HOST ! Turbulent Kinetic Energy Dissipation [m2/s3] qx(:,:,itked)
	241	real, dimension(klon,klev) :: dpkt___dt ! Reduced Potential Temperature TENDENCY, ALL Contribut.[KX/s] a recalculer en fonction de d_t de LMDZ
	242	real, dimension(klon,klev) :: dua____dt ! Wind Speed (x-direc.) TENDENCY, ALL Contribut.[m/s2] output commune MAR-LMD
	243	real, dimension(klon,klev) :: dva____dt ! Wind Speed (y-direc.) TENDENCY, ALL Contribut.[m/s2] output commune MAR-LMD
	244	real, dimension(klon,klev) :: dqv____dt ! Specific Humidity TENDENCY, ALL Contr. [kg/kg/s] output commune MAR-LMD
	245	real, dimension(klon,klev) :: dqw____dt ! Cloud Droplets Concentration TENDENCY, ALL Contr. [kg/kg/s] dqx
	246	! #cw real, dimension(klon,klev) :: dCw____dt ! CCN Concentration TENDENCY, ALL Contr. [1/kg/s] dqx
	247	real, dimension(klon,klev) :: dqi____dt ! Cloud Crystals Concentration TENDENCY, ALL Contr. [kg/kg/s] dqx
	248	real, dimension(klon,klev) :: dCi____dt ! CIN Concentration TENDENCY, ALL Contr. [1/kg/s] dqx
	249	real, dimension(klon,klev) :: dCF____dt ! Cloud Fraction TENDENCY, ALL Contr. [kg/kg/s] dqx
	250	real, dimension(klon,klev) :: dqs____dt ! Snow Particles Concentration TENDENCY, ALL Contr. [kg/kg/s] dqx
	251	real, dimension(klon,klev) :: dqr____dt ! Rain Drops Concentration TENDENCY, ALL Contr. [kg/kg/s] dqx
	252	! #dT& real, dimension(klon,klev) :: dpktSV_dt ! Reduced Potential Temperature Tendency, SISVAT [KX/s]
	253	! #dT real, dimension(klon,klev) :: dpktAT_dt ! Reduced Potential Temperature Tendency, Atm_AT [KX/s]
	254	! #dT real, dimension(klon,klev) :: dqv_AT_dt ! Specific Humidity TENDENCY, Atm_AT [kg/kg/s] dqx
	255	! #dT real, dimension(klon,klev) :: dqw_AT_dt ! Cloud Droplets Concentration TENDENCY, Atm_AT [kg/kg/s] dqx
	256	! #cw real, dimension(klon,klev) :: dCw_AT_dt ! CCN Concentration TENDENCY, Atm_AT [1/s] dqx
	257	! #dT real, dimension(klon,klev) :: dqi_AT_dt ! Cloud Crystals Concentration TENDENCY, Atm_AT [kg/kg/s] dqx
	258	! #dT real, dimension(klon,klev) :: dCi_AT_dt ! CIN Concentration TENDENCY, Atm_AT [1/s] dqx
	259	! #dT real, dimension(klon,klev) :: dqs_AT_dt ! Snow Particles Concentration TENDENCY, Atm_AT [kg/kg/s] dqx
	260	! #dT real, dimension(klon,klev) :: dqr_AT_dt ! Rain Drops Concentration TENDENCY, Atm_AT [kg/kg/s] dqx
	261	! #dT real, dimension(klon,klev) :: dpktCM_dt ! Reduced Potential Temperature Tendency, CMiPhy [KX/s] dqx
	262	! #dT real, dimension(klon,klev) :: dqv_CM_dt ! Specific Humidity TENDENCY, CMiPhy [kg/kg/s] dqx
	263	! #dT real, dimension(klon,klev) :: dqw_CM_dt ! Cloud Droplets Concentration TENDENCY, CMiPhy [kg/kg/s] dqx
	264	! #cw real, dimension(klon,klev) :: dCw_CM_dt ! CCN Concentration TENDENCY, CMiPhy [1/s] dqx
	265	! #dT real, dimension(klon,klev) :: dqi_CM_dt ! Cloud Crystals Concentration TENDENCY, CMiPhy [kg/kg/s] dqx
	266	! #dT real, dimension(klon,klev) :: dCi_CM_dt ! CIN Concentration TENDENCY, CMiPhy [1/s] dqx
	267	! #dT real, dimension(klon,klev) :: dCF_CM_dt ! Cloud Fraction TENDENCY, CMiPhy [1/s] dqx
	268	! #dT real, dimension(klon,klev) :: dqs_CM_dt ! Snow Particles Concentration TENDENCY, CMiPhy [kg/kg/s] dqx
	269	! #dT real, dimension(klon,klev) :: dqr_CM_dt ! Rain Drops Concentration TENDENCY, CMiPhy [kg/kg/s] dqx
	270	! #dT real, dimension(klon,klev) :: dpktCP_dt ! Reduced Potential Temperature TENDENCY, CVAmnh [KX/s]
	271	! #dT real, dimension(klon,klev) :: dqv_CP_dt ! Specific Humidity TENDENCY, CVAmnh [kg/kg/s] dqx
	272	! #dT real, dimension(klon,klev) :: dqw_CP_dt ! Cloud Droplets Concentration TENDENCY, CVAmnh [kg/kg/s] dqx
	273	! #dT real, dimension(klon,klev) :: dqi_CP_dt ! Cloud Crystals Concentration TENDENCY, CVAmnh [kg/kg/s] dqx
	274	! #dT real, dimension(klon,klev) :: dpktRT_dt ! Reduced Potential Temperature TENDENCY, radCEP [KX/s]
	275
	276	integer,save :: it0EXP ! Incremente en fonction de debut et lafin. NB: MAR fait des restarts dans le Fortran, LMDZ avec des scripts shells. On suit ici la logique LMDZ. it0EXP et it0RUN sont donc identiques.
	277	integer,save :: it0RUN ! Incremente en fonction de debut et lafin.
	278	integer :: Year_H ! Time [year]
	279	integer :: Mon__H ! Time [month]
	280	integer :: Day__H ! Time [Day]
	281	integer :: Hour_H ! Time [hour]
	282	integer :: minu_H ! Time [minute]
	283	integer :: sec__H ! Time [s]
	284	integer :: ixq1,i0x0,mxqq ! Domain Dimension: x [-]
	285	integer :: jyq1,j0y0,myqq ! Domain Dimension: y [-]
	286	!integer :: klev ! Domain Dimension: z [-] Commun à MAR et LMDZ
	287	!integer :: klevp1 ! Domain Dimension: z [-] Commun à MAR et LMDZ
	288	integer :: mwq ! Domain Dimension: mosaic [-]
	289	!integer :: klon ! Domain Dimension: x * y [-] Commun à MAR et LMDZ
	290	integer :: kcolw ! Domain Dimension: x * y * mosaic [-]
	291	integer :: m_azim ! Mountain Mask, nb of directions taken into account [-]
	292	integer,parameter :: n0pt=1 ! nombre de points qui peuvent être choisis pour le fichier PHY____.OUT
	293	integer, dimension(n0pt) :: IOi0SV !
	294	integer, dimension(n0pt) :: IOj0SV !
	295	real, dimension(klon) :: sst__HOST ! Ocean FORCING (SST) [K]
	296	! #IP real, dimension(klon) :: sif__HOST ! Ocean FORCING (Sea-Ice Fraction ) [-]
	297	! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq) :: s_T__HOST ! A - O COUPLING n=1: Open Ocean [K]
	298	! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq) :: Alb__HOST ! A - O COUPLING (Surface Albedo ) n=2: Sea Ice [-]
	299	! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq) :: dSdT2HOST ! A - O COUPLING ( d(SH Flux) / dT ) [W/m2/K]
	300	! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq) :: dLdT2HOST ! A - O COUPLING ( d(SH Flux) / dT ) [W/m2/K]
	301	! #TC integer, parameter :: ntrac = 28 !
	302	! #TC real, dimension(ixq1:mxqq,jyq1:myqq,klev,ntrac) :: qxTC ! Aerosols: Atmospheric Contentration
	303	! #TC real, dimension(ixq1:mxqq,jyq1:myqq ,ntrac) :: qsTC ! Aerosols: Near Surface Contentration
	304	! #TC real, dimension(ixq1:mxqq,jyq1:myqq ,ntrac) :: uqTC ! Aerosols: Surf.Flux
	305
	306	!======================================================================
	307	! Variables locales:
	308	!======================================================================
	309
	310	real :: secondes ! pour le calcul de la date
	311	integer :: i ! Pour les boucles
	312	integer :: daysec ! Nombre de secondes dans ue journée
	313	integer,save :: annee_ref ! Année du début de la simulation
	314	!integer :: count_year ! le nombre d'années depuis le début de la simulation
	315	!real,parameter :: epsilon = 1.e-6 ! pour le calcul du temps
	316
	317	! Constantes définies comme dans MAR et utilisées ici:
	318	real :: RdA = 287. ! Perfect Gas Law Constant for Dry Air [J/kg/K]
	319	real :: CpA = 1004. ! Specific Heat Capacity for Dry Air [J/kg/K]
	320	real :: kap ! RdA/CpA [-]
	321
	322	! Pour l'initialisation:
	323	real, dimension(klon) :: Ts___HOST ! Temperature de surface pour l'initialisation.
	324	real, dimension(klon) :: pkt_DY_surf_tmp ! pkt de surface pour l'initialisation.
	325	real, dimension(klon) :: qv_HOST_surf_tmp ! qv de surface pour l'initialisation.
	326
	327	! Pour ecriture en netcdf
	328	real, dimension(klon,klev) :: cldfra ! CF___HOST vertically inverted
	329
	330	!======================================================================
	331	!! Ecriture des sorties netcdf:
	332	!======================================================================
	333
	334	! NB : La version de MAR utilisée ici ne permet pas d'écrire de sorties netcdf.
	335	! On n'utilise pas encore XIOS
	336	! En attendant, on utilise ici le classique ioipsl, ainsi que l'outil iotd de Frédéric Hourdin.
	337
	338	if (debut) then ! Things to do only for the first call to physics
	339	! load initial conditions for physics (including the grid)
	340	call phys_state_var_init() ! some initializations, required before calling phyetat0
	341	call phyetat0("startphy.nc")
	342
	343	! Initialize outputs:
	344	itau0=0
	345	!$OMP MASTER
	346	iwrite_phys_omp=1 !default: output every physics timestep
	347	! NB: getin() is not threadsafe; only one thread should call it.
	348	call getin("iwrite_phys",iwrite_phys_omp)
	349	!$OMP END MASTER
	350	!$OMP BARRIER
	351	iwrite_phys=iwrite_phys_omp
	352	t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
	353	t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
	354	! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
	355	!CALL ymds2ju(annee0, month, dayref, hour, zjulian)
	356	call getin('anneeref',anneeref)
	357	call ymds2ju(anneeref, 1, 1, 0.0, zjulian)
	358	annee_ref=anneeref ! On la sauve pour le reste de la simulation
	359	dtime=pdtphys
	360
	361	! Initialisation de iophy pour l'utilisation de iotd créé par Frédéric Hourdin
	362	call iophys_ini ! NB: cette initialisation concerne l'ecriture de fichiers phys.nc dans PHY_MAR (rechercher iophys_ecrit)
	363
	364	!Gilles : replace zjulian by jD_cur (if correct init of the year)
	365	! call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
	366	call histbeg_phy("histins.nc",itau0,jD_cur,dtime,nhori,nid_hist)
	367
	368	!$OMP MASTER
	369
	370	! define vertical coordinate
	371	call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
	372	ppresnivs,zvertid,'down')
	373	! define variables which will be written in "histins.nc" file
	374	call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
	375	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	376	'inst(X)',t_ops,t_wrt)
	377	call histdef(nid_hist,'pplay','atmospheric pressure','K', &
	378	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	379	'inst(X)',t_ops,t_wrt)
	380	call histdef(nid_hist,'ps','Surface Pressure','Pa', &
	381	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	382	'inst(X)',t_ops,t_wrt)
	383	call histdef(nid_hist,'qx_vap','specific humidity','kg/kg', &
	384	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	385	'inst(X)',t_ops,t_wrt)
	386	call histdef(nid_hist,'qx_liq','atmospheric liquid water content','kg/kg', &
	387	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	388	'inst(X)',t_ops,t_wrt)
	389	call histdef(nid_hist,'qx_ice','atmospheric ice content','kg/kg', &
	390	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	391	'inst(X)',t_ops,t_wrt)
	392	call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
	393	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	394	'inst(X)',t_ops,t_wrt)
	395	call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
	396	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	397	'inst(X)',t_ops,t_wrt)
	398	call histdef(nid_hist,'lonh','longitude','Hours', &
	399	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	400	'inst(X)',t_ops,t_wrt)
	401	call histdef(nid_hist,'latr','latitude','radian', &
	402	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	403	'inst(X)',t_ops,t_wrt)
	404	call histdef(nid_hist,'sst','Prescribed SST','K', &
	405	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	406	'inst(X)',t_ops,t_wrt)
	407	! Tendances physiques de MAR ; pour le developpement MAR-LMDZ:
	408	call histdef(nid_hist,'d_t','Atmospheric temperature trends after MAR physics','K/s', &
	409	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	410	'inst(X)',t_ops,t_wrt)
	411	call histdef(nid_hist,'d_u','Eastward Zonal Wind trend','m/s/s', &
	412	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	413	'inst(X)',t_ops,t_wrt)
	414	call histdef(nid_hist,'d_v','Northward Meridional Wind trend','m/s/s', &
	415	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	416	'inst(X)',t_ops,t_wrt)
	417	! Définition de variables utilisée dans des sous-routines MAR
	418	call histdef(nid_hist,'snow','snowfall convectif + stratiform','m', &
	419	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	420	'inst(X)',t_ops,t_wrt)
	421	call histdef(nid_hist,'rain','rainfall convectif + stratiform','m', &
	422	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	423	'inst(X)',t_ops,t_wrt)
	424	call histdef(nid_hist,'snowCP','snowfall convective','m', &
	425	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	426	'inst(X)',t_ops,t_wrt)
	427	call histdef(nid_hist,'rainCP','rainfall convective','m', &
	428	iim,jj_nb,nhori,1,1,1,zvertid,32, &
	429	'inst(X)',t_ops,t_wrt)
	430	call histdef(nid_hist,'nebulosity','CF___HOST','-', &
	431	iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
	432	'inst(X)',t_ops,t_wrt)
	433	! end definition sequence
	434	call histend(nid_hist)
	435
	436	!XIOS
	437	#ifdef CPP_XIOS
	438	! Déclaration de l'axe vertical du fichier:
	439	!CALL wxios_add_vaxis("presnivs", "histins", klev, ppresnivs)
	440
	441	!Déclaration du pas de temps:
	442	!CALL wxios_set_timestep(dtime)
	443
	444	!Finalisation du contexte:
	445	!CALL wxios_closedef()
	446	#endif
	447	!$OMP END MASTER
	448	endif ! of if (debut)
	449
	450	!======================================================================
	451	! Calcul de variables environnement nécessaires à la physique de MAR
	452	!======================================================================
	453
	454	! NB : Calculs à faire avant l'initialisation de la physique
	455
	456	! Calcul de constantes
	457
	458	kap = RdA / CpA
	459	rpi = 4.*ATAN(1.)
	460	daysec = 86400.
	461
	462	! increment counter itau
	463	itau=itau+1
	464
	465	PRINT*, 'dans physiq.F90 :'
	466	PRINT*, 'itau=',itau
	467
	468	! set all tendencies to zero
	469	PRINT*, 'On initialise les tendances à zéro:'
	470	d_u(1:klon,1:klev)=0.
	471	d_v(1:klon,1:klev)=0.
	472	d_t(1:klon,1:klev)=0.
	473	d_qx(1:klon,1:klev,1:nqtot)=0.
	474	d_ps(1:klon)=0.
	475
	476	! Allocation de tableaux:
	477	! Pour l'instant :
	478	m_azim=1 ! Azimuths pour le calcul des masques de montagne
	479	ALLOCATE(MMaskHOST(klon,m_azim)) ! Mountain Mask [-]
	480
	481	! Appeller iniaqua et phyetat0 pour l'état initial avant le transfert vers les variables MAR !
	482	! A faire pour lorsque l'on prendra en compte les continents
	483
	484	! Timesteps pour MAR (définis ici à partir des pas de temps LMDZ). A modifier éventuellement (via les fichiers .def)
	485	! dt0_CP & dt0_RT modifiees ici selon Hubert avec daystep=720 ie dtdyn=2mn
	486	! pdtphys modifie via gcm.def avec iphysiq=1 & nsplit_phys=1
	487	dt0DYn=daysec/REAL(day_step) ! Time Step, Dynamics [s]
	488	dt0_AT=pdtphys ! Time Step, SISVAT [s]
	489	dt0_SV=pdtphys ! Time Step, Atm_AT [s]
	490	dt0_CM=pdtphys ! Time Step, CMiPhy [s]
	491	dt0_CP=pdtphys*10 ! Time Step, CVAmnh [s]
	492	dt0_RT=pdtphys*20 ! Time Step, radCEP [s]
	493
	494	! Correspondances de grilles: inutilisé car on utilise directement les variables LMD
	495
	496	! Dans un premier temps, une seule mozaique:
	497	mwq=1
	498
	499	kcolw=mwq*klon
	500
	501	DO i=1,klon
	502	IF (rlond(i) .LT. 0) THEN
	503	lonh_HOST(i)=rlond(i)*12./rpi+24. ! from radians to hours
	504	ELSE
	505	lonh_HOST(i)=rlond(i)*12./rpi ! from radians to hours
	506	ENDIF
	507	ENDDO
	508	latr_HOST(:)=rlatd(:) ! from radians to radians
	509
	510	!PRINT*, 'lonh_HOST(:)=',lonh_HOST(:)
	511	!PRINT*, 'latr_HOST(:)=',latr_HOST(:)
	512
	513	! Niveaux verticaux:
	514	! Attention, il faut faire un test pour vérifier si les coordonnées verticaux définis sont adéquats avec ceux de MAR!
	515	DO k = 1, klev
	516	i=klev+1-k
	517	IF (ap(k) .NE. 0) THEN
	518	CALL abort_gcm('','MAR physic can be applied only with sigma levels, check vertical resolution in disvert.F90',1)
	519	ELSE
	520	s_HOST(i)=0.5*(bp(k)+bp(k+1))
	521	ENDIF
	522	END DO
	523	s_HOST(klev+1)=1 ! Dans MAR, sigma est défini depuis tout en haut jusqu'à la surface. s_HOST doit être égal à 1 en surface.
	524
	525	! En attendant que les routines PHY_____INI.F90, PHY_MAR_DAT.F90 et PHY_SISVAT_INP.F90 soient écrites selon un vecteur et non pas selon un tableau 2D lon-lat.
	526	! On utilise ici la grille physique de LMDZ, qui a la taille klon (en argument de la physique): 2 + (JJM-1)IIM (soit 1682 pour une grille dynamique de 4836 déclarée à la compilation).
	527
	528	mxqq=iim
	529	myqq=jjm+1
	530
	531	! Pour MAR, on a besoin des indices du premier point de la grille. Avec une grille globale, c'esty toujours (1,1):
	532	ixq1=1
	533	jyq1=1
	534
	535
	536	! Temps:
	537	IF (debut) THEN
	538	it0EXP=1
	539	ELSE
	540	it0EXP=it0EXP+1
	541	ENDIF
	542
	543	IF (debut) THEN
	544	it0RUN=1
	545	ELSE
	546	it0RUN=it0RUN+1
	547	ENDIF
	548
	549	! Pour l'instant, on ne peut démarrer une simulation seulement le 1er jour de l'année de référence:
	550	! Attention, ju2ymds ne donne pas l'année courante, mais compte les années de simulations à partir de 0 !!! A checker !!!
	551	CALL ymds2ju(annee_ref, 1, 1, 0.0, zjulian)
	552	!Gilles: zjulian already counted in jD_cur if year correctly initialised
	553	!CALL ju2ymds(jD_cur+jH_cur+zjulian,Year_H,Mon__H,Day__H,secondes)
	554	CALL ju2ymds(jD_cur+jH_cur,Year_H,Mon__H,Day__H,secondes)
	555	Hour_H=INT(secondes)/3600
	556	minu_H=(INT(secondes)-Hour_H*3600)/60
	557	!minu_H=INT((secondes-float(Hour_H)*3600.+epsilon)/60.) ! attention aux arrondis machine
	558	sec__H=INT(secondes)-(Hour_H3600)-(minu_H60)
	559
	560	Print*,'CONTROL du temps avant d appeler la physique'
	561	Print*,'zjulian=',zjulian
	562	Print*,'annee_ref=',annee_ref
	563	Print*, 'jD_cur=',JD_cur
	564	Print*, 'jH_cur=',jH_cur
	565	Print*, 'Year_H=',Year_H
	566	Print*, 'Mon__H',Mon__H
	567	Print*, 'Day__H',Day__H
	568	Print*, 'secondes=',secondes
	569	Print*, 'Hour_H=',Hour_H
	570	Print*, 'minu_H=',minu_H
	571	Print*, 'sec__H=',sec__H
	572	Print*, 'klev=',klev
	573	Print*, 'klevp1=',klevp1
	574	Print*, 'nlev=',nlev
	575	Print*, 'nlon=',nlon
	576	Print*, 'klon=',klon
	577	PRINT*, 'RCp=',RCp
	578
	579	! Coordonnées du point (en indice) pour lequel on imprime les sorties en ASCII dans le fichier PHY___________.OUT
	580	i0x0=23
	581	j0y0=48
	582	! Variable utilisée juste pour les outputs:
	583	IOi0SV(:)=1
	584	IOj0SV(:)=1
	585
	586	! Variables inutilisees:
	587	dy___HOST=0.
	588	!CF___HOST=0. ! Pour l'instant pas utilisé, il faut tout de même l'initialiser?
	589	dCF____dt=0.
	590
	591
	592	!======================================================================
	593	! INITIALISATIONS MAR
	594	!======================================================================
	595	!
	596	! NB : les routines MAR ci-aprés doivent encore être travaillées. Ici, elles ne fonctionnent que pour une aquaplanète !
	597
	598	it_RUN = it0RUN
	599	it_EXP = it0EXP
	600
	601	IF (it0RUN.LE.1) THEN
	602
	603	PRINT*,'it0RUN=0; On appelle les routines d initialisation de la physique de MAR'
	604
	605	! Initialization of Mod_PHY____grd
	606	! --------------------------------
	607
	608	! Anciennes initialisations MAR
	609
	610	dxHOST = dx
	611	DirAxX = DD_AxX
	612	ixp1 = ixq1
	613	i_x0 = i0x0
	614	mxpp = mxqq
	615	jyp1 = jyq1
	616	j_y0 = j0y0
	617	mypp = myqq
	618	mzpp = klevp1
	619	kcolp = klon
	620	kcolv = kcolw
	621
	622	! Initialization of Mod_PHY_S0_grd
	623	! --------------------------------
	624	! Pour l'instant, en attendant du relief :
	625	n_azim = m_azim
	626
	627	! Initialization of Mod_SISVAT_grd (Mosaic Discretization only )
	628	! -------------------------------- ( Discretization from Mod_SISVAT_dim : see PHY_SISVAT_ALLOC)
	629
	630	mwp = mwq
	631
[2104]	632	! Initialization of Mod_SISVAT_dat
	633	! --------------------------------
	634	! pour aquaplanete: FixedSST=1
	635	VarSST = 0.
	636
[2089]	637	! Initialization of Physical Constants and Allocation of main Variables
	638	! ---------------------------------------------------------------------
	639
	640	! CONTROL
	641	PRINT*, 'Appel de PHY INI'
	642	! **************
	643	CALL PHY________INI
	644	! **************
	645
	646
	647	! INPUT DATA surface et Calcul des dependances a des Points Horizontaux autres que le Point courant
	648	! =================================================================================================
	649	!
	650
	651	! Pour l'instant, PHY_MAR_DAT n'est pas vectorisé, donc on initilise tout en dur pour l'aquaplanète.
	652	! Ultérieurement, il faudra lires les données topo de LMDZ pour renseigner ces variables à l'initialisation.
	653
	654	PRINT*,'Initialisation en dur des variables topo pour une aquaplanete'
	655	sh___HOST(:)=0.00
	656	sh_a_HOST(:)=0.00
	657	dx___HOST(:) = 1.e3
	658	dy___HOST(:) = 1.e3
	659	slopxHOST(:) = 0.00
	660	slopyHOST(:) = 0.00
	661	slopeHOST(:) = 0.00
	662	MMaskHOST(:,:) = 0.00
	663
	664	! CONTROL
	665	!PRINT*, 'Appel de PHY MAR DAT'
	666	! ***********
	667	! CALL PHY_MAR_DAT( &
	668	!! ***********
	669	!& sh___HOST &! Topographie [m]
	670	!& ,sh_a_HOST &! Topographie, Anomalie [m]
	671	!& ,dx___HOST &! Discretisation, x [m]
	672	!& ,dy___HOST &! Discretisation, y [m]
	673	!& ,slopxHOST &! Pente selon l'Axe x [-]
	674	!& ,slopyHOST &! Pente selon l'Axe y [-]
	675	!& ,slopeHOST &! Pente [-]
	676	!& ,MMaskHOST &! Masque de Montagnes [-]
	677	!& ,ixq1,mxqq &!
	678	!& ,jyq1,myqq &!
	679	!& ,m_azim &!
	680	!& )
	681
	682	! FIRST Initialization of the Surface Variables
	683	! =============================================
	684
	685
	686	! Initial Surface Temperature ( initialement dans DYN_to_PHY_MAR)
	687	! ---------------------------
	688
	689	! CONTROL
	690	PRINT*, 'Calcul de Ts___HOST dans la physique'
	691	PRINT*, 'Initialisation de la temperature de surface avec les SST aquaplanète'
	692	DO i = 1,klon
	693	Ts___HOST(i)=273.+27.(1-sin(1.5rlatd(i))**2)
	694	IF ((rlatd(i).GT.1.0471975).OR.(rlatd(i).LT.-1.0471975)) THEN
	695	Ts___HOST(i)=273.
	696	ENDIF
	697
	698	!PRINT*, 'Initialisation de la temperature de surface avec la temperature du premier niveau'
	699	!Ts___HOST(i) = pkta_HOST(i,klevp1)(psa__HOST(i)+ptop_HOST)*RCp
	700	!Ts___HOST(i) = t(i,klev)
	701	ENDDO
	702
	703	! CONTROL
	704	PRINT*, 'Avant PHY SISVAT INP'
	705	PRINT*,'MAXVAL(Ts___HOST(:)=)',maxval(Ts___HOST(:))
	706	PRINT*,'MINVAL(Ts___HOST(:)=)',minval(Ts___HOST(:))
	707	PRINT*, 'Appel de PHY SISVAT INP'
	708	! **************
	709	CALL PHY_SISVAT_INP(FlagSV_SBC,Ts___HOST,klon)
	710	! **************
	711
	712	PRINT*, 'Apres PHY SISVAT INP'
	713	PRINT*,'MAXVAL(Ts___HOST(:)=)',maxval(Ts___HOST(:))
	714	PRINT*,'MINVAL(Ts___HOST(:)=)',minval(Ts___HOST(:))
	715
	716	END IF ! (it0RUN.LE.1)
	717
	718
	719	! Horizontal Interactions
	720	! ======================
	721
	722
	723	! Interface: DATA for Stand-Alone Version of the Physical Package
	724	! ===============================================================
	725
	726	IF (.NOT. FlagRT .OR. &
	727	& .NOT. FlagAT .OR. &
	728	& .NOT. FlagCM) THEN
	729	! & .NOT. FlagCP) THEN
	730	! CONTROL
	731	PRINT*, 'Appel de PHY SISVAT UBC'
	732	! **************
	733	CALL PHY_SISVAT_UBC(FlagRT,FlagAT,FlagCM,FlagCP)
	734	! **************
	735
	736	END IF
	737
	738	PRINT*, 'Apres PHY_SISVAT_UBC'
	739	Print*, 'mzp=',mzp
	740	Print*, 'mzpp=',mzpp
	741
	742	!======================================================================
	743	! Calcul des inputs de MAR à partir de LMDZ
	744	!======================================================================
	745
	746	! NB : Calculs à faire aprés l'initialisation de la physique du MAR
	747
	748	PRINT*, 'Calcul des inputs de MAR à partir de LMDZ'
	749
	750	! Correspondance ou transformation de variables:
	751	! RKAPPA=R/Cp dans un module LMD ;
	752	! Probleme, on n'a pas la temperature de surface dans les arguments !!! a voir avec Hubert comment on fait...
	753	! En attendant, on prend la température du premier niveau au lieu de celle de la surface (idem pour le haut de l'atmosphère):
	754	! NB : Attention aux unités !
	755
	756	IF (debut) THEN
[2104]	757	! On initialise avec la temperature du premier niveau au lieu de la temperature de surface à la première itération:
	758	pkt_DY_surf_tmp(:)=t(:,1) / (pplay(:,1)/1000) ** kap
[2089]	759	ELSE
[2104]	760	pkt_DY_surf_tmp(:)=pkt_DY(:,klev+1)
[2089]	761	ENDIF
[2104]	762
[2089]	763	pkta_HOST(:,klev+1)=pkt_DY_surf_tmp(:)
	764	DO k = 1, klev
	765	i=klev+1-k
	766	pkta_HOST(:,i) = t(:,k) / (pplay(:,k)/1000) ** kap
	767	END DO
	768
	769	ptop_HOST=0.
	770	psa__HOST(:)=(paprs(:,1)-paprs(:,klev+1))/1000. ! [kPa]
	771
	772	! Correspondance du géopotentiel entre LMDZ et MAR :
	773
	774	! Aux niveaux:
	775	DO k = 1, klev
	776	i=klev+1-k
	777	gZa__HOST(:,i)=pphi(:,k)
	778	END DO
	779	gZa__HOST(:,klev+1)=pphis(:)
	780
	781	! Aux inter-niveaux:
	782	DO k = 1, klev
	783	gZam_HOST(:,k+1)=0.5*(gZa__HOST(:,k+1)+gZa__HOST(:,k))
	784	END DO
	785	gZam_HOST(:,1)=1.5gZa__HOST(:,1)-0.5gZa__HOST(:,2) ! Extrapolation au sommet de l'atmosphère
	786	DO k = 1, klev
	787	i=klev+1-k
	788	Ua___HOST(:,i)=u(:,k)
	789	Va___HOST(:,i)=v(:,k)
	790	ENDDO
	791
	792	! Calcul de vitesse verticale a partir de flux de masse verticale :
	793	DO k = 1, klev
	794	i=klev+1-k
	795	!omega(i,k) = RG*flxmass_w(i,k) / airephy(i) ! omega en Pa/s
	796	Wa___HOST(:,i)=flxmass_w(:,k) / airephy(:) * (gZam_HOST(:,i+1) - gZam_HOST(:,i))/(paprs(:,k+1)-paprs(:,k)) ! Equilibre hydrostatique
	797	END DO
	798	Wa___HOST(:,klev)=0 ! Vitesse nulle à la surface.
	799	Wa___HOST(:,1)=0 ! Vitesse nulle au sommet.
	800
	801	! Tendances:
	802	DO k = 1, klev
	803	i=klev+1-k
	804	dpkt___dt(:,i)=d_t(:,k) / (pplay(:,k)/1000) ** kap
	805	dua____dt(:,i)=d_u(:,k)
	806	dva____dt(:,i)=d_v(:,k)
	807	ENDDO
	808
	809
	810	DO k = 1, klev
	811	i=klev+1-k
	812
	813	! Traceurs:
	814	qv___HOST(:,i)=qx(:,k,ivap) ! Specific Humidity [kg/kg]
	815	qw___HOST(:,i)=qx(:,k,iliq) ! Cloud Droplets Concentration [kg/kg]
	816	qi___HOST(:,i)=qx(:,k,iice) ! Cloud Crystals Concentration [kg/kg]
	817	CIN__HOST(:,i)=qx(:,k,icin) ! CIN Concentration [-/kg]
	818	qs___HOST(:,i)=qx(:,k,isnow) ! Snow Particles Concentration [kg/kg]
	819	qr___HOST(:,i)=qx(:,k,irain) ! Rain Drops Concentration [kg/kg]
	820	TKE__HOST(:,i)=qx(:,k,itke) ! Turbulent Kinetic Energy [m2/s2]
	821	eps__HOST(:,i)=qx(:,k,itked) ! Turbulent Kinetic Energy Dissipation [m2/s3]
	822	! #cw CCN__HOST(:,:)=qx(:,:,iccn) ! CCN Concentration [-/kg]
	823
	824	END DO
	825
	826	! Dans MAR, qv(klev+1)=qv(klev)
	827	IF (debut) THEN
	828	! On initialise avec le qv du premier niveau au lieu du qv de surface à la première itération:
	829	qv_HOST_surf_tmp(:)=qv___HOST(:,klev)
	830	ELSE
	831	qv_HOST_surf_tmp(:)=qv__DY(:,klev+1)
	832	ENDIF
	833	qv___HOST(:,klev+1)=qv_HOST_surf_tmp(:) ! Specific Humidity [kg/kg]
	834
	835	DO k = 1, klev
	836	i=klev+1-k
	837
	838	! Dérivée des traceurs:
	839	dqv____dt(:,i)=d_qx(:,k,ivap) ! Specific Humidity TENDENCY, ALL Contr. [kg/kg/s]
	840	dqw____dt(:,i)=d_qx(:,k,iliq) ! Cloud Droplets Concentration TENDENCY, ALL Contr. [kg/kg/s]
	841	dqi____dt(:,i)=d_qx(:,k,iice) ! Cloud Crystals Concentration TENDENCY, ALL Contr. [kg/kg/s]
	842	dCi____dt(:,i)=d_qx(:,k,icin) ! CIN Concentration TENDENCY, ALL Contr. [1/kg/s]
	843	dqs____dt(:,i)=d_qx(:,k,isnow) ! Snow Particles Concentration TENDENCY, ALL Contr. [kg/kg/s]
	844	dqr____dt(:,i)=d_qx(:,k,irain) ! Rain Drops Concentration TENDENCY, ALL Contr. [kg/kg/s]
	845	! #cw dCw____dt ! CCN Concentration TENDENCY, ALL Contr. [1/kg/s]
	846
	847	END DO
	848
	849	!======================================================================
	850	! Test Aquaplanète:
	851	!======================================================================
	852
	853	! CONTROL
	854	PRINT*, 'Calcul des SST en aquaplanète pour MAR comme dans le cas n°101 de LMDZ' ! (cf iniaqua de LMDZ)
	855	!IF (debut) THEN
	856	DO i=1,klon
	857	sst__HOST(i)=273.+27.(1-sin(1.5rlatd(i))**2)
	858	IF ((rlatd(i).GT.1.0471975).OR.(rlatd(i).LT.-1.0471975)) THEN
	859	sst__HOST(i)=273.
	860	ENDIF
	861	ENDDO
	862	!ENDIF
	863
	864	!======================================================================
	865	! CONTROL
	866	!======================================================================
	867	PRINT*, 'Control variables LMDZ entree de la physique'
	868	PRINT*,'MAXVAL(pplay(:,klev))=',MAXVAL(pplay(:,klev))
	869	PRINT*,'MINVAL(pplay(:,klev))=',MINVAL(pplay(:,klev))
	870	PRINT*,'MAXVAL(t(:,klev))=',MAXVAL(t(:,klev))
	871	PRINT*,'MINVAL(t(:,klev))=',MINVAL(t(:,klev))
	872	PRINT*,'MAXVAL(pplay(:,1))=',MAXVAL(pplay(:,1))
	873	PRINT*,'MINVAL(pplay(:,1))=',MINVAL(pplay(:,1))
	874	PRINT*,'MAXVAL(t(:,1))=',MAXVAL(t(:,1))
	875	PRINT*,'MINVAL(t(:,1))=',MINVAL(t(:,1))
	876	PRINT*,'MAXVAL(pplay)=',MAXVAL(pplay)
	877	PRINT*,'MINVAL(pplay)=',MINVAL(pplay)
	878	PRINT*,'MAXVAL(t)=',MAXVAL(t)
	879	PRINT*,'MINVAL(t)=',MINVAL(t)
	880	PRINT*,'MAXVAL(pphis)=',MAXVAL(pphis)
	881	PRINT*,'MINVAL(pphis)=',MINVAL(pphis)
	882	PRINT*,'kap=',kap
	883	PRINT*, 'CONTROL variables MAR entree de la physique'
	884	PRINT*,'MAXVAL(pkta_HOST)=',MAXVAL(pkta_HOST)
	885	PRINT*,'MAXVAL(psa__HOST)=',MAXVAL(psa__HOST)
	886	PRINT*,'MAXVAL(gZa__HOST)=',MAXVAL(gZa__HOST)
	887	PRINT*,'MAXVAL(gZam_HOST)=',MAXVAL(gZam_HOST)
	888	PRINT*,'MAXVAL(Ua___HOST)=',MAXVAL(Ua___HOST)
	889	PRINT*,'MAXVAL(Va___HOST)=',MAXVAL(Va___HOST)
	890	PRINT*,'MAXVAL(qv___HOST)=',MAXVAL(qv___HOST)
	891	PRINT*,'MINVAL(pkta_HOST)=',MINVAL(pkta_HOST)
	892	PRINT*,'MINVAL(psa__HOST)=',MINVAL(psa__HOST)
	893	PRINT*,'MINVAL(gZa__HOST)=',MINVAL(gZa__HOST)
	894	PRINT*,'MINVAL(gZam_HOST)=',MINVAL(gZam_HOST)
	895	PRINT*,'MINVAL(Ua___HOST)=',MINVAL(Ua___HOST)
	896	PRINT*,'MINVAL(Va___HOST)=',MINVAL(Va___HOST)
	897	PRINT*,'MINVAL(qv___HOST)=',MINVAL(qv___HOST)
	898	PRINT*,'Control iterations:'
	899	PRINT*,'it0RUN=',it0RUN
	900	PRINT*,'it0EXP=',it0EXP
	901	PRINT*,'it_RUN=',it_RUN
	902	PRINT*,'it_EXP=',it_EXP
	903
	904	! On choisi l'endroit où on fait les diagnostics de verification MAR (fichier ASCII):
	905	! Pour la phase de débug, on cherche des Nan !
	906
	907	DO i=1,klon
	908	DO k=1,klev
	909	IF (isnan(Va___HOST(i,k))) THEN
	910	ikl0 = i
	911	PRINT*,'Attention : NaN at'
	912	PRINT,'longitude=',rlond(i)180/rpi
	913	PRINT,'latitude=',rlatd(i)180/rpi
	914	ENDIF
	915	ENDDO
	916	ENDDO
	917
	918	!======================================================================
	919	! Appel de PHY_MAR
	920	!======================================================================
	921
	922	! NB: Logiquement, il faut conserver cet appel tel quel, car Hubert Gallée s'engage à ne pas le changer lors des mises à jour !
	923
	924	! CONTROL
	925	PRINT*, 'Appel de PHY MAR'
	926	!PRINT*, 'FlagSV =',FlagSV
	927	!PRINT*, 'FlagSV_Veg =',FlagSV_Veg
	928	!PRINT*, 'FlagSV_SNo =',FlagSV_SNo
	929	!PRINT*, 'FlagSV_BSn =',FlagSV_BSn
	930	!PRINT*, 'FlagSV_KzT =',FlagSV_KzT
	931	!PRINT*, 'FlagSV_SWD =',FlagSV_SWD
	932	!PRINT*, 'FlagSV_SBC =',FlagSV_SBC
	933	!PRINT*, 'FlagSV_UBC =',FlagSV_UBC
	934	!PRINT*, 'FlagAT =',FlagAT
	935
	936	CALL PHY_MAR(FlagSV & ! FLAG for SISVAT: (T,F) = (active OR NOT)
	937	& ,FlagSV_Veg & ! FLAG for SISVAT: (T,F) = (Variable Vegetation OR NOT)
	938	& ,FlagSV_SNo & ! FLAG for SISVAT: (T,F) = (Snow Model active OR NOT)
	939	& ,FlagSV_BSn & ! FLAG for SISVAT: (T,F) = (Blowing Snow Model active OR NOT)
	940	& ,FlagSV_KzT & ! FLAG for SISVAT: (T,F) = (pkt Turb.Transfert active OR NOT in SISVAT)
	941	& ,FlagSV_SWD & ! FLAG for SISVAT: (T,F) = (Modify SW INPUT->downward OR NOT)
	942	& ,FlagSV_SBC & ! FLAG for SISVAT: (T,F) = (INPUT of Soil & Vege DATA OR NOT in SISVAT)
	943	& ,FlagSV_UBC & ! FLAG for SISVAT: (T,F) = (pkt UpperBC is Von Neuman OR NOT in SISVAT)
	944	& ,FlagAT & ! FLAG for Atm_AT: (T,F) = (Turbulent Transfer active OR NOT)
	945	& ,TypeAT & ! TYPE of Atm_AT: (e= Ee Duynkerke, K= Ee Kitada, L= EL, H= Ee Huan-R)
	946	& ,FlagAT_TKE & ! FLAG for genTKE: (T,F) = (TKE-e Model active OR NOT)
	947	& ,FlagCM & ! FLAG for CMiPhy: (T,F) = (Cloud Microphysics active OR NOT)
	948	& ,FlagCM_UpD & ! FLAG for CMiPhy: (T,F) = (qv & hydrometeors updated OR NOT IN CMiPhy)
	949	& ,FlagCP & ! FLAG for Convection Paramet.
	950	& ,FlagRT & ! FLAG for Radiative Transfer
	951	& ,FlagS0_SLO & ! FLAG for Insolation, Surfac.Slope Impact included NEW
	952	& ,FlagS0_MtM & ! FLAG for Insolation, Surfac.Slope & Mountain Mask Impacts included NEW
	953	& ,Flag_O & ! FLAG for OUTPUT
	954	& ,FlagVR & ! FLAG for OUTPUT for VERIFICATION
	955	& ,dt0DYn & ! Time STEP between 2 CALLs of PHY_MAR [s] I, fix
	956	& ,dt0_SV & ! Time STEP between 2 CALLs of SISVAT [s] I, fix
	957	& ,dt0_AT & ! Time STEP between 2 CALLs of Atm_AT [s] I, fix
	958	& ,dt0_CM & ! Time STEP between 2 CALLs of CMiPhy [s] I, fix
	959	& ,dt0_CP & ! Time STEP between 2 CALLs of CVamnh [s] I, fix
	960	& ,dt0_RT & ! Time STEP between 2 CALLs of radCEP [s] I, fix
	961	& ,dx & ! Grid Mesh size (Horizontal) [m] I, fix
	962	& ,DD_AxX & ! Grid x-Axis Direction [degree] I, fix
	963	& ,s_HOST & ! Grid (Vertical) of HOST (NORMALIZED PRESSURE assumed) [-] I, fix
	964	& ,sh___HOST & ! Topography [m] I, fix
	965	& ,sh_a_HOST & ! Topography Anomaly [m] I, fix NEW
	966	& ,slopxHOST & ! Slope, x-direction [-] I, fix NEW
	967	& ,slopyHOST & ! Slope, y-direction [-] I, fix NEW
	968	& ,slopeHOST & ! Slope [-] I, fix NEW
	969	& ,MMaskHOST & ! Mountain Mask [-] I, fix NEW
	970	& ,lonh_HOST & ! Longitude [hour] I, fix
	971	& ,latr_HOST & ! Latitude [radian] I, fix
	972	& ,pkta_HOST & ! Reduced Potential Temperature [XK] I, O
	973	& ,ptop_HOST & ! Pressure, Model Top [kPa] I, fix
	974	& ,psa__HOST & ! Pressure Thickness [kPa] I
	975	& ,gZa__HOST & ! Geopotential Height [m2/s2] I
	976	& ,gZam_HOST & ! Geopotential Height, mid-level [m2/s2] I
	977	& ,Ua___HOST & ! Wind , x-Direction [m/s] I
	978	& ,Va___HOST & ! Wind , y-Direction [m/s] I
	979	& ,Wa___HOST & ! Wind , z-Direction [m/s] I
	980	& ,qv___HOST & ! Specific Humidity [kg/kg] I, O
	981	& ,qw___HOST & ! Cloud Droplets Concentration [kg/kg] I, O
	982	! #cw& ,CCN__HOST & ! CCN Concentration [-/kg]
	983	& ,qi___HOST & ! Cloud Crystals Concentration [kg/kg] I, O
	984	& ,CIN__HOST & ! CIN Concentration [-/kg] I, O
	985	& ,CF___HOST & ! Cloud Fraction [-] I, O
	986	& ,qs___HOST & ! Snow Particles Concentration [kg/kg] I, O
	987	& ,qr___HOST & ! Rain Drops Concentration [kg/kg] I, O
	988	& ,TKE__HOST & ! Turbulent Kinetic Energy [m2/s2] I, O
	989	& ,eps__HOST & ! Turbulent Kinetic Energy Dissipation [m2/s3] I, O
	990	& ,dpkt___dt & ! Reduced Potential Temperature TENDENCY, ALL Contribut.[KX/s] O
	991	& ,dua____dt & ! Wind Speed (x-direc.) TENDENCY, ALL Contribut.[m/s2] O
	992	& ,dva____dt & ! Wind Speed (y-direc.) TENDENCY, ALL Contribut.[m/s2] O
	993	& ,dqv____dt & ! Specific Humidity TENDENCY, ALL Contr. [kg/kg/s] O
	994	& ,dqw____dt & ! Cloud Droplets Concentration TENDENCY, ALL Contr. [kg/kg/s] O
	995	! #cw& ,dCw____dt & ! CCN Concentration TENDENCY, ALL Contr. [1/s]
	996	& ,dqi____dt & ! Cloud Crystals Concentration TENDENCY, ALL Contr. [kg/kg/s] O
	997	& ,dCi____dt & ! CIN Concentration TENDENCY, ALL Contr. [1/s] O
	998	& ,dCF____dt & ! Cloud Fraction TENDENCY, ALL Contr. [1/s] O
	999	& ,dqs____dt & ! Snow Particles Concentration TENDENCY, ALL Contr. [kg/kg/s] O
	1000	& ,dqr____dt & ! Rain Drops Concentration TENDENCY, ALL Contr. [kg/kg/s] O
	1001	! #dT& ,dpktSV_dt & ! Reduced Potential Temperature TENDENCY, SISVAT [KX/s] (O)
	1002	! #dT& ,dpktAT_dt & ! Reduced Potential Temperature TENDENCY, Atm_AT [KX/s] (O)
	1003	! #dT& ,dqv_AT_dt & ! Specific Humidity TENDENCY, Atm_AT [kg/kg/s] (O)
	1004	! #dT& ,dqw_AT_dt & ! Cloud Droplets Concentration TENDENCY, Atm_AT [kg/kg/s] (O)
	1005	! #dT& ,dqi_AT_dt & ! Cloud Crystals Concentration TENDENCY, Atm_AT [kg/kg/s] (O)
	1006	! #dT& ,dqs_AT_dt & ! Snow Particles Concentration TENDENCY, Atm_AT [kg/kg/s] (O)
	1007	! #dT& ,dqr_AT_dt & ! Rain Drops Concentration TENDENCY, Atm_AT [kg/kg/s] (O)
	1008	! #cw& ,dCw_AT_dt & ! CCN Concentration TENDENCY, Atm_AT [1/s] (O)
	1009	! #dT& ,dCi_AT_dt & ! CIN Concentration TENDENCY, Atm_AT [1/s] (O)
	1010	! #dT& ,dpktCM_dt & ! Reduced Potential Temperature TENDENCY, CMiPhy [KX/s] (O)
	1011	! #dT& ,dqv_CM_dt & ! Specific Humidity TENDENCY, CMiPhy [kg/kg/s] (O)
	1012	! #dT& ,dqw_CM_dt & ! Cloud Droplets Concentration TENDENCY, CMiPhy [kg/kg/s] (O)
	1013	! #dT& ,dCF_CM_dt & ! Cloud Fraction TENDENCY, CMiPhy [1/s] (O)
	1014	! #dT& ,dqi_CM_dt & ! Cloud Crystals Concentration TENDENCY, CMiPhy [kg/kg/s] (O)
	1015	! #dT& ,dqs_CM_dt & ! Snow Particles Concentration TENDENCY, CMiPhy [kg/kg/s] (O)
	1016	! #dT& ,dqr_CM_dt & ! Rain Drops Concentration TENDENCY, CMiPhy [kg/kg/s] (O)
	1017	! #cw& ,dCw_CM_dt & ! CCN Concentration TENDENCY, CMiPhy [1/s] (O)
	1018	! #dT& ,dCi_CM_dt & ! CIN Concentration TENDENCY, CMiPhy [1/s] (O)
	1019	! #dT& ,dpktCP_dt & ! Reduced Potential Temperature TENDENCY, CVAmnh [KX/s] (O)
	1020	! #dT& ,dqv_CP_dt & ! Specific Humidity TENDENCY, CVAmnh [kg/kg/s] (O)
	1021	! #dT& ,dqw_CP_dt & ! Cloud Droplets Concentration TENDENCY, CVAmnh [kg/kg/s] (O)
	1022	! #dT& ,dqi_CP_dt & ! Cloud Crystals Concentration TENDENCY, CVAmnh [kg/kg/s] (O)
	1023	! #dT& ,dpktRT_dt & ! Reduced Potential Temperature TENDENCY, radCEP [KX/s] (O)
	1024	& ,sst__HOST & ! Ocean FORCING (SST) [K] I
	1025	! #IP& ,sif__HOST & ! Ocean FORCING (Sea-Ice Fraction ) [-] I
	1026	! #AO& ,s_T__HOST & ! Ocean COUPLING (Surface Temperat.) n=1: Open Ocean [-] I,NEMO
	1027	! #AO& ,Alb__HOST & ! Ocean COUPLING (Surface Albedo ) n=2: Sea Ice [-] I,NEMO
	1028	! #AO& ,dSdT2HOST & ! Ocean COUPLING ( d(SH Flux) / dT ) [W/m2/K] O
	1029	! #AO& ,dLdT2HOST & ! Ocean COUPLING ( d(LH Flux) / dT ) [W/m2/K] O
	1030	!dead& ,it0EXP,it0RUN & ! Iteration
	1031	& ,Year_H,Mon__H,Day__H,Hour_H,minu_H,sec__H & ! Time
	1032	& ,ixq1 ,i0x0 ,mxqq & ! Domain Dimension: x
	1033	& ,jyq1 ,j0y0 ,myqq & ! Domain Dimension: y
	1034	& ,klev ,klevp1 & ! Domain Dimension: z
	1035	& ,mwq & ! Domain Dimension: mosaic
	1036	& ,klon & ! Domain Dimension: x * y NEW
	1037	& ,kcolw & ! Domain Dimension: x * y * mosaic NEW
	1038	& ,m_azim & ! Mountain Mask, nb of directions taken into account [-] NEW
	1039	& ,IOi0SV,IOj0SV,n0pt) ! Indices of OUTPUT Grid Point
	1040	! *******
	1041
	1042	!======================================================================
	1043	! On renvoie les variables nécessaires à la dynamique de LMDZ:
	1044	!======================================================================
	1045
	1046	PRINT*, 'On est sorti de PHY_MAR, on peut actualiser les variables pour LMDZ'
	1047
	1048	d_ps(:)=0.0 ! d_ps----output-R-tendance physique de la pression au sol ! égal à zéro sur la terre (car on néglige la perte de poid de la colonne lorsqu'il pleut).
	1049
	1050	DO k = 1, klev
	1051	i=klev+1-k
	1052
	1053	d_u(:,k)=dua____dt(:,i) ! d_u-----output-R-tendance physique de "u" (m/s/s)
	1054	d_v(:,k)=dva____dt(:,i) ! d_v-----output-R-tendance physique de "v" (m/s/s)
	1055	d_t(:,k)=dpkt___dt(:,i)(pplay(:,k)/1000) * kap !d_t(:,:)= ! d_t-----output-R-tendance physique de "t" (K/s). Note: la pression ne change pas au cours de la physique.
	1056	! d_qx----output-R-tendance physique de "qx" (kg/kg/s):
	1057	d_qx(:,k,ivap)=dqv____dt(:,i) ! Specific Humidity TENDENCY, ALL Contr. [kg/kg/s]
	1058	d_qx(:,k,iliq)=dqw____dt(:,i) ! Cloud Droplets Concentration TENDENCY, ALL Contr. [kg/kg/s]
	1059	d_qx(:,k,iice)=dqi____dt(:,i) ! Cloud Crystals Concentration TENDENCY, ALL Contr. [kg/kg/s]
	1060	d_qx(:,k,icin)=dCi____dt(:,i) ! CIN Concentration TENDENCY, ALL Contr. [1/kg/s]
	1061	d_qx(:,k,isnow)=dqs____dt(:,i) ! Snow Particles Concentration TENDENCY, ALL Contr. [kg/kg/s]
	1062	d_qx(:,k,irain)=dqr____dt(:,i) ! Rain Drops Concentration TENDENCY, ALL Contr. [kg/kg/s]
	1063	! #cw dCw____dt ! CCN Concentration TENDENCY, ALL Contr. [1/kg/s]
	1064	cldfra(:,k)=CF___HOST(:,i) ! cloud fraction
	1065
	1066	END DO
	1067
	1068	! Heat flux at surface (diagnostic)
	1069	HsenSV_gpt(:) = -uts_SV_gpt(:) 1.e-3 1005. ! sensible heat flx in W/m2
	1070	HLatSV_gpt(:) = -uqs_SV_gpt(:) 1.e-3 2.5e6 ! latent heat flx in W/m2
	1071
	1072	! Question à Hubert : tendences de TKE non nécessaires ? Réponse: la TKE est mise à jour dans la physique, on n'a pas besoin de transporter la dérivée de TKE dans la dynamique.
	1073
	1074	! On désalloue les tableaux:
	1075	DEALLOCATE(MMaskHOST) ! Est-ce qu'il faut le faire ?
	1076
	1077	PRINT*, 'Control des tendances aprés la physique de MAR'
	1078
	1079	PRINT*, 'maxval(d_u)=',maxval(d_u)
	1080	PRINT*, 'maxval(d_v)=',maxval(d_v)
	1081	PRINT*, 'maxval(d_t)=',maxval(d_t)
	1082	PRINT*, 'maxval(paprs(:,1))=',maxval(paprs(:,1))
	1083	PRINT*, 'minval(d_u)=',minval(d_u)
	1084	PRINT*, 'minval(d_v)=',minval(d_v)
	1085	PRINT*, 'minval(d_t)=',minval(d_t)
	1086	PRINT*, 'minval(paprs(:,1))=',minval(paprs(:,1))
	1087	PRINT*,'Test snowCM'
	1088	PRINT*,'MAXVAL(snowCM)=',MAXVAL(snowCM)
	1089	PRINT*,'MAXVAL(rainCM)=',MAXVAL(rainCM)
	1090	PRINT*,'maxval(cldfra)=',maxval(cldfra)
	1091	PRINT*,'minval(cldfra)=',minval(cldfra)
	1092
	1093	!======================================================================
	1094	! Ecriture des sorties netcdf
	1095	!======================================================================
	1096
	1097	PRINT*, 'On ecrit les variables dans un fichier netcdf instantané'
	1098	! write some outputs:
	1099	if (debut) then ! On imprime les variables de sortie intemporelles seulement au début
	1100	!call iophys_ini
	1101	call histwrite_phy(nid_hist,.false.,"lonh",itau,lonh_HOST)
	1102	call histwrite_phy(nid_hist,.false.,"latr",itau,latr_HOST)
	1103	call histwrite_phy(nid_hist,.false.,"sst",itau,sst__HOST)
	1104	endif
	1105	if (modulo(itau,iwrite_phys)==0) then
	1106	call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
	1107	call histwrite_phy(nid_hist,.false.,"pplay",itau,pplay)
	1108	call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
	1109	call histwrite_phy(nid_hist,.false.,"qx_vap",itau,qx(:,:,ivap))
	1110	call histwrite_phy(nid_hist,.false.,"qx_liq",itau,qx(:,:,iliq))
	1111	call histwrite_phy(nid_hist,.false.,"qx_ice",itau,qx(:,:,iice))
	1112	call histwrite_phy(nid_hist,.false.,"u",itau,u)
	1113	call histwrite_phy(nid_hist,.false.,"v",itau,v)
	1114	call histwrite_phy(nid_hist,.false.,"d_t",itau,d_t)
	1115	call histwrite_phy(nid_hist,.false.,"d_u",itau,d_u)
	1116	call histwrite_phy(nid_hist,.false.,"d_v",itau,d_v)
	1117	call histwrite_phy(nid_hist,.false.,"snow",itau,SnowCM)
	1118	call histwrite_phy(nid_hist,.false.,"rain",itau,RainCM)
	1119	call histwrite_phy(nid_hist,.false.,"snowCP",itau,snowCP)
	1120	call histwrite_phy(nid_hist,.false.,"rainCP",itau,rainCP)
	1121	call histwrite_phy(nid_hist,.false.,"nebulosity",itau,cldfra)
	1122	endif
	1123
	1124	!XIOS
	1125	#ifdef CPP_XIOS
	1126	!$OMP MASTER
	1127	!On incrémente le pas de temps XIOS
	1128	!CALL wxios_update_calendar(itau)
	1129
	1130	!Et on écrit, avec la routine histwrite dédiée:
	1131	!CALL histwrite_phy("temperature",t)
	1132	!CALL histwrite_phy("u",u)
	1133	!CALL histwrite_phy("v",v)
	1134	!CALL histwrite_phy("ps",paprs(:,1))
	1135	!$OMP END MASTER
	1136	#endif
	1137	!
	1138	! On synchronise la fin de l'écriture du fichier à chaque pas de temps pour que le fichier soit OK même si le code plante...
	1139	CALL histsync(nid_hist)
	1140
	1141	PRINT*, 'Fin de physiq.f90'
	1142
	1143	return
	1144	end
	1145
	1146	! sub-routine inutilisée :
	1147
	1148	!SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	1149	! IMPLICIT none
	1150	!
	1151	! Tranformer une variable de la grille physique a
	1152	! la grille d'ecriture
	1153	! INTEGER nfield,nlon,iim,jjmp1, jjm
	1154	! REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	1155	! INTEGER i, n, ig
	1156	! jjm = jjmp1 - 1
	1157	! DO n = 1, nfield
	1158	! DO i=1,iim
	1159	! ecrit(i,n) = fi(1,n)
	1160	!, ecrit(i+jjm*iim,n) = fi(nlon,n)
	1161	! ENDDO
	1162	! DO ig = 1, nlon - 2
	1163	! ecrit(iim+ig,n) = fi(1+ig,n)
	1164	! ENDDO
	1165	! ENDDO
	1166	!RETURN
	1167	!END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format