Context Navigation

physiq_mod.F90 @ 5437

Last change on this file since 5437 was 3406, checked in by jghattas, 6 years ago
Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 165.2 KB

Rev	Line
[2418]	1	!
[1279]	2	! $Id: physiq_mod.F90 3406 2018-10-23 09:16:41Z evignon $
[2418]	3	!
[1862]	4	!#define IO_DEBUG
[2418]	5	MODULE physiq_mod
[766]	6
[2469]	7	IMPLICIT NONE
[2418]	8
	9	CONTAINS
	10
[2469]	11	SUBROUTINE physiq (nlon,nlev, &
	12	debut,lafin,pdtphys_, &
	13	paprs,pplay,pphi,pphis,presnivs, &
	14	u,v,rot,t,qx, &
	15	flxmass_w, &
	16	d_u, d_v, d_t, d_qx, d_ps)
[524]	17
[2769]	18	use assert_m, only: assert
[2469]	19	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
	20	histwrite, ju2ymds, ymds2ju, getin
	21	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
	22	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
[2783]	23	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
[2469]	24	USE write_field_phy
	25	USE dimphy
	26	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
[3406]	27	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
[2469]	28	USE mod_phys_lmdz_para
	29	USE iophy
	30	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
	31	USE phystokenc_mod, ONLY: offline, phystokenc
[2783]	32	USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
[2469]	33	USE vampir
	34	USE pbl_surface_mod, ONLY : pbl_surface
	35	USE change_srf_frac_mod
	36	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
[2788]	37	USE tropopause_m, ONLY: dyn_tropopause
[2630]	38	#ifdef CPP_Dust
	39	USE phytracr_spl_mod, ONLY: phytracr_spl
	40	#endif
[2606]	41	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
	42	! [Variables internes non sauvegardees de la physique]
	43	! Variables locales pour effectuer les appels en serie
	44	t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
	45	! Dynamic tendencies (diagnostics)
	46	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
	47	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
	48	! Physic tendencies
	49	d_t_con,d_q_con,d_u_con,d_v_con, &
	50	d_tr, & !! to be removed?? (jyg)
	51	d_t_wake,d_q_wake, &
	52	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
	53	d_t_ajsb,d_q_ajsb, &
	54	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
	55	d_t_ajs_w,d_q_ajs_w, &
	56	d_t_ajs_x,d_q_ajs_x, &
	57	!
[2705]	58	d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
[2606]	59	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
	60	d_t_lscst,d_q_lscst, &
	61	d_t_lscth,d_q_lscth, &
	62	plul_st,plul_th, &
	63	!
	64	d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
	65	d_t_vdf_w,d_q_vdf_w, &
	66	d_t_vdf_x,d_q_vdf_x, &
	67	d_ts, &
	68	!
	69	d_t_oli,d_u_oli,d_v_oli, &
	70	d_t_oro,d_u_oro,d_v_oro, &
[2897]	71	d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
[2606]	72	d_t_lif,d_u_lif,d_v_lif, &
	73	d_t_ec, &
	74	!
	75	du_gwd_hines,dv_gwd_hines,d_t_hin, &
	76	dv_gwd_rando,dv_gwd_front, &
	77	east_gwstress,west_gwstress, &
	78	d_q_ch4, &
	79	! Special RRTM
	80	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
	81	ZFLUP0,ZFSDN0,ZFSUP0, &
	82	!
	83	topswad_aero,solswad_aero, &
	84	topswai_aero,solswai_aero, &
	85	topswad0_aero,solswad0_aero, &
	86	!LW additional
	87	toplwad_aero,sollwad_aero, &
	88	toplwai_aero,sollwai_aero, &
	89	toplwad0_aero,sollwad0_aero, &
	90	!
	91	topsw_aero,solsw_aero, &
	92	topsw0_aero,solsw0_aero, &
	93	topswcf_aero,solswcf_aero, &
	94	tausum_aero,tau3d_aero, &
[2854]	95	drytausum_aero, &
[2606]	96	!
	97	!variables CFMIP2/CMIP5
	98	topswad_aerop, solswad_aerop, &
	99	topswai_aerop, solswai_aerop, &
	100	topswad0_aerop, solswad0_aerop, &
	101	topsw_aerop, topsw0_aerop, &
	102	solsw_aerop, solsw0_aerop, &
	103	topswcf_aerop, solswcf_aerop, &
	104	!LW diagnostics
	105	toplwad_aerop, sollwad_aerop, &
	106	toplwai_aerop, sollwai_aerop, &
	107	toplwad0_aerop, sollwad0_aerop, &
	108	!
	109	ptstar, pt0, slp, &
	110	!
	111	bils, &
	112	!
	113	cldh, cldl,cldm, cldq, cldt, &
	114	JrNt, &
	115	dthmin, evap, fder, plcl, plfc, &
	116	prw, prlw, prsw, &
	117	s_lcl, s_pblh, s_pblt, s_therm, &
	118	cdragm, cdragh, &
	119	zustar, zu10m, zv10m, rh2m, qsat2m, &
	120	zq2m, zt2m, weak_inversion, &
	121	zt2m_min_mon, zt2m_max_mon, & ! pour calcul_divers.h
	122	t2m_min_mon, t2m_max_mon, & ! pour calcul_divers.h
	123	!
	124	s_pblh_x, s_pblh_w, &
	125	s_lcl_x, s_lcl_w, &
	126	!
	127	slab_wfbils, tpot, tpote, &
	128	ue, uq, ve, vq, zxffonte, &
	129	zxfqcalving, zxfluxlat, &
	130	zxrunofflic, &
	131	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
	132	rain_lsc, rain_num, &
	133	!
	134	sens_x, sens_w, &
	135	zxfluxlat_x, zxfluxlat_w, &
	136	!
	137	dtvdf_x, dtvdf_w, &
	138	dqvdf_x, dqvdf_w, &
	139	pbl_tke_input, &
	140	t_therm, q_therm, u_therm, v_therm, &
	141	cdragh_x, cdragh_w, &
	142	cdragm_x, cdragm_w, &
	143	kh, kh_x, kh_w, &
	144	!
[2730]	145	wake_k, &
[2606]	146	ale_wake, alp_wake, &
[2635]	147	wake_h, wake_omg, &
	148	! tendencies of delta T and delta q:
	149	d_deltat_wk, d_deltaq_wk, & ! due to wakes
	150	d_deltat_wk_gw, d_deltaq_wk_gw, & ! due to wake induced gravity waves
	151	d_deltat_vdf, d_deltaq_vdf, & ! due to vertical diffusion
	152	d_deltat_the, d_deltaq_the, & ! due to thermals
	153	d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
	154	! tendencies of wake fractional area and wake number per unit area:
	155	d_s_wk, d_dens_wk, & ! due to wakes
	156	!!! d_s_vdf, d_dens_vdf, & ! due to vertical diffusion
	157	!!! d_s_the, d_dens_the, & ! due to thermals
	158	!
[2877]	159	ptconv, &
[2824]	160	wbeff, convoccur, zmax_th, &
[2606]	161	sens, flwp, fiwp, &
	162	ale_bl_stat,alp_bl_conv,alp_bl_det, &
	163	alp_bl_fluct_m,alp_bl_fluct_tke, &
	164	alp_bl_stat, n2, s2, &
	165	proba_notrig, random_notrig, &
	166	!
	167	dnwd, dnwd0, &
	168	upwd, omega, &
	169	epmax_diag, &
	170	ep, &
	171	cldemi, &
	172	cldfra, cldtau, fiwc, &
	173	fl, re, flwc, &
	174	ref_liq, ref_ice, theta, &
	175	ref_liq_pi, ref_ice_pi, &
[2635]	176	zphi, zx_rh, &
[2606]	177	pmfd, pmfu, &
	178	!
	179	t2m, fluxlat, &
	180	fsollw, evap_pot, &
	181	fsolsw, wfbils, wfbilo, &
[2670]	182	wfevap, wfrain, wfsnow, &
[2606]	183	pmflxr, pmflxs, prfl, &
	184	psfl, fraca, Vprecip, &
	185	zw2, &
	186
	187	fluxu, fluxv, &
	188	fluxt, &
	189
	190	uwriteSTD, vwriteSTD, & !pour calcul_STDlev.h
	191	wwriteSTD, phiwriteSTD, & !pour calcul_STDlev.h
	192	qwriteSTD, twriteSTD, rhwriteSTD, & !pour calcul_STDlev.h
	193
	194	wdtrainA, wdtrainM, &
	195	beta_prec, &
	196	rneb, &
[2968]	197	zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
[2606]	198	!
[2469]	199	USE phys_state_var_mod ! Variables sauvegardees de la physique
[2630]	200	#ifdef CPP_Dust
	201	USE phys_output_write_spl_mod
	202	#else
[2469]	203	USE phys_output_var_mod ! Variables pour les ecritures des sorties
[2630]	204	#endif
	205
[2469]	206	USE phys_output_write_mod
	207	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
	208	USE phys_output_mod
	209	USE phys_output_ctrlout_mod
	210	use open_climoz_m, only: open_climoz ! ozone climatology from a file
[2788]	211	use regr_pr_time_av_m, only: regr_pr_time_av
[2469]	212	use netcdf95, only: nf95_close
	213	!IM for NMC files
	214	! use netcdf, only: nf90_fill_real
[2997]	215	use netcdf, only: nf90_fill_real
[2469]	216	use mod_phys_lmdz_mpi_data, only: is_mpi_root
	217	USE aero_mod
	218	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
	219	use conf_phys_m, only: conf_phys
	220	use radlwsw_m, only: radlwsw
	221	use phyaqua_mod, only: zenang_an
	222	USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
	223	start_time, pdtphys, day_ini
	224	USE tracinca_mod, ONLY: config_inca
[2271]	225	#ifdef CPP_XIOS
[2469]	226	USE wxios, ONLY: missing_val, missing_val_omp
[2679]	227	USE xios, ONLY: xios_get_field_attr, xios_field_is_active
[2271]	228	#endif
[1565]	229	#ifdef REPROBUS
[2469]	230	USE CHEM_REP, ONLY : Init_chem_rep_xjour
[1565]	231	#endif
[2469]	232	USE indice_sol_mod
	233	USE phytrac_mod, ONLY : phytrac
[782]	234
[2009]	235	#ifdef CPP_RRTM
[2517]	236	USE YOERAD, ONLY : NRADLP
[2524]	237	USE YOESW, ONLY : RSUN
[2009]	238	#endif
[2469]	239	USE ioipsl_getin_p_mod, ONLY : getin_p
[2003]	240
[2651]	241	#ifndef CPP_XIOS
[2590]	242	USE paramLMDZ_phy_mod
[2651]	243	#endif
[2294]	244
[2611]	245	USE cmp_seri_mod
[2902]	246	USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
[2799]	247	& fl_ebil, fl_cor_ebil
[2611]	248
[2469]	249	!IM stations CFMIP
	250	USE CFMIP_point_locations
	251	use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
	252	use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
[2832]	253	USE VERTICAL_LAYERS_MOD, ONLY: aps,bps
[3406]	254	USE etat0_limit_unstruct_mod
	255	#ifdef CPP_XIOS
	256	USE xios, ONLY: xios_update_calendar, xios_context_finalize
	257	#endif
	258	USE climoz_mod
	259	USE limit_read_mod, ONLY : init_limit_read
	260	USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
[1938]	261
[2832]	262
[2469]	263	IMPLICIT none
	264	!>======================================================================
	265	!!
	266	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	267	!!
	268	!! Objet: Moniteur general de la physique du modele
	269	!!AA Modifications quant aux traceurs :
	270	!!AA - uniformisation des parametrisations ds phytrac
	271	!!AA - stockage des moyennes des champs necessaires
	272	!!AA en mode traceur off-line
	273	!!======================================================================
	274	!! CLEFS CPP POUR LES IO
	275	!! =====================
[1352]	276	#define histNMC
[2469]	277	!!======================================================================
	278	!! modif ( P. Le Van , 12/10/98 )
	279	!!
	280	!! Arguments:
	281	!!
	282	!! nlon----input-I-nombre de points horizontaux
	283	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
	284	!! debut---input-L-variable logique indiquant le premier passage
	285	!! lafin---input-L-variable logique indiquant le dernier passage
	286	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
	287	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
	288	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
	289	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
	290	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	291	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	292	!! pphis---input-R-geopotentiel du sol
	293	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	294	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
	295	!! v-------input-R-vitesse Y (de S a N) en m/s
	296	!! t-------input-R-temperature (K)
	297	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	298	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	299	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
[2496]	300	!! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
	301	!! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
[2469]	302	!! flxmass_w -input-R- flux de masse verticale
	303	!! d_u-----output-R-tendance physique de "u" (m/s/s)
	304	!! d_v-----output-R-tendance physique de "v" (m/s/s)
	305	!! d_t-----output-R-tendance physique de "t" (K/s)
	306	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	307	!! d_ps----output-R-tendance physique de la pression au sol
	308	!!======================================================================
	309	integer jjmp1
	310	! parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
	311	! integer iip1
	312	! parameter (iip1=iim+1)
[782]	313
[2469]	314	include "regdim.h"
	315	include "dimsoil.h"
	316	include "clesphys.h"
	317	include "thermcell.h"
[3011]	318	include "dimpft.h"
[2469]	319	!======================================================================
	320	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	321	PARAMETER (ok_cvl=.TRUE.)
	322	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	323	PARAMETER (ok_gust=.FALSE.)
	324	integer iflag_radia ! active ou non le rayonnement (MPL)
	325	save iflag_radia
	326	!$OMP THREADPRIVATE(iflag_radia)
	327	!======================================================================
	328	LOGICAL check ! Verifier la conservation du modele en eau
	329	PARAMETER (check=.FALSE.)
	330	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	331	PARAMETER (ok_stratus=.FALSE.)
	332	!======================================================================
	333	REAL amn, amx
	334	INTEGER igout
	335	!======================================================================
	336	! Clef controlant l'activation du cycle diurne:
	337	! en attente du codage des cles par Fred
	338	INTEGER iflag_cycle_diurne
	339	PARAMETER (iflag_cycle_diurne=1)
	340	!======================================================================
	341	! Modele thermique du sol, a activer pour le cycle diurne:
	342	!cc LOGICAL soil_model
	343	!cc PARAMETER (soil_model=.FALSE.)
	344	!======================================================================
	345	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	346	! le calcul du rayonnement est celle apres la precipitation des nuages.
	347	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	348	! la condensation et la precipitation. Cette cle augmente les impacts
	349	! radiatifs des nuages.
	350	!cc LOGICAL new_oliq
	351	!cc PARAMETER (new_oliq=.FALSE.)
	352	!======================================================================
	353	! Clefs controlant deux parametrisations de l'orographie:
	354	!c LOGICAL ok_orodr
	355	!cc PARAMETER (ok_orodr=.FALSE.)
	356	!cc LOGICAL ok_orolf
	357	!cc PARAMETER (ok_orolf=.FALSE.)
	358	!======================================================================
	359	LOGICAL ok_journe ! sortir le fichier journalier
	360	save ok_journe
	361	!$OMP THREADPRIVATE(ok_journe)
	362	!
	363	LOGICAL ok_mensuel ! sortir le fichier mensuel
	364	save ok_mensuel
	365	!$OMP THREADPRIVATE(ok_mensuel)
	366	!
	367	LOGICAL ok_instan ! sortir le fichier instantane
	368	save ok_instan
	369	!$OMP THREADPRIVATE(ok_instan)
	370	!
	371	LOGICAL ok_LES ! sortir le fichier LES
	372	save ok_LES
	373	!$OMP THREADPRIVATE(ok_LES)
	374	!
	375	LOGICAL callstats ! sortir le fichier stats
	376	save callstats
	377	!$OMP THREADPRIVATE(callstats)
	378	!
	379	LOGICAL ok_region ! sortir le fichier regional
	380	PARAMETER (ok_region=.FALSE.)
	381	!======================================================================
	382	real seuil_inversion
	383	save seuil_inversion
	384	!$OMP THREADPRIVATE(seuil_inversion)
	385	integer iflag_ratqs
	386	save iflag_ratqs
	387	!$OMP THREADPRIVATE(iflag_ratqs)
	388	real facteur
[1507]	389
[2469]	390	REAL wmax_th(klon)
	391	REAL tau_overturning_th(klon)
[878]	392
[2469]	393	integer lmax_th(klon)
	394	integer limbas(klon)
	395	real ratqscth(klon,klev)
	396	real ratqsdiff(klon,klev)
	397	real zqsatth(klon,klev)
[878]	398
[2469]	399	!======================================================================
	400	!
	401	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	402	PARAMETER (ivap=1)
	403	INTEGER iliq ! indice de traceurs pour eau liquide
	404	PARAMETER (iliq=2)
	405	!CR: on ajoute la phase glace
	406	INTEGER isol ! indice de traceurs pour eau glace
	407	PARAMETER (isol=3)
	408	!
	409	!
	410	! Variables argument:
	411	!
	412	INTEGER nlon
	413	INTEGER nlev
	414	REAL,INTENT(IN) :: pdtphys_
	415	! NB: pdtphys to be used in physics is in time_phylmdz_mod
	416	LOGICAL debut, lafin
	417	REAL paprs(klon,klev+1)
	418	REAL pplay(klon,klev)
	419	REAL pphi(klon,klev)
	420	REAL pphis(klon)
	421	REAL presnivs(klev)
[2799]	422	!JLD REAL znivsig(klev)
	423	!JLD real pir
[719]	424
[2469]	425	REAL u(klon,klev)
	426	REAL v(klon,klev)
[2333]	427
[2469]	428	REAL, intent(in):: rot(klon, klev)
	429	! relative vorticity, in s-1, needed for frontal waves
[2333]	430
[2469]	431	REAL t(klon,klev),thetal(klon,klev)
	432	! thetal: ligne suivante a decommenter si vous avez les fichiers
	433	! MPL 20130625
	434	! fth_fonctions.F90 et parkind1.F90
	435	! sinon thetal=theta
	436	! REAL fth_thetae,fth_thetav,fth_thetal
	437	REAL qx(klon,klev,nqtot)
	438	REAL flxmass_w(klon,klev)
	439	REAL d_u(klon,klev)
	440	REAL d_v(klon,klev)
	441	REAL d_t(klon,klev)
	442	REAL d_qx(klon,klev,nqtot)
	443	REAL d_ps(klon)
[2897]	444	! variables pour tend_to_tke
	445	REAL duadd(klon,klev)
	446	REAL dvadd(klon,klev)
	447	REAL dtadd(klon,klev)
	448
[2469]	449	! Variables pour le transport convectif
	450	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
	451	real wght_cvfd(klon,klev)
[2271]	452	#ifndef CPP_XIOS
[2997]	453	REAL, SAVE :: missing_val=nf90_fill_real
[2271]	454	#endif
[2469]	455	! Variables pour le lessivage convectif
	456	! RomP >>>
	457	real phi2(klon,klev,klev)
	458	real d1a(klon,klev),dam(klon,klev)
[2481]	459	real ev(klon,klev)
[2469]	460	real clw(klon,klev),elij(klon,klev,klev)
	461	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
	462	! RomP <<<
	463	!IM definition dynamique o_trac dans phys_output_open
	464	! type(ctrl_out) :: o_trac(nqtot)
[524]	465
[2469]	466	! variables a une pression donnee
	467	!
	468	include "declare_STDlev.h"
	469	!
	470	!
	471	include "radopt.h"
	472	!
	473	!
	474	INTEGER debug
	475	INTEGER n
	476	!ym INTEGER npoints
	477	!ym PARAMETER(npoints=klon)
	478	!
	479	INTEGER nregISCtot
	480	PARAMETER(nregISCtot=1)
	481	!
	482	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
	483	! sur 1 region rectangulaire y compris pour 1 point
	484	! imin_debut : indice minimum de i; nbpti : nombre de points en
	485	! direction i (longitude)
	486	! jmin_debut : indice minimum de j; nbptj : nombre de points en
	487	! direction j (latitude)
[2799]	488	!JLD INTEGER imin_debut, nbpti
	489	!JLD INTEGER jmin_debut, nbptj
[2469]	490	!IM: region='3d' <==> sorties en global
	491	CHARACTER*3 region
	492	PARAMETER(region='3d')
	493	logical ok_hf
	494	!
	495	save ok_hf
	496	!$OMP THREADPRIVATE(ok_hf)
[524]	497
[2469]	498	INTEGER,PARAMETER :: longcles=20
	499	REAL,SAVE :: clesphy0(longcles)
	500	!$OMP THREADPRIVATE(clesphy0)
	501	!
	502	! Variables propres a la physique
	503	INTEGER itap
	504	SAVE itap ! compteur pour la physique
	505	!$OMP THREADPRIVATE(itap)
[2235]	506
[2469]	507	INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys
	508	!$OMP THREADPRIVATE(abortphy)
	509	!
	510	REAL,save :: solarlong0
	511	!$OMP THREADPRIVATE(solarlong0)
[987]	512
[2469]	513	!
	514	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	515	!
	516	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	517	REAL zulow(klon),zvlow(klon)
	518	!
	519	INTEGER igwd,idx(klon),itest(klon)
	520	!
	521	! REAL,allocatable,save :: run_off_lic_0(:)
	522	! !$OMP THREADPRIVATE(run_off_lic_0)
	523	!ym SAVE run_off_lic_0
	524	!KE43
	525	! Variables liees a la convection de K. Emanuel (sb):
	526	!
	527	REAL bas, top ! cloud base and top levels
	528	SAVE bas
	529	SAVE top
	530	!$OMP THREADPRIVATE(bas, top)
	531	!------------------------------------------------------------------
	532	! Upmost level reached by deep convection and related variable (jyg)
	533	!
	534	INTEGER izero
	535	INTEGER k_upper_cv
	536	!------------------------------------------------------------------
	537	!
	538	!==========================================================================
	539	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema
	540	!de convection avec poches froides
	541	! Variables li\'ees \`a la poche froide (jyg)
[879]	542
[2469]	543	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
	544	!
	545	REAL wape_prescr, fip_prescr
	546	INTEGER it_wape_prescr
	547	SAVE wape_prescr, fip_prescr, it_wape_prescr
	548	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
	549	!
	550	! variables supplementaires de concvl
	551	REAL Tconv(klon,klev)
	552	REAL sij(klon,klev,klev)
[2812]	553	!! !
	554	!! ! variables pour tester la conservation de l'energie dans concvl
	555	!! REAL, DIMENSION(klon,klev) :: d_t_con_sat
	556	!! REAL, DIMENSION(klon,klev) :: d_q_con_sat
	557	!! REAL, DIMENSION(klon,klev) :: dql_sat
[970]	558
[2469]	559	real, save :: alp_bl_prescr=0.
	560	real, save :: ale_bl_prescr=0.
[979]	561
[2469]	562	real, save :: wake_s_min_lsp=0.1
[1516]	563
[2469]	564	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
	565	!$OMP THREADPRIVATE(wake_s_min_lsp)
[970]	566
[1516]	567
[2469]	568	real ok_wk_lsp(klon)
[1516]	569
[2469]	570	!RC
	571	! Variables li\'ees \`a la poche froide (jyg et rr)
[879]	572
[2635]	573	INTEGER, SAVE :: iflag_wake_tend ! wake: if =0, then wake state variables are
	574	! updated within calwake
	575	!$OMP THREADPRIVATE(iflag_wake_tend)
[3000]	576	INTEGER, SAVE :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
	577	! power provided by the wakes; else, Alp_wk is the
	578	! lifting power conditionned on the presence of a
	579	! gust-front in the grid cell.
	580	!$OMP THREADPRIVATE(iflag_alp_wk_cond)
[2635]	581	REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
	582	REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
[879]	583
[2469]	584	REAL wake_dth(klon,klev) ! wake : temp pot difference
[879]	585
[2469]	586	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta
	587	! transported by LS omega
	588	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of
	589	! large scale omega
	590	REAL wake_dtKE(klon,klev) ! Wake : differential heating
	591	! (wake - unpertubed) CONV
	592	REAL wake_dqKE(klon,klev) ! Wake : differential moistening
	593	! (wake - unpertubed) CONV
	594	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
	595	REAL wake_spread(klon,klev) ! spreading term in wake_delt
	596	!
	597	!pourquoi y'a pas de save??
	598	!
[2730]	599	!!! INTEGER, SAVE, DIMENSION(klon) :: wake_k
	600	!!! !$OMP THREADPRIVATE(wake_k)
[2469]	601	!
	602	!jyg<
	603	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
	604	!>jyg
[879]	605
[3000]	606	REAL wake_fip_0(klon) ! Average Front Incoming Power (unconditionned)
[2469]	607	REAL wake_gfl(klon) ! Gust Front Length
[2635]	608	!!! REAL wake_dens(klon) ! moved to phys_state_var_mod
[2469]	609	!
	610	!
	611	REAL dt_dwn(klon,klev)
	612	REAL dq_dwn(klon,klev)
	613	REAL M_dwn(klon,klev)
	614	REAL M_up(klon,klev)
	615	REAL dt_a(klon,klev)
	616	REAL dq_a(klon,klev)
	617	REAL d_t_adjwk(klon,klev) !jyg
	618	REAL d_q_adjwk(klon,klev) !jyg
	619	LOGICAL,SAVE :: ok_adjwk=.FALSE.
	620	!$OMP THREADPRIVATE(ok_adjwk)
[2882]	621	INTEGER,SAVE :: iflag_adjwk=0 !jyg
	622	!$OMP THREADPRIVATE(iflag_adjwk) !jyg
[2657]	623	REAL,SAVE :: oliqmax=999.,oicemax=999.
	624	!$OMP THREADPRIVATE(oliqmax,oicemax)
[2469]	625	REAL, SAVE :: alp_offset
	626	!$OMP THREADPRIVATE(alp_offset)
[2897]	627
[2469]	628	!
	629	!RR:fin declarations poches froides
	630	!==========================================================================
[1032]	631
[2469]	632	REAL ztv(klon,klev),ztva(klon,klev)
	633	REAL zpspsk(klon,klev)
	634	REAL ztla(klon,klev),zqla(klon,klev)
	635	REAL zthl(klon,klev)
[1638]	636
[2469]	637	!cc nrlmd le 10/04/2012
[1638]	638
[2469]	639	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
	640	!---Propri\'et\'es du thermiques au LCL
	641	real zlcl_th(klon) ! Altitude du LCL calcul\'e
	642	! continument (pcon dans
	643	! thermcell_main.F90)
	644	real fraca0(klon) ! Fraction des thermiques au LCL
	645	real w0(klon) ! Vitesse des thermiques au LCL
	646	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
	647	real tke0(klon,klev+1) ! TKE au d\'ebut du pas de temps
	648	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
	649	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
[1638]	650
[2799]	651	!JLD !---D\'eclenchement stochastique
	652	!JLD integer :: tau_trig(klon)
[1638]	653
[2469]	654	REAL,SAVE :: random_notrig_max=1.
	655	!$OMP THREADPRIVATE(random_notrig_max)
[2294]	656
[2469]	657	!--------Statistical Boundary Layer Closure: ALP_BL--------
	658	!---Profils de TKE dans et hors du thermique
	659	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
	660	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
[1638]	661
[2897]	662	!-------Activer les tendances de TKE due a l'orograp??ie---------
	663	INTEGER, SAVE :: addtkeoro
	664	!$OMP THREADPRIVATE(addtkeoro)
	665	REAL, SAVE :: alphatkeoro
	666	!$OMP THREADPRIVATE(alphatkeoro)
	667	LOGICAL, SAVE :: smallscales_tkeoro
	668	!$OMP THREADPRIVATE(smallscales_tkeoro)
[1638]	669
[2897]	670
	671
[2469]	672	!cc fin nrlmd le 10/04/2012
[782]	673
[2469]	674	! Variables locales pour la couche limite (al1):
	675	!
	676	!Al1 REAL pblh(klon) ! Hauteur de couche limite
	677	!Al1 SAVE pblh
	678	!34EK
	679	!
	680	! Variables locales:
	681	!
	682	!AA
	683	!AA Pour phytrac
	684	REAL u1(klon) ! vents dans la premiere couche U
	685	REAL v1(klon) ! vents dans la premiere couche V
[524]	686
[2469]	687	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
	688	!@$$ PARAMETER (offline=.false.)
	689	!@$$ INTEGER physid
	690	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	691	REAL frac_nucl(klon,klev) ! idem (nucleation)
	692	! RomP >>>
	693	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
	694	! RomP <<<
[2832]	695	REAL :: calday
[2068]	696
[2469]	697	!IM cf FH pour Tiedtke 080604
	698	REAL rain_tiedtke(klon),snow_tiedtke(klon)
	699	!
	700	!IM 050204 END
	701	REAL devap(klon) ! evaporation et sa derivee
	702	REAL dsens(klon) ! chaleur sensible et sa derivee
[1279]	703
[2469]	704	!
	705	! Conditions aux limites
	706	!
	707	!
	708	REAL :: day_since_equinox
	709	! Date de l'equinoxe de printemps
	710	INTEGER, parameter :: mth_eq=3, day_eq=21
	711	REAL :: jD_eq
[1279]	712
[2469]	713	LOGICAL, parameter :: new_orbit = .true.
[524]	714
[2469]	715	!
	716	INTEGER lmt_pas
	717	SAVE lmt_pas ! frequence de mise a jour
	718	!$OMP THREADPRIVATE(lmt_pas)
	719	real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
	720	! (column-density of mass of air in a cell, in kg m-2)
	721	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
[1797]	722
[2469]	723	!IM sorties
	724	REAL un_jour
	725	PARAMETER(un_jour=86400.)
	726	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
	727	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
	728	!$OMP THREADPRIVATE(itapm1)
	729	!======================================================================
	730	!
	731	! Declaration des procedures appelees
	732	!
	733	EXTERNAL angle ! calculer angle zenithal du soleil
	734	EXTERNAL alboc ! calculer l'albedo sur ocean
	735	EXTERNAL ajsec ! ajustement sec
	736	EXTERNAL conlmd ! convection (schema LMD)
	737	!KE43
	738	EXTERNAL conema3 ! convect4.3
	739	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	740	!AA
	741	! JBM (3/14) fisrtilp_tr not loaded
	742	! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	743	! ! stockage des coefficients necessaires au
	744	! ! lessivage OFF-LINE et ON-LINE
	745	EXTERNAL hgardfou ! verifier les temperatures
	746	EXTERNAL nuage ! calculer les proprietes radiatives
	747	!C EXTERNAL o3cm ! initialiser l'ozone
	748	EXTERNAL orbite ! calculer l'orbite terrestre
	749	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	750	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	751	EXTERNAL suphel ! initialiser certaines constantes
	752	EXTERNAL transp ! transport total de l'eau et de l'energie
	753	!IM
	754	EXTERNAL haut2bas !variables de haut en bas
	755	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
	756	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
	757	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
	758	! EXTERNAL moyglo_aire
	759	! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
	760	! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
	761	!
	762	!
	763	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	764	! Local variables
	765	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	766	!
	767	REAL rhcl(klon,klev) ! humiditi relative ciel clair
	768	REAL dialiq(klon,klev) ! eau liquide nuageuse
	769	REAL diafra(klon,klev) ! fraction nuageuse
	770	REAL cldliq(klon,klev) ! eau liquide nuageuse
	771	!
	772	!XXX PB
	773	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	774	!
	775	REAL zxfluxt(klon, klev)
	776	REAL zxfluxq(klon, klev)
	777	REAL zxfluxu(klon, klev)
	778	REAL zxfluxv(klon, klev)
[1797]	779
[2469]	780	! Le rayonnement n'est pas calcule tous les pas, il faut donc
	781	! sauvegarder les sorties du rayonnement
	782	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
	783	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
	784	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
	785	!
	786	INTEGER itaprad
	787	SAVE itaprad
	788	!$OMP THREADPRIVATE(itaprad)
	789	!
	790	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	791	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
	792	!
[2799]	793	#ifdef INCA
[2469]	794	REAL zxsnow_dummy(klon)
[2799]	795	#endif
[2469]	796	REAL zsav_tsol(klon)
	797	!
	798	REAL dist, rmu0(klon), fract(klon)
	799	REAL zrmu0(klon), zfract(klon)
	800	REAL zdtime, zdtime1, zdtime2, zlongi
	801	!
	802	REAL qcheck
	803	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	804	LOGICAL zx_ajustq
	805	!
[2799]	806	REAL za
	807	REAL zx_t, zx_qs, zdelta, zcor
[2469]	808	real zqsat(klon,klev)
	809	!
[2897]	810	INTEGER i, k, iq, j, nsrf, ll, l
[2469]	811	!
	812	REAL t_coup
	813	PARAMETER (t_coup=234.0)
[1797]	814
[2469]	815	!ym A voir plus tard !!
	816	!ym REAL zx_relief(iim,jjmp1)
	817	!ym REAL zx_aire(iim,jjmp1)
	818	!
	819	! Grandeurs de sorties
	820	REAL s_capCL(klon)
	821	REAL s_oliqCL(klon), s_cteiCL(klon)
	822	REAL s_trmb1(klon), s_trmb2(klon)
	823	REAL s_trmb3(klon)
[2707]	824
	825	! La convection n'est pas calculee tous les pas, il faut donc
	826	! sauvegarder les sorties de la convection
	827	!ym SAVE
	828	!ym SAVE
	829	!ym SAVE
	830	!
[2730]	831	INTEGER itapcv, itapwk
	832	SAVE itapcv, itapwk
	833	!$OMP THREADPRIVATE(itapcv, itapwk)
[2707]	834
[2469]	835	!KE43
	836	! Variables locales pour la convection de K. Emanuel (sb):
[524]	837
[2469]	838	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	839	CHARACTER*40 capemaxcels !max(CAPE)
[1412]	840
[2469]	841	REAL rflag(klon) ! flag fonctionnement de convect
	842	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
[904]	843
[2469]	844	! -- convect43:
	845	INTEGER ntra ! nb traceurs pour convect4.3
	846	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	847	REAL dplcldt(klon), dplcldr(klon)
	848	!? . condm_con(klon,klev),conda_con(klon,klev),
	849	!? . mr_con(klon,klev),ep_con(klon,klev)
	850	!? . ,sadiab(klon,klev),wadiab(klon,klev)
	851	! --
	852	!34EK
	853	!
	854	! Variables du changement
	855	!
	856	! con: convection
	857	! lsc: condensation a grande echelle (Large-Scale-Condensation)
	858	! ajs: ajustement sec
	859	! eva: evaporation de l'eau liquide nuageuse
	860	! vdf: couche limite (Vertical DiFfusion)
[2611]	861	!
[2469]	862	! tendance nulles
[2812]	863	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
	864	REAL, dimension(klon) :: dsig0, ddens0
	865	INTEGER, dimension(klon) :: wkoccur1
[2801]	866	! tendance buffer pour appel de add_phys_tend
	867	REAL, DIMENSION(klon,klev) :: d_q_ch4_dtime
[2611]	868	!
	869	! Flag pour pouvoir ne pas ajouter les tendances.
	870	! Par defaut, les tendances doivente etre ajoutees et
	871	! flag_inhib_tend = 0
	872	! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
	873	! croissant de print quand la valeur du flag augmente
	874	!!! attention, ce flag doit etre change avec prudence !!!
	875	INTEGER :: flag_inhib_tend = 0 ! 0 is the default value
	876	!! INTEGER :: flag_inhib_tend = 2
[524]	877
[2469]	878	!
	879	!********************************************************
	880	! declarations
[524]	881
[2469]	882	!********************************************************
	883	!IM 081204 END
	884	!
	885	REAL pen_u(klon,klev), pen_d(klon,klev)
	886	REAL pde_u(klon,klev), pde_d(klon,klev)
	887	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	888	!
	889	REAL ratqsc(klon,klev)
	890	real ratqsbas,ratqshaut,tau_ratqs
	891	save ratqsbas,ratqshaut,tau_ratqs
	892	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
[2534]	893	REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
	894	!$OMP THREADPRIVATE(ratqsp0, ratqsdp)
[644]	895
[2469]	896	! Parametres lies au nouveau schema de nuages (SB, PDF)
	897	real fact_cldcon
	898	real facttemps
	899	logical ok_newmicro
	900	save ok_newmicro
	901	!$OMP THREADPRIVATE(ok_newmicro)
	902	!real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
	903	save fact_cldcon,facttemps
	904	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
[524]	905
[2469]	906	integer iflag_cld_th
	907	save iflag_cld_th
	908	!$OMP THREADPRIVATE(iflag_cld_th)
[2877]	909	!IM logical ptconv(klon,klev) !passe dans phys_local_var_mod
[2469]	910	!IM cf. AM 081204 BEG
	911	logical ptconvth(klon,klev)
	912	!IM cf. AM 081204 END
	913	!
	914	! Variables liees a l'ecriture de la bande histoire physique
	915	!
	916	!======================================================================
	917	!
[2068]	918
[2469]	919	!
[2799]	920	!JLD integer itau_w ! pas de temps ecriture = itap + itau_phy
[2469]	921	!
	922	!
	923	! Variables locales pour effectuer les appels en serie
	924	!
	925	!IM RH a 2m (la surface)
	926	REAL Lheat
[524]	927
[2469]	928	INTEGER length
	929	PARAMETER ( length = 100 )
	930	REAL tabcntr0( length )
	931	!
[2799]	932	!JLD INTEGER ndex2d(nbp_lon*nbp_lat)
[2469]	933	!IM
	934	!
	935	!IM AMIP2 BEG
[2799]	936	!JLD REAL moyglo, mountor
[2469]	937	!IM 141004 BEG
	938	REAL zustrdr(klon), zvstrdr(klon)
	939	REAL zustrli(klon), zvstrli(klon)
	940	REAL zustrph(klon), zvstrph(klon)
	941	REAL aam, torsfc
	942	!IM 141004 END
	943	!IM 190504 BEG
	944	! INTEGER imp1jmp1
	945	! PARAMETER(imp1jmp1=(iim+1)*jjmp1)
	946	!ym A voir plus tard
	947	! REAL zx_tmp((nbp_lon+1)*nbp_lat)
	948	! REAL airedyn(nbp_lon+1,nbp_lat)
	949	!IM 190504 END
[2799]	950	!JLD LOGICAL ok_msk
	951	!JLD REAL msk(klon)
[2469]	952	!ym A voir plus tard
	953	!ym REAL zm_wo(jjmp1, klev)
	954	!IM AMIP2 END
	955	!
	956	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	957	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
[2799]	958	!JLD REAL zx_tmp_2d(nbp_lon,nbp_lat)
	959	!JLD REAL zx_lon(nbp_lon,nbp_lat)
	960	!JLD REAL zx_lat(nbp_lon,nbp_lat)
[2469]	961	!
[2630]	962	INTEGER nid_ctesGCM
	963	SAVE nid_ctesGCM
	964	!$OMP THREADPRIVATE(nid_ctesGCM)
[2469]	965	!
	966	!IM 280405 BEG
	967	! INTEGER nid_bilKPins, nid_bilKPave
	968	! SAVE nid_bilKPins, nid_bilKPave
	969	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
	970	!
	971	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
	972	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
	973	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
	974	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
	975	!
[2799]	976	!JLD REAL zjulian
	977	!JLD SAVE zjulian
	978	!JLD!$OMP THREADPRIVATE(zjulian)
[2590]	979
[2799]	980	!JLD INTEGER nhori, nvert
	981	!JLD REAL zsto
	982	!JLD REAL zstophy, zout
[2068]	983
[2469]	984	character*20 modname
	985	character*80 abort_message
	986	logical, save :: ok_sync, ok_sync_omp
	987	!$OMP THREADPRIVATE(ok_sync)
	988	real date0
[524]	989
[2469]	990	! essai writephys
	991	integer fid_day, fid_mth, fid_ins
	992	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	993	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	994	parameter (prof2d_on = 1, prof3d_on = 2, &
	995	prof2d_av = 3, prof3d_av = 4)
	996	REAL ztsol(klon)
	997	REAL q2m(klon,nbsrf) ! humidite a 2m
[524]	998
[2469]	999	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
	1000	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
[2799]	1001	CHARACTER*40 tinst, tave
[2469]	1002	REAL cldtaupi(klon,klev) ! Cloud optical thickness for
	1003	! pre-industrial (pi) aerosols
[524]	1004
[2863]	1005	INTEGER :: naero
[2469]	1006	! Aerosol optical properties
	1007	CHARACTER*4, DIMENSION(naero_grp) :: rfname
	1008	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass
	1009	! concentration
	1010	! for all soluble
	1011	! aerosols[ug/m3]
	1012	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi
	1013	! - " - (pre-industrial value)
[1279]	1014
[2469]	1015	! Parameters
	1016	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
[2738]	1017	LOGICAL ok_alw ! Apply aerosol LW effect or not
[2469]	1018	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
	1019	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
[2738]	1020	SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
	1021	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
[2469]	1022	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
	1023	! false : lecture des aerosol dans un fichier
	1024	!$OMP THREADPRIVATE(aerosol_couple)
	1025	INTEGER, SAVE :: flag_aerosol
	1026	!$OMP THREADPRIVATE(flag_aerosol)
[2644]	1027	LOGICAL, SAVE :: flag_bc_internal_mixture
	1028	!$OMP THREADPRIVATE(flag_bc_internal_mixture)
[2469]	1029	LOGICAL, SAVE :: new_aod
	1030	!$OMP THREADPRIVATE(new_aod)
	1031	!
	1032	!--STRAT AEROSOL
[2530]	1033	INTEGER, SAVE :: flag_aerosol_strat
[2469]	1034	!$OMP THREADPRIVATE(flag_aerosol_strat)
	1035	!c-fin STRAT AEROSOL
	1036	!
	1037	! Declaration des constantes et des fonctions thermodynamiques
	1038	!
	1039	LOGICAL,SAVE :: first=.true.
	1040	!$OMP THREADPRIVATE(first)
[1279]	1041
[2788]	1042	! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
	1043	! Note that pressure vectors are in Pa and in stricly ascending order
	1044	INTEGER,SAVE :: read_climoz ! Read ozone climatology
[2469]	1045	! (let it keep the default OpenMP shared attribute)
	1046	! Allowed values are 0, 1 and 2
	1047	! 0: do not read an ozone climatology
	1048	! 1: read a single ozone climatology that will be used day and night
	1049	! 2: read two ozone climatologies, the average day and night
	1050	! climatology and the daylight climatology
[2788]	1051	INTEGER,SAVE :: ncid_climoz ! NetCDF file identifier
	1052	REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
	1053	REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
	1054	REAL, POINTER, SAVE :: time_climoz(:) ! Time vector
	1055	CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
	1056	= ["tro3 ","tro3_daylight"]
	1057	! vars_climoz(1:read_climoz): variables names in climoz file.
[2819]	1058	! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
	1059	REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
	1060	! the ozone fields, old method.
[1279]	1061
[2469]	1062	include "YOMCST.h"
	1063	include "YOETHF.h"
	1064	include "FCTTRE.h"
	1065	!IM 100106 BEG : pouvoir sortir les ctes de la physique
	1066	include "conema3.h"
	1067	include "fisrtilp.h"
	1068	include "nuage.h"
	1069	include "compbl.h"
	1070	!IM 100106 END : pouvoir sortir les ctes de la physique
	1071	!
	1072	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1073	! Declarations pour Simulateur COSP
	1074	!============================================================
[2799]	1075	#ifdef CPP_COSP
[2469]	1076	real :: mr_ozone(klon,klev)
[2799]	1077	#endif
[2469]	1078	!IM stations CFMIP
	1079	INTEGER, SAVE :: nCFMIP
	1080	!$OMP THREADPRIVATE(nCFMIP)
	1081	INTEGER, PARAMETER :: npCFMIP=120
	1082	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
	1083	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
	1084	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
	1085	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
	1086	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
	1087	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
	1088	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
	1089	!$OMP THREADPRIVATE(iGCM, jGCM)
	1090	logical, dimension(nfiles) :: phys_out_filestations
	1091	logical, parameter :: lNMC=.FALSE.
[1539]	1092
[2469]	1093	!IM betaCRF
	1094	REAL, SAVE :: pfree, beta_pbl, beta_free
	1095	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
	1096	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
	1097	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
	1098	LOGICAL, SAVE :: mskocean_beta
	1099	!$OMP THREADPRIVATE(mskocean_beta)
	1100	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
	1101	! cldemirad pour evaluer les
	1102	! retros liees aux CRF
	1103	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
	1104	! pour radlwsw pour
	1105	! tester "CRF off"
	1106	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
	1107	! pour radlwsw pour
	1108	! tester "CRF off"
	1109	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
	1110	! radlwsw pour tester
	1111	! "CRF off"
	1112	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
[1735]	1113
[2469]	1114	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
	1115	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
[2784]	1116	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
[2469]	1117	integer iostat
[1539]	1118
[2469]	1119	REAL zzz
	1120	!albedo SB >>>
	1121	real,dimension(6),save :: SFRWL
[3065]	1122	!$OMP THREADPRIVATE(SFRWL)
[2469]	1123	!albedo SB <<<
[1955]	1124
[2485]	1125	!--OB variables for mass fixer (hard coded for now)
[2477]	1126	logical, parameter :: mass_fixer=.false.
[2799]	1127	real qql1(klon),qql2(klon),corrqql
[2476]	1128
[2469]	1129	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
	1130	jjmp1=nbp_lat
[2344]	1131
[2469]	1132	!======================================================================
	1133	! Gestion calendrier : mise a jour du module phys_cal_mod
	1134	!
	1135	pdtphys=pdtphys_
	1136	CALL update_time(pdtphys)
[3320]	1137	phys_tstep=NINT(pdtphys)
[3406]	1138	#ifdef CPP_XIOS
	1139	IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
	1140	#endif
[1355]	1141
[2469]	1142	!======================================================================
	1143	! Ecriture eventuelle d'un profil verticale en entree de la physique.
	1144	! Utilise notamment en 1D mais peut etre active egalement en 3D
	1145	! en imposant la valeur de igout.
	1146	!======================================================================d
[2692]	1147	IF (prt_level.ge.1) THEN
[2469]	1148	igout=klon/2+1/klon
	1149	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1150	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
	1151	longitude_deg(igout)
	1152	write(lunout,*) &
	1153	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
	1154	write(lunout,*) &
	1155	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
[879]	1156
[2469]	1157	write(lunout,*) 'paprs, play, phi, u, v, t'
[2692]	1158	DO k=1,klev
[2469]	1159	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
	1160	u(igout,k),v(igout,k),t(igout,k)
[2692]	1161	ENDDO
[2469]	1162	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
[2692]	1163	DO k=1,klev
[2469]	1164	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
[2692]	1165	ENDDO
	1166	ENDIF
[879]	1167
[2769]	1168	! Quick check on pressure levels:
[2774]	1169	call assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
[2769]	1170	"physiq_mod paprs bad order")
[879]	1171
[2692]	1172	IF (first) THEN
[3406]	1173	CALL init_etat0_limit_unstruct
	1174	IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
[2469]	1175	!CR:nvelles variables convection/poches froides
	1176
	1177	print*, '================================================='
	1178	print*, 'Allocation des variables locales et sauvegardees'
[2692]	1179	CALL phys_local_var_init
[2469]	1180	!
	1181	! appel a la lecture du run.def physique
[2692]	1182	CALL conf_phys(ok_journe, ok_mensuel, &
[2469]	1183	ok_instan, ok_hf, &
	1184	ok_LES, &
	1185	callstats, &
	1186	solarlong0,seuil_inversion, &
	1187	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
	1188	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
[2738]	1189	ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, &
[2469]	1190	flag_aerosol, flag_aerosol_strat, new_aod, &
[2644]	1191	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
[2469]	1192	! nv flags pour la convection et les
	1193	! poches froides
	1194	read_climoz, &
	1195	alp_offset)
[2692]	1196	CALL phys_state_var_init(read_climoz)
	1197	CALL phys_output_var_init
[3406]	1198	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
	1199
[2469]	1200	print*, '================================================='
	1201	!
	1202	!CR: check sur le nb de traceurs de l eau
[2692]	1203	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
[2469]	1204	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
	1205	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
[2224]	1206	STOP
[2692]	1207	ENDIF
[2224]	1208
[2469]	1209	dnwd0=0.0
	1210	ftd=0.0
	1211	fqd=0.0
	1212	cin=0.
	1213	!ym Attention pbase pas initialise dans concvl !!!!
	1214	pbase=0
	1215	!IM 180608
[904]	1216
[2999]	1217
[2469]	1218	itau_con=0
	1219	first=.false.
[1797]	1220
[2692]	1221	ENDIF ! first
[1797]	1222
[2469]	1223	!ym => necessaire pour iflag_con != 2
	1224	pmfd(:,:) = 0.
	1225	pen_u(:,:) = 0.
	1226	pen_d(:,:) = 0.
	1227	pde_d(:,:) = 0.
	1228	pde_u(:,:) = 0.
	1229	aam=0.
	1230	d_t_adjwk(:,:)=0
	1231	d_q_adjwk(:,:)=0
[1797]	1232
[2469]	1233	alp_bl_conv(:)=0.
[2245]	1234
[2469]	1235	torsfc=0.
	1236	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
[1797]	1237
[2469]	1238	modname = 'physiq'
[644]	1239
[2469]	1240	IF (debut) THEN
	1241	CALL suphel ! initialiser constantes et parametres phys.
[3000]	1242	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
[2469]	1243	CALL getin_p('random_notrig_max',random_notrig_max)
[2882]	1244	CALL getin_p('ok_adjwk',ok_adjwk)
	1245	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
	1246	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
	1247	! 1 => convective adjustment but state variables are unchanged
	1248	! 2 => convective adjustment and state variables are changed
	1249	CALL getin_p('iflag_adjwk',iflag_adjwk)
[2613]	1250	CALL getin_p('oliqmax',oliqmax)
[2657]	1251	CALL getin_p('oicemax',oicemax)
[2534]	1252	CALL getin_p('ratqsp0',ratqsp0)
	1253	CALL getin_p('ratqsdp',ratqsdp)
[2635]	1254	iflag_wake_tend = 0
	1255	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
[2799]	1256	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
	1257	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
	1258	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
	1259	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
	1260	CALL getin_p('fl_ebil',fl_ebil)
	1261	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
	1262	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
[2984]	1263	iflag_phytrac = 1 ! by default we do want to call phytrac
[2973]	1264	CALL getin_p('iflag_phytrac',iflag_phytrac)
[3011]	1265	nvm_lmdz = 13
	1266	CALL getin_p('NVM',nvm_lmdz)
[2469]	1267	ENDIF
[878]	1268
[2692]	1269	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
[1279]	1270
[959]	1271
[2469]	1272	!======================================================================
	1273	! Gestion calendrier : mise a jour du module phys_cal_mod
	1274	!
	1275	! CALL phys_cal_update(jD_cur,jH_cur)
[1279]	1276
[2469]	1277	!
	1278	! Si c'est le debut, il faut initialiser plusieurs choses
	1279	! ********
	1280	!
	1281	IF (debut) THEN
	1282	!rv CRinitialisation de wght_th et lalim_conv pour la
	1283	!definition de la couche alimentation de la convection a partir
	1284	!des caracteristiques du thermique
	1285	wght_th(:,:)=1.
	1286	lalim_conv(:)=1
	1287	!RC
	1288	ustar(:,:)=0.
[2569]	1289	! u10m(:,:)=0.
	1290	! v10m(:,:)=0.
[2469]	1291	rain_con(:)=0.
	1292	snow_con(:)=0.
	1293	topswai(:)=0.
	1294	topswad(:)=0.
	1295	solswai(:)=0.
	1296	solswad(:)=0.
[959]	1297
[2469]	1298	wmax_th(:)=0.
	1299	tau_overturning_th(:)=0.
[645]	1300
[2469]	1301	IF (type_trac == 'inca') THEN
	1302	! jg : initialisation jusqu'au ces variables sont dans restart
	1303	ccm(:,:,:) = 0.
	1304	tau_aero(:,:,:,:) = 0.
	1305	piz_aero(:,:,:,:) = 0.
	1306	cg_aero(:,:,:,:) = 0.
[2372]	1307
[2469]	1308	config_inca='none' ! default
	1309	CALL getin_p('config_inca',config_inca)
[2372]	1310
[2469]	1311	ELSE
	1312	config_inca='none' ! default
[2692]	1313	ENDIF
[782]	1314
[2469]	1315	IF (aerosol_couple .AND. (config_inca /= "aero" &
	1316	.AND. config_inca /= "aeNP ")) THEN
	1317	abort_message &
	1318	= 'if aerosol_couple is activated, config_inca need to be ' &
	1319	// 'aero or aeNP'
	1320	CALL abort_physic (modname,abort_message,1)
	1321	ENDIF
[2443]	1322
	1323
[1863]	1324
[2469]	1325	rnebcon0(:,:) = 0.0
	1326	clwcon0(:,:) = 0.0
	1327	rnebcon(:,:) = 0.0
	1328	clwcon(:,:) = 0.0
[1863]	1329
[2469]	1330	!
	1331	print*,'iflag_coupl,iflag_clos,iflag_wake', &
	1332	iflag_coupl,iflag_clos,iflag_wake
	1333	print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
	1334	!
	1335	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
	1336	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
	1337	CALL abort_physic (modname,abort_message,1)
	1338	ENDIF
	1339	!
	1340	!
	1341	! Initialiser les compteurs:
	1342	!
	1343	itap = 0
	1344	itaprad = 0
[2707]	1345	itapcv = 0
[2730]	1346	itapwk = 0
[878]	1347
[2469]	1348	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1349	!! Un petit travail \`a faire ici.
	1350	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[878]	1351
[2692]	1352	IF (iflag_pbl>1) THEN
[2469]	1353	PRINT*, "Using method MELLOR&YAMADA"
[2692]	1354	ENDIF
[956]	1355
[2469]	1356	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1357	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
	1358	! phylmd plutot que dyn3d
	1359	! Attention : la version precedente n'etait pas tres propre.
	1360	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
	1361	! pour obtenir le meme resultat.
[2731]	1362	!jyg for fh<
[3065]	1363	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
	1364	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
[2731]	1365	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
	1366	CALL abort_physic(modname,abort_message,1)
	1367	ENDIF
	1368	!>jyg
[3065]	1369	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
	1370	radpas = NINT( 86400./phys_tstep)/nbapp_rad
[2469]	1371	ELSE
	1372	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
	1373	'multiple de nbapp_rad'
	1374	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
	1375	'mais 1+1<>2'
	1376	abort_message='nbre de pas de temps physique n est pas multiple ' &
	1377	// 'de nbapp_rad'
[2692]	1378	CALL abort_physic(modname,abort_message,1)
[2469]	1379	ENDIF
[3065]	1380	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
	1381	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
[2730]	1382	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
[3065]	1383	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
	1384	cvpas = NINT( 86400./phys_tstep)/nbapp_cv
[2707]	1385	print *,'physiq, cvpas ',cvpas
	1386	ELSE
	1387	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
	1388	'multiple de nbapp_cv'
	1389	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
	1390	'mais 1+1<>2'
	1391	abort_message='nbre de pas de temps physique n est pas multiple ' &
	1392	// 'de nbapp_cv'
	1393	call abort_physic(modname,abort_message,1)
	1394	ENDIF
[3065]	1395	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
	1396	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
[2730]	1397	print *,'physiq, wkpas ',wkpas
	1398	ELSE
	1399	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
	1400	'multiple de nbapp_wk'
	1401	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
	1402	'mais 1+1<>2'
	1403	abort_message='nbre de pas de temps physique n est pas multiple ' &
	1404	// 'de nbapp_wk'
	1405	call abort_physic(modname,abort_message,1)
	1406	ENDIF
[2469]	1407	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3322]	1408	CALL init_iophy_new(latitude_deg,longitude_deg)
[524]	1409
[3406]	1410	!===================================================================
	1411	!IM stations CFMIP
	1412	nCFMIP=npCFMIP
	1413	OPEN(98,file='npCFMIP_param.data',status='old', &
	1414	form='formatted',iostat=iostat)
	1415	IF (iostat == 0) THEN
	1416	READ(98,*,end=998) nCFMIP
	1417	998 CONTINUE
	1418	CLOSE(98)
	1419	CONTINUE
	1420	IF(nCFMIP.GT.npCFMIP) THEN
	1421	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
	1422	CALL abort_physic("physiq", "", 1)
	1423	ELSE
	1424	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
	1425	ENDIF
[3322]	1426
[3406]	1427	!
	1428	ALLOCATE(tabCFMIP(nCFMIP))
	1429	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
	1430	ALLOCATE(tabijGCM(nCFMIP))
	1431	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
	1432	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
	1433	!
	1434	! lecture des nCFMIP stations CFMIP, de leur numero
	1435	! et des coordonnees geographiques lonCFMIP, latCFMIP
	1436	!
	1437	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
	1438	lonCFMIP, latCFMIP)
	1439	!
	1440	! identification des
	1441	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
	1442	! grille de LMDZ
	1443	! 2) indices points tabijGCM de la grille physique 1d sur
	1444	! klon points
	1445	! 3) indices iGCM, jGCM de la grille physique 2d
	1446	!
	1447	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
	1448	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
	1449	!
	1450	ELSE
	1451	ALLOCATE(tabijGCM(0))
	1452	ALLOCATE(lonGCM(0), latGCM(0))
	1453	ALLOCATE(iGCM(0), jGCM(0))
	1454	ENDIF
	1455
	1456	#ifdef CPP_IOIPSL
	1457
	1458	!$OMP MASTER
	1459	! FH : if ok_sync=.true. , the time axis is written at each time step
	1460	! in the output files. Only at the end in the opposite case
	1461	ok_sync_omp=.false.
	1462	CALL getin('ok_sync',ok_sync_omp)
	1463	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
	1464	iGCM,jGCM,lonGCM,latGCM, &
	1465	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
	1466	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
	1467	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
	1468	read_climoz, phys_out_filestations, &
	1469	new_aod, aerosol_couple, &
	1470	flag_aerosol_strat, pdtphys, paprs, pphis, &
	1471	pplay, lmax_th, ptconv, ptconvth, ivap, &
	1472	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
	1473	!$OMP END MASTER
	1474	!$OMP BARRIER
	1475	ok_sync=ok_sync_omp
	1476
	1477	freq_outNMC(1) = ecrit_files(7)
	1478	freq_outNMC(2) = ecrit_files(8)
	1479	freq_outNMC(3) = ecrit_files(9)
	1480	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
	1481	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
	1482	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
	1483
	1484	#ifndef CPP_XIOS
	1485	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
	1486	#endif
	1487
	1488	#endif
	1489	ecrit_reg = ecrit_reg * un_jour
	1490	ecrit_tra = ecrit_tra * un_jour
	1491
	1492	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
	1493	date0 = jD_ref
	1494	WRITE(,) 'physiq date0 : ',date0
	1495	!
	1496
	1497	! CALL create_climoz(read_climoz)
	1498
	1499	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
	1500	pplay, lmax_th, aerosol_couple, &
	1501	ok_ade, ok_aie, ivap, iliq, isol, new_aod, ok_sync,&
	1502	ptconv, read_climoz, clevSTD, &
	1503	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
	1504	flag_aerosol, flag_aerosol_strat, ok_cdnc)
	1505
	1506	#ifdef CPP_XIOS
	1507	IF (is_omp_master) CALL xios_update_calendar(1)
	1508	#endif
	1509	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
	1510	CALL create_etat0_limit_unstruct
[2469]	1511	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
[3406]	1512
[2565]	1513	!jyg<
[2469]	1514	IF (klon_glo==1) THEN
[2565]	1515	pbl_tke(:,:,is_ave) = 0.
	1516	DO nsrf=1,nbsrf
	1517	DO k = 1,klev+1
	1518	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
	1519	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
	1520	ENDDO
	1521	ENDDO
[3406]	1522	ELSE
	1523	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
[2565]	1524	!>jyg
[2469]	1525	ENDIF
	1526	!IM begin
	1527	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
	1528	,ratqs(1,1)
	1529	!IM end
[878]	1530
	1531
[2469]	1532	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1533	!
	1534	! on remet le calendrier a zero
	1535	!
	1536	IF (raz_date .eq. 1) THEN
	1537	itau_phy = 0
	1538	ENDIF
[524]	1539
[3065]	1540	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
	1541	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
	1542	! pdtphys
	1543	! abort_message='Pas physique n est pas correct '
	1544	! ! call abort_physic(modname,abort_message,1)
	1545	! phys_tstep=pdtphys
	1546	! ENDIF
[2469]	1547	IF (nlon .NE. klon) THEN
	1548	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
	1549	klon
	1550	abort_message='nlon et klon ne sont pas coherents'
[2692]	1551	CALL abort_physic(modname,abort_message,1)
[2469]	1552	ENDIF
	1553	IF (nlev .NE. klev) THEN
	1554	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
	1555	klev
	1556	abort_message='nlev et klev ne sont pas coherents'
[2692]	1557	CALL abort_physic(modname,abort_message,1)
[2469]	1558	ENDIF
	1559	!
[3065]	1560	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
[2469]	1561	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1562	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1563	abort_message='Nbre d appels au rayonnement insuffisant'
[2692]	1564	CALL abort_physic(modname,abort_message,1)
[2469]	1565	ENDIF
	1566	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1567	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
	1568	ok_cvl
	1569	!
	1570	!KE43
	1571	! Initialisation pour la convection de K.E. (sb):
	1572	IF (iflag_con.GE.3) THEN
[524]	1573
[2469]	1574	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
	1575	WRITE(lunout,*) &
	1576	"On va utiliser le melange convectif des traceurs qui"
	1577	WRITE(lunout,*)"est calcule dans convect4.3"
	1578	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]	1579
[2469]	1580	DO i = 1, klon
	1581	ema_cbmf(i) = 0.
	1582	ema_pcb(i) = 0.
	1583	ema_pct(i) = 0.
	1584	! ema_workcbmf(i) = 0.
	1585	ENDDO
	1586	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1587	DO i = 1, klon
	1588	ibas_con(i) = 1
	1589	itop_con(i) = 1
	1590	ENDDO
	1591	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	1592	!================================================================
	1593	!CR:04.12.07: initialisations poches froides
	1594	! Controle de ALE et ALP pour la fermeture convective (jyg)
[2692]	1595	IF (iflag_wake>=1) THEN
[2469]	1596	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
	1597	,alp_bl_prescr, ale_bl_prescr)
	1598	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1599	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
[2638]	1600	!
	1601	! Initialize tendencies of wake state variables (for some flag values
	1602	! they are not computed).
	1603	d_deltat_wk(:,:) = 0.
	1604	d_deltaq_wk(:,:) = 0.
	1605	d_deltat_wk_gw(:,:) = 0.
	1606	d_deltaq_wk_gw(:,:) = 0.
	1607	d_deltat_vdf(:,:) = 0.
	1608	d_deltaq_vdf(:,:) = 0.
	1609	d_deltat_the(:,:) = 0.
	1610	d_deltaq_the(:,:) = 0.
	1611	d_deltat_ajs_cv(:,:) = 0.
	1612	d_deltaq_ajs_cv(:,:) = 0.
	1613	d_s_wk(:) = 0.
	1614	d_dens_wk(:) = 0.
[2692]	1615	ENDIF
[973]	1616
[2469]	1617	! do i = 1,klon
	1618	! Ale_bl(i)=0.
	1619	! Alp_bl(i)=0.
	1620	! enddo
[1638]	1621
[3406]	1622	!ELSE
	1623	! ALLOCATE(tabijGCM(0))
	1624	! ALLOCATE(lonGCM(0), latGCM(0))
	1625	! ALLOCATE(iGCM(0), jGCM(0))
[2469]	1626	ENDIF
[878]	1627
[2469]	1628	DO i=1,klon
	1629	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1630	ENDDO
[1863]	1631
[2469]	1632	!34EK
	1633	IF (ok_orodr) THEN
[524]	1634
[2469]	1635	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1636	! FH sans doute a enlever de finitivement ou, si on le
	1637	! garde, l'activer justement quand ok_orodr = false.
	1638	! ce rugoro est utilise par la couche limite et fait double emploi
	1639	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
	1640	! DO i=1,klon
	1641	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1642	! ENDDO
	1643	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1644	IF (ok_strato) THEN
	1645	CALL SUGWD_strato(klon,klev,paprs,pplay)
	1646	ELSE
	1647	CALL SUGWD(klon,klev,paprs,pplay)
	1648	ENDIF
[1863]	1649
[2469]	1650	DO i=1,klon
	1651	zuthe(i)=0.
	1652	zvthe(i)=0.
[2692]	1653	IF (zstd(i).gt.10.) THEN
[2469]	1654	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1655	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
[2692]	1656	ENDIF
[2469]	1657	ENDDO
	1658	ENDIF
	1659	!
	1660	!
[3065]	1661	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
[2469]	1662	WRITE(lunout,*)'La frequence de lecture surface est de ', &
	1663	lmt_pas
	1664	!
	1665	capemaxcels = 't_max(X)'
	1666	t2mincels = 't_min(X)'
	1667	t2maxcels = 't_max(X)'
	1668	tinst = 'inst(X)'
	1669	tave = 'ave(X)'
	1670	!IM cf. AM 081204 BEG
	1671	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1672	!IM cf. AM 081204 END
	1673	!
	1674	!=============================================================
	1675	! Initialisation des sorties
	1676	!=============================================================
	1677
[3322]	1678	#ifdef CPP_XIOS
	1679	! Get "missing_val" value from XML files (from temperature variable)
	1680	!$OMP MASTER
	1681	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
	1682	!$OMP END MASTER
	1683	!$OMP BARRIER
	1684	missing_val=missing_val_omp
	1685	#endif
	1686
[2679]	1687	#ifdef CPP_XIOS
	1688	!--setting up swaero_diag to TRUE in XIOS case
	1689	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
	1690	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
	1691	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
	1692	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
	1693	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
	1694	!!!--for now these fields are not in the XML files so they are omitted
	1695	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
	1696	swaero_diag=.TRUE.
[2854]	1697
	1698	!--setting up dryaod_diag to TRUE in XIOS case
	1699	DO naero = 1, naero_tot-1
[2863]	1700	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
[2854]	1701	ENDDO
[2989]	1702	!
	1703	!--setting up ok_4xCO2atm to TRUE in XIOS case
	1704	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
	1705	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
	1706	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
	1707	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
	1708	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
	1709	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
	1710	ok_4xCO2atm=.TRUE.
[2679]	1711	#endif
	1712
[524]	1713
[3406]	1714	CALL printflag( tabcntr0,radpas,ok_journe, &
	1715	ok_instan, ok_region )
[2469]	1716	!
	1717	!
	1718	!
	1719	! Prescrire l'ozone dans l'atmosphere
	1720	!
	1721	!
	1722	!c DO i = 1, klon
	1723	!c DO k = 1, klev
	1724	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1725	!c ENDDO
	1726	!c ENDDO
	1727	!
	1728	IF (type_trac == 'inca') THEN
[524]	1729	#ifdef INCA
[2469]	1730	CALL VTe(VTphysiq)
	1731	CALL VTb(VTinca)
	1732	calday = REAL(days_elapsed) + jH_cur
	1733	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
[959]	1734
[2469]	1735	CALL chemini( &
	1736	rg, &
	1737	ra, &
	1738	cell_area, &
	1739	latitude_deg, &
	1740	longitude_deg, &
	1741	presnivs, &
	1742	calday, &
	1743	klon, &
	1744	nqtot, &
[2566]	1745	nqo, &
[2469]	1746	pdtphys, &
	1747	annee_ref, &
[2906]	1748	year_cur, &
[2469]	1749	day_ref, &
	1750	day_ini, &
	1751	start_time, &
	1752	itau_phy, &
[2831]	1753	date0, &
[2469]	1754	io_lon, &
	1755	io_lat)
[959]	1756
[2469]	1757	CALL VTe(VTinca)
	1758	CALL VTb(VTphysiq)
[524]	1759	#endif
[2692]	1760	ENDIF
[2469]	1761	!
	1762	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1763	! Nouvelle initialisation pour le rayonnement RRTM
	1764	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[998]	1765
[2692]	1766	CALL iniradia(klon,klev,paprs(1,1:klev+1))
[998]	1767
[2469]	1768	!$omp single
[2788]	1769	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
[2820]	1770	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
[2469]	1771	!$omp end single
	1772	!
	1773	!IM betaCRF
	1774	pfree=70000. !Pa
	1775	beta_pbl=1.
	1776	beta_free=1.
	1777	lon1_beta=-180.
	1778	lon2_beta=+180.
	1779	lat1_beta=90.
	1780	lat2_beta=-90.
	1781	mskocean_beta=.FALSE.
[1539]	1782
[2469]	1783	!albedo SB >>>
	1784	select case(nsw)
	1785	case(2)
	1786	SFRWL(1)=0.45538747
	1787	SFRWL(2)=0.54461211
	1788	case(4)
	1789	SFRWL(1)=0.45538747
	1790	SFRWL(2)=0.32870591
	1791	SFRWL(3)=0.18568763
	1792	SFRWL(4)=3.02191470E-02
	1793	case(6)
	1794	SFRWL(1)=1.28432794E-03
	1795	SFRWL(2)=0.12304168
	1796	SFRWL(3)=0.33106142
	1797	SFRWL(4)=0.32870591
	1798	SFRWL(5)=0.18568763
	1799	SFRWL(6)=3.02191470E-02
	1800	end select
[2227]	1801
	1802
[2469]	1803	!albedo SB <<<
[2227]	1804
[2469]	1805	OPEN(99,file='beta_crf.data',status='old', &
	1806	form='formatted',err=9999)
	1807	READ(99,*,end=9998) pfree
	1808	READ(99,*,end=9998) beta_pbl
	1809	READ(99,*,end=9998) beta_free
	1810	READ(99,*,end=9998) lon1_beta
	1811	READ(99,*,end=9998) lon2_beta
	1812	READ(99,*,end=9998) lat1_beta
	1813	READ(99,*,end=9998) lat2_beta
	1814	READ(99,*,end=9998) mskocean_beta
	1815	9998 Continue
	1816	CLOSE(99)
	1817	9999 Continue
	1818	WRITE(,)'pfree=',pfree
	1819	WRITE(,)'beta_pbl=',beta_pbl
	1820	WRITE(,)'beta_free=',beta_free
	1821	WRITE(,)'lon1_beta=',lon1_beta
	1822	WRITE(,)'lon2_beta=',lon2_beta
	1823	WRITE(,)'lat1_beta=',lat1_beta
	1824	WRITE(,)'lat2_beta=',lat2_beta
	1825	WRITE(,)'mskocean_beta=',mskocean_beta
	1826	ENDIF
	1827	!
	1828	! ************** Fin de IF ( debut ) *************
	1829	!
	1830	!
	1831	! Incrementer le compteur de la physique
	1832	!
	1833	itap = itap + 1
[2795]	1834	IF (is_master .OR. prt_level > 9) THEN
[2783]	1835	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
[2795]	1836	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
	1837	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
	1838	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
[2783]	1839	ENDIF
	1840	ENDIF
[2469]	1841	!
	1842	!
	1843	! Update fraction of the sub-surfaces (pctsrf) and
	1844	! initialize, where a new fraction has appeared, all variables depending
	1845	! on the surface fraction.
	1846	!
[3065]	1847	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
[2469]	1848	pctsrf, fevap, z0m, z0h, agesno, &
	1849	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
[996]	1850
[2469]	1851	! Update time and other variables in Reprobus
	1852	IF (type_trac == 'repr') THEN
[1565]	1853	#ifdef REPROBUS
[2469]	1854	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
	1855	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
	1856	CALL Rtime(debut)
[1565]	1857	#endif
[2692]	1858	ENDIF
[1565]	1859
	1860
[2469]	1861	! Tendances bidons pour les processus qui n'affectent pas certaines
	1862	! variables.
	1863	du0(:,:)=0.
	1864	dv0(:,:)=0.
	1865	dt0 = 0.
	1866	dq0(:,:)=0.
	1867	dql0(:,:)=0.
	1868	dqi0(:,:)=0.
[2635]	1869	dsig0(:) = 0.
	1870	ddens0(:) = 0.
	1871	wkoccur1(:)=1
[2469]	1872	!
	1873	! Mettre a zero des variables de sortie (pour securite)
	1874	!
	1875	DO i = 1, klon
	1876	d_ps(i) = 0.0
	1877	ENDDO
	1878	DO k = 1, klev
	1879	DO i = 1, klon
	1880	d_t(i,k) = 0.0
	1881	d_u(i,k) = 0.0
	1882	d_v(i,k) = 0.0
	1883	ENDDO
	1884	ENDDO
	1885	DO iq = 1, nqtot
	1886	DO k = 1, klev
	1887	DO i = 1, klon
	1888	d_qx(i,k,iq) = 0.0
	1889	ENDDO
	1890	ENDDO
	1891	ENDDO
	1892	da(:,:)=0.
	1893	mp(:,:)=0.
	1894	phi(:,:,:)=0.
	1895	! RomP >>>
	1896	phi2(:,:,:)=0.
	1897	beta_prec_fisrt(:,:)=0.
	1898	beta_prec(:,:)=0.
	1899	epmlmMm(:,:,:)=0.
	1900	eplaMm(:,:)=0.
	1901	d1a(:,:)=0.
	1902	dam(:,:)=0.
	1903	pmflxr=0.
	1904	pmflxs=0.
	1905	! RomP <<<
[1742]	1906
[2469]	1907	!
	1908	! Ne pas affecter les valeurs entrees de u, v, h, et q
	1909	!
	1910	DO k = 1, klev
	1911	DO i = 1, klon
	1912	t_seri(i,k) = t(i,k)
	1913	u_seri(i,k) = u(i,k)
	1914	v_seri(i,k) = v(i,k)
	1915	q_seri(i,k) = qx(i,k,ivap)
	1916	ql_seri(i,k) = qx(i,k,iliq)
	1917	!CR: ATTENTION, on rajoute la variable glace
[2692]	1918	IF (nqo.eq.2) THEN
[2469]	1919	qs_seri(i,k) = 0.
[2692]	1920	ELSE IF (nqo.eq.3) THEN
[2469]	1921	qs_seri(i,k) = qx(i,k,isol)
[2692]	1922	ENDIF
[2469]	1923	ENDDO
	1924	ENDDO
[2476]	1925	!
	1926	!--OB mass fixer
	1927	IF (mass_fixer) THEN
	1928	!--store initial water burden
	1929	qql1(:)=0.0
[2499]	1930	DO k = 1, klev
	1931	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2476]	1932	ENDDO
	1933	ENDIF
	1934	!--fin mass fixer
	1935
[2469]	1936	tke0(:,:)=pbl_tke(:,:,is_ave)
	1937	!CR:Nombre de traceurs de l'eau: nqo
	1938	! IF (nqtot.GE.3) THEN
	1939	IF (nqtot.GE.(nqo+1)) THEN
	1940	! DO iq = 3, nqtot
	1941	DO iq = nqo+1, nqtot
	1942	DO k = 1, klev
	1943	DO i = 1, klon
	1944	! tr_seri(i,k,iq-2) = qx(i,k,iq)
	1945	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
	1946	ENDDO
	1947	ENDDO
	1948	ENDDO
	1949	ELSE
	1950	DO k = 1, klev
	1951	DO i = 1, klon
	1952	tr_seri(i,k,1) = 0.0
	1953	ENDDO
	1954	ENDDO
	1955	ENDIF
	1956	!
	1957	DO i = 1, klon
	1958	ztsol(i) = 0.
	1959	ENDDO
	1960	DO nsrf = 1, nbsrf
	1961	DO i = 1, klon
	1962	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1963	ENDDO
	1964	ENDDO
[2611]	1965	! Initialize variables used for diagnostic purpose
[2692]	1966	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
[524]	1967
[2469]	1968	! Diagnostiquer la tendance dynamique
	1969	!
	1970	IF (ancien_ok) THEN
[2499]	1971	!
[3065]	1972	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
	1973	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
	1974	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
	1975	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
	1976	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
	1977	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
[2499]	1978	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
[3065]	1979	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
[2499]	1980	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
[3065]	1981	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
[2499]	1982	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
[3065]	1983	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
[2469]	1984	! !! RomP >>> td dyn traceur
[2499]	1985	IF (nqtot.GT.nqo) THEN ! jyg
[2469]	1986	DO iq = nqo+1, nqtot ! jyg
[3065]	1987	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
[2469]	1988	ENDDO
	1989	ENDIF
	1990	! !! RomP <<<
	1991	ELSE
[2499]	1992	d_u_dyn(:,:) = 0.0
	1993	d_v_dyn(:,:) = 0.0
	1994	d_t_dyn(:,:) = 0.0
	1995	d_q_dyn(:,:) = 0.0
	1996	d_ql_dyn(:,:) = 0.0
	1997	d_qs_dyn(:,:) = 0.0
	1998	d_q_dyn2d(:) = 0.0
	1999	d_ql_dyn2d(:) = 0.0
	2000	d_qs_dyn2d(:) = 0.0
[2469]	2001	! !! RomP >>> td dyn traceur
[2499]	2002	IF (nqtot.GT.nqo) THEN ! jyg
	2003	DO iq = nqo+1, nqtot ! jyg
	2004	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
[2469]	2005	ENDDO
	2006	ENDIF
	2007	! !! RomP <<<
	2008	ancien_ok = .TRUE.
	2009	ENDIF
	2010	!
	2011	! Ajouter le geopotentiel du sol:
	2012	!
	2013	DO k = 1, klev
	2014	DO i = 1, klon
	2015	zphi(i,k) = pphi(i,k) + pphis(i)
	2016	ENDDO
	2017	ENDDO
	2018	!
	2019	! Verifier les temperatures
	2020	!
	2021	!IM BEG
	2022	IF (check) THEN
	2023	amn=MIN(ftsol(1,is_ter),1000.)
	2024	amx=MAX(ftsol(1,is_ter),-1000.)
	2025	DO i=2, klon
	2026	amn=MIN(ftsol(i,is_ter),amn)
	2027	amx=MAX(ftsol(i,is_ter),amx)
	2028	ENDDO
	2029	!
	2030	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	2031	ENDIF !(check) THEN
	2032	!IM END
	2033	!
	2034	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
	2035	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
[2235]	2036
[2469]	2037	!
	2038	!IM BEG
	2039	IF (check) THEN
	2040	amn=MIN(ftsol(1,is_ter),1000.)
	2041	amx=MAX(ftsol(1,is_ter),-1000.)
	2042	DO i=2, klon
	2043	amn=MIN(ftsol(i,is_ter),amn)
	2044	amx=MAX(ftsol(i,is_ter),amx)
	2045	ENDDO
	2046	!
	2047	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	2048	ENDIF !(check) THEN
	2049	!IM END
	2050	!
	2051	! Mettre en action les conditions aux limites (albedo, sst, etc.).
	2052	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
	2053	!
[2661]	2054	! Update ozone if day change
	2055	IF (MOD(itap-1,lmt_pas) == 0) THEN
[2774]	2056	IF (read_climoz <= 0) THEN
	2057	! Once per day, update ozone from Royer:
	2058	IF (solarlong0<-999.) then
	2059	! Generic case with evolvoing season
	2060	zzz=real(days_elapsed+1)
	2061	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
	2062	! Particular case with annual mean insolation
	2063	zzz=real(90) ! could be revisited
	2064	IF (read_climoz/=-1) THEN
	2065	abort_message ='read_climoz=-1 is recommended when ' &
	2066	// 'solarlong0=1000.'
	2067	CALL abort_physic (modname,abort_message,1)
	2068	ENDIF
	2069	ELSE
	2070	! Case where the season is imposed with solarlong0
	2071	zzz=real(90) ! could be revisited
	2072	ENDIF
[2661]	2073
[2774]	2074	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
	2075	ELSE
[2820]	2076	!--- ro3i = elapsed days number since current year 1st january, 0h
	2077	ro3i=days_elapsed+jh_cur-jh_1jan
	2078	!--- scaling for old style files (360 records)
	2079	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
[2788]	2080	IF(adjust_tropopause) THEN
[2820]	2081	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
	2082	ro3i, press_edg_climoz, paprs, wo, time_climoz, &
[2968]	2083	longitude_deg, latitude_deg, press_cen_climoz, &
	2084	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
[2774]	2085	ELSE
[2820]	2086	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
	2087	ro3i, press_edg_climoz, paprs, wo, time_climoz)
[2788]	2088	END IF
[2774]	2089	! Convert from mole fraction of ozone to column density of ozone in a
	2090	! cell, in kDU:
	2091	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
	2092	* zmasse / dobson_u / 1e3
[2788]	2093	! (By regridding ozone values for LMDZ only once a day, we
[2774]	2094	! have already neglected the variation of pressure in one
	2095	! day. So do not recompute "wo" at each time step even if
	2096	! "zmasse" changes a little.)
	2097	ENDIF
[2469]	2098	ENDIF
	2099	!
	2100	! Re-evaporer l'eau liquide nuageuse
	2101	!
[2705]	2102	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
	2103	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
[1849]	2104
[2705]	2105	CALL add_phys_tend &
	2106	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
[2812]	2107	'eva',abortphy,flag_inhib_tend,itap,0)
[2799]	2108	call prt_enerbil('eva',itap)
[2086]	2109
[2469]	2110	!=========================================================================
	2111	! Calculs de l'orbite.
	2112	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	2113	! doit donc etre plac\'e avant radlwsw et pbl_surface
[883]	2114
[2469]	2115	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2692]	2116	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
[2469]	2117	day_since_equinox = (jD_cur + jH_cur) - jD_eq
	2118	!
	2119	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	2120	! solarlong0
[2692]	2121	IF (solarlong0<-999.) THEN
	2122	IF (new_orbit) THEN
[2469]	2123	! calcul selon la routine utilisee pour les planetes
[2692]	2124	CALL solarlong(day_since_equinox, zlongi, dist)
	2125	ELSE
[2469]	2126	! calcul selon la routine utilisee pour l'AR4
	2127	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
[2692]	2128	ENDIF
	2129	ELSE
[2469]	2130	zlongi=solarlong0 ! longitude solaire vraie
	2131	dist=1. ! distance au soleil / moyenne
[2692]	2132	ENDIF
[1529]	2133
[2692]	2134	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
[524]	2135
[2692]	2136
[2469]	2137	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2138	! Calcul de l'ensoleillement :
	2139	! ============================
	2140	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
	2141	! l'annee a partir d'une formule analytique.
	2142	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
	2143	! non nul aux poles.
[2692]	2144	IF (abs(solarlong0-1000.)<1.e-4) THEN
	2145	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
[2469]	2146	latitude_deg,longitude_deg,rmu0,fract)
[2979]	2147	swradcorr(:) = 1.0
	2148	JrNt(:) = 1.0
	2149	zrmu0(:) = rmu0(:)
[2469]	2150	ELSE
	2151	! recode par Olivier Boucher en sept 2015
	2152	SELECT CASE (iflag_cycle_diurne)
	2153	CASE(0)
	2154	! Sans cycle diurne
	2155	CALL angle(zlongi, latitude_deg, fract, rmu0)
	2156	swradcorr = 1.0
	2157	JrNt = 1.0
	2158	zrmu0 = rmu0
	2159	CASE(1)
	2160	! Avec cycle diurne sans application des poids
	2161	! bit comparable a l ancienne formulation cycle_diurne=true
	2162	! on integre entre gmtime et gmtime+radpas
[3065]	2163	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
[2469]	2164	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
	2165	latitude_deg,longitude_deg,rmu0,fract)
	2166	zrmu0 = rmu0
	2167	swradcorr = 1.0
	2168	! Calcul du flag jour-nuit
	2169	JrNt = 0.0
	2170	WHERE (fract.GT.0.0) JrNt = 1.0
	2171	CASE(2)
	2172	! Avec cycle diurne sans application des poids
	2173	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
	2174	! Comme cette routine est appele a tous les pas de temps de
	2175	! la physique meme si le rayonnement n'est pas appele je
	2176	! remonte en arriere les radpas-1 pas de temps
	2177	! suivant. Petite ruse avec MOD pour prendre en compte le
	2178	! premier pas de temps de la physique pendant lequel
	2179	! itaprad=0
[3065]	2180	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
	2181	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
[2469]	2182	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
	2183	latitude_deg,longitude_deg,rmu0,fract)
	2184	!
	2185	! Calcul des poids
	2186	!
[3065]	2187	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
[2469]	2188	zdtime2=0.0 !--pas de temps de la physique qui se termine
	2189	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
	2190	latitude_deg,longitude_deg,zrmu0,zfract)
	2191	swradcorr = 0.0
	2192	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
	2193	swradcorr=zfract/fract*zrmu0/rmu0
	2194	! Calcul du flag jour-nuit
	2195	JrNt = 0.0
	2196	WHERE (zfract.GT.0.0) JrNt = 1.0
	2197	END SELECT
	2198	ENDIF
[782]	2199
[2692]	2200	IF (mydebug) THEN
	2201	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	2202	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	2203	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	2204	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	2205	ENDIF
[883]	2206
[2469]	2207	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2208	! Appel au pbl_surface : Planetary Boudary Layer et Surface
	2209	! Cela implique tous les interactions des sous-surfaces et la
	2210	! partie diffusion turbulent du couche limit.
	2211	!
	2212	! Certains varibales de sorties de pbl_surface sont utiliser que pour
	2213	! ecriture des fihiers hist_XXXX.nc, ces sont :
	2214	! qsol, zq2m, s_pblh, s_lcl,
	2215	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	2216	! s_therm, s_trmb1, s_trmb2, s_trmb3,
	2217	! zu10m, zv10m, fder,
	2218	! zxqsurf, rh2m, zxfluxu, zxfluxv,
	2219	! frugs, agesno, fsollw, fsolsw,
	2220	! d_ts, fevap, fluxlat, t2m,
	2221	! wfbils, wfbilo, fluxt, fluxu, fluxv,
	2222	!
	2223	! Certains ne sont pas utiliser du tout :
	2224	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
	2225	!
[1724]	2226
[2469]	2227	! Calcul de l'humidite de saturation au niveau du sol
[1724]	2228
	2229
[996]	2230
[2692]	2231	IF (iflag_pbl/=0) THEN
[2240]	2232
[2469]	2233	!jyg+nrlmd<
[2852]	2234	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
	2235	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
[2469]	2236	print *,'debut du splitting de la PBL'
	2237	ENDIF
	2238	! !!
	2239	!>jyg+nrlmd
	2240	!
	2241	!-------gustiness calculation-------!
[3406]	2242	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
	2243	gustiness=0 !ym missing init
	2244
[2469]	2245	IF (iflag_gusts==0) THEN
	2246	gustiness(1:klon)=0
	2247	ELSE IF (iflag_gusts==1) THEN
	2248	do i = 1, klon
	2249	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
	2250	enddo
	2251	! ELSE IF (iflag_gusts==2) THEN
	2252	! do i = 1, klon
	2253	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
	2254	! *ale_wake(i) !! need to make sigma_wk accessible here
	2255	! enddo
	2256	! ELSE IF (iflag_gusts==3) THEN
	2257	! do i = 1, klon
	2258	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
	2259	! enddo
	2260	ENDIF
[2278]	2261
	2262
[1067]	2263
[2469]	2264	CALL pbl_surface( &
[3065]	2265	phys_tstep, date0, itap, days_elapsed+1, &
[2469]	2266	debut, lafin, &
	2267	longitude_deg, latitude_deg, rugoro, zrmu0, &
	2268	zsig, sollwdown, pphi, cldt, &
	2269	rain_fall, snow_fall, solsw, sollw, &
	2270	gustiness, &
	2271	t_seri, q_seri, u_seri, v_seri, &
	2272	!nrlmd+jyg<
	2273	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
	2274	!>nrlmd+jyg
	2275	pplay, paprs, pctsrf, &
	2276	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
	2277	!albedo SB <<<
	2278	cdragh, cdragm, u1, v1, &
	2279	!albedo SB >>>
	2280	! albsol1, albsol2, sens, evap, &
	2281	albsol_dir, albsol_dif, sens, evap, &
	2282	!albedo SB <<<
	2283	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
	2284	zxtsol, zxfluxlat, zt2m, qsat2m, &
	2285	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
	2286	!nrlmd<
	2287	!jyg<
	2288	d_t_vdf_w, d_q_vdf_w, &
	2289	d_t_vdf_x, d_q_vdf_x, &
	2290	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
	2291	!>jyg
	2292	delta_tsurf,wake_dens, &
	2293	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
	2294	kh,kh_x,kh_w, &
	2295	!>nrlmd
	2296	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
	2297	slab_wfbils, &
	2298	qsol, zq2m, s_pblh, s_lcl, &
	2299	!jyg<
	2300	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
	2301	!>jyg
	2302	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
	2303	s_therm, s_trmb1, s_trmb2, s_trmb3, &
	2304	zustar, zu10m, zv10m, fder, &
	2305	zxqsurf, rh2m, zxfluxu, zxfluxv, &
	2306	z0m, z0h, agesno, fsollw, fsolsw, &
	2307	d_ts, fevap, fluxlat, t2m, &
[2670]	2308	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
	2309	fluxt, fluxu, fluxv, &
[2469]	2310	dsens, devap, zxsnow, &
	2311	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
	2312	!nrlmd+jyg<
[2952]	2313	wake_delta_pbl_TKE, &
[2469]	2314	!>nrlmd+jyg
[2952]	2315	treedrg )
	2316	!FC
[2469]	2317	!
	2318	! Add turbulent diffusion tendency to the wake difference variables
[2852]	2319	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
	2320	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
[2635]	2321	!jyg<
	2322	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
	2323	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
	2324	CALL add_wake_tend &
	2325	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy)
	2326	ELSE
	2327	d_deltat_vdf(:,:) = 0.
	2328	d_deltaq_vdf(:,:) = 0.
	2329	!>jyg
[2469]	2330	ENDIF
[1624]	2331
[766]	2332
[2897]	2333
	2334
	2335
[2469]	2336	!---------------------------------------------------------------------
	2337	! ajout des tendances de la diffusion turbulente
	2338	IF (klon_glo==1) THEN
	2339	CALL add_pbl_tend &
	2340	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2799]	2341	'vdf',abortphy,flag_inhib_tend,itap)
[2469]	2342	ELSE
	2343	CALL add_phys_tend &
	2344	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2812]	2345	'vdf',abortphy,flag_inhib_tend,itap,0)
[2469]	2346	ENDIF
[2799]	2347	call prt_enerbil('vdf',itap)
[2469]	2348	!--------------------------------------------------------------------
[766]	2349
[2692]	2350	IF (mydebug) THEN
	2351	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	2352	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	2353	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	2354	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	2355	ENDIF
[2227]	2356
[2469]	2357	!albedo SB >>>
	2358	albsol1=0.
	2359	albsol2=0.
	2360	falb1=0.
	2361	falb2=0.
[2692]	2362	SELECT CASE(nsw)
	2363	CASE(2)
[2469]	2364	albsol1=albsol_dir(:,1)
	2365	albsol2=albsol_dir(:,2)
	2366	falb1=falb_dir(:,1,:)
	2367	falb2=falb_dir(:,2,:)
[2692]	2368	CASE(4)
[2469]	2369	albsol1=albsol_dir(:,1)
	2370	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
	2371	+albsol_dir(:,4)*SFRWL(4)
	2372	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
	2373	falb1=falb_dir(:,1,:)
	2374	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
	2375	+falb_dir(:,4,:)*SFRWL(4)
	2376	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
[2692]	2377	CASE(6)
[2469]	2378	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
	2379	+albsol_dir(:,3)*SFRWL(3)
	2380	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
	2381	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
	2382	+albsol_dir(:,6)*SFRWL(6)
	2383	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
	2384	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
	2385	+falb_dir(:,3,:)*SFRWL(3)
	2386	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
	2387	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
	2388	+falb_dir(:,6,:)*SFRWL(6)
	2389	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
[2692]	2390	END SELECt
[2469]	2391	!albedo SB <<<
[2227]	2392
[766]	2393
[2469]	2394	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
	2395	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
[1724]	2396
[2469]	2397	ENDIF
	2398	! =================================================================== c
	2399	! Calcul de Qsat
[881]	2400
[2469]	2401	DO k = 1, klev
	2402	DO i = 1, klon
	2403	zx_t = t_seri(i,k)
	2404	IF (thermcep) THEN
	2405	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2406	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2407	zx_qs = MIN(0.5,zx_qs)
	2408	zcor = 1./(1.-retv*zx_qs)
	2409	zx_qs = zx_qs*zcor
	2410	ELSE
	2411	!! IF (zx_t.LT.t_coup) THEN !jyg
	2412	IF (zx_t.LT.rtt) THEN !jyg
	2413	zx_qs = qsats(zx_t)/pplay(i,k)
	2414	ELSE
	2415	zx_qs = qsatl(zx_t)/pplay(i,k)
	2416	ENDIF
	2417	ENDIF
	2418	zqsat(i,k)=zx_qs
	2419	ENDDO
	2420	ENDDO
[959]	2421
[2692]	2422	IF (prt_level.ge.1) THEN
[2469]	2423	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	2424	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
[2692]	2425	ENDIF
[2469]	2426	!
	2427	! Appeler la convection (au choix)
	2428	!
	2429	DO k = 1, klev
	2430	DO i = 1, klon
	2431	conv_q(i,k) = d_q_dyn(i,k) &
[3065]	2432	+ d_q_vdf(i,k)/phys_tstep
[2469]	2433	conv_t(i,k) = d_t_dyn(i,k) &
[3065]	2434	+ d_t_vdf(i,k)/phys_tstep
[2469]	2435	ENDDO
	2436	ENDDO
	2437	IF (check) THEN
	2438	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	2439	WRITE(lunout,*) "avantcon=", za
	2440	ENDIF
	2441	zx_ajustq = .FALSE.
	2442	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	2443	IF (zx_ajustq) THEN
	2444	DO i = 1, klon
	2445	z_avant(i) = 0.0
	2446	ENDDO
	2447	DO k = 1, klev
	2448	DO i = 1, klon
	2449	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2450	*(paprs(i,k)-paprs(i,k+1))/RG
	2451	ENDDO
	2452	ENDDO
	2453	ENDIF
[959]	2454
[2469]	2455	! Calcule de vitesse verticale a partir de flux de masse verticale
	2456	DO k = 1, klev
	2457	DO i = 1, klon
	2458	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
[2692]	2459	ENDDO
	2460	ENDDO
	2461
	2462	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
[2469]	2463	omega(igout, :)
[2707]	2464	!
	2465	! Appel de la convection tous les "cvpas"
	2466	!
	2467	IF (MOD(itapcv,cvpas).EQ.0) THEN
	2468
[2469]	2469	IF (iflag_con.EQ.1) THEN
	2470	abort_message ='reactiver le call conlmd dans physiq.F'
	2471	CALL abort_physic (modname,abort_message,1)
[3065]	2472	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
[2469]	2473	! . d_t_con, d_q_con,
	2474	! . rain_con, snow_con, ibas_con, itop_con)
	2475	ELSE IF (iflag_con.EQ.2) THEN
[3065]	2476	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
[2469]	2477	conv_t, conv_q, -evap, omega, &
	2478	d_t_con, d_q_con, rain_con, snow_con, &
	2479	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	2480	kcbot, kctop, kdtop, pmflxr, pmflxs)
	2481	d_u_con = 0.
	2482	d_v_con = 0.
[879]	2483
[2469]	2484	WHERE (rain_con < 0.) rain_con = 0.
	2485	WHERE (snow_con < 0.) snow_con = 0.
	2486	DO i = 1, klon
	2487	ibas_con(i) = klev+1 - kcbot(i)
	2488	itop_con(i) = klev+1 - kctop(i)
	2489	ENDDO
	2490	ELSE IF (iflag_con.GE.3) THEN
	2491	! nb of tracers for the KE convection:
	2492	! MAF la partie traceurs est faite dans phytrac
	2493	! on met ntra=1 pour limiter les appels mais on peut
	2494	! supprimer les calculs / ftra.
	2495	ntra = 1
	2496
	2497	!=======================================================================
	2498	!ajout pour la parametrisation des poches froides: calcul de
[2635]	2499	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
[2692]	2500	IF (iflag_wake>=1) THEN
	2501	DO k=1,klev
	2502	DO i=1,klon
	2503	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
	2504	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
	2505	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
	2506	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
	2507	ENDDO
	2508	ENDDO
	2509	ELSE
	2510	t_w(:,:) = t_seri(:,:)
[2635]	2511	q_w(:,:) = q_seri(:,:)
	2512	t_x(:,:) = t_seri(:,:)
	2513	q_x(:,:) = q_seri(:,:)
[2692]	2514	ENDIF
[2469]	2515	!
	2516	!jyg<
	2517	! Perform dry adiabatic adjustment on wake profile
	2518	! The corresponding tendencies are added to the convective tendencies
	2519	! after the call to the convective scheme.
	2520	IF (iflag_wake>=1) then
[2882]	2521	IF (iflag_adjwk >= 1) THEN
[2469]	2522	limbas(:) = 1
[2635]	2523	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
[2309]	2524	d_t_adjwk, d_q_adjwk)
[2638]	2525	!
	2526	DO k=1,klev
	2527	DO i=1,klon
	2528	IF (wake_s(i) .GT. 1.e-3) THEN
	2529	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
	2530	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
	2531	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
	2532	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
	2533	ELSE
	2534	d_deltat_ajs_cv(i,k) = 0.
	2535	d_deltaq_ajs_cv(i,k) = 0.
	2536	ENDIF
	2537	ENDDO
[2469]	2538	ENDDO
[2882]	2539	IF (iflag_adjwk == 2) THEN
	2540	CALL add_wake_tend &
[2638]	2541	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
[2882]	2542	ENDIF ! (iflag_adjwk == 2)
	2543	ENDIF ! (iflag_adjwk >= 1)
[2469]	2544	ENDIF ! (iflag_wake>=1)
	2545	!>jyg
	2546	!
[2638]	2547
	2548	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
	2549	!! (k, q_w(1,k), q_x(1,k),k=1,25)
	2550
[2513]	2551	!jyg<
[3065]	2552	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
[2513]	2553	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
	2554	ale_bl_prescr, alp_bl_prescr, &
	2555	wake_pe, wake_fip, &
	2556	Ale_bl, Ale_bl_trig, Alp_bl, &
[2554]	2557	Ale, Alp , Ale_wake, Alp_wake)
[2513]	2558	!>jyg
	2559	!
[2469]	2560	! sb, oct02:
	2561	! Schema de convection modularise et vectorise:
	2562	! (driver commun aux versions 3 et 4)
	2563	!
	2564	IF (ok_cvl) THEN ! new driver for convectL
	2565	!
	2566	!jyg<
	2567	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2568	! Calculate the upmost level of deep convection loops: k_upper_cv
	2569	! (near 22 km)
	2570	izero = klon/2+1/klon
	2571	k_upper_cv = klev
	2572	DO k = klev,1,-1
	2573	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
	2574	ENDDO
	2575	IF (prt_level .ge. 5) THEN
	2576	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
	2577	k_upper_cv
	2578	ENDIF
	2579	!
	2580	!>jyg
	2581	IF (type_trac == 'repr') THEN
	2582	nbtr_tmp=ntra
	2583	ELSE
	2584	nbtr_tmp=nbtr
[2692]	2585	ENDIF
[2469]	2586	!jyg iflag_con est dans clesphys
	2587	!c CALL concvl (iflag_con,iflag_clos,
	2588	CALL concvl (iflag_clos, &
[3065]	2589	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
[2635]	2590	t_w,q_w,wake_s, &
[2469]	2591	u_seri,v_seri,tr_seri,nbtr_tmp, &
	2592	ALE,ALP, &
	2593	sig1,w01, &
	2594	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2595	rain_con, snow_con, ibas_con, itop_con, sigd, &
[2824]	2596	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
[2469]	2597	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
	2598	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
	2599	! RomP >>>
	2600	!! . pmflxr,pmflxs,da,phi,mp,
	2601	!! . ftd,fqd,lalim_conv,wght_th)
	2602	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
	2603	ftd,fqd,lalim_conv,wght_th, &
	2604	ev, ep,epmlmMm,eplaMm, &
	2605	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
[2481]	2606	tau_cld_cv,coefw_cld_cv,epmax_diag)
[2630]	2607
[2469]	2608	! RomP <<<
[619]	2609
[2469]	2610	!IM begin
	2611	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2612	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
	2613	!IM end
	2614	!IM cf. FH
	2615	clwcon0=qcondc
	2616	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[524]	2617
[2692]	2618	DO i = 1, klon
	2619	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2620	ENDDO
[2469]	2621	!
	2622	!jyg<
	2623	! Add the tendency due to the dry adjustment of the wake profile
	2624	IF (iflag_wake>=1) THEN
[2882]	2625	IF (iflag_adjwk == 2) THEN
	2626	DO k=1,klev
	2627	DO i=1,klon
[3065]	2628	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
	2629	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
[2882]	2630	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
	2631	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
	2632	ENDDO
	2633	ENDDO
	2634	ENDIF ! (iflag_adjwk = 2)
	2635	ENDIF ! (iflag_wake>=1)
[2469]	2636	!>jyg
	2637	!
	2638	ELSE ! ok_cvl
[1412]	2639
[2469]	2640	! MAF conema3 ne contient pas les traceurs
[3065]	2641	CALL conema3 (phys_tstep, &
[2469]	2642	paprs,pplay,t_seri,q_seri, &
	2643	u_seri,v_seri,tr_seri,ntra, &
	2644	sig1,w01, &
	2645	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2646	rain_con, snow_con, ibas_con, itop_con, &
	2647	upwd,dnwd,dnwd0,bas,top, &
	2648	Ma,cape,tvp,rflag, &
	2649	pbase &
	2650	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
	2651	,clwcon0)
[524]	2652
[2469]	2653	ENDIF ! ok_cvl
[524]	2654
[2469]	2655	!
	2656	! Correction precip
	2657	rain_con = rain_con * cvl_corr
	2658	snow_con = snow_con * cvl_corr
	2659	!
[766]	2660
[2469]	2661	IF (.NOT. ok_gust) THEN
	2662	do i = 1, klon
	2663	wd(i)=0.0
	2664	enddo
	2665	ENDIF
[524]	2666
[2469]	2667	! =================================================================== c
	2668	! Calcul des proprietes des nuages convectifs
	2669	!
[524]	2670
[2469]	2671	! calcul des proprietes des nuages convectifs
	2672	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2673	IF (iflag_cld_cv == 0) THEN
[2692]	2674	CALL clouds_gno &
[2469]	2675	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2676	ELSE
[2692]	2677	CALL clouds_bigauss &
[2469]	2678	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
	2679	ENDIF
[524]	2680
[2205]	2681
[2469]	2682	! =================================================================== c
[524]	2683
[2469]	2684	DO i = 1, klon
	2685	itop_con(i) = min(max(itop_con(i),1),klev)
	2686	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
	2687	ENDDO
[1428]	2688
[2469]	2689	DO i = 1, klon
	2690	ema_pcb(i) = paprs(i,ibas_con(i))
	2691	ENDDO
	2692	DO i = 1, klon
	2693	! L'idicage de itop_con peut cacher un pb potentiel
	2694	! FH sous la dictee de JYG, CR
	2695	ema_pct(i) = paprs(i,itop_con(i)+1)
[879]	2696
[2692]	2697	IF (itop_con(i).gt.klev-3) THEN
	2698	IF (prt_level >= 9) THEN
[2469]	2699	write(lunout,*)'La convection monte trop haut '
	2700	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
[2692]	2701	ENDIF
	2702	ENDIF
[2469]	2703	ENDDO
	2704	ELSE IF (iflag_con.eq.0) THEN
	2705	write(lunout,*) 'On n appelle pas la convection'
	2706	clwcon0=0.
	2707	rnebcon0=0.
	2708	d_t_con=0.
	2709	d_q_con=0.
	2710	d_u_con=0.
	2711	d_v_con=0.
	2712	rain_con=0.
	2713	snow_con=0.
	2714	bas=1
	2715	top=1
	2716	ELSE
	2717	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
[2692]	2718	CALL abort_physic("physiq", "", 1)
[2469]	2719	ENDIF
[524]	2720
[2469]	2721	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2722	! . d_u_con, d_v_con)
[524]	2723
[2730]	2724	!jyg Reinitialize proba_notrig and itapcv when convection has been called
	2725	proba_notrig(:) = 1.
[2707]	2726	itapcv = 0
	2727	ENDIF ! (MOD(itapcv,cvpas).EQ.0)
[2730]	2728	!
[2707]	2729	itapcv = itapcv+1
	2730
[2812]	2731	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
	2732	!!! l'energie dans les courants satures.
	2733	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
	2734	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
	2735	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
	2736	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
	2737	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
	2738	!! itap, 1)
	2739	!! call prt_enerbil('convection_sat',itap)
	2740	!!
	2741	!!
[2469]	2742	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
[2812]	2743	'convection',abortphy,flag_inhib_tend,itap,0)
[2799]	2744	call prt_enerbil('convection',itap)
[2235]	2745
[2469]	2746	!-------------------------------------------------------------------------
[766]	2747
[2692]	2748	IF (mydebug) THEN
	2749	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	2750	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	2751	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	2752	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	2753	ENDIF
[766]	2754
[2469]	2755	IF (check) THEN
	2756	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	2757	WRITE(lunout,*)"aprescon=", za
	2758	zx_t = 0.0
	2759	za = 0.0
	2760	DO i = 1, klon
	2761	za = za + cell_area(i)/REAL(klon)
	2762	zx_t = zx_t + (rain_con(i)+ &
	2763	snow_con(i))*cell_area(i)/REAL(klon)
	2764	ENDDO
[3065]	2765	zx_t = zx_t/za*phys_tstep
[2469]	2766	WRITE(lunout,*)"Precip=", zx_t
	2767	ENDIF
	2768	IF (zx_ajustq) THEN
	2769	DO i = 1, klon
	2770	z_apres(i) = 0.0
	2771	ENDDO
	2772	DO k = 1, klev
	2773	DO i = 1, klon
	2774	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2775	*(paprs(i,k)-paprs(i,k+1))/RG
	2776	ENDDO
	2777	ENDDO
	2778	DO i = 1, klon
[3065]	2779	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
[2469]	2780	/z_apres(i)
	2781	ENDDO
	2782	DO k = 1, klev
	2783	DO i = 1, klon
	2784	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
	2785	z_factor(i).LT.(1.0-1.0E-08)) THEN
	2786	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2787	ENDIF
	2788	ENDDO
	2789	ENDDO
	2790	ENDIF
	2791	zx_ajustq=.FALSE.
[879]	2792
[2469]	2793	!
	2794	!==========================================================================
	2795	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2796	!pour la couche limite diffuse pour l instant
	2797	!
	2798	!
	2799	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
	2800	! il faut rajouter cette tendance calcul\'ee hors des poches
	2801	! froides
	2802	!
[2692]	2803	IF (iflag_wake>=1) THEN
[2707]	2804	!
	2805	!
[2730]	2806	! Call wakes every "wkpas" step
	2807	!
	2808	IF (MOD(itapwk,wkpas).EQ.0) THEN
	2809	!
	2810	DO k=1,klev
[2469]	2811	DO i=1,klon
[2730]	2812	dt_dwn(i,k) = ftd(i,k)
	2813	dq_dwn(i,k) = fqd(i,k)
	2814	M_dwn(i,k) = dnwd0(i,k)
	2815	M_up(i,k) = upwd(i,k)
[3065]	2816	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
	2817	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
[2469]	2818	ENDDO
	2819	ENDDO
[2730]	2820
	2821	IF (iflag_wake==2) THEN
	2822	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2823	DO k = 1,klev
	2824	dt_dwn(:,k)= dt_dwn(:,k)+ &
[3065]	2825	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
[2730]	2826	dq_dwn(:,k)= dq_dwn(:,k)+ &
[3065]	2827	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
[2730]	2828	ENDDO
	2829	ELSEIF (iflag_wake==3) THEN
	2830	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2831	DO k = 1,klev
	2832	DO i=1,klon
	2833	IF (rneb(i,k)==0.) THEN
	2834	! On ne tient compte des tendances qu'en dehors des
	2835	! nuages (c'est-\`a-dire a priri dans une region ou
	2836	! l'eau se reevapore).
	2837	dt_dwn(i,k)= dt_dwn(i,k)+ &
[3065]	2838	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
[2730]	2839	dq_dwn(i,k)= dq_dwn(i,k)+ &
[3065]	2840	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
[2730]	2841	ENDIF
	2842	ENDDO
	2843	ENDDO
	2844	ENDIF
	2845
	2846	!
	2847	!calcul caracteristiques de la poche froide
[3065]	2848	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
[2730]	2849	t_seri, q_seri, omega, &
	2850	dt_dwn, dq_dwn, M_dwn, M_up, &
	2851	dt_a, dq_a, &
	2852	sigd, &
	2853	wake_deltat, wake_deltaq, wake_s, wake_dens, &
	2854	wake_dth, wake_h, &
[3000]	2855	!! wake_pe, wake_fip, wake_gfl, &
	2856	wake_pe, wake_fip_0, wake_gfl, & !! jyg
[2730]	2857	d_t_wake, d_q_wake, &
	2858	wake_k, t_x, q_x, &
	2859	wake_omgbdth, wake_dp_omgb, &
	2860	wake_dtKE, wake_dqKE, &
	2861	wake_omg, wake_dp_deltomg, &
	2862	wake_spread, wake_Cstar, d_deltat_wk_gw, &
	2863	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
	2864	!
	2865	!jyg Reinitialize itapwk when wakes have been called
	2866	itapwk = 0
	2867	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
[2469]	2868	!
[2730]	2869	itapwk = itapwk+1
[2469]	2870	!
	2871	!-----------------------------------------------------------------------
	2872	! ajout des tendances des poches froides
	2873	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
[2812]	2874	abortphy,flag_inhib_tend,itap,0)
[2799]	2875	call prt_enerbil('wake',itap)
[2469]	2876	!------------------------------------------------------------------------
[879]	2877
[2730]	2878	! Increment Wake state variables
[2635]	2879	IF (iflag_wake_tend .GT. 0.) THEN
	2880
	2881	CALL add_wake_tend &
	2882	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
	2883	'wake', abortphy)
[2799]	2884	call prt_enerbil('wake',itap)
[2635]	2885	ENDIF ! (iflag_wake_tend .GT. 0.)
	2886
[3000]	2887	IF (iflag_alp_wk_cond .GT. 0.) THEN
	2888
[3065]	2889	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
[3000]	2890	wake_fip)
	2891	ELSE
	2892	wake_fip(:) = wake_fip_0(:)
	2893	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
	2894
[2692]	2895	ENDIF ! (iflag_wake>=1)
[2469]	2896	!
	2897	!===================================================================
	2898	! Convection seche (thermiques ou ajustement)
	2899	!===================================================================
	2900	!
[2692]	2901	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
[2469]	2902	,seuil_inversion,weak_inversion,dthmin)
[878]	2903
	2904
	2905
[2469]	2906	d_t_ajsb(:,:)=0.
	2907	d_q_ajsb(:,:)=0.
	2908	d_t_ajs(:,:)=0.
	2909	d_u_ajs(:,:)=0.
	2910	d_v_ajs(:,:)=0.
	2911	d_q_ajs(:,:)=0.
	2912	clwcon0th(:,:)=0.
	2913	!
	2914	! fm_therm(:,:)=0.
	2915	! entr_therm(:,:)=0.
	2916	! detr_therm(:,:)=0.
	2917	!
[2692]	2918	IF (prt_level>9) WRITE(lunout,*) &
[2469]	2919	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
	2920	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[2692]	2921	IF (iflag_thermals<0) THEN
[2469]	2922	! Rien
	2923	! ====
[2692]	2924	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
[541]	2925
[878]	2926
[2692]	2927	ELSE
[878]	2928
[2469]	2929	! Thermiques
	2930	! ==========
[2692]	2931	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
[2469]	2932	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[878]	2933
	2934
[2469]	2935	!cc nrlmd le 10/04/2012
	2936	DO k=1,klev+1
	2937	DO i=1,klon
	2938	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
	2939	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
	2940	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
	2941	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
[2159]	2942	ENDDO
[2469]	2943	ENDDO
	2944	!cc fin nrlmd le 10/04/2012
[1403]	2945
[2692]	2946	IF (iflag_thermals>=1) THEN
[2469]	2947	!jyg<
[2852]	2948	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
	2949	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]	2950	! Appel des thermiques avec les profils exterieurs aux poches
	2951	DO k=1,klev
	2952	DO i=1,klon
	2953	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
	2954	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
[2606]	2955	u_therm(i,k) = u_seri(i,k)
	2956	v_therm(i,k) = v_seri(i,k)
[2469]	2957	ENDDO
	2958	ENDDO
	2959	ELSE
	2960	! Appel des thermiques avec les profils moyens
	2961	DO k=1,klev
	2962	DO i=1,klon
	2963	t_therm(i,k) = t_seri(i,k)
	2964	q_therm(i,k) = q_seri(i,k)
[2606]	2965	u_therm(i,k) = u_seri(i,k)
	2966	v_therm(i,k) = v_seri(i,k)
[2469]	2967	ENDDO
	2968	ENDDO
	2969	ENDIF
	2970	!>jyg
[2692]	2971	CALL calltherm(pdtphys &
[2469]	2972	,pplay,paprs,pphi,weak_inversion &
[2606]	2973	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
	2974	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
[2469]	2975	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
	2976	,fm_therm,entr_therm,detr_therm &
	2977	,zqasc,clwcon0th,lmax_th,ratqscth &
	2978	,ratqsdiff,zqsatth &
	2979	!on rajoute ale et alp, et les
	2980	!caracteristiques de la couche alim
	2981	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
	2982	,ztv,zpspsk,ztla,zthl &
	2983	!cc nrlmd le 10/04/2012
	2984	,pbl_tke_input,pctsrf,omega,cell_area &
	2985	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	2986	,n2,s2,ale_bl_stat &
	2987	,therm_tke_max,env_tke_max &
	2988	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	2989	,alp_bl_conv,alp_bl_stat &
	2990	!cc fin nrlmd le 10/04/2012
	2991	,zqla,ztva )
	2992	!
	2993	!jyg<
[2852]	2994	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
	2995	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]	2996	! Si les thermiques ne sont presents que hors des
	2997	! poches, la tendance moyenne associ\'ee doit etre
	2998	! multipliee par la fraction surfacique qu'ils couvrent.
	2999	DO k=1,klev
	3000	DO i=1,klon
	3001	!
[2635]	3002	d_deltat_the(i,k) = - d_t_ajs(i,k)
	3003	d_deltaq_the(i,k) = - d_q_ajs(i,k)
[2469]	3004	!
	3005	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
	3006	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
	3007	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
	3008	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
	3009	!
	3010	ENDDO
	3011	ENDDO
[2606]	3012	!
[2638]	3013	CALL add_wake_tend &
	3014	(d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
[2799]	3015	call prt_enerbil('the',itap)
[2638]	3016	!
[2852]	3017	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
[2638]	3018	!
[2606]	3019	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
[2812]	3020	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
[2799]	3021	call prt_enerbil('thermals',itap)
[2606]	3022	!
[2513]	3023	!
[3065]	3024	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
[2513]	3025	cin, s2, n2, &
	3026	ale_bl_trig, ale_bl_stat, ale_bl, &
[2556]	3027	alp_bl, alp_bl_stat, &
	3028	proba_notrig, random_notrig)
[2635]	3029	!>jyg
[1638]	3030
[2554]	3031	! ------------------------------------------------------------------
	3032	! Transport de la TKE par les panaches thermiques.
	3033	! FH : 2010/02/01
	3034	! if (iflag_pbl.eq.10) then
	3035	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
	3036	! s rg,paprs,pbl_tke)
	3037	! endif
	3038	! -------------------------------------------------------------------
	3039
[2692]	3040	DO i=1,klon
[2469]	3041	! zmax_th(i)=pphi(i,lmax_th(i))/rg
	3042	!CR:04/05/12:correction calcul zmax
	3043	zmax_th(i)=zmax0(i)
[2692]	3044	ENDDO
[1507]	3045
[2692]	3046	ENDIF
[878]	3047
[2469]	3048	! Ajustement sec
	3049	! ==============
[878]	3050
[2469]	3051	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
	3052	! a partir du sommet des thermiques.
	3053	! Dans le cas contraire, on demarre au niveau 1.
[878]	3054
[2692]	3055	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
[878]	3056
[2692]	3057	IF (iflag_thermals.eq.0) THEN
	3058	IF (prt_level>9) WRITE(lunout,*)'ajsec'
[2469]	3059	limbas(:)=1
[2692]	3060	ELSE
[2469]	3061	limbas(:)=lmax_th(:)
[2692]	3062	ENDIF
[878]	3063
[2469]	3064	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	3065	! pour des test de convergence numerique.
	3066	! Le nouveau ajsec est a priori mieux, meme pour le cas
	3067	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	3068	! non nulles numeriquement pour des mailles non concernees.
[878]	3069
[2692]	3070	IF (iflag_thermals==0) THEN
[2469]	3071	! Calling adjustment alone (but not the thermal plume model)
	3072	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
	3073	, d_t_ajsb, d_q_ajsb)
[2692]	3074	ELSE IF (iflag_thermals>0) THEN
[2469]	3075	! Calling adjustment above the top of thermal plumes
	3076	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
	3077	, d_t_ajsb, d_q_ajsb)
[2692]	3078	ENDIF
[878]	3079
[2469]	3080	!--------------------------------------------------------------------
	3081	! ajout des tendances de l'ajustement sec ou des thermiques
	3082	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
[2812]	3083	'ajsb',abortphy,flag_inhib_tend,itap,0)
[2799]	3084	call prt_enerbil('ajsb',itap)
[2469]	3085	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	3086	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
[904]	3087
[2469]	3088	!---------------------------------------------------------------------
[878]	3089
[2692]	3090	ENDIF
[524]	3091
[2692]	3092	ENDIF
[2469]	3093	!
	3094	!===================================================================
	3095	! Computation of ratqs, the width (normalized) of the subrid scale
	3096	! water distribution
	3097	CALL calcratqs(klon,klev,prt_level,lunout, &
	3098	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
[2534]	3099	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
	3100	tau_ratqs,fact_cldcon, &
[2469]	3101	ptconv,ptconvth,clwcon0th, rnebcon0th, &
	3102	paprs,pplay,q_seri,zqsat,fm_therm, &
	3103	ratqs,ratqsc)
[1032]	3104
[2100]	3105
[2469]	3106	!
	3107	! Appeler le processus de condensation a grande echelle
	3108	! et le processus de precipitation
	3109	!-------------------------------------------------------------------------
	3110	IF (prt_level .GE.10) THEN
	3111	print *,'itap, ->fisrtilp ',itap
	3112	ENDIF
	3113	!
[3065]	3114	CALL fisrtilp(phys_tstep,paprs,pplay, &
[2469]	3115	t_seri, q_seri,ptconv,ratqs, &
	3116	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
	3117	rain_lsc, snow_lsc, &
	3118	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	3119	frac_impa, frac_nucl, beta_prec_fisrt, &
	3120	prfl, psfl, rhcl, &
	3121	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
	3122	iflag_ice_thermo)
	3123	!
	3124	WHERE (rain_lsc < 0) rain_lsc = 0.
	3125	WHERE (snow_lsc < 0) snow_lsc = 0.
[766]	3126
[2799]	3127	!+JLD
	3128	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
	3129	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
	3130	! & ,rain_lsc,snow_lsc
	3131	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
	3132	!-JLD
[2469]	3133	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
[2812]	3134	'lsc',abortphy,flag_inhib_tend,itap,0)
[2799]	3135	call prt_enerbil('lsc',itap)
[2613]	3136	rain_num(:)=0.
[2657]	3137	DO k = 1, klev
[2613]	3138	DO i = 1, klon
	3139	IF (ql_seri(i,k)>oliqmax) THEN
	3140	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
	3141	ql_seri(i,k)=oliqmax
	3142	ENDIF
	3143	ENDDO
	3144	ENDDO
[2657]	3145	IF (nqo==3) THEN
	3146	DO k = 1, klev
	3147	DO i = 1, klon
	3148	IF (qs_seri(i,k)>oicemax) THEN
	3149	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
	3150	qs_seri(i,k)=oicemax
	3151	ENDIF
	3152	ENDDO
	3153	ENDDO
	3154	ENDIF
[2613]	3155
[2524]	3156	!---------------------------------------------------------------------------
[2469]	3157	DO k = 1, klev
	3158	DO i = 1, klon
	3159	cldfra(i,k) = rneb(i,k)
	3160	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
	3161	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	3162	ENDDO
	3163	ENDDO
	3164	IF (check) THEN
	3165	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	3166	WRITE(lunout,*)"apresilp=", za
	3167	zx_t = 0.0
	3168	za = 0.0
	3169	DO i = 1, klon
	3170	za = za + cell_area(i)/REAL(klon)
	3171	zx_t = zx_t + (rain_lsc(i) &
	3172	+ snow_lsc(i))*cell_area(i)/REAL(klon)
	3173	ENDDO
[3065]	3174	zx_t = zx_t/za*phys_tstep
[2469]	3175	WRITE(lunout,*)"Precip=", zx_t
	3176	ENDIF
[766]	3177
[2692]	3178	IF (mydebug) THEN
	3179	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	3180	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	3181	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	3182	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	3183	ENDIF
[524]	3184
[2469]	3185	!
	3186	!-------------------------------------------------------------------
	3187	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	3188	!-------------------------------------------------------------------
[524]	3189
[2469]	3190	! 1. NUAGES CONVECTIFS
	3191	!
	3192	!IM cf FH
	3193	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
	3194	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
	3195	snow_tiedtke=0.
	3196	! print*,'avant calcul de la pseudo precip '
	3197	! print*,'iflag_cld_th',iflag_cld_th
[2692]	3198	IF (iflag_cld_th.eq.-1) THEN
[2469]	3199	rain_tiedtke=rain_con
[2692]	3200	ELSE
[2469]	3201	! print*,'calcul de la pseudo precip '
	3202	rain_tiedtke=0.
	3203	! print*,'calcul de la pseudo precip 0'
[2692]	3204	DO k=1,klev
	3205	DO i=1,klon
	3206	IF (d_q_con(i,k).lt.0.) THEN
[2469]	3207	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
	3208	*(paprs(i,k)-paprs(i,k+1))/rg
[2692]	3209	ENDIF
	3210	ENDDO
	3211	ENDDO
	3212	ENDIF
[2469]	3213	!
	3214	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	3215	!
[524]	3216
[2469]	3217	! Nuages diagnostiques pour Tiedtke
	3218	CALL diagcld1(paprs,pplay, &
	3219	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
	3220	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
	3221	diafra,dialiq)
	3222	DO k = 1, klev
	3223	DO i = 1, klon
	3224	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3225	cldliq(i,k) = dialiq(i,k)
	3226	cldfra(i,k) = diafra(i,k)
	3227	ENDIF
	3228	ENDDO
	3229	ENDDO
[524]	3230
[2469]	3231	ELSE IF (iflag_cld_th.ge.3) THEN
	3232	! On prend pour les nuages convectifs le max du calcul de la
	3233	! convection et du calcul du pas de temps precedent diminue d'un facteur
	3234	! facttemps
	3235	facteur = pdtphys *facttemps
[2692]	3236	DO k=1,klev
	3237	DO i=1,klon
[2469]	3238	rnebcon(i,k)=rnebcon(i,k)*facteur
[2692]	3239	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
[2469]	3240	rnebcon(i,k)=rnebcon0(i,k)
	3241	clwcon(i,k)=clwcon0(i,k)
[2692]	3242	ENDIF
	3243	ENDDO
	3244	ENDDO
[2469]	3245
	3246	! On prend la somme des fractions nuageuses et des contenus en eau
[524]	3247
[2692]	3248	IF (iflag_cld_th>=5) THEN
[1411]	3249
[2692]	3250	DO k=1,klev
[2469]	3251	ptconvth(:,k)=fm_therm(:,k+1)>0.
[2692]	3252	ENDDO
[1496]	3253
[2692]	3254	IF (iflag_coupl==4) THEN
[1496]	3255
[2469]	3256	! Dans le cas iflag_coupl==4, on prend la somme des convertures
	3257	! convectives et lsc dans la partie des thermiques
	3258	! Le controle par iflag_coupl est peut etre provisoire.
[2692]	3259	DO k=1,klev
	3260	DO i=1,klon
	3261	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
[2469]	3262	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3263	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
[2692]	3264	ELSE IF (ptconv(i,k)) THEN
[2469]	3265	cldfra(i,k)=rnebcon(i,k)
	3266	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]	3267	ENDIF
	3268	ENDDO
	3269	ENDDO
[1496]	3270
[2692]	3271	ELSE IF (iflag_coupl==5) THEN
	3272	DO k=1,klev
	3273	DO i=1,klon
[2469]	3274	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3275	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
[2692]	3276	ENDDO
	3277	ENDDO
[1525]	3278
[2692]	3279	ELSE
[1525]	3280
[2469]	3281	! Si on est sur un point touche par la convection
	3282	! profonde et pas par les thermiques, on prend la
	3283	! couverture nuageuse et l'eau nuageuse de la convection
	3284	! profonde.
[1411]	3285
[2469]	3286	!IM/FH: 2011/02/23
	3287	! definition des points sur lesquels ls thermiques sont actifs
[1496]	3288
[2692]	3289	DO k=1,klev
	3290	DO i=1,klon
	3291	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
[2469]	3292	cldfra(i,k)=rnebcon(i,k)
	3293	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]	3294	ENDIF
	3295	ENDDO
	3296	ENDDO
[1496]	3297
[2692]	3298	ENDIF
[1496]	3299
[2692]	3300	ELSE
[1496]	3301
[2469]	3302	! Ancienne version
	3303	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	3304	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
[2692]	3305	ENDIF
[1411]	3306
[2469]	3307	ENDIF
[1507]	3308
[2469]	3309	! plulsc(:)=0.
	3310	! do k=1,klev,-1
	3311	! do i=1,klon
	3312	! zzz=prfl(:,k)+psfl(:,k)
	3313	! if (.not.ptconvth.zzz.gt.0.)
	3314	! enddo prfl, psfl,
	3315	! enddo
	3316	!
	3317	! 2. NUAGES STARTIFORMES
	3318	!
	3319	IF (ok_stratus) THEN
	3320	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	3321	DO k = 1, klev
	3322	DO i = 1, klon
	3323	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3324	cldliq(i,k) = dialiq(i,k)
	3325	cldfra(i,k) = diafra(i,k)
	3326	ENDIF
	3327	ENDDO
	3328	ENDDO
	3329	ENDIF
	3330	!
	3331	! Precipitation totale
	3332	!
	3333	DO i = 1, klon
	3334	rain_fall(i) = rain_con(i) + rain_lsc(i)
	3335	snow_fall(i) = snow_con(i) + snow_lsc(i)
	3336	ENDDO
	3337	!
	3338	! Calculer l'humidite relative pour diagnostique
	3339	!
	3340	DO k = 1, klev
	3341	DO i = 1, klon
	3342	zx_t = t_seri(i,k)
	3343	IF (thermcep) THEN
	3344	!! if (iflag_ice_thermo.eq.0) then !jyg
	3345	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	3346	!! else !jyg
	3347	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
	3348	!! endif !jyg
	3349	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	3350	zx_qs = MIN(0.5,zx_qs)
	3351	zcor = 1./(1.-retv*zx_qs)
	3352	zx_qs = zx_qs*zcor
	3353	ELSE
	3354	!! IF (zx_t.LT.t_coup) THEN !jyg
	3355	IF (zx_t.LT.rtt) THEN !jyg
	3356	zx_qs = qsats(zx_t)/pplay(i,k)
	3357	ELSE
	3358	zx_qs = qsatl(zx_t)/pplay(i,k)
	3359	ENDIF
	3360	ENDIF
	3361	zx_rh(i,k) = q_seri(i,k)/zx_qs
	3362	zqsat(i,k)=zx_qs
	3363	ENDDO
	3364	ENDDO
[782]	3365
[2469]	3366	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	3367	! equivalente a 2m (tpote) pour diagnostique
	3368	!
	3369	DO i = 1, klon
	3370	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	3371	IF (thermcep) THEN
	3372	IF(zt2m(i).LT.RTT) then
	3373	Lheat=RLSTT
	3374	ELSE
	3375	Lheat=RLVTT
	3376	ENDIF
	3377	ELSE
	3378	IF (zt2m(i).LT.RTT) THEN
	3379	Lheat=RLSTT
	3380	ELSE
	3381	Lheat=RLVTT
	3382	ENDIF
	3383	ENDIF
	3384	tpote(i) = tpot(i)* &
	3385	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	3386	ENDDO
[524]	3387
[2469]	3388	IF (type_trac == 'inca') THEN
[524]	3389	#ifdef INCA
[2469]	3390	CALL VTe(VTphysiq)
	3391	CALL VTb(VTinca)
	3392	calday = REAL(days_elapsed + 1) + jH_cur
[524]	3393
[3065]	3394	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
[2469]	3395	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
	3396	CALL AEROSOL_METEO_CALC( &
	3397	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
	3398	prfl,psfl,pctsrf,cell_area, &
	3399	latitude_deg,longitude_deg,u10m,v10m)
[2692]	3400	ENDIF
[524]	3401
[2469]	3402	zxsnow_dummy(:) = 0.0
[625]	3403
[2469]	3404	CALL chemhook_begin (calday, &
	3405	days_elapsed+1, &
	3406	jH_cur, &
	3407	pctsrf(1,1), &
	3408	latitude_deg, &
	3409	longitude_deg, &
	3410	cell_area, &
	3411	paprs, &
	3412	pplay, &
	3413	coefh(1:klon,1:klev,is_ave), &
	3414	pphi, &
	3415	t_seri, &
	3416	u, &
	3417	v, &
	3418	wo(:, :, 1), &
	3419	q_seri, &
	3420	zxtsol, &
	3421	zxsnow_dummy, &
	3422	solsw, &
	3423	albsol1, &
	3424	rain_fall, &
	3425	snow_fall, &
	3426	itop_con, &
	3427	ibas_con, &
	3428	cldfra, &
	3429	nbp_lon, &
	3430	nbp_lat-1, &
	3431	tr_seri, &
	3432	ftsol, &
	3433	paprs, &
	3434	cdragh, &
	3435	cdragm, &
	3436	pctsrf, &
	3437	pdtphys, &
	3438	itap)
[616]	3439
[2469]	3440	CALL VTe(VTinca)
	3441	CALL VTb(VTphysiq)
[959]	3442	#endif
[2692]	3443	ENDIF !type_trac = inca
[2618]	3444
	3445
[2469]	3446	!
[2618]	3447	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	3448	!
	3449	IF (MOD(itaprad,radpas).EQ.0) THEN
[959]	3450
[2618]	3451	!
	3452	!jq - introduce the aerosol direct and first indirect radiative forcings
	3453	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
[2738]	3454	IF (flag_aerosol .GT. 0) THEN
[2618]	3455	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3456	IF (.NOT. aerosol_couple) THEN
	3457	!
	3458	CALL readaerosol_optic( &
	3459	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3460	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3461	mass_solu_aero, mass_solu_aero_pi, &
	3462	tau_aero, piz_aero, cg_aero, &
	3463	tausum_aero, tau3d_aero)
	3464	ENDIF
	3465	ELSE ! RRTM radiation
	3466	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
	3467	abort_message='config_inca=aero et rrtm=1 impossible'
[2692]	3468	CALL abort_physic(modname,abort_message,1)
[2618]	3469	ELSE
	3470	!
	3471	#ifdef CPP_RRTM
	3472	IF (NSW.EQ.6) THEN
[2738]	3473	!--new aerosol properties SW and LW
[2618]	3474	!
[2753]	3475	#ifdef CPP_Dust
	3476	!--SPL aerosol model
	3477	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
	3478	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
	3479	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3480	tausum_aero, tau3d_aero)
	3481	#else
	3482	!--climatologies or INCA aerosols
[2738]	3483	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
[2644]	3484	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
[2618]	3485	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3486	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
	3487	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
[2854]	3488	tausum_aero, drytausum_aero, tau3d_aero)
[2753]	3489	#endif
[2738]	3490	!
[2618]	3491	ELSE IF (NSW.EQ.2) THEN
	3492	!--for now we use the old aerosol properties
	3493	!
	3494	CALL readaerosol_optic( &
	3495	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3496	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3497	mass_solu_aero, mass_solu_aero_pi, &
	3498	tau_aero, piz_aero, cg_aero, &
	3499	tausum_aero, tau3d_aero)
	3500	!
	3501	!--natural aerosols
	3502	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
	3503	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
	3504	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
	3505	!--all aerosols
	3506	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
	3507	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
	3508	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
[2738]	3509	!
	3510	!--no LW optics
	3511	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
	3512	!
[2618]	3513	ELSE
	3514	abort_message='Only NSW=2 or 6 are possible with ' &
	3515	// 'aerosols and iflag_rrtm=1'
[2692]	3516	CALL abort_physic(modname,abort_message,1)
[2618]	3517	ENDIF
	3518	#else
	3519	abort_message='You should compile with -rrtm if running ' &
	3520	// 'with iflag_rrtm=1'
[2692]	3521	CALL abort_physic(modname,abort_message,1)
[2618]	3522	#endif
	3523	!
	3524	ENDIF
	3525	ENDIF
[2738]	3526	ELSE !--flag_aerosol = 0
[2618]	3527	tausum_aero(:,:,:) = 0.
[2854]	3528	drytausum_aero(:,:) = 0.
[2640]	3529	mass_solu_aero(:,:) = 0.
	3530	mass_solu_aero_pi(:,:) = 0.
[2618]	3531	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3532	tau_aero(:,:,:,:) = 1.e-15
	3533	piz_aero(:,:,:,:) = 1.
	3534	cg_aero(:,:,:,:) = 0.
	3535	ELSE
	3536	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
	3537	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
	3538	piz_aero_sw_rrtm(:,:,:,:) = 1.0
	3539	cg_aero_sw_rrtm(:,:,:,:) = 0.0
	3540	ENDIF
	3541	ENDIF
	3542	!
[2994]	3543	!--WMO criterion to determine tropopause
	3544	CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
	3545	!
[2618]	3546	!--STRAT AEROSOL
	3547	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
	3548	IF (flag_aerosol_strat.GT.0) THEN
	3549	IF (prt_level .GE.10) THEN
	3550	PRINT *,'appel a readaerosolstrat', mth_cur
	3551	ENDIF
	3552	IF (iflag_rrtm.EQ.0) THEN
	3553	IF (flag_aerosol_strat.EQ.1) THEN
	3554	CALL readaerosolstrato(debut)
	3555	ELSE
	3556	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
	3557	CALL abort_physic(modname,abort_message,1)
	3558	ENDIF
	3559	ELSE
[2009]	3560	#ifdef CPP_RRTM
[2690]	3561	#ifndef CPP_StratAer
	3562	!--prescribed strat aerosols
	3563	!--only in the case of non-interactive strat aerosols
[2618]	3564	IF (flag_aerosol_strat.EQ.1) THEN
	3565	CALL readaerosolstrato1_rrtm(debut)
	3566	ELSEIF (flag_aerosol_strat.EQ.2) THEN
	3567	CALL readaerosolstrato2_rrtm(debut)
	3568	ELSE
	3569	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
	3570	CALL abort_physic(modname,abort_message,1)
	3571	ENDIF
[2690]	3572	#endif
[2618]	3573	#else
	3574	abort_message='You should compile with -rrtm if running ' &
	3575	// 'with iflag_rrtm=1'
	3576	CALL abort_physic(modname,abort_message,1)
	3577	#endif
	3578	ENDIF
	3579	ENDIF
[2690]	3580	!
	3581	#ifdef CPP_RRTM
	3582	#ifdef CPP_StratAer
[2692]	3583	!--compute stratospheric mask
[2992]	3584	CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
[2690]	3585	!--interactive strat aerosols
	3586	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
	3587	#endif
	3588	#endif
[2618]	3589	!--fin STRAT AEROSOL
	3590	!
	3591
	3592	! Calculer les parametres optiques des nuages et quelques
	3593	! parametres pour diagnostiques:
	3594	!
	3595	IF (aerosol_couple.AND.config_inca=='aero') THEN
	3596	mass_solu_aero(:,:) = ccm(:,:,1)
	3597	mass_solu_aero_pi(:,:) = ccm(:,:,2)
[2692]	3598	ENDIF
[2618]	3599
	3600	IF (ok_newmicro) then
	3601	IF (iflag_rrtm.NE.0) THEN
	3602	#ifdef CPP_RRTM
	3603	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
[2469]	3604	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
	3605	// 'pour ok_cdnc'
[2618]	3606	CALL abort_physic(modname,abort_message,1)
	3607	ENDIF
[2009]	3608	#else
	3609
[2618]	3610	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
	3611	CALL abort_physic(modname,abort_message,1)
[2009]	3612	#endif
[2618]	3613	ENDIF
	3614	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
	3615	paprs, pplay, t_seri, cldliq, cldfra, &
	3616	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
	3617	flwp, fiwp, flwc, fiwc, &
	3618	mass_solu_aero, mass_solu_aero_pi, &
	3619	cldtaupi, re, fl, ref_liq, ref_ice, &
	3620	ref_liq_pi, ref_ice_pi)
	3621	ELSE
	3622	CALL nuage (paprs, pplay, &
	3623	t_seri, cldliq, cldfra, cldtau, cldemi, &
	3624	cldh, cldl, cldm, cldt, cldq, &
	3625	ok_aie, &
	3626	mass_solu_aero, mass_solu_aero_pi, &
	3627	bl95_b0, bl95_b1, &
	3628	cldtaupi, re, fl)
[2469]	3629	ENDIF
	3630	!
[2618]	3631	!IM betaCRF
[2469]	3632	!
[2618]	3633	cldtaurad = cldtau
	3634	cldtaupirad = cldtaupi
	3635	cldemirad = cldemi
	3636	cldfrarad = cldfra
	3637
[2469]	3638	!
[2618]	3639	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
	3640	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
	3641	!
	3642	! global
	3643	!
	3644	DO k=1, klev
	3645	DO i=1, klon
	3646	IF (pplay(i,k).GE.pfree) THEN
[2469]	3647	beta(i,k) = beta_pbl
[2618]	3648	ELSE
[2469]	3649	beta(i,k) = beta_free
[2618]	3650	ENDIF
	3651	IF (mskocean_beta) THEN
[2469]	3652	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
[2618]	3653	ENDIF
[2469]	3654	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3655	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3656	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3657	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
[2618]	3658	ENDDO
	3659	ENDDO
	3660	!
	3661	ELSE
	3662	!
	3663	! regional
	3664	!
	3665	DO k=1, klev
	3666	DO i=1,klon
	3667	!
	3668	IF (longitude_deg(i).ge.lon1_beta.AND. &
	3669	longitude_deg(i).le.lon2_beta.AND. &
	3670	latitude_deg(i).le.lat1_beta.AND. &
	3671	latitude_deg(i).ge.lat2_beta) THEN
	3672	IF (pplay(i,k).GE.pfree) THEN
	3673	beta(i,k) = beta_pbl
	3674	ELSE
	3675	beta(i,k) = beta_free
	3676	ENDIF
	3677	IF (mskocean_beta) THEN
	3678	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3679	ENDIF
	3680	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3681	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3682	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3683	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3684	ENDIF
[2469]	3685	!
[2618]	3686	ENDDO
[2469]	3687	ENDDO
	3688	!
[2618]	3689	ENDIF
[766]	3690
[2618]	3691	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
	3692	IF (ok_chlorophyll) THEN
[2469]	3693	print*,"-- reading chlorophyll"
[2618]	3694	CALL readchlorophyll(debut)
	3695	ENDIF
[1863]	3696
[2524]	3697	!--if ok_suntime_rrtm we use ancillay data for RSUN
	3698	!--previous values are therefore overwritten
	3699	!--this is needed for CMIP6 runs
	3700	!--and only possible for new radiation scheme
	3701	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
[2525]	3702	#ifdef CPP_RRTM
[2524]	3703	CALL read_rsun_rrtm(debut)
[2525]	3704	#endif
[2524]	3705	ENDIF
	3706
[2692]	3707	IF (mydebug) THEN
	3708	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	3709	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	3710	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	3711	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	3712	ENDIF
[2524]	3713
[2469]	3714	!
	3715	!sonia : If Iflag_radia >=2, pertubation of some variables
	3716	!input to radiation (DICE)
	3717	!
	3718	IF (iflag_radia .ge. 2) THEN
	3719	zsav_tsol (:) = zxtsol(:)
[2692]	3720	CALL perturb_radlwsw(zxtsol,iflag_radia)
[2469]	3721	ENDIF
[2328]	3722
[2469]	3723	IF (aerosol_couple.AND.config_inca=='aero') THEN
[959]	3724	#ifdef INCA
[2469]	3725	CALL radlwsw_inca &
	3726	(kdlon,kflev,dist, rmu0, fract, solaire, &
	3727	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
[2684]	3728	size(wo,3), wo, &
[2469]	3729	cldfrarad, cldemirad, cldtaurad, &
	3730	heat,heat0,cool,cool0,albpla, &
	3731	topsw,toplw,solsw,sollw, &
	3732	sollwdown, &
	3733	topsw0,toplw0,solsw0,sollw0, &
	3734	lwdn0, lwdn, lwup0, lwup, &
	3735	swdn0, swdn, swup0, swup, &
	3736	ok_ade, ok_aie, &
	3737	tau_aero, piz_aero, cg_aero, &
	3738	topswad_aero, solswad_aero, &
	3739	topswad0_aero, solswad0_aero, &
	3740	topsw_aero, topsw0_aero, &
	3741	solsw_aero, solsw0_aero, &
	3742	cldtaupirad, &
	3743	topswai_aero, solswai_aero)
[955]	3744	#endif
[2469]	3745	ELSE
	3746	!
	3747	!IM calcul radiatif pour le cas actuel
	3748	!
	3749	RCO2 = RCO2_act
	3750	RCH4 = RCH4_act
	3751	RN2O = RN2O_act
	3752	RCFC11 = RCFC11_act
	3753	RCFC12 = RCFC12_act
	3754	!
	3755	IF (prt_level .GE.10) THEN
	3756	print *,' ->radlwsw, number 1 '
	3757	ENDIF
	3758	!
	3759	CALL radlwsw &
	3760	(dist, rmu0, fract, &
	3761	!albedo SB >>>
	3762	! paprs, pplay,zxtsol,albsol1, albsol2, &
	3763	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
	3764	!albedo SB <<<
	3765	t_seri,q_seri,wo, &
	3766	cldfrarad, cldemirad, cldtaurad, &
[2530]	3767	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
[2469]	3768	flag_aerosol_strat, &
	3769	tau_aero, piz_aero, cg_aero, &
	3770	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3771	! Rajoute par OB pour RRTM
	3772	tau_aero_lw_rrtm, &
	3773	cldtaupirad,new_aod, &
	3774	zqsat, flwc, fiwc, &
	3775	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
	3776	heat,heat0,cool,cool0,albpla, &
	3777	topsw,toplw,solsw,sollw, &
	3778	sollwdown, &
	3779	topsw0,toplw0,solsw0,sollw0, &
	3780	lwdn0, lwdn, lwup0, lwup, &
	3781	swdn0, swdn, swup0, swup, &
	3782	topswad_aero, solswad_aero, &
	3783	topswai_aero, solswai_aero, &
	3784	topswad0_aero, solswad0_aero, &
	3785	topsw_aero, topsw0_aero, &
	3786	solsw_aero, solsw0_aero, &
	3787	topswcf_aero, solswcf_aero, &
	3788	!-C. Kleinschmitt for LW diagnostics
	3789	toplwad_aero, sollwad_aero,&
	3790	toplwai_aero, sollwai_aero, &
	3791	toplwad0_aero, sollwad0_aero,&
	3792	!-end
	3793	ZLWFT0_i, ZFLDN0, ZFLUP0, &
	3794	ZSWFT0_i, ZFSDN0, ZFSUP0)
[879]	3795
[2679]	3796	#ifndef CPP_XIOS
	3797	!--OB 30/05/2016 modified 21/10/2016
[2854]	3798	!--here we return swaero_diag and dryaod_diag to FALSE
[2529]	3799	!--and histdef will switch it back to TRUE if necessary
	3800	!--this is necessary to get the right swaero at first step
[2679]	3801	!--but only in the case of no XIOS as XIOS is covered elsewhere
[2529]	3802	IF (debut) swaero_diag = .FALSE.
[2854]	3803	IF (debut) dryaod_diag = .FALSE.
[2989]	3804	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
	3805	!--as for swaero_diag, see above
	3806	IF (debut) ok_4xCO2atm = .FALSE.
[2679]	3807	#endif
[2469]	3808	!
	3809	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
	3810	!IM des taux doit etre different du taux actuel
	3811	!IM Par defaut on a les taux perturbes egaux aux taux actuels
	3812	!
[2989]	3813	IF (RCO2_per.NE.RCO2_act.OR. &
	3814	RCH4_per.NE.RCH4_act.OR. &
	3815	RN2O_per.NE.RN2O_act.OR. &
	3816	RCFC11_per.NE.RCFC11_act.OR. &
	3817	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
	3818	!
[2692]	3819	IF (ok_4xCO2atm) THEN
[2469]	3820	!
	3821	RCO2 = RCO2_per
	3822	RCH4 = RCH4_per
	3823	RN2O = RN2O_per
	3824	RCFC11 = RCFC11_per
	3825	RCFC12 = RCFC12_per
	3826	!
	3827	IF (prt_level .GE.10) THEN
	3828	print *,' ->radlwsw, number 2 '
	3829	ENDIF
	3830	!
	3831	CALL radlwsw &
	3832	(dist, rmu0, fract, &
	3833	!albedo SB >>>
	3834	! paprs, pplay,zxtsol,albsol1, albsol2, &
	3835	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
	3836	!albedo SB <<<
	3837	t_seri,q_seri,wo, &
[2640]	3838	cldfrarad, cldemirad, cldtaurad, &
[2530]	3839	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
[2469]	3840	flag_aerosol_strat, &
	3841	tau_aero, piz_aero, cg_aero, &
	3842	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3843	! Rajoute par OB pour RRTM
	3844	tau_aero_lw_rrtm, &
	3845	cldtaupi,new_aod, &
	3846	zqsat, flwc, fiwc, &
	3847	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
	3848	heatp,heat0p,coolp,cool0p,albplap, &
	3849	topswp,toplwp,solswp,sollwp, &
	3850	sollwdownp, &
	3851	topsw0p,toplw0p,solsw0p,sollw0p, &
	3852	lwdn0p, lwdnp, lwup0p, lwupp, &
	3853	swdn0p, swdnp, swup0p, swupp, &
	3854	topswad_aerop, solswad_aerop, &
	3855	topswai_aerop, solswai_aerop, &
	3856	topswad0_aerop, solswad0_aerop, &
	3857	topsw_aerop, topsw0_aerop, &
	3858	solsw_aerop, solsw0_aerop, &
	3859	topswcf_aerop, solswcf_aerop, &
	3860	!-C. Kleinschmitt for LW diagnostics
	3861	toplwad_aerop, sollwad_aerop,&
	3862	toplwai_aerop, sollwai_aerop, &
	3863	toplwad0_aerop, sollwad0_aerop,&
	3864	!-end
	3865	ZLWFT0_i, ZFLDN0, ZFLUP0, &
	3866	ZSWFT0_i, ZFSDN0, ZFSUP0)
	3867	endif
	3868	!
	3869	ENDIF ! aerosol_couple
	3870	itaprad = 0
	3871	!
	3872	! If Iflag_radia >=2, reset pertubed variables
	3873	!
	3874	IF (iflag_radia .ge. 2) THEN
	3875	zxtsol(:) = zsav_tsol (:)
	3876	ENDIF
	3877	ENDIF ! MOD(itaprad,radpas)
	3878	itaprad = itaprad + 1
[879]	3879
[2469]	3880	IF (iflag_radia.eq.0) THEN
	3881	IF (prt_level.ge.9) THEN
	3882	PRINT *,'--------------------------------------------------'
	3883	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	3884	PRINT *,'>>>> heat et cool mis a zero '
	3885	PRINT *,'--------------------------------------------------'
[2692]	3886	ENDIF
[2469]	3887	heat=0.
	3888	cool=0.
	3889	sollw=0. ! MPL 01032011
	3890	solsw=0.
	3891	radsol=0.
	3892	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
	3893	swup0=0.
	3894	lwup=0.
	3895	lwup0=0.
	3896	lwdn=0.
	3897	lwdn0=0.
[2692]	3898	ENDIF
[782]	3899
[2469]	3900	!
	3901	! Calculer radsol a l'exterieur de radlwsw
	3902	! pour prendre en compte le cycle diurne
	3903	! recode par Olivier Boucher en sept 2015
	3904	!
	3905	radsol=solsw*swradcorr+sollw
[2618]	3906
[2692]	3907	IF (ok_4xCO2atm) THEN
[2469]	3908	radsolp=solswp*swradcorr+sollwp
[2692]	3909	ENDIF
[2359]	3910
[2469]	3911	!
	3912	! Ajouter la tendance des rayonnements (tous les pas)
	3913	! avec une correction pour le cycle diurne dans le SW
	3914	!
[2359]	3915
[2469]	3916	DO k=1, klev
[3065]	3917	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
	3918	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
	3919	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
	3920	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
[2469]	3921	ENDDO
[2194]	3922
[2812]	3923	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
[2799]	3924	call prt_enerbil('SW',itap)
[2812]	3925	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
[2799]	3926	call prt_enerbil('LW',itap)
[1863]	3927
[2469]	3928	!
[2692]	3929	IF (mydebug) THEN
	3930	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	3931	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	3932	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	3933	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	3934	ENDIF
[1863]	3935
[2469]	3936	! Calculer l'hydrologie de la surface
	3937	!
	3938	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	3939	! . agesno, ftsol,fqsurf,fsnow, ruis)
	3940	!
[1001]	3941
[2469]	3942	!
	3943	! Calculer le bilan du sol et la derive de temperature (couplage)
	3944	!
	3945	DO i = 1, klon
	3946	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	3947	! a la demande de JLD
	3948	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	3949	ENDDO
	3950	!
	3951	!moddeblott(jan95)
	3952	! Appeler le programme de parametrisation de l'orographie
	3953	! a l'echelle sous-maille:
	3954	!
	3955	IF (prt_level .GE.10) THEN
	3956	print *,' call orography ? ', ok_orodr
	3957	ENDIF
	3958	!
	3959	IF (ok_orodr) THEN
	3960	!
	3961	! selection des points pour lesquels le shema est actif:
	3962	igwd=0
	3963	DO i=1,klon
	3964	itest(i)=0
	3965	! IF ((zstd(i).gt.10.0)) THEN
	3966	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	3967	itest(i)=1
	3968	igwd=igwd+1
	3969	idx(igwd)=i
	3970	ENDIF
	3971	ENDDO
	3972	! igwdim=MAX(1,igwd)
	3973	!
	3974	IF (ok_strato) THEN
[1863]	3975
[3065]	3976	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
[2469]	3977	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3978	igwd,idx,itest, &
	3979	t_seri, u_seri, v_seri, &
	3980	zulow, zvlow, zustrdr, zvstrdr, &
	3981	d_t_oro, d_u_oro, d_v_oro)
[1863]	3982
[2469]	3983	ELSE
[3065]	3984	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]	3985	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3986	igwd,idx,itest, &
	3987	t_seri, u_seri, v_seri, &
	3988	zulow, zvlow, zustrdr, zvstrdr, &
	3989	d_t_oro, d_u_oro, d_v_oro)
	3990	ENDIF
	3991	!
	3992	! ajout des tendances
	3993	!-----------------------------------------------------------------------
	3994	! ajout des tendances de la trainee de l'orographie
	3995	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
[2812]	3996	abortphy,flag_inhib_tend,itap,0)
[2799]	3997	call prt_enerbil('oro',itap)
[2469]	3998	!----------------------------------------------------------------------
	3999	!
	4000	ENDIF ! fin de test sur ok_orodr
	4001	!
[2692]	4002	IF (mydebug) THEN
	4003	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	4004	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	4005	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	4006	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	4007	ENDIF
[1001]	4008
[2469]	4009	IF (ok_orolf) THEN
	4010	!
	4011	! selection des points pour lesquels le shema est actif:
	4012	igwd=0
	4013	DO i=1,klon
	4014	itest(i)=0
	4015	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	4016	itest(i)=1
	4017	igwd=igwd+1
	4018	idx(igwd)=i
	4019	ENDIF
	4020	ENDDO
	4021	! igwdim=MAX(1,igwd)
	4022	!
	4023	IF (ok_strato) THEN
[1001]	4024
[3065]	4025	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
[2469]	4026	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
	4027	igwd,idx,itest, &
	4028	t_seri, u_seri, v_seri, &
	4029	zulow, zvlow, zustrli, zvstrli, &
	4030	d_t_lif, d_u_lif, d_v_lif )
[2333]	4031
[2469]	4032	ELSE
[3065]	4033	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]	4034	latitude_deg,zmea,zstd,zpic, &
	4035	itest, &
	4036	t_seri, u_seri, v_seri, &
	4037	zulow, zvlow, zustrli, zvstrli, &
	4038	d_t_lif, d_u_lif, d_v_lif)
	4039	ENDIF
[1638]	4040
[2469]	4041	! ajout des tendances de la portance de l'orographie
	4042	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
[2812]	4043	'lif', abortphy,flag_inhib_tend,itap,0)
[2799]	4044	call prt_enerbil('lif',itap)
[2469]	4045	ENDIF ! fin de test sur ok_orolf
[1638]	4046
[2469]	4047	IF (ok_hines) then
	4048	! HINES GWD PARAMETRIZATION
	4049	east_gwstress=0.
	4050	west_gwstress=0.
	4051	du_gwd_hines=0.
	4052	dv_gwd_hines=0.
[3065]	4053	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
[2469]	4054	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
	4055	du_gwd_hines, dv_gwd_hines)
	4056	zustr_gwd_hines=0.
	4057	zvstr_gwd_hines=0.
	4058	DO k = 1, klev
[3065]	4059	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
[2469]	4060	* (paprs(:, k)-paprs(:, k+1))/rg
[3065]	4061	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
[2469]	4062	* (paprs(:, k)-paprs(:, k+1))/rg
	4063	ENDDO
[1001]	4064
[2469]	4065	d_t_hin(:, :)=0.
	4066	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
[2812]	4067	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
[2799]	4068	call prt_enerbil('hin',itap)
[2469]	4069	ENDIF
[2333]	4070
[2469]	4071	IF (.not. ok_hines .and. ok_gwd_rando) then
[3406]	4072	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
[3065]	4073	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
[2469]	4074	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
	4075	dv_gwd_front, east_gwstress, west_gwstress)
	4076	zustr_gwd_front=0.
	4077	zvstr_gwd_front=0.
	4078	DO k = 1, klev
[3065]	4079	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
[2469]	4080	* (paprs(:, k)-paprs(:, k+1))/rg
[3065]	4081	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
[2469]	4082	* (paprs(:, k)-paprs(:, k+1))/rg
	4083	ENDDO
[644]	4084
[2469]	4085	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
[2812]	4086	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[2799]	4087	call prt_enerbil('front_gwd_rando',itap)
[2469]	4088	ENDIF
[1938]	4089
[2692]	4090	IF (ok_gwd_rando) THEN
[3065]	4091	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
[2469]	4092	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
	4093	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
	4094	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
[2812]	4095	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[2799]	4096	call prt_enerbil('flott_gwd_rando',itap)
[2469]	4097	zustr_gwd_rando=0.
	4098	zvstr_gwd_rando=0.
	4099	DO k = 1, klev
[3065]	4100	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
[2469]	4101	* (paprs(:, k)-paprs(:, k+1))/rg
[3065]	4102	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
[2469]	4103	* (paprs(:, k)-paprs(:, k+1))/rg
	4104	ENDDO
[2692]	4105	ENDIF
[766]	4106
[2469]	4107	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
[1279]	4108
[2692]	4109	IF (mydebug) THEN
	4110	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	4111	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	4112	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	4113	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	4114	ENDIF
[2136]	4115
[2469]	4116	DO i = 1, klon
	4117	zustrph(i)=0.
	4118	zvstrph(i)=0.
	4119	ENDDO
	4120	DO k = 1, klev
	4121	DO i = 1, klon
[3065]	4122	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
[2469]	4123	(paprs(i,k)-paprs(i,k+1))/rg
[3065]	4124	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
[2469]	4125	(paprs(i,k)-paprs(i,k+1))/rg
	4126	ENDDO
	4127	ENDDO
	4128	!
	4129	!IM calcul composantes axiales du moment angulaire et couple des montagnes
	4130	!
	4131	IF (is_sequential .and. ok_orodr) THEN
	4132	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
	4133	ra,rg,romega, &
	4134	latitude_deg,longitude_deg,pphis, &
	4135	zustrdr,zustrli,zustrph, &
	4136	zvstrdr,zvstrli,zvstrph, &
	4137	paprs,u,v, &
	4138	aam, torsfc)
	4139	ENDIF
	4140	!IM cf. FLott END
	4141	!DC Calcul de la tendance due au methane
	4142	IF(ok_qch4) THEN
	4143	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
	4144	! ajout de la tendance d'humidite due au methane
[3065]	4145	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
[2801]	4146	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
[2812]	4147	'q_ch4', abortphy,flag_inhib_tend,itap,0)
[3065]	4148	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
[2692]	4149	ENDIF
[2469]	4150	!
	4151	!
[2897]	4152
	4153	!===============================================================
	4154	! Additional tendency of TKE due to orography
	4155	!===============================================================
	4156	!
	4157	! Inititialization
	4158	!------------------
	4159
	4160
	4161
	4162	addtkeoro=0
	4163	CALL getin_p('addtkeoro',addtkeoro)
	4164
	4165	IF (prt_level.ge.5) &
	4166	print*,'addtkeoro', addtkeoro
	4167
	4168	alphatkeoro=1.
	4169	CALL getin_p('alphatkeoro',alphatkeoro)
	4170	alphatkeoro=min(max(0.,alphatkeoro),1.)
	4171
	4172	smallscales_tkeoro=.false.
	4173	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
	4174
	4175
	4176	dtadd(:,:)=0.
	4177	duadd(:,:)=0.
	4178	dvadd(:,:)=0.
	4179
	4180
	4181
	4182	! Choices for addtkeoro:
	4183	! ** 0 no TKE tendency from orography
	4184	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
	4185	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
	4186	!
	4187
	4188	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
	4189	! -------------------------------------------
	4190
	4191
	4192	! selection des points pour lesquels le schema est actif:
	4193
	4194
	4195
	4196	IF (addtkeoro .EQ. 1 ) THEN
	4197
	4198	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
	4199	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
	4200
	4201	ELSE IF (addtkeoro .EQ. 2) THEN
	4202
	4203
	4204
	4205	IF (smallscales_tkeoro) THEN
	4206	igwd=0
	4207	DO i=1,klon
	4208	itest(i)=0
	4209	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
	4210	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
	4211	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
	4212	IF (zstd(i).GT.1.0) THEN
	4213	itest(i)=1
	4214	igwd=igwd+1
	4215	idx(igwd)=i
	4216	ENDIF
	4217	ENDDO
	4218
	4219	ELSE
	4220
	4221	igwd=0
	4222	DO i=1,klon
	4223	itest(i)=0
	4224	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	4225	itest(i)=1
	4226	igwd=igwd+1
	4227	idx(igwd)=i
	4228	ENDIF
	4229	ENDDO
	4230
	4231	END IF
	4232
	4233
	4234
	4235
[3065]	4236	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
[2897]	4237	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	4238	igwd,idx,itest, &
	4239	t_seri, u_seri, v_seri, &
	4240	zulow, zvlow, zustrdr, zvstrdr, &
	4241	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
	4242
	4243	zustrdr(:)=0.
	4244	zvstrdr(:)=0.
	4245	zulow(:)=0.
	4246	zvlow(:)=0.
	4247
	4248	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
	4249	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
	4250	END IF
	4251
	4252
	4253
	4254	! TKE update from subgrid temperature and wind tendencies
	4255	!----------------------------------------------------------
	4256	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
	4257
	4258
	4259	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pbl_tke)
	4260
	4261
	4262
	4263	ENDIF
	4264	! -----
	4265	!===============================================================
	4266
	4267
	4268
[2469]	4269	!====================================================================
	4270	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
	4271	!====================================================================
	4272	! Abderrahmane 24.08.09
	4273
	4274	IF (ok_cosp) THEN
	4275	! adeclarer
[1279]	4276	#ifdef CPP_COSP
[3065]	4277	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[1279]	4278
[2469]	4279	IF (prt_level .GE.10) THEN
	4280	print*,'freq_cosp',freq_cosp
	4281	ENDIF
	4282	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
	4283	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
	4284	! s ref_liq,ref_ice
[3065]	4285	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
[2469]	4286	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
[2794]	4287	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
[2469]	4288	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
	4289	JrNt,ref_liq,ref_ice, &
	4290	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
	4291	zu10m,zv10m,pphis, &
	4292	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
	4293	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
	4294	prfl(:,1:klev),psfl(:,1:klev), &
	4295	pmflxr(:,1:klev),pmflxs(:,1:klev), &
	4296	mr_ozone,cldtau, cldemi)
[1412]	4297
[2469]	4298	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
	4299	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
	4300	! M clMISR,
	4301	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
	4302	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
[1279]	4303
[2469]	4304	ENDIF
[1279]	4305
	4306	#endif
[2469]	4307	ENDIF !ok_cosp
[2580]	4308
	4309
	4310	! Marine
	4311
	4312	IF (ok_airs) then
	4313
[3065]	4314	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
[2692]	4315	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
	4316	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
	4317	& map_prop_hc,map_prop_hist,&
	4318	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
	4319	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
	4320	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
	4321	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
	4322	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
	4323	& map_ntot,map_hc,map_hist,&
	4324	& map_Cb,map_ThCi,map_Anv,&
	4325	& alt_tropo )
[2580]	4326	ENDIF
	4327
	4328	ENDIF ! ok_airs
	4329
	4330
[2469]	4331	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4332	!AA
	4333	!AA Installation de l'interface online-offline pour traceurs
	4334	!AA
	4335	!====================================================================
	4336	! Calcul des tendances traceurs
	4337	!====================================================================
	4338	!
[959]	4339
[2469]	4340	IF (type_trac=='repr') THEN
	4341	sh_in(:,:) = q_seri(:,:)
	4342	ELSE
	4343	sh_in(:,:) = qx(:,:,ivap)
[2784]	4344	ch_in(:,:) = qx(:,:,iliq)
[2692]	4345	ENDIF
[1565]	4346
[2973]	4347	IF (iflag_phytrac == 1 ) THEN
	4348
[2630]	4349	#ifdef CPP_Dust
	4350	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
	4351	pdtphys,ftsol, & ! I
	4352	t,q_seri,paprs,pplay,RHcl, & ! I
	4353	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
	4354	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
	4355	u_seri, v_seri, latitude_deg, longitude_deg, &
	4356	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
	4357	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
	4358	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
	4359	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
	4360	fm_therm, entr_therm, rneb, & ! I
	4361	beta_prec_fisrt,beta_prec, & !I
	4362	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
	4363	d_tr_dyn,tr_seri)
	4364
	4365	#else
	4366
[2692]	4367	CALL phytrac ( &
[2469]	4368	itap, days_elapsed+1, jH_cur, debut, &
[3065]	4369	lafin, phys_tstep, u, v, t, &
[2469]	4370	paprs, pplay, pmfu, pmfd, &
	4371	pen_u, pde_u, pen_d, pde_d, &
	4372	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
	4373	u1, v1, ftsol, pctsrf, &
	4374	zustar, zu10m, zv10m, &
	4375	wstar(:,is_ave), ale_bl, ale_wake, &
	4376	latitude_deg, longitude_deg, &
	4377	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
	4378	presnivs, pphis, pphi, albsol1, &
[2784]	4379	sh_in, ch_in, rhcl, cldfra, rneb, &
[2469]	4380	diafra, cldliq, itop_con, ibas_con, &
	4381	pmflxr, pmflxs, prfl, psfl, &
	4382	da, phi, mp, upwd, &
	4383	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
	4384	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
	4385	ev, ep, epmlmMm, eplaMm, & !<<RomP
	4386	dnwd, aerosol_couple, flxmass_w, &
	4387	tau_aero, piz_aero, cg_aero, ccm, &
	4388	rfname, &
	4389	d_tr_dyn, & !<<RomP
	4390	tr_seri)
[2630]	4391	#endif
[2973]	4392	ENDIF ! (iflag_phytrac=1)
[524]	4393
[2469]	4394	IF (offline) THEN
[524]	4395
[2469]	4396	IF (prt_level.ge.9) &
	4397	print*,'Attention on met a 0 les thermiques pour phystoke'
[2692]	4398	CALL phystokenc ( &
[2469]	4399	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
	4400	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	4401	fm_therm,entr_therm, &
	4402	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
	4403	frac_impa, frac_nucl, &
[3065]	4404	pphis,cell_area,phys_tstep,itap, &
[2469]	4405	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
[524]	4406
	4407
[2469]	4408	ENDIF
[524]	4409
[2469]	4410	!
	4411	! Calculer le transport de l'eau et de l'energie (diagnostique)
	4412	!
	4413	CALL transp (paprs,zxtsol, &
	4414	t_seri, q_seri, u_seri, v_seri, zphi, &
	4415	ve, vq, ue, uq)
	4416	!
	4417	!IM global posePB BEG
	4418	IF(1.EQ.0) THEN
	4419	!
	4420	CALL transp_lay (paprs,zxtsol, &
	4421	t_seri, q_seri, u_seri, v_seri, zphi, &
	4422	ve_lay, vq_lay, ue_lay, uq_lay)
	4423	!
	4424	ENDIF !(1.EQ.0) THEN
	4425	!IM global posePB END
	4426	! Accumuler les variables a stocker dans les fichiers histoire:
	4427	!
[1279]	4428
[2469]	4429	!================================================================
	4430	! Conversion of kinetic and potential energy into heat, for
	4431	! parameterisation of subgrid-scale motions
	4432	!================================================================
[1753]	4433
[2469]	4434	d_t_ec(:,:)=0.
	4435	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
[2851]	4436	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
	4437	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
[2469]	4438	zmasse,exner,d_t_ec)
	4439	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
[1753]	4440
[2469]	4441	!=======================================================================
	4442	! SORTIES
	4443	!=======================================================================
	4444	!
	4445	!IM initialisation + calculs divers diag AMIP2
	4446	!
	4447	include "calcul_divers.h"
	4448	!
	4449	!IM Interpolation sur les niveaux de pression du NMC
	4450	! -------------------------------------------------
	4451	!
	4452	include "calcul_STDlev.h"
	4453	!
	4454	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
	4455	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
	4456	!
[2496]	4457	!cc prw = eau precipitable
	4458	! prlw = colonne eau liquide
	4459	! prlw = colonne eau solide
[2499]	4460	prw(:) = 0.
	4461	prlw(:) = 0.
	4462	prsw(:) = 0.
	4463	DO k = 1, klev
	4464	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
	4465	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
	4466	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
[2469]	4467	ENDDO
	4468	!
	4469	IF (type_trac == 'inca') THEN
[655]	4470	#ifdef INCA
[2469]	4471	CALL VTe(VTphysiq)
	4472	CALL VTb(VTinca)
[959]	4473
[2469]	4474	CALL chemhook_end ( &
[3065]	4475	phys_tstep, &
[2469]	4476	pplay, &
	4477	t_seri, &
	4478	tr_seri, &
	4479	nbtr, &
	4480	paprs, &
	4481	q_seri, &
	4482	cell_area, &
	4483	pphi, &
	4484	pphis, &
[2832]	4485	zx_rh, &
	4486	aps, bps)
[959]	4487
[2469]	4488	CALL VTe(VTinca)
	4489	CALL VTb(VTphysiq)
[655]	4490	#endif
[2692]	4491	ENDIF
[655]	4492
[1753]	4493
[2469]	4494	!
	4495	! Convertir les incrementations en tendances
	4496	!
	4497	IF (prt_level .GE.10) THEN
	4498	print *,'Convertir les incrementations en tendances '
	4499	ENDIF
	4500	!
[2692]	4501	IF (mydebug) THEN
	4502	CALL writefield_phy('u_seri',u_seri,nbp_lev)
	4503	CALL writefield_phy('v_seri',v_seri,nbp_lev)
	4504	CALL writefield_phy('t_seri',t_seri,nbp_lev)
	4505	CALL writefield_phy('q_seri',q_seri,nbp_lev)
	4506	ENDIF
[766]	4507
[2469]	4508	DO k = 1, klev
	4509	DO i = 1, klon
[3065]	4510	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
	4511	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
	4512	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
	4513	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
	4514	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
[2469]	4515	!CR: on ajoute le contenu en glace
[2692]	4516	IF (nqo.eq.3) THEN
[3065]	4517	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
[2692]	4518	ENDIF
[2469]	4519	ENDDO
	4520	ENDDO
	4521	!
	4522	!CR: nb de traceurs eau: nqo
	4523	! IF (nqtot.GE.3) THEN
	4524	IF (nqtot.GE.(nqo+1)) THEN
	4525	! DO iq = 3, nqtot
	4526	DO iq = nqo+1, nqtot
	4527	DO k = 1, klev
	4528	DO i = 1, klon
[3065]	4529	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
	4530	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
[2469]	4531	ENDDO
	4532	ENDDO
	4533	ENDDO
	4534	ENDIF
	4535	!
	4536	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	4537	!IM global posePB include "write_bilKP_ins.h"
	4538	!IM global posePB include "write_bilKP_ave.h"
	4539	!
[1412]	4540
[2489]	4541	!--OB mass fixer
	4542	!--profile is corrected to force mass conservation of water
	4543	IF (mass_fixer) THEN
	4544	qql2(:)=0.0
[2499]	4545	DO k = 1, klev
	4546	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2489]	4547	ENDDO
	4548	DO i = 1, klon
	4549	!--compute ratio of what q+ql should be with conservation to what it is
	4550	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
	4551	DO k = 1, klev
	4552	q_seri(i,k) =q_seri(i,k)*corrqql
	4553	ql_seri(i,k)=ql_seri(i,k)*corrqql
	4554	ENDDO
	4555	ENDDO
	4556	ENDIF
	4557	!--fin mass fixer
	4558
[2469]	4559	! Sauvegarder les valeurs de t et q a la fin de la physique:
	4560	!
[2499]	4561	u_ancien(:,:) = u_seri(:,:)
	4562	v_ancien(:,:) = v_seri(:,:)
	4563	t_ancien(:,:) = t_seri(:,:)
	4564	q_ancien(:,:) = q_seri(:,:)
	4565	ql_ancien(:,:) = ql_seri(:,:)
	4566	qs_ancien(:,:) = qs_seri(:,:)
	4567	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
	4568	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
	4569	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
[2469]	4570	! !! RomP >>>
	4571	!CR: nb de traceurs eau: nqo
[2499]	4572	IF (nqtot.GT.nqo) THEN
[2469]	4573	DO iq = nqo+1, nqtot
[2499]	4574	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
[2469]	4575	ENDDO
	4576	ENDIF
	4577	! !! RomP <<<
	4578	!==========================================================================
	4579	! Sorties des tendances pour un point particulier
	4580	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	4581	! pour le debug
	4582	! La valeur de igout est attribuee plus haut dans le programme
	4583	!==========================================================================
[879]	4584
[2692]	4585	IF (prt_level.ge.1) THEN
[2469]	4586	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	4587	write(lunout,*) &
	4588	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
	4589	write(lunout,*) &
	4590	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
	4591	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
	4592	pctsrf(igout,is_sic)
	4593	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
[2692]	4594	DO k=1,klev
[2469]	4595	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
	4596	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
	4597	d_t_eva(igout,k)
[2692]	4598	ENDDO
[2469]	4599	write(lunout,*) 'cool,heat'
[2692]	4600	DO k=1,klev
[2469]	4601	write(lunout,*) cool(igout,k),heat(igout,k)
[2692]	4602	ENDDO
[879]	4603
[2469]	4604	!jyg< (En attendant de statuer sur le sort de d_t_oli)
	4605	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	4606	!jyg! do k=1,klev
	4607	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
	4608	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	4609	!jyg! enddo
	4610	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
[2692]	4611	DO k=1,klev
[2469]	4612	write(lunout,*) d_t_vdf(igout,k), &
	4613	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
[2692]	4614	ENDDO
[2469]	4615	!>jyg
[879]	4616
[2469]	4617	write(lunout,*) 'd_ps ',d_ps(igout)
	4618	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
[2692]	4619	DO k=1,klev
[2469]	4620	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
	4621	d_qx(igout,k,1),d_qx(igout,k,2)
[2692]	4622	ENDDO
	4623	ENDIF
[879]	4624
[2469]	4625	!============================================================
	4626	! Calcul de la temperature potentielle
	4627	!============================================================
	4628	DO k = 1, klev
	4629	DO i = 1, klon
	4630	!JYG/IM theta en debut du pas de temps
	4631	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4632	!JYG/IM theta en fin de pas de temps de physique
	4633	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4634	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
	4635	! MPL 20130625
	4636	! fth_fonctions.F90 et parkind1.F90
	4637	! sinon thetal=theta
	4638	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
	4639	! : ql_seri(i,k))
	4640	thetal(i,k)=theta(i,k)
	4641	ENDDO
	4642	ENDDO
	4643	!
[879]	4644
[2469]	4645	! 22.03.04 BEG
	4646	!=============================================================
	4647	! Ecriture des sorties
	4648	!=============================================================
[524]	4649	#ifdef CPP_IOIPSL
	4650
[2469]	4651	! Recupere des varibles calcule dans differents modules
	4652	! pour ecriture dans histxxx.nc
[782]	4653
[2469]	4654	! Get some variables from module fonte_neige_mod
	4655	CALL fonte_neige_get_vars(pctsrf, &
[2517]	4656	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
[782]	4657
[1507]	4658
[2469]	4659	!=============================================================
	4660	! Separation entre thermiques et non thermiques dans les sorties
	4661	! de fisrtilp
	4662	!=============================================================
[1507]	4663
[2692]	4664	IF (iflag_thermals>=1) THEN
[2469]	4665	d_t_lscth=0.
	4666	d_t_lscst=0.
	4667	d_q_lscth=0.
	4668	d_q_lscst=0.
[2692]	4669	DO k=1,klev
	4670	DO i=1,klon
	4671	IF (ptconvth(i,k)) THEN
[2469]	4672	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4673	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]	4674	ELSE
[2469]	4675	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4676	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]	4677	ENDIF
	4678	ENDDO
	4679	ENDDO
[1507]	4680
[2692]	4681	DO i=1,klon
[2469]	4682	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
	4683	plul_th(i)=prfl(i,1)+psfl(i,1)
[2692]	4684	ENDDO
	4685	ENDIF
[909]	4686
[2469]	4687	!On effectue les sorties:
[1791]	4688
[2630]	4689	#ifdef CPP_Dust
	4690	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
	4691	pplay, lmax_th, aerosol_couple, &
	4692	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
	4693	ptconv, read_climoz, clevSTD, &
	4694	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
	4695	flag_aerosol, flag_aerosol_strat, ok_cdnc)
	4696	#else
[2469]	4697	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
	4698	pplay, lmax_th, aerosol_couple, &
[2496]	4699	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
	4700	ok_sync, ptconv, read_climoz, clevSTD, &
[2665]	4701	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
[2469]	4702	flag_aerosol, flag_aerosol_strat, ok_cdnc)
[2630]	4703	#endif
[1791]	4704
[2651]	4705	#ifndef CPP_XIOS
[2590]	4706	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
[2651]	4707	#endif
[687]	4708
[524]	4709	#endif
	4710
[2235]	4711
[2469]	4712	!====================================================================
	4713	! Arret du modele apres hgardfou en cas de detection d'un
	4714	! plantage par hgardfou
	4715	!====================================================================
[2235]	4716
	4717	IF (abortphy==1) THEN
	4718	abort_message ='Plantage hgardfou'
[2311]	4719	CALL abort_physic (modname,abort_message,1)
[2235]	4720	ENDIF
	4721
[2469]	4722	! 22.03.04 END
	4723	!
	4724	!====================================================================
	4725	! Si c'est la fin, il faut conserver l'etat de redemarrage
	4726	!====================================================================
	4727	!
[782]	4728
[2469]	4729	IF (lafin) THEN
	4730	itau_phy = itau_phy + itap
	4731	CALL phyredem ("restartphy.nc")
	4732	! open(97,form="unformatted",file="finbin")
	4733	! write(97) u_seri,v_seri,t_seri,q_seri
	4734	! close(97)
[3406]	4735	! !$OMP MASTER
[2692]	4736	IF (read_climoz >= 1) THEN
	4737	IF (is_mpi_root) THEN
	4738	CALL nf95_close(ncid_climoz)
	4739	ENDIF
[2788]	4740	DEALLOCATE(press_edg_climoz) ! pointer
	4741	DEALLOCATE(press_cen_climoz) ! pointer
[2692]	4742	ENDIF
[3406]	4743	! !$OMP END MASTER
	4744	#ifdef CPP_XIOS
	4745	IF (is_omp_master) CALL xios_context_finalize
	4746	#endif
[2469]	4747	ENDIF
[1863]	4748
[2469]	4749	! first=.false.
[1863]	4750
[2418]	4751
[2469]	4752	END SUBROUTINE physiq
[2418]	4753
[2902]	4754	END MODULE physiq_mod

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

Download in other formats:

Original Format