Context Navigation

physiq_mod.F90 @ 5452

Last change on this file since 5452 was 5452, checked in by aborella, 39 hours ago
First implementation of the contrails parameterisation Lacks the emission of H2O + the impact on radiative transfer
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: 216.5 KB

Line
1
2	! $Id: physiq_mod.F90 5452 2024-12-23 17:12:29Z aborella $
3	!
4	!#define IO_DEBUG
5	MODULE physiq_mod
6
7	IMPLICIT NONE
8
9	CONTAINS
10
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)
17
18	! For clarity, the "USE" section is now arranged in alphabetical order,
19	! with a separate section for CPP keys
20	! PLEASE try to follow this rule
21
22	USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
23	USE aero_mod
24	USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
25	& fl_ebil, fl_cor_ebil
26	USE assert_m, only: assert
27	USE change_srf_frac_mod
28	USE conf_phys_m, only: conf_phys
29	USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
30	USE CFMIP_point_locations ! IM stations CFMIP
31	USE cmp_seri_mod
32	USE dimphy
33	USE etat0_limit_unstruct_mod
34	USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
35	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
36	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
37	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
38	histwrite, ju2ymds, ymds2ju, getin
39	USE ioipsl_getin_p_mod, ONLY : getin_p
40	USE indice_sol_mod
41	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, addPhase
42	USE strings_mod, ONLY: strIdx
43	USE iophy
44	USE limit_read_mod, ONLY : init_limit_read
45	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
46	USE mod_phys_lmdz_mpi_data, only: is_mpi_root
47	USE mod_phys_lmdz_para
48	USE netcdf95, only: nf95_close
49	USE netcdf, only: nf90_fill_real ! IM for NMC files
50	USE open_climoz_m, only: open_climoz ! ozone climatology from a file
51	USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
52	USE pbl_surface_mod, ONLY : pbl_surface
53	USE phyaqua_mod, only: zenang_an
54	USE phyetat0_mod, only: phyetat0
55	USE phystokenc_mod, ONLY: offline, phystokenc
56	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
57	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend
58	!! USE phys_local_var_mod, ONLY : a long list of variables
59	!! ==> see below, after "CPP Keys" section
60	USE phys_state_var_mod ! Variables sauvegardees de la physique
61	USE phys_output_mod
62	USE phys_output_ctrlout_mod
63	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level, &
64	alert_first_call, call_alert, prt_alerte
65	USE readaerosol_mod, ONLY : init_aero_fromfile
66	USE readaerosolstrato_m, ONLY : init_readaerosolstrato
67	USE radlwsw_m, only: radlwsw
68	USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
69	USE regr_pr_time_av_m, only: regr_pr_time_av
70	USE surface_data, ONLY : type_ocean, ok_veget
71	USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time
72	USE tracinca_mod, ONLY: config_inca
73	USE tropopause_m, ONLY: dyn_tropopause
74	USE vampir
75	USE write_field_phy
76	use wxios_mod, ONLY: g_ctx, wxios_set_context
77	USE lmdz_lscp, ONLY : lscp
78	USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
79	USE lmdz_lscp_old, ONLY : fisrtilp
80	USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
81	USE lmdz_wake_ini, ONLY : wake_ini
82	USE lmdz_surf_wind_ini, ONLY : surf_wind_ini, iflag_surf_wind
83	USE lmdz_surf_wind, ONLY : surf_wind
84	USE yamada_ini_mod, ONLY : yamada_ini
85	USE lmdz_atke_turbulence_ini, ONLY : atke_ini
86	USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
87	USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
88	USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
89	USE lmdz_lscp_ini, ONLY : lscp_ini
90	USE lmdz_ratqs_main, ONLY : ratqs_main
91	USE lmdz_ratqs_ini, ONLY : ratqs_ini
92	USE lmdz_cloud_optics_prop_ini, ONLY : cloud_optics_prop_ini
93	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
94	USE phys_output_var_mod, ONLY : cloud_cover_sw, cloud_cover_sw_s2
95
96
97	!USE cmp_seri_mod
98	! USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
99	! & fl_ebil, fl_cor_ebil
100
101
102	USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
103	USE phys_output_write_spl_mod, ONLY: phys_output_write_spl
104	USE phytrac_mod, ONLY: phytrac_init, phytrac
105	USE phys_output_write_mod, ONLY: phys_output_write
106
107
108	USE geometry_mod, ONLY: longitude, latitude, boundslon, boundslat, ind_cell_glo
109	USE time_phylmdz_mod, ONLY: ndays
110	USE infotrac_phy, ONLY: nqCO2
111	USE lmdz_reprobus_wrappers, ONLY: init_chem_rep_xjour, d_q_rep, d_ql_rep, d_qi_rep, &
112	ptrop, ttrop, ztrop, gravit, itroprep, Z1, Z2, fac, B, chemini_rep, chemtime_rep, coord_hyb_rep, &
113	rtime
114	USE strataer_local_var_mod
115	USE strataer_emiss_mod, ONLY: strataer_emiss_init
116	USE time_phylmdz_mod, ONLY: annee_ref, day_ini, day_ref, start_time
117	USE vertical_layers_mod, ONLY: aps, bps, ap, bp
118	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA, CPPKEY_REPROBUS, CPPKEY_STRATAER, CPPKEY_DUST, CPPKEY_COSP, CPPKEY_COSP2, CPPKEY_COSPV2
119
120	#ifdef CPP_RRTM
121	USE YOERAD, ONLY : NRADLP
122	! USE YOESW, ONLY : RSUN
123	#endif
124
125	USE phys_local_var_mod, ONLY: d_q_emiss
126	USE strataer_local_var_mod
127	USE strataer_nuc_mod, ONLY: strataer_nuc_init
128	USE strataer_emiss_mod, ONLY: strataer_emiss_init
129
130	USE lmdz_xios, ONLY: xios_update_calendar, xios_context_finalize
131	USE lmdz_xios, ONLY: xios_get_field_attr, xios_field_is_active, xios_context
132	USE lmdz_xios, ONLY: xios_set_current_context
133	use wxios_mod, ONLY: missing_val, using_xios
134
135	#ifndef CPP_XIOS
136	USE paramLMDZ_phy_mod
137	#endif
138	!
139	!
140	!!!!!!!!!!!!!!!!!! END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
141
142	USE physiqex_mod, ONLY : physiqex
143	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
144	! [Variables internes non sauvegardees de la physique]
145	! Variables locales pour effectuer les appels en serie
146	t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
147	u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
148	rhcl, &
149	! Dynamic tendencies (diagnostics)
150	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
151	d_u_dyn,d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
152	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
153	! Physic tendencies
154	d_t_con,d_q_con,d_q_con_zmasse,d_u_con,d_v_con, &
155	d_tr, & !! to be removed?? (jyg)
156	d_t_wake,d_q_wake, &
157	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
158	d_t_ajsb,d_q_ajsb, &
159	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
160	! d_t_ajs_w,d_q_ajs_w, &
161	! d_t_ajs_x,d_q_ajs_x, &
162	!
163	d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
164	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
165	d_t_lscst,d_q_lscst, &
166	d_t_lscth,d_q_lscth, &
167	plul_st,plul_th, &
168	!
169	d_t_vdf,d_q_vdf, d_qbs_vdf, d_u_vdf,d_v_vdf,d_t_diss, &
170	d_t_vdf_x, d_t_vdf_w, &
171	d_q_vdf_x, d_q_vdf_w, &
172	d_ts, &
173	!
174	d_t_bsss,d_q_bsss,d_qbs_bsss, &
175	!
176	! d_t_oli,d_u_oli,d_v_oli, &
177	d_t_oro,d_u_oro,d_v_oro, &
178	d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
179	d_t_lif,d_u_lif,d_v_lif, &
180	d_t_ec, &
181	!
182	du_gwd_hines,dv_gwd_hines,d_t_hin, &
183	dv_gwd_rando,dv_gwd_front, &
184	east_gwstress,west_gwstress, &
185	d_q_ch4, &
186	! proprecip
187	qraindiag, qsnowdiag, &
188	dqreva, dqssub, &
189	dqrauto,dqrcol,dqrmelt,dqrfreez, &
190	dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
191	! Special RRTM
192	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
193	ZFLUP0,ZFSDN0,ZFSUP0, &
194	!
195	topswad_aero,solswad_aero, &
196	topswai_aero,solswai_aero, &
197	topswad0_aero,solswad0_aero, &
198	!LW additional
199	toplwad_aero,sollwad_aero, &
200	toplwai_aero,sollwai_aero, &
201	toplwad0_aero,sollwad0_aero, &
202	!pour Ecrad
203	topswad_aero_s2, solswad_aero_s2, &
204	topswai_aero_s2, solswai_aero_s2, &
205	topswad0_aero_s2, solswad0_aero_s2, &
206	topsw_aero_s2, topsw0_aero_s2, &
207	solsw_aero_s2, solsw0_aero_s2, &
208	topswcf_aero_s2, solswcf_aero_s2, &
209	!LW diagnostics
210	toplwad_aero_s2, sollwad_aero_s2, &
211	toplwai_aero_s2, sollwai_aero_s2, &
212	toplwad0_aero_s2, sollwad0_aero_s2, &
213	!
214	topsw_aero,solsw_aero, &
215	topsw0_aero,solsw0_aero, &
216	topswcf_aero,solswcf_aero, &
217	tausum_aero,tau3d_aero, &
218	drytausum_aero, &
219	!
220	!variables CFMIP2/CMIP5
221	topswad_aerop, solswad_aerop, &
222	topswai_aerop, solswai_aerop, &
223	topswad0_aerop, solswad0_aerop, &
224	topsw_aerop, topsw0_aerop, &
225	solsw_aerop, solsw0_aerop, &
226	topswcf_aerop, solswcf_aerop, &
227	!LW diagnostics
228	toplwad_aerop, sollwad_aerop, &
229	toplwai_aerop, sollwai_aerop, &
230	toplwad0_aerop, sollwad0_aerop, &
231	!pour Ecrad
232	topswad_aero_s2, solswad_aero_s2, &
233	topswai_aero_s2, solswai_aero_s2, &
234	topswad0_aero_s2, solswad0_aero_s2, &
235	topsw_aero_s2, topsw0_aero_s2, &
236	solsw_aero_s2, solsw0_aero_s2, &
237	topswcf_aero_s2, solswcf_aero_s2, &
238	!LW diagnostics
239	toplwad_aero_s2, sollwad_aero_s2, &
240	toplwai_aero_s2, sollwai_aero_s2, &
241	toplwad0_aero_s2, sollwad0_aero_s2, &
242	!
243	ptstar, pt0, slp, &
244	!
245	bils, &
246	!
247	cldh, cldl,cldm, cldq, cldt, &
248	JrNt, &
249	dthmin, evap, snowerosion,fder, plcl, plfc, &
250	prw, prlw, prsw, prbsw, water_budget, &
251	s_lcl, s_pblh, s_pblt, s_therm, &
252	cdragm, cdragh, &
253	zustar, zu10m, zv10m, rh2m, qsat2m, &
254	zq2m, zt2m, zn2mout, weak_inversion, &
255	!
256	s_pblh_x, s_pblh_w, &
257	s_lcl_x, s_lcl_w, &
258	!
259	slab_wfbils, tpot, tpote, &
260	ue, uq, ve, vq, zxffonte, &
261	uwat, vwat, &
262	zxfqcalving, zxfluxlat, &
263	zxrunofflic, &
264	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
265	delta_qsurf, &
266	rain_lsc, rain_num, &
267	!
268	sens_x, sens_w, &
269	zxfluxlat_x, zxfluxlat_w, &
270	!
271	pbl_tke_input, pbl_eps, l_mix, wprime,&
272	t_therm, q_therm, u_therm, v_therm, &
273	cdragh_x, cdragh_w, &
274	cdragm_x, cdragm_w, &
275	kh, kh_x, kh_w, &
276	!
277	wake_k, &
278	alp_wake, &
279	wake_h, wake_omg, &
280	! tendencies of delta T and delta q:
281	d_deltat_wk, d_deltaq_wk, & ! due to wakes
282	d_deltat_wk_gw, d_deltaq_wk_gw, & ! due to wake induced gravity waves
283	d_deltat_vdf, d_deltaq_vdf, & ! due to vertical diffusion
284	d_deltat_the, d_deltaq_the, & ! due to thermals
285	d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
286	! tendencies of wake fractional area and wake number per unit area:
287	d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, & ! due to wakes
288	!!! d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
289	!!! d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
290	!
291	ptconv, ratqsc, &
292	wbeff, convoccur, zmax_th, &
293	sens, flwp, fiwp, &
294	alp_bl_conv,alp_bl_det, &
295	alp_bl_fluct_m,alp_bl_fluct_tke, &
296	alp_bl_stat, n2, s2, strig, zcong, zlcl_th, &
297	proba_notrig, random_notrig, &
298	!! cv_gen, & !moved to phys_state_var_mod
299	!
300	dnwd0, &
301	omega, &
302	epmax_diag, &
303	! Deep convective variables used in phytrac
304	pmflxr, pmflxs, &
305	wdtrainA, wdtrainS, wdtrainM, wdtrainAS, &
306	upwd, dnwd, &
307	ep, &
308	da, mp, &
309	phi, &
310	wght_cvfd, &
311	phi2, &
312	d1a, dam, &
313	ev, &
314	elij, &
315	qtaa, &
316	clw, &
317	epmlmMm, eplaMm, &
318	sij, &
319	!
320	rneblsvol, &
321	pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
322	distcltop, temp_cltop, &
323	!-- LSCP - condensation and ice supersaturation variables
324	qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
325	ql_seri_lscp, ratio_ql_qtot, qi_seri_lscp, ratio_qi_qtot, &
326	dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
327	dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
328	!-- LSCP - aviation and contrails variables
329	rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
330	Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
331	dcf_avi, dqi_avi, dqvc_avi, &
332	!
333	cldemi, &
334	cldfra, cldtau, fiwc, &
335	fl, re, flwc, &
336	ref_liq, ref_ice, theta, &
337	ref_liq_pi, ref_ice_pi, &
338	zphi, zx_rh, zx_rhl, zx_rhi, &
339	pmfd, pmfu, &
340	!
341	t2m, fluxlat, &
342	fsollw, evap_pot, &
343	fsolsw, wfbils, wfevap, &
344	prfl, psfl,bsfl, fraca, Vprecip, &
345	zw2, &
346	!
347	fluxu, fluxv, &
348	fluxt, &
349	!
350	uwriteSTD, vwriteSTD, & !pour calcul_STDlev.h
351	wwriteSTD, phiwriteSTD, & !pour calcul_STDlev.h
352	qwriteSTD, twriteSTD, rhwriteSTD, & !pour calcul_STDlev.h
353	!
354	beta_prec, &
355	rneb, &
356	zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
357	zxfluxt,zxfluxq
358	!
359	USE phys_local_var_mod, ONLY: zfice, dNovrN, ptconv
360	USE phys_output_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
361	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra
362	USE output_physiqex_mod, ONLY: output_physiqex
363	USE yomcst_mod_h
364	USE clesphys_mod_h
365	USE conema3_mod_h
366	USE alpale_mod
367	USE yoethf_mod_h
368	USE calcul_divers_mod_h, ONLY: calcul_divers
369	USE compbl_mod_h
370	USE nuage_params_mod_h
371	USE dimpft_mod_h, ONLY: nvm_lmdz
372	USE radepsi_mod_h
373	USE radopt_mod_h
374	USE regdim_mod_h
375	USE phys_output_write_spl_mod, ONLY: phys_output_write_spl
376	USE phytracr_spl_mod, ONLY: phytracr_spl_out_init, phytracr_spl
377	IMPLICIT NONE
378	!>======================================================================
379	!!
380	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
381	!!
382	!! Objet: Moniteur general de la physique du modele
383	!!AA Modifications quant aux traceurs :
384	!!AA - uniformisation des parametrisations ds phytrac
385	!!AA - stockage des moyennes des champs necessaires
386	!!AA en mode traceur off-line
387	!!======================================================================
388	!! CLEFS CPP POUR LES IO
389	!! =====================
390	#define histNMC
391	!!======================================================================
392	!! modif ( P. Le Van , 12/10/98 )
393	!!
394	!! Arguments:
395	!!
396	!! nlon----input-I-nombre de points horizontaux
397	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
398	!! debut---input-L-variable logique indiquant le premier passage
399	!! lafin---input-L-variable logique indiquant le dernier passage
400	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
401	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
402	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
403	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
404	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
405	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
406	!! pphis---input-R-geopotentiel du sol
407	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
408	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
409	!! v-------input-R-vitesse Y (de S a N) en m/s
410	!! t-------input-R-temperature (K)
411	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
412	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
413	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
414	!! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
415	!! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
416	!! flxmass_w -input-R- flux de masse verticale
417	!! d_u-----output-R-tendance physique de "u" (m/s/s)
418	!! d_v-----output-R-tendance physique de "v" (m/s/s)
419	!! d_t-----output-R-tendance physique de "t" (K/s)
420	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
421	!! d_ps----output-R-tendance physique de la pression au sol
422	!!======================================================================
423	integer jjmp1
424	! parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
425	! integer iip1
426	! parameter (iip1=iim+1)
427
428	!======================================================================
429	LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
430	!$OMP THREADPRIVATE(ok_volcan)
431	INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato
432	!$OMP THREADPRIVATE(flag_volc_surfstrat)
433	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
434	PARAMETER (ok_cvl=.TRUE.)
435	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
436	PARAMETER (ok_gust=.FALSE.)
437	INTEGER, SAVE :: iflag_radia ! active ou non le rayonnement (MPL)
438	!$OMP THREADPRIVATE(iflag_radia)
439	!======================================================================
440	LOGICAL check ! Verifier la conservation du modele en eau
441	PARAMETER (check=.FALSE.)
442	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
443	PARAMETER (ok_stratus=.FALSE.)
444	!======================================================================
445	REAL amn, amx
446	INTEGER igout
447	!======================================================================
448	! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
449	! en attente du codage des cles par Fred
450	! iflag_cycle_diurne est initialise par conf_phys et se trouve
451	! dans clesphys.h (IM)
452	!======================================================================
453	! Modele thermique du sol, a activer pour le cycle diurne:
454	!cc LOGICAL soil_model
455	!cc PARAMETER (soil_model=.FALSE.)
456	!======================================================================
457	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
458	! le calcul du rayonnement est celle apres la precipitation des nuages.
459	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
460	! la condensation et la precipitation. Cette cle augmente les impacts
461	! radiatifs des nuages.
462	!cc LOGICAL new_oliq
463	!cc PARAMETER (new_oliq=.FALSE.)
464	!======================================================================
465	! Clefs controlant deux parametrisations de l'orographie:
466	!c LOGICAL ok_orodr
467	!cc PARAMETER (ok_orodr=.FALSE.)
468	!cc LOGICAL ok_orolf
469	!cc PARAMETER (ok_orolf=.FALSE.)
470	!======================================================================
471	LOGICAL ok_journe ! sortir le fichier journalier
472	SAVE ok_journe
473	!$OMP THREADPRIVATE(ok_journe)
474	!
475	LOGICAL ok_mensuel ! sortir le fichier mensuel
476	SAVE ok_mensuel
477	!$OMP THREADPRIVATE(ok_mensuel)
478	!
479	LOGICAL ok_instan ! sortir le fichier instantane
480	SAVE ok_instan
481	!$OMP THREADPRIVATE(ok_instan)
482	!
483	LOGICAL ok_LES ! sortir le fichier LES
484	SAVE ok_LES
485	!$OMP THREADPRIVATE(ok_LES)
486	!
487	LOGICAL callstats ! sortir le fichier stats
488	SAVE callstats
489	!$OMP THREADPRIVATE(callstats)
490	!
491	LOGICAL ok_region ! sortir le fichier regional
492	PARAMETER (ok_region=.FALSE.)
493	!======================================================================
494	REAL seuil_inversion
495	SAVE seuil_inversion
496	!$OMP THREADPRIVATE(seuil_inversion)
497
498
499
500	real facteur
501
502	REAL wmax_th(klon)
503	REAL tau_overturning_th(klon)
504
505	INTEGER lmax_th(klon)
506	INTEGER limbas(klon)
507	REAL ratqscth(klon,klev)
508	REAL ratqsdiff(klon,klev)
509	REAL zqsatth(klon,klev)
510
511	!======================================================================
512	!
513	! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
514	INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc, ircont
515	!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc, ircont)
516	!
517	!
518	! Variables argument:
519	!
520	INTEGER nlon
521	INTEGER nlev
522	REAL,INTENT(IN) :: pdtphys_
523	! NB: pdtphys to be used in physics is in time_phylmdz_mod
524	LOGICAL debut, lafin
525	REAL paprs(klon,klev+1)
526	REAL pplay(klon,klev)
527	REAL pphi(klon,klev)
528	REAL pphis(klon)
529	REAL presnivs(klev)
530	!JLD REAL znivsig(klev)
531	!JLD real pir
532
533	REAL u(klon,klev)
534	REAL v(klon,klev)
535
536	REAL, intent(in):: rot(klon, klev)
537	! relative vorticity, in s-1, needed for frontal waves
538
539	REAL t(klon,klev),thetal(klon,klev)
540	! thetal: ligne suivante a decommenter si vous avez les fichiers
541	! MPL 20130625
542	! fth_fonctions.F90 et parkind1.F90
543	! sinon thetal=theta
544	! REAL fth_thetae,fth_thetav,fth_thetal
545	REAL qx(klon,klev,nqtot)
546	REAL flxmass_w(klon,klev)
547	REAL d_u(klon,klev)
548	REAL d_v(klon,klev)
549	REAL d_t(klon,klev)
550	REAL d_qx(klon,klev,nqtot)
551	REAL d_ps(klon)
552	! variables pour tend_to_tke
553	REAL duadd(klon,klev)
554	REAL dvadd(klon,klev)
555	REAL dtadd(klon,klev)
556
557	!! Variables moved to phys_local_var_mod
558	!! ! Variables pour le transport convectif
559	!! real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
560	!! real wght_cvfd(klon,klev)
561	!! ! Variables pour le lessivage convectif
562	!! ! RomP >>>
563	!! real phi2(klon,klev,klev)
564	!! real d1a(klon,klev),dam(klon,klev)
565	!! real ev(klon,klev)
566	!! real clw(klon,klev),elij(klon,klev,klev)
567	!! real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
568	!! ! RomP <<<
569	!IM definition dynamique o_trac dans phys_output_open
570	! type(ctrl_out) :: o_trac(nqtot)
571
572	! variables a une pression donnee
573	!
574	include "declare_STDlev.h"
575
576	INTEGER n
577	!ym INTEGER npoints
578	!ym PARAMETER(npoints=klon)
579	!
580	INTEGER nregISCtot
581	PARAMETER(nregISCtot=1)
582	!
583	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
584	! sur 1 region rectangulaire y compris pour 1 point
585	! imin_debut : indice minimum de i; nbpti : nombre de points en
586	! direction i (longitude)
587	! jmin_debut : indice minimum de j; nbptj : nombre de points en
588	! direction j (latitude)
589	!JLD INTEGER imin_debut, nbpti
590	!JLD INTEGER jmin_debut, nbptj
591	!IM: region='3d' <==> sorties en global
592	CHARACTER*3 region
593	PARAMETER(region='3d')
594	LOGICAL ok_hf
595	!
596	SAVE ok_hf
597	!$OMP THREADPRIVATE(ok_hf)
598
599	INTEGER, PARAMETER :: longcles=20
600	REAL, SAVE :: clesphy0(longcles)
601	!$OMP THREADPRIVATE(clesphy0)
602	!
603	! Variables propres a la physique
604	INTEGER, SAVE :: itap ! compteur pour la physique
605	!$OMP THREADPRIVATE(itap)
606
607	INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys
608	!$OMP THREADPRIVATE(abortphy)
609	!
610	REAL,SAVE :: solarlong0
611	!$OMP THREADPRIVATE(solarlong0)
612
613	!
614	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
615	!
616	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
617	REAL zulow(klon),zvlow(klon)
618	!
619	INTEGER igwd,idx(klon),itest(klon)
620	!
621	! REAL,allocatable,save :: run_off_lic_0(:)
622	! !$OMP THREADPRIVATE(run_off_lic_0)
623	!ym SAVE run_off_lic_0
624	!KE43
625	! Variables liees a la convection de K. Emanuel (sb):
626	!
627	REAL, SAVE :: bas, top ! cloud base and top levels
628	!$OMP THREADPRIVATE(bas, top)
629	!------------------------------------------------------------------
630	! Upmost level reached by deep convection and related variable (jyg)
631	!
632	! INTEGER izero
633	INTEGER k_upper_cv
634	!------------------------------------------------------------------
635	! Compteur de l'occurence de cvpas=1
636	INTEGER Ncvpaseq1
637	SAVE Ncvpaseq1
638	!$OMP THREADPRIVATE(Ncvpaseq1)
639	!
640	!==========================================================================
641	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema
642	!de convection avec poches froides
643	! Variables li\'ees \`a la poche froide (jyg)
644
645	!! REAL mipsh(klon,klev) ! mass flux shed by the adiab ascent at each level
646	!! Moved to phys_state_var_mod
647	!
648	REAL wape_prescr, fip_prescr
649	INTEGER it_wape_prescr
650	SAVE wape_prescr, fip_prescr, it_wape_prescr
651	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
652	!
653	! variables supplementaires de concvl
654	REAL Tconv(klon,klev)
655	!! variable moved to phys_local_var_mod
656	!! REAL sij(klon,klev,klev)
657	!! !
658	!! ! variables pour tester la conservation de l'energie dans concvl
659	!! REAL, DIMENSION(klon,klev) :: d_t_con_sat
660	!! REAL, DIMENSION(klon,klev) :: d_q_con_sat
661	!! REAL, DIMENSION(klon,klev) :: dql_sat
662
663	REAL, SAVE :: alp_bl_prescr=0.
664	REAL, SAVE :: ale_bl_prescr=0.
665	REAL, SAVE :: wake_s_min_lsp=0.1
666	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
667	!$OMP THREADPRIVATE(wake_s_min_lsp)
668
669	REAL ok_wk_lsp(klon)
670
671	!RC
672	! Variables li\'ees \`a la poche froide (jyg et rr)
673
674	INTEGER, SAVE :: iflag_wake_tend ! wake: if =0, then wake state variables are
675	! updated within calwake
676	!$OMP THREADPRIVATE(iflag_wake_tend)
677	INTEGER, SAVE :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
678	! power provided by the wakes; else, Alp_wk is the
679	! lifting power conditionned on the presence of a
680	! gust-front in the grid cell.
681	!$OMP THREADPRIVATE(iflag_alp_wk_cond)
682
683	REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
684	REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
685
686	REAL wake_dth(klon,klev) ! wake : temp pot difference
687
688	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta
689	! transported by LS omega
690	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of
691	! large scale omega
692	REAL wake_dtKE(klon,klev) ! Wake : differential heating
693	! (wake - unpertubed) CONV
694	REAL wake_dqKE(klon,klev) ! Wake : differential moistening
695	! (wake - unpertubed) CONV
696	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
697	REAL wake_spread(klon,klev) ! spreading term in wake_delt
698	!
699	!pourquoi y'a pas de save??
700	!
701	!!! INTEGER, SAVE, DIMENSION(klon) :: wake_k
702	!!! !$OMP THREADPRIVATE(wake_k)
703	!
704	!jyg<
705	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
706	!>jyg
707
708	REAL wake_fip_0(klon) ! Average Front Incoming Power (unconditionned)
709	REAL wake_gfl(klon) ! Gust Front Length
710	!!! REAL wake_dens(klon) ! moved to phys_state_var_mod
711	!
712	!
713	REAL dt_dwn(klon,klev)
714	REAL dq_dwn(klon,klev)
715	REAL M_dwn(klon,klev)
716	REAL M_up(klon,klev)
717	REAL dt_a(klon,klev)
718	REAL dq_a(klon,klev)
719	REAL d_t_adjwk(klon,klev) !jyg
720	REAL d_q_adjwk(klon,klev) !jyg
721	LOGICAL,SAVE :: ok_adjwk=.FALSE.
722	!$OMP THREADPRIVATE(ok_adjwk)
723	INTEGER,SAVE :: iflag_adjwk=0 !jyg
724	!$OMP THREADPRIVATE(iflag_adjwk) !jyg
725	REAL,SAVE :: oliqmax=999.,oicemax=999.
726	!$OMP THREADPRIVATE(oliqmax,oicemax)
727	REAL, SAVE :: alp_offset
728	!$OMP THREADPRIVATE(alp_offset)
729	REAL, SAVE :: dtcon_multistep_max=1.e6
730	!$OMP THREADPRIVATE(dtcon_multistep_max)
731	REAL, SAVE :: dqcon_multistep_max=1.e6
732	!$OMP THREADPRIVATE(dqcon_multistep_max)
733
734
735	!
736	!RR:fin declarations poches froides
737	!==========================================================================
738
739	REAL ztv(klon,klev),ztva(klon,klev)
740	REAL zpspsk(klon,klev)
741	REAL ztla(klon,klev),zqla(klon,klev)
742	REAL zthl(klon,klev)
743
744	!cc nrlmd le 10/04/2012
745
746	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
747	!---Propri\'et\'es du thermiques au LCL
748	! real zlcl_th(klon) ! Altitude du LCL calcul\'e
749	! continument (pcon dans
750	! thermcell_main.F90)
751	real fraca0(klon) ! Fraction des thermiques au LCL
752	real w0(klon) ! Vitesse des thermiques au LCL
753	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
754	real tke0(klon,klev+1) ! TKE au d\'ebut du pas de temps
755	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
756	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
757	INTEGER, SAVE :: iflag_thermcell_tke ! transtport TKE by thermals
758	!$OMP THREADPRIVATE(iflag_thermcell_tke)
759
760	!JLD !---D\'eclenchement stochastique
761	!JLD integer :: tau_trig(klon)
762
763	REAL,SAVE :: random_notrig_max=1.
764	!$OMP THREADPRIVATE(random_notrig_max)
765
766	!--------Statistical Boundary Layer Closure: ALP_BL--------
767	!---Profils de TKE dans et hors du thermique
768	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
769	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
770
771	!-------Activer les tendances de TKE due a l'orograp??ie---------
772	INTEGER, SAVE :: addtkeoro
773	!$OMP THREADPRIVATE(addtkeoro)
774	REAL, SAVE :: alphatkeoro
775	!$OMP THREADPRIVATE(alphatkeoro)
776	LOGICAL, SAVE :: smallscales_tkeoro
777	!$OMP THREADPRIVATE(smallscales_tkeoro)
778
779
780
781	!cc fin nrlmd le 10/04/2012
782
783	! Variables locales pour la couche limite (al1):
784	!
785	!Al1 REAL pblh(klon) ! Hauteur de couche limite
786	!Al1 SAVE pblh
787	!34EK
788	!
789	! Variables locales:
790	!
791	!AA
792	!AA Pour phytrac
793	REAL u1(klon) ! vents dans la premiere couche U
794	REAL v1(klon) ! vents dans la premiere couche V
795
796	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
797	!@$$ PARAMETER (offline=.false.)
798	!@$$ INTEGER physid
799	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
800	REAL frac_nucl(klon,klev) ! idem (nucleation)
801	! RomP >>>
802	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
803	! RomP <<<
804
805	!IM cf FH pour Tiedtke 080604
806	REAL rain_tiedtke(klon),snow_tiedtke(klon)
807	!
808	!IM 050204 END
809	REAL devap(klon) ! evaporation et sa derivee
810	REAL dsens(klon) ! chaleur sensible et sa derivee
811
812	!
813	! Conditions aux limites
814	!
815	!
816	REAL :: day_since_equinox
817	! Date de l'equinoxe de printemps
818	INTEGER, parameter :: mth_eq=3, day_eq=21
819	REAL :: jD_eq
820
821	LOGICAL, parameter :: new_orbit = .TRUE.
822
823	!
824	INTEGER lmt_pas
825	SAVE lmt_pas ! frequence de mise a jour
826	!$OMP THREADPRIVATE(lmt_pas)
827	real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
828	! (column-density of mass of air in a cell, in kg m-2)
829	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
830
831	!IM sorties
832	REAL un_jour
833	PARAMETER(un_jour=86400.)
834	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
835	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
836	!$OMP THREADPRIVATE(itapm1)
837	!======================================================================
838	!
839	! Declaration des procedures appelees
840	!
841	EXTERNAL angle ! calculer angle zenithal du soleil
842	EXTERNAL alboc ! calculer l'albedo sur ocean
843	EXTERNAL ajsec ! ajustement sec
844	EXTERNAL conlmd ! convection (schema LMD)
845	EXTERNAL conema3 ! convect4.3
846	EXTERNAL hgardfou ! verifier les temperatures
847	EXTERNAL nuage ! calculer les proprietes radiatives
848	!C EXTERNAL o3cm ! initialiser l'ozone
849	EXTERNAL orbite ! calculer l'orbite terrestre
850	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
851	EXTERNAL suphel ! initialiser certaines constantes
852	EXTERNAL transp ! transport total de l'eau et de l'energie
853	!IM
854	EXTERNAL haut2bas !variables de haut en bas
855	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
856	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
857	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
858	! EXTERNAL moyglo_aire
859	! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
860	! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
861	!
862	!
863	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
864	! Local variables
865	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
866	!
867	! REAL rhcl(klon,klev) ! humiditi relative ciel clair
868	REAL dialiq(klon,klev) ! eau liquide nuageuse
869	REAL diafra(klon,klev) ! fraction nuageuse
870	REAL radocond(klon,klev) ! eau condensee nuageuse
871	!
872	!XXX PB
873	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
874	REAL fluxqbs(klon,klev, nbsrf) ! flux turbulent de neige soufflee
875	!
876	!FC REAL zxfluxt(klon, klev)
877	!FC REAL zxfluxq(klon, klev)
878	REAL zxfluxqbs(klon,klev)
879	REAL zxfluxu(klon, klev)
880	REAL zxfluxv(klon, klev)
881
882	! Le rayonnement n'est pas calcule tous les pas, il faut donc
883	! sauvegarder les sorties du rayonnement
884	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
885	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
886	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
887	!
888	INTEGER itaprad
889	SAVE itaprad
890	!$OMP THREADPRIVATE(itaprad)
891	!
892	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
893	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
894	!
895	REAL zsav_tsol(klon)
896	!
897	REAL dist, rmu0(klon), fract(klon)
898	REAL zrmu0(klon), zfract(klon)
899	REAL zdtime, zdtime1, zdtime2, zlongi
900	!
901	REAL z_avant(klon), z_apres(klon), z_factor(klon)
902	LOGICAL zx_ajustq
903	!
904	REAL za
905	REAL zx_t, zx_qs, zdelta, zcor
906	real zqsat(klon,klev)
907	!
908	INTEGER i, k, iq, nsrf, l, itr
909	!
910	REAL t_coup
911	PARAMETER (t_coup=234.0)
912
913	!ym A voir plus tard !!
914	!ym REAL zx_relief(iim,jjmp1)
915	!ym REAL zx_aire(iim,jjmp1)
916	!
917	! Grandeurs de sorties
918	REAL s_capCL(klon)
919	REAL s_oliqCL(klon), s_cteiCL(klon)
920	REAL s_trmb1(klon), s_trmb2(klon)
921	REAL s_trmb3(klon)
922
923	! La convection n'est pas calculee tous les pas, il faut donc
924	! sauvegarder les sorties de la convection
925	!ym SAVE
926	!ym SAVE
927	!ym SAVE
928	!
929	INTEGER itapcv, itapwk
930	SAVE itapcv, itapwk
931	!$OMP THREADPRIVATE(itapcv, itapwk)
932
933	!KE43
934	! Variables locales pour la convection de K. Emanuel (sb):
935
936	REAL tvp(klon,klev) ! virtual temp of lifted parcel
937	CHARACTER*40 capemaxcels !max(CAPE)
938
939	REAL rflag(klon) ! flag fonctionnement de convect
940	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
941
942	! -- convect43:
943	INTEGER ntra ! nb traceurs pour convect4.3
944	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
945	REAL dplcldt(klon), dplcldr(klon)
946	!? . condm_con(klon,klev),conda_con(klon,klev),
947	!? . mr_con(klon,klev),ep_con(klon,klev)
948	!? . ,sadiab(klon,klev),wadiab(klon,klev)
949	! --
950	!34EK
951	!
952	! Variables du changement
953	!
954	! con: convection
955	! lsc: condensation a grande echelle (Large-Scale-Condensation)
956	! ajs: ajustement sec
957	! eva: evaporation de l'eau liquide nuageuse
958	! vdf: couche limite (Vertical DiFfusion)
959	!
960	! tendance nulles
961	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
962	REAL, dimension(klon) :: dsig0, ddens0
963	INTEGER, dimension(klon) :: wkoccur1
964	! tendance buffer pour appel de add_phys_tend
965	REAL, DIMENSION(klon,klev) :: d_q_ch4_dtime
966	!
967	! Flag pour pouvoir ne pas ajouter les tendances.
968	! Par defaut, les tendances doivente etre ajoutees et
969	! flag_inhib_tend = 0
970	! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
971	! croissant de print quand la valeur du flag augmente
972	!!! attention, ce flag doit etre change avec prudence !!!
973	INTEGER :: flag_inhib_tend = 0 ! 0 is the default value
974	!! INTEGER :: flag_inhib_tend = 2
975	!
976	! Logical switch to a bug : reseting to 0 convective variables at the
977	! begining of physiq.
978	LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
979	!$OMP THREADPRIVATE(ok_bug_cv_trac)
980	!
981	! Logical switch to a bug : changing wake_deltat when thermals are active
982	! even when there are no wakes.
983	LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
984	!$OMP THREADPRIVATE(ok_bug_split_th)
985
986	! Logical switch to a bug : modifying directly wake_deltat by adding
987	! the (w) dry adjustment tendency to wake_deltat
988	LOGICAL, SAVE :: ok_bug_ajs_cv = .TRUE.
989	!$OMP THREADPRIVATE(ok_bug_ajs_cv)
990
991	!
992	!********************************************************
993	! declarations
994
995	!********************************************************
996	!IM 081204 END
997	!
998	REAL pen_u(klon,klev), pen_d(klon,klev)
999	REAL pde_u(klon,klev), pde_d(klon,klev)
1000	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
1001	!
1002	REAL ratqsbas,ratqshaut,tau_ratqs
1003	SAVE ratqsbas,ratqshaut,tau_ratqs
1004	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
1005	REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
1006	!$OMP THREADPRIVATE(ratqsp0, ratqsdp)
1007
1008	! Parametres lies au nouveau schema de nuages (SB, PDF)
1009	REAL, SAVE :: fact_cldcon
1010	REAL, SAVE :: facttemps
1011	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
1012	LOGICAL, SAVE :: ok_newmicro
1013	!$OMP THREADPRIVATE(ok_newmicro)
1014
1015	INTEGER, SAVE :: iflag_cld_th
1016	!$OMP THREADPRIVATE(iflag_cld_th)
1017	!IM logical ptconv(klon,klev) !passe dans phys_local_var_mod
1018	!IM cf. AM 081204 BEG
1019	LOGICAL ptconvth(klon,klev)
1020
1021	REAL picefra(klon,klev)
1022	REAL zrel_oro(klon)
1023	!IM cf. AM 081204 END
1024	!
1025	! Variables liees a l'ecriture de la bande histoire physique
1026	!
1027	!======================================================================
1028	!
1029	!
1030	!JLD integer itau_w ! pas de temps ecriture = itap + itau_phy
1031	!
1032	!
1033	! Variables locales pour effectuer les appels en serie
1034	!
1035	!IM RH a 2m (la surface)
1036	REAL Lheat
1037
1038	INTEGER length
1039	PARAMETER ( length = 100 )
1040	REAL tabcntr0( length )
1041	!
1042	!JLD INTEGER ndex2d(nbp_lon*nbp_lat)
1043	!IM
1044	!
1045	!IM AMIP2 BEG
1046	!JLD REAL moyglo, mountor
1047	!IM 141004 BEG
1048	REAL zustrdr(klon), zvstrdr(klon)
1049	REAL zustrli(klon), zvstrli(klon)
1050	REAL zustrph(klon), zvstrph(klon)
1051	REAL aam, torsfc
1052	!IM 141004 END
1053	!IM 190504 BEG
1054	! INTEGER imp1jmp1
1055	! PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1056	!ym A voir plus tard
1057	! REAL zx_tmp((nbp_lon+1)*nbp_lat)
1058	! REAL airedyn(nbp_lon+1,nbp_lat)
1059	!IM 190504 END
1060	!JLD LOGICAL ok_msk
1061	!JLD REAL msk(klon)
1062	!ym A voir plus tard
1063	!ym REAL zm_wo(jjmp1, klev)
1064	!IM AMIP2 END
1065	!
1066	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
1067	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
1068	!JLD REAL zx_tmp_2d(nbp_lon,nbp_lat)
1069	!JLD REAL zx_lon(nbp_lon,nbp_lat)
1070	!JLD REAL zx_lat(nbp_lon,nbp_lat)
1071	!
1072	INTEGER nid_ctesGCM
1073	SAVE nid_ctesGCM
1074	!$OMP THREADPRIVATE(nid_ctesGCM)
1075	!
1076	!IM 280405 BEG
1077	! INTEGER nid_bilKPins, nid_bilKPave
1078	! SAVE nid_bilKPins, nid_bilKPave
1079	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1080	!
1081	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1082	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1083	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1084	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1085	!
1086	!JLD REAL zjulian
1087	!JLD SAVE zjulian
1088	!JLD!$OMP THREADPRIVATE(zjulian)
1089
1090	!JLD INTEGER nhori, nvert
1091	!JLD REAL zsto
1092	!JLD REAL zstophy, zout
1093
1094	CHARACTER (LEN=20) :: modname='physiq_mod'
1095	CHARACTER*80 abort_message
1096	LOGICAL, SAVE :: ok_sync, ok_sync_omp
1097	!$OMP THREADPRIVATE(ok_sync)
1098	REAL date0
1099
1100	! essai writephys
1101	INTEGER fid_day, fid_mth, fid_ins
1102	PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3)
1103	INTEGER prof2d_on, prof3d_on, prof2d_av, prof3d_av
1104	PARAMETER (prof2d_on = 1, prof3d_on = 2, prof2d_av = 3, prof3d_av = 4)
1105	REAL ztsol(klon)
1106	REAL q2m(klon,nbsrf) ! humidite a 2m
1107	REAL fsnowerosion(klon,nbsrf) ! blowing snow flux at surface
1108	REAL qbsfra ! blowing snow fraction
1109	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1110	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1111	CHARACTER*40 tinst, tave
1112	REAL cldtaupi(klon,klev) ! Cloud optical thickness for
1113	! pre-industrial (pi) aerosols
1114
1115	INTEGER :: naero
1116	! Aerosol optical properties
1117	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1118	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass
1119	! concentration
1120	! for all soluble
1121	! aerosols[ug/m3]
1122	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi
1123	! - " - (pre-industrial value)
1124	REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer
1125
1126	! Parameters
1127	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1128	LOGICAL ok_alw ! Apply aerosol LW effect or not
1129	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1130	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1131	SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1132	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
1133	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1134	! false : lecture des aerosol dans un fichier
1135	!$OMP THREADPRIVATE(aerosol_couple)
1136	LOGICAL, SAVE :: chemistry_couple ! true : use INCA chemistry O3
1137	! false : use offline chemistry O3
1138	!$OMP THREADPRIVATE(chemistry_couple)
1139	INTEGER, SAVE :: flag_aerosol
1140	!$OMP THREADPRIVATE(flag_aerosol)
1141	LOGICAL, SAVE :: flag_bc_internal_mixture
1142	!$OMP THREADPRIVATE(flag_bc_internal_mixture)
1143	!
1144	!--STRAT AEROSOL
1145	INTEGER, SAVE :: flag_aerosol_strat
1146	!$OMP THREADPRIVATE(flag_aerosol_strat)
1147	!
1148	!--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1149	LOGICAL, SAVE :: flag_aer_feedback
1150	!$OMP THREADPRIVATE(flag_aer_feedback)
1151
1152	!c-fin STRAT AEROSOL
1153	!
1154	! Declaration des constantes et des fonctions thermodynamiques
1155	!
1156	LOGICAL,SAVE :: first=.TRUE.
1157	!$OMP THREADPRIVATE(first)
1158
1159	! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1160	! Note that pressure vectors are in Pa and in stricly ascending order
1161	INTEGER,SAVE :: read_climoz ! Read ozone climatology
1162	! (let it keep the default OpenMP shared attribute)
1163	! Allowed values are 0, 1 and 2
1164	! 0: do not read an ozone climatology
1165	! 1: read a single ozone climatology that will be used day and night
1166	! 2: read two ozone climatologies, the average day and night
1167	! climatology and the daylight climatology
1168	INTEGER,SAVE :: ncid_climoz ! NetCDF file identifier
1169	REAL, ALLOCATABLE, SAVE :: press_cen_climoz(:) ! Pressure levels
1170	REAL, ALLOCATABLE, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1171	REAL, ALLOCATABLE, SAVE :: time_climoz(:) ! Time vector
1172
1173	CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1174	= ["tro3 ","tro3_daylight"]
1175	! vars_climoz(1:read_climoz): variables names in climoz file.
1176	! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1177	REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1178	! the ozone fields, old method.
1179
1180	include "FCTTRE.h"
1181	!IM 100106 END : pouvoir sortir les ctes de la physique
1182	!
1183	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1184	! Declarations pour Simulateur COSP
1185	!============================================================
1186	! AI 10-22
1187	include "ini_COSP.h"
1188	real :: mr_ozone(klon,klev), phicosp(klon,klev)
1189
1190	!IM stations CFMIP
1191	INTEGER, SAVE :: nCFMIP
1192	!$OMP THREADPRIVATE(nCFMIP)
1193	INTEGER, PARAMETER :: npCFMIP=120
1194	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1195	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1196	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1197	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1198	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1199	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1200	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1201	!$OMP THREADPRIVATE(iGCM, jGCM)
1202	logical, dimension(nfiles) :: phys_out_filestations
1203	logical, parameter :: lNMC=.FALSE.
1204
1205	!IM betaCRF
1206	REAL, SAVE :: pfree, beta_pbl, beta_free
1207	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1208	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
1209	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
1210	LOGICAL, SAVE :: mskocean_beta
1211	!$OMP THREADPRIVATE(mskocean_beta)
1212	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1213	! cldemirad pour evaluer les
1214	! retros liees aux CRF
1215	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
1216	! pour radlwsw pour
1217	! tester "CRF off"
1218	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1219	! pour radlwsw pour
1220	! tester "CRF off"
1221	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
1222	! radlwsw pour tester
1223	! "CRF off"
1224	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
1225
1226	REAL :: calday, zxsnow_dummy(klon)
1227	! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1228	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1229	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1230	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1231	REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1232	REAL, DIMENSION(klon,nbtr) :: init_source
1233
1234	!lwoff=y : offset LW CRE for radiation code and other schemes
1235	REAL, SAVE :: betalwoff
1236	!$OMP THREADPRIVATE(betalwoff)
1237	!
1238	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1239	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1240	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1241	integer iostat
1242
1243	REAL, dimension(klon,klev+1) :: l_mix_ave, wprime_ave
1244	REAL zzz
1245	!albedo SB >>>
1246	REAL,DIMENSION(6), SAVE :: SFRWL
1247	!$OMP THREADPRIVATE(SFRWL)
1248	!albedo SB <<<
1249
1250	!--OB variables for mass fixer (hard coded for now)
1251	REAL qql1(klon),qql2(klon),corrqql
1252
1253	!--OB flag to activate better conservation of water tendency when convection is not called every timestep
1254	LOGICAL, PARAMETER :: ok_conserv_d_q_con=.FALSE.
1255
1256	REAL, dimension(klon,klev) :: t_env,q_env
1257
1258	REAL, dimension(klon) :: pr_et
1259	REAL, dimension(klon) :: w_et, jlr_g_c, jlr_g_s
1260
1261	REAL pi
1262	REAL viscom, viscoh
1263	INTEGER ieru
1264
1265	!AI namelist pour gerer le double appel de Ecrad
1266	CHARACTER(len=512) :: namelist_ecrad_file
1267
1268
1269	! Subgrid scale wind :
1270	! Need to be allocatable/save because the number of bin is not known (provided by surf_wind_ini)
1271	integer, save :: nsrfwnd=1
1272	real, dimension(:,:), allocatable, save :: surf_wind_value, surf_wind_proba ! module and probability of sugrdi wind wind sample
1273	!$OMP THREADPRIVATE(nsrfwnd,surf_wind_value, surf_wind_proba)
1274
1275
1276
1277	!======================================================================!
1278	! Bifurcation vers un nouveau moniteur physique pour experimenter !
1279	! des solutions et préparer le couplage avec la physique de MesoNH !
1280	! 14 mai 2023 !
1281	!======================================================================!
1282	if (debut) then !
1283	iflag_physiq=0
1284	call getin_p('iflag_physiq', iflag_physiq) !
1285	endif !
1286	if ( iflag_physiq == 2 ) then !
1287	call physiqex (nlon,nlev, & !
1288	debut,lafin,pdtphys_, & !
1289	paprs,pplay,pphi,pphis,presnivs, & !
1290	u,v,rot,t,qx, & !
1291	flxmass_w, & !
1292	d_u, d_v, d_t, d_qx, d_ps) !
1293	return !
1294	endif !
1295	!======================================================================!
1296
1297
1298	pi = 4. * ATAN(1.)
1299
1300	! set-up call to alerte function
1301	call_alert = (alert_first_call .AND. is_master)
1302
1303	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1304	jjmp1=nbp_lat
1305
1306	!======================================================================
1307	! Gestion calendrier : mise a jour du module phys_cal_mod
1308	!
1309	pdtphys=pdtphys_
1310	CALL update_time(pdtphys)
1311	phys_tstep=NINT(pdtphys)
1312	IF (.NOT. using_xios) missing_val=nf90_fill_real
1313
1314	IF (using_xios) THEN
1315	! switch to XIOS LMDZ physics context
1316	IF (.NOT. debut .AND. is_omp_master) THEN
1317	CALL wxios_set_context()
1318	CALL xios_update_calendar(itap+1)
1319	ENDIF
1320	ENDIF
1321
1322	!======================================================================
1323	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1324	! Utilise notamment en 1D mais peut etre active egalement en 3D
1325	! en imposant la valeur de igout.
1326	!======================================================================
1327	IF (prt_level.ge.1) THEN
1328	igout=klon/2+1/klon
1329	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1330	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1331	longitude_deg(igout)
1332	write(lunout,*) &
1333	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1334	write(lunout,*) &
1335	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1336
1337	write(lunout,*) 'paprs, play, phi, u, v, t'
1338	DO k=1,klev
1339	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1340	u(igout,k),v(igout,k),t(igout,k)
1341	ENDDO
1342	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1343	DO k=1,klev
1344	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1345	ENDDO
1346	ENDIF
1347
1348	! Quick check on pressure levels:
1349	CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1350	"physiq_mod paprs bad order")
1351
1352	IF (first) THEN
1353	ivap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
1354	iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
1355	isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
1356	ibs = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
1357	icf = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
1358	irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
1359	ircont = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
1360	! CALL init_etat0_limit_unstruct
1361	! IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
1362	!CR:nvelles variables convection/poches froides
1363
1364	WRITE(lunout,*) '================================================='
1365	WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
1366	WRITE(lunout,*) '================================================='
1367	CALL phys_local_var_init
1368	!
1369	! appel a la lecture du run.def physique
1370	CALL conf_phys(ok_journe, ok_mensuel, &
1371	ok_instan, ok_hf, &
1372	ok_LES, &
1373	callstats, &
1374	solarlong0,seuil_inversion, &
1375	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1376	iflag_cld_th,ratqsbas,ratqshaut,tau_ratqs, &
1377	ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
1378	chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
1379	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1380	! nv flags pour la convection et les
1381	! poches froides
1382	read_climoz, &
1383	alp_offset)
1384	CALL init_etat0_limit_unstruct
1385	IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
1386	CALL phys_state_var_init(read_climoz)
1387	CALL phys_output_var_init
1388	IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
1389	CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1390
1391	print*, '================================================='
1392	!
1393	!CR: check sur le nb de traceurs de l eau
1394	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1395	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1396	'(H2O_g, H2O_l, H2O_s) but nqo=', nqo, '. Might as well stop here.'
1397	abort_message='see above'
1398	CALL abort_physic(modname,abort_message,1)
1399	ENDIF
1400
1401	IF (ok_ice_supersat.AND.(iflag_ice_thermo.EQ.0)) THEN
1402	WRITE (lunout, *) ' ok_ice_supersat=y requires iflag_ice_thermo=1 as well'
1403	abort_message='see above'
1404	CALL abort_physic(modname,abort_message,1)
1405	ENDIF
1406
1407	IF (ok_ice_supersat.AND.(nqo.LT.5)) THEN
1408	WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
1409	'(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
1410	abort_message='see above'
1411	CALL abort_physic(modname,abort_message,1)
1412	ENDIF
1413
1414	IF (ok_plane_h2o.AND..NOT.ok_ice_supersat) THEN
1415	WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_supersat=y '
1416	abort_message='see above'
1417	CALL abort_physic(modname,abort_message,1)
1418	ENDIF
1419
1420	IF (ok_plane_contrail.AND..NOT.ok_ice_supersat) THEN
1421	WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_supersat=y '
1422	abort_message='see above'
1423	CALL abort_physic(modname,abort_message,1)
1424	ENDIF
1425
1426	IF (ok_plane_contrail.AND.(nqo.LT.6)) THEN
1427	WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
1428	'(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
1429	abort_message='see above'
1430	CALL abort_physic(modname,abort_message,1)
1431	ENDIF
1432
1433	IF (ok_bs) THEN
1434	IF ((ok_ice_supersat.AND.nqo .LT.6).OR.(.NOT.ok_ice_supersat.AND.nqo.LT.4)) THEN
1435	WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
1436	'but nqo=', nqo
1437	abort_message='see above'
1438	CALL abort_physic(modname,abort_message, 1)
1439	ENDIF
1440	ENDIF
1441
1442	Ncvpaseq1 = 0
1443	dnwd0=0.0
1444	ftd=0.0
1445	fqd=0.0
1446	cin=0.
1447	!ym Attention pbase pas initialise dans concvl !!!!
1448	pbase=0
1449	!IM 180608
1450
1451	itau_con=0
1452	first=.FALSE.
1453
1454	ENDIF ! first
1455
1456	!ym => necessaire pour iflag_con != 2
1457	pmfd(:,:) = 0.
1458	pen_u(:,:) = 0.
1459	pen_d(:,:) = 0.
1460	pde_d(:,:) = 0.
1461	pde_u(:,:) = 0.
1462	aam=0.
1463	d_t_adjwk(:,:)=0
1464	d_q_adjwk(:,:)=0
1465
1466	alp_bl_conv(:)=0.
1467
1468	torsfc=0.
1469	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1470
1471
1472	IF (debut) THEN
1473	CALL suphel ! initialiser constantes et parametres phys.
1474	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1475	tau_gl=5.
1476	CALL getin_p('tau_gl', tau_gl)
1477	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1478	! secondes
1479	tau_gl=86400.*tau_gl
1480	WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1481	iflag_thermcell_tke=0
1482	call getin_p('iflag_thermcell_tke', iflag_thermcell_tke) !
1483
1484	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1485	CALL getin_p('random_notrig_max',random_notrig_max)
1486	CALL getin_p('ok_adjwk',ok_adjwk)
1487	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1488	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1489	! 1 => convective adjustment but state variables are unchanged
1490	! 2 => convective adjustment and state variables are changed
1491	CALL getin_p('iflag_adjwk',iflag_adjwk)
1492	CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1493	CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
1494	CALL getin_p('oliqmax',oliqmax)
1495	CALL getin_p('oicemax',oicemax)
1496	CALL getin_p('ratqsp0',ratqsp0)
1497	CALL getin_p('ratqsdp',ratqsdp)
1498	iflag_wake_tend = 0
1499	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1500	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1501	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1502	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1503	CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
1504	CALL getin_p('ok_bug_split_th',ok_bug_split_th)
1505	CALL getin_p('ok_bug_ajs_cv',ok_bug_ajs_cv)
1506	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1507	CALL getin_p('fl_ebil',fl_ebil)
1508	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1509	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1510	iflag_phytrac = 1 ! by default we do want to call phytrac
1511	CALL getin_p('iflag_phytrac',iflag_phytrac)
1512
1513	ok_water_mass_fixer=.FALSE. ! OB: by default we do not apply the mass fixer
1514	CALL getin_p('ok_water_mass_fixer',ok_water_mass_fixer)
1515	IF (CPPKEY_DUST) THEN
1516	IF (iflag_phytrac.EQ.0) THEN
1517	WRITE(lunout,*) 'In order to run with SPLA, iflag_phytrac will be forced to 1'
1518	iflag_phytrac = 1
1519	ENDIF
1520	END IF
1521	nvm_lmdz = 13
1522	CALL getin_p('NVM',nvm_lmdz)
1523
1524	WRITE(lunout,*) 'iflag_alp_wk_cond=', iflag_alp_wk_cond
1525	WRITE(lunout,*) 'random_ntrig_max=', random_notrig_max
1526	WRITE(lunout,*) 'ok_adjwk=', ok_adjwk
1527	WRITE(lunout,*) 'iflag_adjwk=', iflag_adjwk
1528	WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1529	WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1530	WRITE(lunout,*) 'ratqsp0=', ratqsp0
1531	WRITE(lunout,*) 'ratqsdp=', ratqsdp
1532	WRITE(lunout,*) 'iflag_wake_tend=', iflag_wake_tend
1533	WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1534	WRITE(lunout,*) 'ok_bug_cv_trac=', ok_bug_cv_trac
1535	WRITE(lunout,*) 'ok_bug_split_th=', ok_bug_split_th
1536	WRITE(lunout,*) 'fl_ebil=', fl_ebil
1537	WRITE(lunout,*) 'fl_cor_ebil=', fl_cor_ebil
1538	WRITE(lunout,*) 'iflag_phytrac=', iflag_phytrac
1539	WRITE(lunout,*) 'ok_water_mass_fixer=',ok_water_mass_fixer
1540	WRITE(lunout,*) 'NVM=', nvm_lmdz
1541
1542	!--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1543	WRITE(lunout,*) 'Call to infocfields from physiq'
1544	CALL infocfields_init
1545
1546	!AI 08 2023
1547	#ifdef CPP_ECRAD
1548	ok_3Deffect=.false.
1549	CALL getin_p('ok_3Deffect',ok_3Deffect)
1550	namelist_ecrad_file='namelist_ecrad'
1551	#endif
1552
1553	ENDIF
1554
1555	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1556
1557	!======================================================================
1558	! Gestion calendrier : mise a jour du module phys_cal_mod
1559	!
1560	! CALL phys_cal_update(jD_cur,jH_cur)
1561
1562	!
1563	! Si c'est le debut, il faut initialiser plusieurs choses
1564	! ********
1565	!
1566	IF (debut) THEN
1567	!rv CRinitialisation de wght_th et lalim_conv pour la
1568	!definition de la couche alimentation de la convection a partir
1569	!des caracteristiques du thermique
1570	wght_th(:,:)=1.
1571	lalim_conv(:)=1
1572	!RC
1573	ustar(:,:)=0.
1574	! u10m(:,:)=0.
1575	! v10m(:,:)=0.
1576	rain_con(:)=0.
1577	snow_con(:)=0.
1578	topswai(:)=0.
1579	topswad(:)=0.
1580	solswai(:)=0.
1581	solswad(:)=0.
1582
1583	wmax_th(:)=0.
1584	tau_overturning_th(:)=0.
1585
1586	IF (ANY(type_trac == ['inca','inco'])) THEN
1587	! jg : initialisation jusqu'au ces variables sont dans restart
1588	ccm(:,:,:) = 0.
1589	tau_aero(:,:,:,:) = 0.
1590	piz_aero(:,:,:,:) = 0.
1591	cg_aero(:,:,:,:) = 0.
1592	d_q_ch4(:,:) = 0.
1593
1594	config_inca='none' ! default
1595	CALL getin_p('config_inca',config_inca)
1596
1597	ELSE
1598	config_inca='none' ! default
1599	ENDIF
1600
1601	tau_aero(:,:,:,:) = 1.e-15
1602	piz_aero(:,:,:,:) = 1.
1603	cg_aero(:,:,:,:) = 0.
1604	d_q_ch4(:,:) = 0.
1605
1606	IF (aerosol_couple .AND. (config_inca /= "aero" &
1607	.AND. config_inca /= "aeNP ")) THEN
1608	abort_message &
1609	= 'if aerosol_couple is activated, config_inca need to be ' &
1610	// 'aero or aeNP'
1611	CALL abort_physic (modname,abort_message,1)
1612	ENDIF
1613
1614	rnebcon0(:,:) = 0.0
1615	clwcon0(:,:) = 0.0
1616	rnebcon(:,:) = 0.0
1617	clwcon(:,:) = 0.0
1618
1619	!
1620	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1621	iflag_coupl,iflag_clos,iflag_wake
1622	print*,'iflag_cycle_diurne', iflag_cycle_diurne
1623	!
1624	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1625	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1626	CALL abort_physic (modname,abort_message,1)
1627	ENDIF
1628	!
1629	!
1630	! Initialiser les compteurs:
1631	!
1632	itap = 0
1633	itaprad = 0
1634	itapcv = 0
1635	itapwk = 0
1636
1637	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1638	!! Un petit travail \`a faire ici.
1639	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1640
1641	IF (iflag_pbl>1) THEN
1642	PRINT*, "Using method MELLOR&YAMADA"
1643	ENDIF
1644
1645	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1646	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1647	! phylmd plutot que dyn3d
1648	! Attention : la version precedente n'etait pas tres propre.
1649	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1650	! pour obtenir le meme resultat.
1651	!jyg for fh<
1652	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1653	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
1654	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1655	CALL abort_physic(modname,abort_message,1)
1656	ENDIF
1657	!>jyg
1658	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1659	radpas = NINT( 86400./phys_tstep)/nbapp_rad
1660	ELSE
1661	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1662	'multiple de nbapp_rad'
1663	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1664	'mais 1+1<>2'
1665	abort_message='nbre de pas de temps physique n est pas multiple ' &
1666	// 'de nbapp_rad'
1667	CALL abort_physic(modname,abort_message,1)
1668	ENDIF
1669	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1670	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
1671	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1672	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1673	cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
1674	cvpas = cvpas_0
1675	print *,'physiq, cvpas ',cvpas
1676	ELSE
1677	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1678	'multiple de nbapp_cv'
1679	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1680	'mais 1+1<>2'
1681	abort_message='nbre de pas de temps physique n est pas multiple ' &
1682	// 'de nbapp_cv'
1683	CALL abort_physic(modname,abort_message,1)
1684	ENDIF
1685	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1686	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
1687	! print *,'physiq, wkpas ',wkpas
1688	ELSE
1689	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1690	'multiple de nbapp_wk'
1691	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1692	'mais 1+1<>2'
1693	abort_message='nbre de pas de temps physique n est pas multiple ' &
1694	// 'de nbapp_wk'
1695	CALL abort_physic(modname,abort_message,1)
1696	ENDIF
1697	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1698	CALL init_iophy_new(latitude_deg,longitude_deg)
1699
1700	!===================================================================
1701	!IM stations CFMIP
1702	nCFMIP=npCFMIP
1703	OPEN(98,file='npCFMIP_param.data',status='old', &
1704	form='formatted',iostat=iostat)
1705	IF (iostat == 0) THEN
1706	READ(98,*,end=998) nCFMIP
1707	998 CONTINUE
1708	CLOSE(98)
1709	CONTINUE
1710	IF(nCFMIP.GT.npCFMIP) THEN
1711	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1712	CALL abort_physic("physiq", "", 1)
1713	ELSE
1714	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1715	ENDIF
1716
1717	!
1718	ALLOCATE(tabCFMIP(nCFMIP))
1719	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1720	ALLOCATE(tabijGCM(nCFMIP))
1721	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1722	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1723	!
1724	! lecture des nCFMIP stations CFMIP, de leur numero
1725	! et des coordonnees geographiques lonCFMIP, latCFMIP
1726	!
1727	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1728	lonCFMIP, latCFMIP)
1729	!
1730	! identification des
1731	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1732	! grille de LMDZ
1733	! 2) indices points tabijGCM de la grille physique 1d sur
1734	! klon points
1735	! 3) indices iGCM, jGCM de la grille physique 2d
1736	!
1737	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1738	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1739	!
1740	ELSE
1741	ALLOCATE(tabijGCM(0))
1742	ALLOCATE(lonGCM(0), latGCM(0))
1743	ALLOCATE(iGCM(0), jGCM(0))
1744	ENDIF
1745
1746	!$OMP MASTER
1747	! FH : if ok_sync=.true. , the time axis is written at each time step
1748	! in the output files. Only at the end in the opposite case
1749	ok_sync_omp=.FALSE.
1750	CALL getin('ok_sync',ok_sync_omp)
1751	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1752	iGCM,jGCM,lonGCM,latGCM, &
1753	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1754	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1755	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1756	read_climoz, phys_out_filestations, &
1757	aerosol_couple, &
1758	flag_aerosol_strat, pdtphys, paprs, pphis, &
1759	pplay, lmax_th, ptconv, ptconvth, ivap, &
1760	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1761	!$OMP END MASTER
1762	!$OMP BARRIER
1763	ok_sync=ok_sync_omp
1764
1765	freq_outNMC(1) = ecrit_files(7)
1766	freq_outNMC(2) = ecrit_files(8)
1767	freq_outNMC(3) = ecrit_files(9)
1768	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1769	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1770	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1771
1772	#ifndef CPP_XIOS
1773	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1774	#endif
1775
1776	ecrit_reg = ecrit_reg * un_jour
1777	ecrit_tra = ecrit_tra * un_jour
1778
1779	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1780	date0 = jD_ref
1781	WRITE(,) 'physiq date0 : ',date0
1782	!
1783
1784	! CALL create_climoz(read_climoz)
1785	IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol, aerosol_couple) !! initialise aero from file for XIOS interpolation (unstructured_grid)
1786	IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat) !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
1787
1788	! A.I : Initialisations pour le 1er passage a Cosp
1789	if (ok_cosp) then
1790
1791	IF (CPPKEY_COSP) THEN
1792	CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
1793	zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
1794	fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
1795	mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
1796
1797	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
1798	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1799	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1800	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1801	JrNt_cosp0,ref_liq_cosp0,ref_ice_cosp0, &
1802	pctsrf_cosp0, &
1803	zu10m_cosp0,zv10m_cosp0,pphis, &
1804	pphi,paprs(:,1:klev),pplay,zxtsol_cosp0,t, &
1805	qx(:,:,ivap),zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0,fiwc_cosp0, &
1806	prfl_cosp0(:,1:klev),psfl_cosp0(:,1:klev), &
1807	pmflxr_cosp0(:,1:klev),pmflxs_cosp0(:,1:klev), &
1808	mr_ozone_cosp0,cldtau_cosp0, cldemi_cosp0)
1809	END IF
1810
1811	IF (CPPKEY_COSPV2) THEN
1812	CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
1813	zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
1814	fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
1815	mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
1816
1817	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1818	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1819	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1820	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1821	JrNt_cosp0,ref_liq_cosp0,ref_ice_cosp0, &
1822	pctsrf_cosp0, &
1823	zu10m_cosp0,zv10m_cosp0,pphis, &
1824	pphi,paprs(:,1:klev),pplay,zxtsol_cosp0,t, &
1825	qx(:,:,ivap),zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0,fiwc_cosp0, &
1826	prfl_cosp0(:,1:klev),psfl_cosp0(:,1:klev), &
1827	pmflxr_cosp0(:,1:klev),pmflxs_cosp0(:,1:klev), &
1828	mr_ozone_cosp0,cldtau_cosp0, cldemi_cosp0)
1829	END IF
1830	ENDIF
1831
1832	!
1833	!
1834	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1835	! Nouvelle initialisation pour le rayonnement RRTM
1836	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1837
1838	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1839
1840	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1841	CALL surf_wind_ini(klon,lunout)
1842	CALL getin_p('nsrfwnd',nsrfwnd)
1843	allocate(surf_wind_value(klon,nsrfwnd),surf_wind_proba(klon,nsrfwnd))
1844
1845	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1846	CALL wake_ini(rg,rd,rv,prt_level)
1847	CALL yamada_ini(klon,lunout,prt_level)
1848	viscom=1.46E-5
1849	viscoh=2.06E-5
1850	CALL atke_ini(RG, RD, RPI, RCPD, RV, viscom, viscoh)
1851	CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
1852	& RG,RD,RCPD,RKAPPA,RLVTT,RETV)
1853	CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
1854	CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_supersat,&
1855	ok_plane_contrail,iflag_ratqs,fl_cor_ebil, &
1856	RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
1857	CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
1858	RVTMP2, RTT,RD,RG, RV, RPI)
1859	! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
1860	IF (ok_newmicro) then
1861	IF (iflag_rrtm.EQ.1) THEN
1862	#ifdef CPP_RRTM
1863	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
1864	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
1865	// 'pour ok_cdnc'
1866	CALL abort_physic(modname,abort_message,1)
1867	ENDIF
1868	#else
1869
1870	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
1871	CALL abort_physic(modname,abort_message,1)
1872	#endif
1873	ENDIF
1874	ENDIF
1875	CALL cloud_optics_prop_ini(klon, prt_level, lunout, flag_aerosol, &
1876	& ok_cdnc, bl95_b0, &
1877	& bl95_b1, latitude_deg, rpi, rg, rd, &
1878	& zepsec, novlp, iflag_ice_thermo, ok_new_lscp)
1879	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1880
1881	!
1882	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1883	! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
1884	!
1885	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1886	IF (CPPKEY_REPROBUS) THEN
1887	CALL strataer_init
1888	CALL strataer_emiss_init
1889	END IF
1890
1891	IF (CPPKEY_STRATAER) THEN
1892	CALL strataer_init
1893	CALL strataer_nuc_init
1894	CALL strataer_emiss_init
1895	END IF
1896
1897	IF (CPPKEY_DUST) THEN
1898	! Quand on utilise SPLA, on force iflag_phytrac=1
1899	CALL phytracr_spl_out_init()
1900	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
1901	pplay, lmax_th, aerosol_couple, &
1902	ok_ade, ok_aie, ivap, ok_sync, &
1903	ptconv, read_climoz, clevSTD, &
1904	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
1905	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1906	ELSE
1907	! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1
1908	! donc seulement dans ce cas on doit appeler phytrac_init()
1909	IF (iflag_phytrac == 1 ) THEN
1910	CALL phytrac_init()
1911	ENDIF
1912	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
1913	pplay, lmax_th, aerosol_couple, &
1914	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs, ok_sync,&
1915	ptconv, read_climoz, clevSTD, &
1916	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
1917	flag_aerosol, flag_aerosol_strat, ok_cdnc, t, u1, v1)
1918	END IF
1919
1920
1921	IF (using_xios) THEN
1922	IF (is_omp_master) CALL xios_update_calendar(1)
1923	ENDIF
1924
1925	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1926	CALL create_etat0_limit_unstruct
1927	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1928
1929	!jyg<
1930	IF (iflag_pbl<=1) THEN
1931	! No TKE for Standard Physics
1932	pbl_tke(:,:,:)=0.
1933
1934	ELSE IF (klon_glo==1) THEN
1935	pbl_tke(:,:,is_ave) = 0.
1936	pbl_eps(:,:,is_ave) = 0.
1937	DO nsrf=1,nbsrf
1938	DO k = 1,klev+1
1939	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1940	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1941	pbl_eps(:,k,is_ave) = pbl_eps(:,k,is_ave) &
1942	+pctsrf(:,nsrf)*pbl_eps(:,k,nsrf)
1943	ENDDO
1944	ENDDO
1945	ELSE
1946	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1947	!>jyg
1948	pbl_eps(:,:,is_ave) = 0.
1949	ENDIF
1950	!IM begin
1951	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1952	,ratqs(1,1)
1953	!IM end
1954
1955
1956	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1957	!
1958	! on remet le calendrier a zero
1959	!
1960	IF (raz_date .eq. 1) THEN
1961	itau_phy = 0
1962	ENDIF
1963
1964	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1965	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1966	! pdtphys
1967	! abort_message='Pas physique n est pas correct '
1968	! ! call abort_physic(modname,abort_message,1)
1969	! phys_tstep=pdtphys
1970	! ENDIF
1971	IF (nlon .NE. klon) THEN
1972	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1973	klon
1974	abort_message='nlon et klon ne sont pas coherents'
1975	CALL abort_physic(modname,abort_message,1)
1976	ENDIF
1977	IF (nlev .NE. klev) THEN
1978	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1979	klev
1980	abort_message='nlev et klev ne sont pas coherents'
1981	CALL abort_physic(modname,abort_message,1)
1982	ENDIF
1983	!
1984	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1985	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1986	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1987	abort_message='Nbre d appels au rayonnement insuffisant'
1988	CALL abort_physic(modname,abort_message,1)
1989	ENDIF
1990
1991	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1992	! Initialisation pour la convection de K.E. et pour les poches froides
1993	!
1994	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1995
1996	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1997	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", ok_cvl
1998	!
1999	!KE43
2000	! Initialisation pour la convection de K.E. (sb):
2001	IF (iflag_con.GE.3) THEN
2002
2003	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
2004	WRITE(lunout,*) &
2005	"On va utiliser le melange convectif des traceurs qui"
2006	WRITE(lunout,*)"est calcule dans convect4.3"
2007	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
2008
2009	DO i = 1, klon
2010	ema_cbmf(i) = 0.
2011	ema_pcb(i) = 0.
2012	ema_pct(i) = 0.
2013	! ema_workcbmf(i) = 0.
2014	ENDDO
2015	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
2016	DO i = 1, klon
2017	ibas_con(i) = 1
2018	itop_con(i) = 1
2019	ENDDO
2020	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
2021	!================================================================
2022	!CR:04.12.07: initialisations poches froides
2023	! Controle de ALE et ALP pour la fermeture convective (jyg)
2024	IF (iflag_wake>=1) THEN
2025	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
2026	,alp_bl_prescr, ale_bl_prescr)
2027	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
2028	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
2029	!
2030	! Initialize tendencies of wake state variables (for some flag values
2031	! they are not computed).
2032	d_deltat_wk(:,:) = 0.
2033	d_deltaq_wk(:,:) = 0.
2034	d_deltat_wk_gw(:,:) = 0.
2035	d_deltaq_wk_gw(:,:) = 0.
2036	d_deltat_vdf(:,:) = 0.
2037	d_deltaq_vdf(:,:) = 0.
2038	d_deltat_the(:,:) = 0.
2039	d_deltaq_the(:,:) = 0.
2040	d_deltat_ajs_cv(:,:) = 0.
2041	d_deltaq_ajs_cv(:,:) = 0.
2042	d_s_wk(:) = 0.
2043	d_s_a_wk(:) = 0.
2044	d_dens_wk(:) = 0.
2045	d_dens_a_wk(:) = 0.
2046	ENDIF ! (iflag_wake>=1)
2047
2048	! do i = 1,klon
2049	! Ale_bl(i)=0.
2050	! Alp_bl(i)=0.
2051	! enddo
2052
2053	!ELSE
2054	! ALLOCATE(tabijGCM(0))
2055	! ALLOCATE(lonGCM(0), latGCM(0))
2056	! ALLOCATE(iGCM(0), jGCM(0))
2057	ENDIF ! (iflag_con.GE.3)
2058	!
2059	DO i=1,klon
2060	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
2061	ENDDO
2062
2063	!34EK
2064	IF (ok_orodr) THEN
2065
2066	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2067	! FH sans doute a enlever de finitivement ou, si on le
2068	! garde, l'activer justement quand ok_orodr = false.
2069	! ce rugoro est utilise par la couche limite et fait double emploi
2070	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
2071	! DO i=1,klon
2072	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
2073	! ENDDO
2074	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2075	IF (ok_strato) THEN
2076	CALL SUGWD_strato(klon,klev,paprs,pplay)
2077	ELSE
2078	CALL SUGWD(klon,klev,paprs,pplay)
2079	ENDIF
2080
2081	DO i=1,klon
2082	zuthe(i)=0.
2083	zvthe(i)=0.
2084	IF (zstd(i).gt.10.) THEN
2085	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
2086	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
2087	ENDIF
2088	ENDDO
2089	ENDIF
2090	!
2091	!
2092	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
2093	WRITE(lunout,*)'La frequence de lecture surface est de ', &
2094	lmt_pas
2095	!
2096	capemaxcels = 't_max(X)'
2097	t2mincels = 't_min(X)'
2098	t2maxcels = 't_max(X)'
2099	tinst = 'inst(X)'
2100	tave = 'ave(X)'
2101	!IM cf. AM 081204 BEG
2102	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
2103	!IM cf. AM 081204 END
2104	!
2105	!=============================================================
2106	! Initialisation des sorties
2107	!=============================================================
2108
2109	IF (using_xios) THEN
2110	! Get "missing_val" value from XML files (from temperature variable)
2111	IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
2112	CALL bcast_omp(missing_val)
2113	ENDIF
2114
2115	IF (using_xios) THEN
2116	! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
2117	! initialised at that moment
2118	! Get "missing_val" value from XML files (from temperature variable)
2119	IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
2120	CALL bcast_omp(missing_val)
2121	!
2122	! Now we activate some double radiation call flags only if some
2123	! diagnostics are requested, otherwise there is no point in doing this
2124	IF (is_master) THEN
2125	!--setting up swaero_diag to TRUE in XIOS case
2126	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
2127	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
2128	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
2129	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
2130	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
2131	!!!--for now these fields are not in the XML files so they are omitted
2132	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
2133	swaero_diag=.TRUE.
2134
2135	!--setting up swaerofree_diag to TRUE in XIOS case
2136	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
2137	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
2138	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
2139	xios_field_is_active("LWupTOAcleanclr")) &
2140	swaerofree_diag=.TRUE.
2141
2142	!--setting up dryaod_diag to TRUE in XIOS case
2143	DO naero = 1, naero_tot-1
2144	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
2145	ENDDO
2146	!
2147	!--setting up ok_4xCO2atm to TRUE in XIOS case
2148	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
2149	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
2150	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
2151	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
2152	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
2153	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
2154	ok_4xCO2atm=.TRUE.
2155	ENDIF
2156	!$OMP BARRIER
2157	CALL bcast(swaero_diag)
2158	CALL bcast(swaerofree_diag)
2159	CALL bcast(dryaod_diag)
2160	CALL bcast(ok_4xCO2atm)
2161	ENDIF !using_xios
2162	!
2163	CALL printflag( tabcntr0,radpas,ok_journe, &
2164	ok_instan, ok_region )
2165	!
2166	!
2167	! Prescrire l'ozone dans l'atmosphere
2168	!
2169	!c DO i = 1, klon
2170	!c DO k = 1, klev
2171	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
2172	!c ENDDO
2173	!c ENDDO
2174	!
2175	IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
2176	IF (CPPKEY_INCA) THEN
2177	CALL VTe(VTphysiq)
2178	CALL VTb(VTinca)
2179	calday = REAL(days_elapsed) + jH_cur
2180	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
2181
2182	call init_const_lmdz( &
2183	ndays, nbsrf, is_oce,is_sic, is_ter,is_lic, calend, &
2184	config_inca)
2185
2186	CALL init_inca_geometry( &
2187	longitude, latitude, &
2188	boundslon, boundslat, &
2189	cell_area, ind_cell_glo)
2190
2191	if (grid_type==unstructured) THEN
2192	CALL chemini( pplay, &
2193	nbp_lon, nbp_lat, &
2194	latitude_deg, &
2195	longitude_deg, &
2196	presnivs, &
2197	calday, &
2198	klon, &
2199	nqtot, &
2200	nqo+nqCO2, &
2201	pdtphys, &
2202	annee_ref, &
2203	year_cur, &
2204	day_ref, &
2205	day_ini, &
2206	start_time, &
2207	itau_phy, &
2208	date0, &
2209	chemistry_couple, &
2210	init_source, &
2211	init_tauinca, &
2212	init_pizinca, &
2213	init_cginca, &
2214	init_ccminca)
2215	ELSE
2216	CALL chemini( pplay, &
2217	nbp_lon, nbp_lat, &
2218	latitude_deg, &
2219	longitude_deg, &
2220	presnivs, &
2221	calday, &
2222	klon, &
2223	nqtot, &
2224	nqo+nqCO2, &
2225	pdtphys, &
2226	annee_ref, &
2227	year_cur, &
2228	day_ref, &
2229	day_ini, &
2230	start_time, &
2231	itau_phy, &
2232	date0, &
2233	chemistry_couple, &
2234	init_source, &
2235	init_tauinca, &
2236	init_pizinca, &
2237	init_cginca, &
2238	init_ccminca, &
2239	io_lon, &
2240	io_lat)
2241	ENDIF
2242
2243
2244	! initialisation des variables depuis le restart de inca
2245	ccm(:,:,:) = init_ccminca
2246	tau_aero(:,:,:,:) = init_tauinca
2247	piz_aero(:,:,:,:) = init_pizinca
2248	cg_aero(:,:,:,:) = init_cginca
2249	!
2250
2251
2252	CALL VTe(VTinca)
2253	CALL VTb(VTphysiq)
2254	END IF
2255	ENDIF
2256	!
2257	IF (type_trac == 'repr') THEN
2258	IF (CPPKEY_REPROBUS) THEN
2259	CALL chemini_rep( &
2260	presnivs, &
2261	pdtphys, &
2262	annee_ref, &
2263	day_ref, &
2264	day_ini, &
2265	start_time, &
2266	itau_phy, &
2267	io_lon, &
2268	io_lat)
2269	END IF
2270	ENDIF
2271
2272	!$omp single
2273	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
2274	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
2275	!$omp end single
2276	!
2277	!IM betaCRF
2278	pfree=70000. !Pa
2279	beta_pbl=1.
2280	beta_free=1.
2281	lon1_beta=-180.
2282	lon2_beta=+180.
2283	lat1_beta=90.
2284	lat2_beta=-90.
2285	mskocean_beta=.FALSE.
2286
2287	!albedo SB >>>
2288	SELECT CASE(nsw)
2289	CASE(2)
2290	SFRWL(1)=0.45538747
2291	SFRWL(2)=0.54461211
2292	CASE(4)
2293	SFRWL(1)=0.45538747
2294	SFRWL(2)=0.32870591
2295	SFRWL(3)=0.18568763
2296	SFRWL(4)=3.02191470E-02
2297	CASE(6)
2298	SFRWL(1)=1.28432794E-03
2299	SFRWL(2)=0.12304168
2300	SFRWL(3)=0.33106142
2301	SFRWL(4)=0.32870591
2302	SFRWL(5)=0.18568763
2303	SFRWL(6)=3.02191470E-02
2304	END SELECT
2305	!albedo SB <<<
2306
2307	OPEN(99,file='beta_crf.data',status='old', &
2308	form='formatted',err=9999)
2309	READ(99,*,end=9998) pfree
2310	READ(99,*,end=9998) beta_pbl
2311	READ(99,*,end=9998) beta_free
2312	READ(99,*,end=9998) lon1_beta
2313	READ(99,*,end=9998) lon2_beta
2314	READ(99,*,end=9998) lat1_beta
2315	READ(99,*,end=9998) lat2_beta
2316	READ(99,*,end=9998) mskocean_beta
2317	9998 Continue
2318	CLOSE(99)
2319	9999 Continue
2320	WRITE(,)'pfree=',pfree
2321	WRITE(,)'beta_pbl=',beta_pbl
2322	WRITE(,)'beta_free=',beta_free
2323	WRITE(,)'lon1_beta=',lon1_beta
2324	WRITE(,)'lon2_beta=',lon2_beta
2325	WRITE(,)'lat1_beta=',lat1_beta
2326	WRITE(,)'lat2_beta=',lat2_beta
2327	WRITE(,)'mskocean_beta=',mskocean_beta
2328
2329	!lwoff=y : offset LW CRE for radiation code and other schemes
2330	!lwoff=y : betalwoff=1.
2331	betalwoff=0.
2332	IF (ok_lwoff) THEN
2333	betalwoff=1.
2334	ENDIF
2335	WRITE(,)'ok_lwoff=',ok_lwoff
2336	!
2337	!lwoff=y to begin only sollw and sollwdown are set up to CS values
2338	sollw = sollw + betalwoff * (sollw0 - sollw)
2339	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
2340	sollwdown(:))
2341
2342	!--Init for LSCP - condensation
2343	ratio_ql_qtot(:,:) = 0.
2344	ratio_qi_qtot(:,:) = 0.
2345
2346
2347	ENDIF
2348	!
2349	! ************** Fin de IF ( debut ) *************
2350	!
2351	!
2352	! Incrementer le compteur de la physique
2353	!
2354	itap = itap + 1
2355	IF (is_master .OR. prt_level > 9) THEN
2356	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
2357	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2358	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2359	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
2360	ENDIF
2361	ENDIF
2362	!
2363	!
2364	! Update fraction of the sub-surfaces (pctsrf) and
2365	! initialize, where a new fraction has appeared, all variables depending
2366	! on the surface fraction.
2367	!
2368	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
2369	pctsrf, fevap, z0m, z0h, agesno, &
2370	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
2371
2372	! Update time and other variables in Reprobus
2373	IF (type_trac == 'repr') THEN
2374	IF (CPPKEY_REPROBUS) THEN
2375	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2376	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2377	CALL Rtime(debut)
2378	END IF
2379	ENDIF
2380
2381	! Tendances bidons pour les processus qui n'affectent pas certaines
2382	! variables.
2383	du0(:,:)=0.
2384	dv0(:,:)=0.
2385	dt0 = 0.
2386	dq0(:,:)=0.
2387	dql0(:,:)=0.
2388	dqi0(:,:)=0.
2389	dqbs0(:,:)=0.
2390	dsig0(:) = 0.
2391	ddens0(:) = 0.
2392	wkoccur1(:)=1
2393	!
2394	! Mettre a zero des variables de sortie (pour securite)
2395	!
2396	DO i = 1, klon
2397	d_ps(i) = 0.0
2398	ENDDO
2399	DO k = 1, klev
2400	DO i = 1, klon
2401	d_t(i,k) = 0.0
2402	d_u(i,k) = 0.0
2403	d_v(i,k) = 0.0
2404	ENDDO
2405	ENDDO
2406	DO iq = 1, nqtot
2407	DO k = 1, klev
2408	DO i = 1, klon
2409	d_qx(i,k,iq) = 0.0
2410	ENDDO
2411	ENDDO
2412	ENDDO
2413	beta_prec_fisrt(:,:)=0.
2414	beta_prec(:,:)=0.
2415	!
2416	! Output variables from the convective scheme should not be set to 0
2417	! since convection is not always called at every time step.
2418	IF (ok_bug_cv_trac) THEN
2419	da(:,:)=0.
2420	mp(:,:)=0.
2421	phi(:,:,:)=0.
2422	! RomP >>>
2423	phi2(:,:,:)=0.
2424	epmlmMm(:,:,:)=0.
2425	eplaMm(:,:)=0.
2426	d1a(:,:)=0.
2427	dam(:,:)=0.
2428	pmflxr(:,:)=0.
2429	pmflxs(:,:)=0.
2430	! RomP <<<
2431	ENDIF
2432	!
2433	! Ne pas affecter les valeurs entrees de u, v, h, et q
2434	!
2435	DO k = 1, klev
2436	DO i = 1, klon
2437	t_seri(i,k) = t(i,k)
2438	u_seri(i,k) = u(i,k)
2439	v_seri(i,k) = v(i,k)
2440	q_seri(i,k) = qx(i,k,ivap)
2441	ql_seri(i,k) = qx(i,k,iliq)
2442	qbs_seri(i,k)= 0.
2443	cf_seri(i,k) = 0.
2444	rvc_seri(i,k)= 0.
2445	rcont_seri(i,k)= 0.
2446	!CR: ATTENTION, on rajoute la variable glace
2447	IF (nqo.EQ.2) THEN !--vapour and liquid only
2448	qs_seri(i,k) = 0.
2449	ELSE IF (nqo.EQ.3) THEN !--vapour, liquid and ice
2450	qs_seri(i,k) = qx(i,k,isol)
2451	ELSE IF (nqo.GE.4) THEN !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
2452	qs_seri(i,k) = qx(i,k,isol)
2453	IF (ok_ice_supersat) THEN
2454	cf_seri(i,k) = qx(i,k,icf)
2455	rvc_seri(i,k) = qx(i,k,irvc)
2456	ENDIF
2457	IF (ok_plane_contrail) THEN
2458	rcont_seri(i,k) = qx(i,k,ircont)
2459	ENDIF
2460	IF (ok_bs) THEN
2461	qbs_seri(i,k)= qx(i,k,ibs)
2462	ENDIF
2463	ENDIF
2464	ENDDO
2465	ENDDO
2466	!
2467	!--OB water mass fixer
2468	IF (ok_water_mass_fixer) THEN
2469	!--store initial water burden
2470	qql1(:)=0.0
2471	DO k = 1, klev
2472	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
2473	IF (nqo >= 3) THEN
2474	qql1(:)=qql1(:)+qs_seri(:,k)*zmasse(:,k)
2475	ENDIF
2476	IF (ok_bs) THEN
2477	qql1(:)=qql1(:)+qbs_seri(:,k)*zmasse(:,k)
2478	ENDIF
2479	ENDDO
2480	ENDIF
2481	!--fin mass fixer
2482
2483	tke0(:,:)=pbl_tke(:,:,is_ave)
2484	IF (nqtot > nqo) THEN
2485	! water isotopes are not included in tr_seri
2486	itr = 0
2487	DO iq = 1, nqtot
2488	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
2489	itr = itr+1
2490	DO k = 1, klev
2491	DO i = 1, klon
2492	tr_seri(i,k,itr) = qx(i,k,iq)
2493	ENDDO
2494	ENDDO
2495	ENDDO
2496	ELSE
2497	! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
2498	tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
2499	ENDIF
2500	!
2501	! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2502	! LF
2503	IF (debut) THEN
2504	WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
2505	itr = 0
2506	do iq = 1, nqtot
2507	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
2508	itr = itr+1
2509	tr_ancien(:,:,itr)=tr_seri(:,:,itr)
2510	enddo
2511	ENDIF
2512	!
2513	DO i = 1, klon
2514	ztsol(i) = 0.
2515	ENDDO
2516	DO nsrf = 1, nbsrf
2517	DO i = 1, klon
2518	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2519	ENDDO
2520	ENDDO
2521	! Initialize variables used for diagnostic purpose
2522	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2523
2524	! Diagnostiquer la tendance dynamique
2525	!
2526	IF (ancien_ok) THEN
2527	!
2528	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2529	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2530	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2531	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2532	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2533	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2534	d_qbs_dyn(:,:)= (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
2535	d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
2536	d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
2537	d_rcont_dyn(:,:)= (rcont_seri(:,:)-rcont_ancien(:,:))/phys_tstep
2538	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2539	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2540	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2541	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2542	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2543	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2544	CALL water_int(klon,klev,qbs_seri,zmasse,zx_tmp_fi2d)
2545	d_qbs_dyn2d(:)=(zx_tmp_fi2d(:)-prbsw_ancien(:))/phys_tstep
2546	! !! RomP >>> td dyn traceur
2547	IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
2548	! !! RomP <<<
2549	ELSE
2550	d_u_dyn(:,:) = 0.0
2551	d_v_dyn(:,:) = 0.0
2552	d_t_dyn(:,:) = 0.0
2553	d_q_dyn(:,:) = 0.0
2554	d_ql_dyn(:,:) = 0.0
2555	d_qs_dyn(:,:) = 0.0
2556	d_qbs_dyn(:,:)= 0.0
2557	d_cf_dyn(:,:) = 0.0
2558	d_rvc_dyn(:,:)= 0.0
2559	d_rcont_dyn(:,:)= 0.0
2560	d_q_dyn2d(:) = 0.0
2561	d_ql_dyn2d(:) = 0.0
2562	d_qs_dyn2d(:) = 0.0
2563	d_qbs_dyn2d(:)= 0.0
2564	! !! RomP >>> td dyn traceur
2565	IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
2566	! !! RomP <<<
2567	ancien_ok = .TRUE.
2568	ENDIF
2569	!
2570	! Ajouter le geopotentiel du sol:
2571	!
2572	DO k = 1, klev
2573	DO i = 1, klon
2574	zphi(i,k) = pphi(i,k) + pphis(i)
2575	ENDDO
2576	ENDDO
2577	!
2578	! Verifier les temperatures
2579	!
2580	!IM BEG
2581	IF (check) THEN
2582	amn=MIN(ftsol(1,is_ter),1000.)
2583	amx=MAX(ftsol(1,is_ter),-1000.)
2584	DO i=2, klon
2585	amn=MIN(ftsol(i,is_ter),amn)
2586	amx=MAX(ftsol(i,is_ter),amx)
2587	ENDDO
2588	!
2589	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2590	ENDIF !(check) THEN
2591	!IM END
2592	!
2593	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2594	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2595
2596	!
2597	!IM BEG
2598	IF (check) THEN
2599	amn=MIN(ftsol(1,is_ter),1000.)
2600	amx=MAX(ftsol(1,is_ter),-1000.)
2601	DO i=2, klon
2602	amn=MIN(ftsol(i,is_ter),amn)
2603	amx=MAX(ftsol(i,is_ter),amx)
2604	ENDDO
2605	!
2606	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2607	ENDIF !(check) THEN
2608	!IM END
2609	!
2610	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2611	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2612	!
2613	! Update ozone if day change
2614	IF (MOD(itap-1,lmt_pas) == 0) THEN
2615	IF (read_climoz <= 0) THEN
2616	! Once per day, update ozone from Royer:
2617	IF (solarlong0<-999.) then
2618	! Generic case with evolvoing season
2619	zzz=real(days_elapsed+1)
2620	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2621	! Particular case with annual mean insolation
2622	zzz=real(90) ! could be revisited
2623	IF (read_climoz/=-1) THEN
2624	abort_message ='read_climoz=-1 is recommended when ' &
2625	// 'solarlong0=1000.'
2626	CALL abort_physic (modname,abort_message,1)
2627	ENDIF
2628	ELSE
2629	! Case where the season is imposed with solarlong0
2630	zzz=real(90) ! could be revisited
2631	ENDIF
2632
2633	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2634	IF (CPPKEY_REPROBUS) THEN
2635	ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2636	DO i = 1, klon
2637	Z1=t_seri(i,itroprep(i)+1)
2638	Z2=t_seri(i,itroprep(i))
2639	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2640	B=Z2-fac*alog(pplay(i,itroprep(i)))
2641	ttrop(i)= fac*alog(ptrop(i))+B
2642	!
2643	Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2644	Z2= 1.e-3 * ( pphi(i,itroprep(i)) +pphis(i) ) / gravit
2645	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2646	B=Z2-fac*alog(pplay(i,itroprep(i)))
2647	ztrop(i)=fac*alog(ptrop(i))+B
2648	ENDDO
2649	END IF
2650	ELSE
2651	!--- ro3i = elapsed days number since current year 1st january, 0h
2652	ro3i=days_elapsed+jh_cur-jh_1jan
2653	!--- scaling for old style files (360 records)
2654	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2655	IF(adjust_tropopause) THEN
2656	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2657	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2658	time_climoz , longitude_deg, latitude_deg, &
2659	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2660	ELSE
2661	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2662	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2663	time_climoz )
2664	ENDIF
2665	! Convert from mole fraction of ozone to column density of ozone in a
2666	! cell, in kDU:
2667	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2668	* zmasse / dobson_u / 1e3
2669	! (By regridding ozone values for LMDZ only once a day, we
2670	! have already neglected the variation of pressure in one
2671	! day. So do not recompute "wo" at each time step even if
2672	! "zmasse" changes a little.)
2673	ENDIF
2674	ENDIF
2675
2676	!-- Needed for LSCP - condensation and ice supersaturation
2677	IF (ok_ice_supersat) THEN
2678	DO k = 1, klev
2679	DO i = 1, klon
2680	IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
2681	ratio_ql_qtot(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
2682	ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
2683	rvc_seri(i,k) = rvc_seri(i,k) * q_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
2684	ELSE
2685	ratio_ql_qtot(i,k) = 0.
2686	ratio_qi_qtot(i,k) = 0.
2687	rvc_seri(i,k) = 0.
2688	ENDIF
2689	ENDDO
2690	ENDDO
2691	ELSE
2692	DO k = 1, klev
2693	DO i = 1, klon
2694	IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
2695	ratio_ql_qtot(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
2696	ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
2697	ELSE
2698	ratio_ql_qtot(i,k) = 0.
2699	ratio_qi_qtot(i,k) = 0.
2700	ENDIF
2701	ENDDO
2702	ENDDO
2703	ENDIF
2704
2705	!
2706	! Re-evaporer l'eau liquide nuageuse
2707	!
2708	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2709	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2710
2711	CALL add_phys_tend &
2712	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,dqbs0,paprs,&
2713	'eva',abortphy,flag_inhib_tend,itap,0)
2714	CALL prt_enerbil('eva',itap)
2715
2716	!=========================================================================
2717	! Calculs de l'orbite.
2718	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2719	! doit donc etre plac\'e avant radlwsw et pbl_surface
2720
2721	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2722	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2723	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2724	!
2725	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2726	! solarlong0
2727	IF (solarlong0<-999.) THEN
2728	IF (new_orbit) THEN
2729	! calcul selon la routine utilisee pour les planetes
2730	CALL solarlong(day_since_equinox, zlongi, dist)
2731	ELSE
2732	! calcul selon la routine utilisee pour l'AR4
2733	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2734	ENDIF
2735	ELSE
2736	zlongi=solarlong0 ! longitude solaire vraie
2737	dist=1. ! distance au soleil / moyenne
2738	ENDIF
2739
2740	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2741
2742
2743	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2744	! Calcul de l'ensoleillement :
2745	! ============================
2746	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2747	! l'annee a partir d'une formule analytique.
2748	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2749	! non nul aux poles.
2750	IF (abs(solarlong0-1000.)<1.e-4) THEN
2751	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2752	latitude_deg,longitude_deg,rmu0,fract)
2753	swradcorr(:) = 1.0
2754	JrNt(:) = 1.0
2755	zrmu0(:) = rmu0(:)
2756	ELSE
2757	! recode par Olivier Boucher en sept 2015
2758	SELECT CASE (iflag_cycle_diurne)
2759	CASE(0)
2760	! Sans cycle diurne
2761	CALL angle(zlongi, latitude_deg, fract, rmu0)
2762	swradcorr = 1.0
2763	JrNt = 1.0
2764	zrmu0 = rmu0
2765	CASE(1)
2766	! Avec cycle diurne sans application des poids
2767	! bit comparable a l ancienne formulation cycle_diurne=true
2768	! on integre entre gmtime et gmtime+radpas
2769	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2770	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2771	latitude_deg,longitude_deg,rmu0,fract)
2772	zrmu0 = rmu0
2773	swradcorr = 1.0
2774	! Calcul du flag jour-nuit
2775	JrNt = 0.0
2776	WHERE (fract.GT.0.0) JrNt = 1.0
2777	CASE(2)
2778	! Avec cycle diurne sans application des poids
2779	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2780	! Comme cette routine est appele a tous les pas de temps de
2781	! la physique meme si le rayonnement n'est pas appele je
2782	! remonte en arriere les radpas-1 pas de temps
2783	! suivant. Petite ruse avec MOD pour prendre en compte le
2784	! premier pas de temps de la physique pendant lequel
2785	! itaprad=0
2786	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2787	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2788	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2789	latitude_deg,longitude_deg,rmu0,fract)
2790	!
2791	! Calcul des poids
2792	!
2793	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2794	zdtime2=0.0 !--pas de temps de la physique qui se termine
2795	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2796	latitude_deg,longitude_deg,zrmu0,zfract)
2797	swradcorr = 0.0
2798	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2799	swradcorr=zfract/fract*zrmu0/rmu0
2800	! Calcul du flag jour-nuit
2801	JrNt = 0.0
2802	WHERE (zfract.GT.0.0) JrNt = 1.0
2803	END SELECT
2804	ENDIF
2805	sza_o = ACOS (rmu0) *180./pi
2806
2807	IF (mydebug) THEN
2808	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2809	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2810	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2811	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2812	ENDIF
2813
2814	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2815	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2816	! Cela implique tous les interactions des sous-surfaces et la
2817	! partie diffusion turbulent du couche limit.
2818	!
2819	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2820	! ecriture des fihiers hist_XXXX.nc, ces sont :
2821	! qsol, zq2m, s_pblh, s_lcl,
2822	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2823	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2824	! zu10m, zv10m, fder,
2825	! zxqsurf, delta_qsurf,
2826	! rh2m, zxfluxu, zxfluxv,
2827	! frugs, agesno, fsollw, fsolsw,
2828	! d_ts, fevap, fluxlat, t2m,
2829	! wfbils, fluxt, fluxu, fluxv,
2830	!
2831	! Certains ne sont pas utiliser du tout :
2832	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2833	!
2834
2835	! Calcul de l'humidite de saturation au niveau du sol
2836
2837	! Tests Fredho, instensibilite au pas de temps -------------------------------
2838	! A detruire en 2024 une fois les tests documentes et les choix faits !
2839	! Conservation des variables avant l'appel à l a diffusion pour les tehrmic !
2840	if (iflag_thermals_tenv / 10 == 1 ) then !
2841	do k=1,klev !
2842	do i=1,klon !
2843	t_env(i,k)=t_seri(i,k) !
2844	q_env(i,k)=q_seri(i,k) !
2845	enddo !
2846	enddo !
2847	else if (iflag_thermals_tenv / 10 == 2 ) then !
2848	do k=1,klev !
2849	do i=1,klon !
2850	t_env(i,k)=t_seri(i,k) !
2851	enddo !
2852	enddo !
2853	endif !
2854	! Tests Fredho, instensibilite au pas de temps -------------------------------
2855
2856
2857	IF (iflag_pbl/=0) THEN
2858
2859	!jyg+nrlmd<
2860	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2861	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2862	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2863	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2864	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2865	ENDIF
2866	! !!
2867	!>jyg+nrlmd
2868	!
2869	!-------gustiness calculation-------!
2870	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2871	gustiness=0 !ym missing init
2872
2873	IF (iflag_gusts==0) THEN
2874	gustiness(1:klon)=0
2875	ELSE IF (iflag_gusts==1) THEN
2876	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2877	ELSE IF (iflag_gusts==2) THEN
2878	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2879	!!!! modif olivier torres
2880	ELSE IF (iflag_gusts==3) THEN
2881	w_et=wstar(1,3)
2882	jlr_g_s=(0.65w_et)*2
2883	pr_et=rain_con*8640
2884	jlr_g_c = (((19.8(pr_et(1:klon)2))/(1.5+pr_et(1:klon)+pr_et(1:klon)2))(0.4))*2
2885	gustiness(1:klon)=jlr_g_c+jlr_g_s
2886	!! write(,) "rain ",pr_et
2887	!! write(,) "jlr_g_c",jlr_g_c
2888	!! write(,) "wstar",wstar(1,3)
2889	!! write(,) "jlr_g_s",jlr_g_s
2890	! ELSE IF (iflag_gusts==2) THEN
2891	! do i = 1, klon
2892	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2893	! *ale_wake(i) !! need to make sigma_wk accessible here
2894	! enddo
2895	! ELSE IF (iflag_gusts==3) THEN
2896	! do i = 1, klon
2897	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2898	! enddo
2899	ENDIF
2900
2901	CALL pbl_surface( &
2902	phys_tstep, date0, itap, days_elapsed+1, &
2903	debut, lafin, &
2904	longitude_deg, latitude_deg, rugoro, zrmu0, &
2905	sollwdown, cldt, &
2906	rain_fall, snow_fall, bs_fall, solsw, solswfdiff, sollw, &
2907	gustiness, &
2908	t_seri, q_seri, qbs_seri, u_seri, v_seri, &
2909	!nrlmd+jyg<
2910	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2911	!>nrlmd+jyg
2912	pplay, paprs, pctsrf, &
2913	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2914	!albedo SB <<<
2915	cdragh, cdragm, u1, v1, &
2916	beta_aridity, &
2917	!albedo SB >>>
2918	! albsol1, albsol2, sens, evap, &
2919	albsol_dir, albsol_dif, sens, evap, snowerosion, &
2920	!albedo SB <<<
2921	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2922	zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, &
2923	d_t_vdf, d_q_vdf, d_qbs_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2924	!nrlmd<
2925	!jyg<
2926	d_t_vdf_w, d_q_vdf_w, &
2927	d_t_vdf_x, d_q_vdf_x, &
2928	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2929	!>jyg
2930	delta_tsurf,wake_dens, &
2931	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2932	kh,kh_x,kh_w, &
2933	!>nrlmd
2934	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2935	slab_wfbils, &
2936	qsol, zq2m, s_pblh, s_lcl, &
2937	!jyg<
2938	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2939	!>jyg
2940	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2941	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2942	zustar, zu10m, zv10m, fder, &
2943	zxqsurf, delta_qsurf, rh2m, zxfluxu, zxfluxv, &
2944	z0m, z0h, agesno, fsollw, fsolsw, &
2945	d_ts, fevap, fluxlat, t2m, &
2946	wfbils, wfevap, &
2947	fluxt, fluxu, fluxv, &
2948	dsens, devap, zxsnow, &
2949	zxfluxt, zxfluxq, zxfluxqbs, q2m, fluxq, fluxqbs, pbl_tke, pbl_eps, &
2950	!nrlmd+jyg<
2951	wake_delta_pbl_TKE, &
2952	!>nrlmd+jyg
2953	treedrg )
2954	!FC
2955	!
2956	! Add turbulent diffusion tendency to the wake difference variables
2957	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2958	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2959	!jyg<
2960	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2961	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2962	CALL add_wake_tend &
2963	(d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2964	ELSE
2965	d_deltat_vdf(:,:) = 0.
2966	d_deltaq_vdf(:,:) = 0.
2967	!>jyg
2968	ENDIF
2969
2970	!---------------------------------------------------------------------
2971	! ajout des tendances de la diffusion turbulente
2972	IF (klon_glo==1) THEN
2973	CALL add_pbl_tend &
2974	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
2975	'vdf',abortphy,flag_inhib_tend,itap)
2976	ELSE
2977	CALL add_phys_tend &
2978	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
2979	'vdf',abortphy,flag_inhib_tend,itap,0)
2980	ENDIF
2981	CALL prt_enerbil('vdf',itap)
2982
2983	!--------------------------------------------------------------------
2984
2985	IF (mydebug) THEN
2986	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2987	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2988	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2989	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2990	ENDIF
2991
2992	!albedo SB >>>
2993	albsol1=0.
2994	albsol2=0.
2995	falb1=0.
2996	falb2=0.
2997	SELECT CASE(nsw)
2998	CASE(2)
2999	albsol1=albsol_dir(:,1)
3000	albsol2=albsol_dir(:,2)
3001	falb1=falb_dir(:,1,:)
3002	falb2=falb_dir(:,2,:)
3003	CASE(4)
3004	albsol1=albsol_dir(:,1)
3005	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
3006	+albsol_dir(:,4)*SFRWL(4)
3007	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
3008	falb1=falb_dir(:,1,:)
3009	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
3010	+falb_dir(:,4,:)*SFRWL(4)
3011	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
3012	CASE(6)
3013	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
3014	+albsol_dir(:,3)*SFRWL(3)
3015	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
3016	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
3017	+albsol_dir(:,6)*SFRWL(6)
3018	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
3019	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
3020	+falb_dir(:,3,:)*SFRWL(3)
3021	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
3022	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
3023	+falb_dir(:,6,:)*SFRWL(6)
3024	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
3025	END SELECt
3026	!albedo SB <<<
3027
3028
3029	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
3030	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
3031
3032	ENDIF
3033
3034	! ==================================================================
3035	! Blowing snow sublimation and sedimentation
3036
3037	d_t_bsss(:,:)=0.
3038	d_q_bsss(:,:)=0.
3039	d_qbs_bsss(:,:)=0.
3040	bsfl(:,:)=0.
3041	bs_fall(:)=0.
3042	IF (ok_bs) THEN
3043
3044	CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
3045	d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
3046
3047	CALL add_phys_tend &
3048	(du0,dv0,d_t_bsss,d_q_bsss,dql0,dqi0,d_qbs_bsss,paprs,&
3049	'bsss',abortphy,flag_inhib_tend,itap,0)
3050
3051	ENDIF
3052
3053	! =================================================================== c
3054	! Calcul de Qsat
3055
3056	DO k = 1, klev
3057	DO i = 1, klon
3058	zx_t = t_seri(i,k)
3059	IF (thermcep) THEN
3060	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3061	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3062	zx_qs = MIN(0.5,zx_qs)
3063	zcor = 1./(1.-retv*zx_qs)
3064	zx_qs = zx_qs*zcor
3065	ELSE
3066	!! IF (zx_t.LT.t_coup) THEN !jyg
3067	IF (zx_t.LT.rtt) THEN !jyg
3068	zx_qs = qsats(zx_t)/pplay(i,k)
3069	ELSE
3070	zx_qs = qsatl(zx_t)/pplay(i,k)
3071	ENDIF
3072	ENDIF
3073	zqsat(i,k)=zx_qs
3074	ENDDO
3075	ENDDO
3076
3077	IF (prt_level.ge.1) THEN
3078	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
3079	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
3080	ENDIF
3081	!
3082	! Appeler la convection (au choix)
3083	!
3084	DO k = 1, klev
3085	DO i = 1, klon
3086	conv_q(i,k) = d_q_dyn(i,k) &
3087	+ d_q_vdf(i,k)/phys_tstep
3088	conv_t(i,k) = d_t_dyn(i,k) &
3089	+ d_t_vdf(i,k)/phys_tstep
3090	ENDDO
3091	ENDDO
3092
3093	! Calcule de vitesse verticale a partir de flux de masse verticale
3094	DO k = 1, klev
3095	DO i = 1, klon
3096	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
3097	ENDDO
3098	ENDDO
3099
3100	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
3101	omega(igout, :)
3102	!
3103	! Appel de la convection tous les "cvpas"
3104	!
3105	!!jyg IF (MOD(itapcv,cvpas).EQ.0) THEN
3106	!! print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
3107	!! itapcv, cvpas, itap-1, cvpas_0
3108	IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
3109
3110	!
3111	! Mettre a zero des variables de sortie (pour securite)
3112	!
3113	pmflxr(:,:) = 0.
3114	pmflxs(:,:) = 0.
3115	wdtrainA(:,:) = 0.
3116	wdtrainS(:,:) = 0.
3117	wdtrainM(:,:) = 0.
3118	wdtrainAS(:,:) = 0.
3119	upwd(:,:) = 0.
3120	dnwd(:,:) = 0.
3121	ep(:,:) = 0.
3122	da(:,:)=0.
3123	mp(:,:)=0.
3124	wght_cvfd(:,:)=0.
3125	phi(:,:,:)=0.
3126	phi2(:,:,:)=0.
3127	epmlmMm(:,:,:)=0.
3128	eplaMm(:,:)=0.
3129	d1a(:,:)=0.
3130	dam(:,:)=0.
3131	elij(:,:,:)=0.
3132	ev(:,:)=0.
3133	qtaa(:,:)=0.
3134	clw(:,:)=0.
3135	sij(:,:,:)=0.
3136	!
3137	IF (iflag_con.EQ.1) THEN
3138	abort_message ='reactiver le call conlmd dans physiq.F'
3139	CALL abort_physic (modname,abort_message,1)
3140	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
3141	! . d_t_con, d_q_con,
3142	! . rain_con, snow_con, ibas_con, itop_con)
3143	ELSE IF (iflag_con.EQ.2) THEN
3144	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
3145	conv_t, conv_q, -evap, omega, &
3146	d_t_con, d_q_con, rain_con, snow_con, &
3147	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3148	kcbot, kctop, kdtop, pmflxr, pmflxs)
3149	d_u_con = 0.
3150	d_v_con = 0.
3151
3152	WHERE (rain_con < 0.) rain_con = 0.
3153	WHERE (snow_con < 0.) snow_con = 0.
3154	DO i = 1, klon
3155	ibas_con(i) = klev+1 - kcbot(i)
3156	itop_con(i) = klev+1 - kctop(i)
3157	ENDDO
3158	ELSE IF (iflag_con.GE.3) THEN
3159	! nb of tracers for the KE convection:
3160	! MAF la partie traceurs est faite dans phytrac
3161	! on met ntra=1 pour limiter les appels mais on peut
3162	! supprimer les calculs / ftra.
3163	ntra = 1
3164
3165	!=======================================================================
3166	!ajout pour la parametrisation des poches froides: calcul de
3167	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
3168	IF (iflag_wake>=1) THEN
3169	DO k=1,klev
3170	DO i=1,klon
3171	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
3172	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
3173	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3174	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3175	ENDDO
3176	ENDDO
3177	ELSE
3178	t_w(:,:) = t_seri(:,:)
3179	q_w(:,:) = q_seri(:,:)
3180	t_x(:,:) = t_seri(:,:)
3181	q_x(:,:) = q_seri(:,:)
3182	ENDIF
3183	!
3184	!jyg<
3185	! Perform dry adiabatic adjustment on wake profile
3186	! The corresponding tendencies are added to the convective tendencies
3187	! after the call to the convective scheme.
3188	IF (iflag_wake>=1) then
3189	IF (iflag_adjwk >= 1) THEN
3190	limbas(:) = 1
3191	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
3192	d_t_adjwk, d_q_adjwk)
3193	!
3194	DO k=1,klev
3195	DO i=1,klon
3196	IF (wake_s(i) .GT. 1.e-3) THEN
3197	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
3198	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
3199	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
3200	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
3201	ELSE
3202	d_deltat_ajs_cv(i,k) = 0.
3203	d_deltaq_ajs_cv(i,k) = 0.
3204	ENDIF
3205	ENDDO
3206	ENDDO
3207	IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
3208	CALL add_wake_tend &
3209	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
3210	ENDIF ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
3211	ENDIF ! (iflag_adjwk >= 1)
3212	ENDIF ! (iflag_wake>=1)
3213	!>jyg
3214	!
3215
3216	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
3217	!! (k, q_w(1,k), q_x(1,k),k=1,25)
3218
3219	!jyg<
3220	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
3221	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
3222	ale_bl_prescr, alp_bl_prescr, &
3223	wake_pe, wake_fip, &
3224	Ale_bl, Ale_bl_trig, Alp_bl, &
3225	Ale, Alp , Ale_wake, Alp_wake)
3226	!>jyg
3227	!
3228	! sb, oct02:
3229	! Schema de convection modularise et vectorise:
3230	! (driver commun aux versions 3 et 4)
3231	!
3232	IF (ok_cvl) THEN ! new driver for convectL
3233	!
3234	!jyg<
3235	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3236	! Calculate the upmost level of deep convection loops: k_upper_cv
3237	! (near 22 km)
3238	k_upper_cv = klev
3239	!izero = klon/2+1/klon
3240	!DO k = klev,1,-1
3241	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
3242	!ENDDO
3243	! FH : nouveau calcul base sur un profil global sans quoi
3244	! le modele etait sensible au decoupage de domaines
3245	DO k = klev,1,-1
3246	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
3247	ENDDO
3248	IF (prt_level .ge. 5) THEN
3249	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
3250	k_upper_cv
3251	ENDIF
3252	!
3253	!>jyg
3254	IF (type_trac == 'repr') THEN
3255	nbtr_tmp=ntra
3256	ELSE
3257	nbtr_tmp=nbtr
3258	ENDIF
3259	!jyg iflag_con est dans clesphys
3260	!c CALL concvl (iflag_con,iflag_clos,
3261	CALL concvl (iflag_clos, &
3262	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
3263	t_w,q_w,wake_s, &
3264	u_seri,v_seri,tr_seri,nbtr_tmp, &
3265	ALE,ALP, &
3266	sig1,w01, &
3267	d_t_con,d_q_con,fqcomp,d_u_con,d_v_con,d_tr, &
3268	rain_con, snow_con, ibas_con, itop_con, sigd, &
3269	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
3270	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
3271	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
3272	! RomP >>>
3273	!! . pmflxr,pmflxs,da,phi,mp,
3274	!! . ftd,fqd,lalim_conv,wght_th)
3275	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
3276	ftd,fqd,lalim_conv,wght_th, &
3277	ev, ep,epmlmMm,eplaMm, &
3278	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv,detrain_cv, &
3279	tau_cld_cv,coefw_cld_cv,epmax_diag)
3280
3281	! RomP <<<
3282
3283	!IM begin
3284	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
3285	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
3286	!IM end
3287	!IM cf. FH
3288	clwcon0=qcondc
3289	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
3290	fm_cv(:,:)=upwd(:,:)+dnwd(:,:)+dnwd0(:,:)
3291	!
3292	!jyg<
3293	! If convective tendencies are too large, then call convection
3294	! every time step
3295	cvpas = cvpas_0
3296	DO k=1,k_upper_cv
3297	DO i=1,klon
3298	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3299	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3300	dtcon_multistep_max = 3.
3301	dqcon_multistep_max = 0.02
3302	ENDIF
3303	ENDDO
3304	ENDDO
3305	!
3306	DO k=1,k_upper_cv
3307	DO i=1,klon
3308	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3309	!! abs(d_q_con(i,k)) > 2.e-2) THEN
3310	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3311	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3312	cvpas = 1
3313	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3314	!! i,k,d_t_con(i,k),d_q_con(i,k)
3315	ENDIF
3316	ENDDO
3317	ENDDO
3318	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3319	!!! call bcast(cvpas)
3320	!!! ------------------------------------------------------------
3321	!>jyg
3322	!
3323	DO i = 1, klon
3324	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
3325	ENDDO
3326	!
3327	!jyg<
3328	! Add the tendency due to the dry adjustment of the wake profile
3329	IF (iflag_wake>=1) THEN
3330	IF (iflag_adjwk == 2) THEN
3331	DO k=1,klev
3332	DO i=1,klon
3333	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3334	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
3335	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3336	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3337	ENDDO
3338	ENDDO
3339	ENDIF ! (iflag_adjwk = 2)
3340	ENDIF ! (iflag_wake>=1)
3341	!>jyg
3342	!
3343	ELSE ! ok_cvl
3344
3345	! MAF conema3 ne contient pas les traceurs
3346	CALL conema3 (phys_tstep, &
3347	paprs,pplay,t_seri,q_seri, &
3348	u_seri,v_seri,tr_seri,ntra, &
3349	sig1,w01, &
3350	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3351	rain_con, snow_con, ibas_con, itop_con, &
3352	upwd,dnwd,dnwd0,bas,top, &
3353	Ma,cape,tvp,rflag, &
3354	pbase &
3355	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3356	,clwcon0)
3357
3358	ENDIF ! ok_cvl
3359
3360	!
3361	! Correction precip
3362	rain_con = rain_con * cvl_corr
3363	snow_con = snow_con * cvl_corr
3364	!
3365
3366	IF (.NOT. ok_gust) THEN
3367	do i = 1, klon
3368	wd(i)=0.0
3369	enddo
3370	ENDIF
3371
3372	! =================================================================== c
3373	! Calcul des proprietes des nuages convectifs
3374	!
3375
3376	! calcul des proprietes des nuages convectifs
3377	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3378	IF (iflag_cld_cv == 0) THEN
3379	CALL clouds_gno &
3380	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3381	ELSE
3382	CALL clouds_bigauss &
3383	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3384	ENDIF
3385
3386
3387	! =================================================================== c
3388
3389	DO i = 1, klon
3390	itop_con(i) = min(max(itop_con(i),1),klev)
3391	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3392	ENDDO
3393
3394	DO i = 1, klon
3395	! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3396	! où i n'est pas un point convectif et donc ibas_con(i)=0
3397	! c'est un pb indépendant des isotopes
3398	if (ibas_con(i) > 0) then
3399	ema_pcb(i) = paprs(i,ibas_con(i))
3400	else
3401	ema_pcb(i) = 0.0
3402	endif
3403	ENDDO
3404	DO i = 1, klon
3405	! L'idicage de itop_con peut cacher un pb potentiel
3406	! FH sous la dictee de JYG, CR
3407	ema_pct(i) = paprs(i,itop_con(i)+1)
3408
3409	IF (itop_con(i).gt.klev-3) THEN
3410	IF (prt_level >= 9) THEN
3411	write(lunout,*)'La convection monte trop haut '
3412	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
3413	ENDIF
3414	ENDIF
3415	ENDDO
3416	ELSE IF (iflag_con.eq.0) THEN
3417	write(lunout,*) 'On n appelle pas la convection'
3418	clwcon0=0.
3419	rnebcon0=0.
3420	d_t_con=0.
3421	d_q_con=0.
3422	d_u_con=0.
3423	d_v_con=0.
3424	rain_con=0.
3425	snow_con=0.
3426	bas=1
3427	top=1
3428	ELSE
3429	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
3430	CALL abort_physic("physiq", "", 1)
3431	ENDIF
3432
3433	!--saving d_q_con * zmass for next timestep if convection is not called every timestep
3434	IF (ok_conserv_d_q_con) THEN
3435	d_q_con_zmasse(:,:) = d_q_con(:,:) * zmasse(:,:)
3436	ENDIF
3437
3438	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3439	! . d_u_con, d_v_con)
3440
3441	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3442	proba_notrig(:) = 1.
3443	itapcv = 0
3444	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3445	!
3446	itapcv = itapcv+1
3447	!
3448	! Compter les steps ou cvpas=1
3449	IF (cvpas == 1) THEN
3450	Ncvpaseq1 = Ncvpaseq1+1
3451	ENDIF
3452	IF (mod(itap,1000) == 0) THEN
3453	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3454	ENDIF
3455
3456	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3457	!!! l'energie dans les courants satures.
3458	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3459	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3460	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3461	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3462	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3463	!! itap, 1)
3464	!! call prt_enerbil('convection_sat',itap)
3465	!!
3466	!!
3467
3468	!--recompute d_q_con with zmasse from new timestep
3469	IF (ok_conserv_d_q_con) THEN
3470	d_q_con(:,:)=d_q_con_zmasse(:,:)/zmasse(:,:)
3471	ENDIF
3472
3473	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, dqbs0, paprs, &
3474	'convection',abortphy,flag_inhib_tend,itap,0)
3475	CALL prt_enerbil('convection',itap)
3476
3477	!-------------------------------------------------------------------------
3478
3479	IF (mydebug) THEN
3480	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3481	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3482	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3483	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3484	ENDIF
3485
3486	!
3487	!==========================================================================
3488	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3489	!pour la couche limite diffuse pour l instant
3490	!
3491	!
3492	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3493	! il faut rajouter cette tendance calcul\'ee hors des poches
3494	! froides
3495	!
3496	IF (iflag_wake>=1) THEN
3497	!
3498	!
3499	! Call wakes every "wkpas" step
3500	!
3501	IF (MOD(itapwk,wkpas).EQ.0) THEN
3502	!
3503	DO k=1,klev
3504	DO i=1,klon
3505	dt_dwn(i,k) = ftd(i,k)
3506	dq_dwn(i,k) = fqd(i,k)
3507	M_dwn(i,k) = dnwd0(i,k)
3508	M_up(i,k) = upwd(i,k)
3509	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3510	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3511	ENDDO
3512	ENDDO
3513
3514	IF (iflag_wake==2) THEN
3515	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3516	DO k = 1,klev
3517	dt_dwn(:,k)= dt_dwn(:,k)+ &
3518	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3519	dq_dwn(:,k)= dq_dwn(:,k)+ &
3520	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3521	ENDDO
3522	ELSEIF (iflag_wake==3) THEN
3523	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3524	DO k = 1,klev
3525	DO i=1,klon
3526	IF (rneb(i,k)==0.) THEN
3527	! On ne tient compte des tendances qu'en dehors des
3528	! nuages (c'est-\`a-dire a priri dans une region ou
3529	! l'eau se reevapore).
3530	dt_dwn(i,k)= dt_dwn(i,k)+ &
3531	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3532	dq_dwn(i,k)= dq_dwn(i,k)+ &
3533	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3534	ENDIF
3535	ENDDO
3536	ENDDO
3537	ENDIF
3538
3539	!
3540	!calcul caracteristiques de la poche froide
3541	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3542	t_seri, q_seri, omega, &
3543	dt_dwn, dq_dwn, M_dwn, M_up, &
3544	dt_a, dq_a, cv_gen, &
3545	sigd, cin, &
3546	wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, &
3547	wake_dth, wake_h, &
3548	!! wake_pe, wake_fip, wake_gfl, &
3549	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3550	d_t_wake, d_q_wake, &
3551	wake_k, t_x, q_x, &
3552	wake_omgbdth, wake_dp_omgb, &
3553	wake_dtKE, wake_dqKE, &
3554	wake_omg, wake_dp_deltomg, &
3555	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3556	d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk)
3557	!
3558	!jyg Reinitialize itapwk when wakes have been called
3559	itapwk = 0
3560	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3561	!
3562	itapwk = itapwk+1
3563	!
3564	!-----------------------------------------------------------------------
3565	! ajout des tendances des poches froides
3566	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,dqbs0,paprs,'wake', &
3567	abortphy,flag_inhib_tend,itap,0)
3568	CALL prt_enerbil('wake',itap)
3569	!------------------------------------------------------------------------
3570
3571	! Increment Wake state variables
3572	IF (iflag_wake_tend .GT. 0.) THEN
3573
3574	CALL add_wake_tend &
3575	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &
3576	'wake', abortphy)
3577	CALL prt_enerbil('wake',itap)
3578	ENDIF ! (iflag_wake_tend .GT. 0.)
3579	!
3580	IF (prt_level .GE. 10) THEN
3581	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3582	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3583	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3584	ENDIF
3585
3586	IF (iflag_alp_wk_cond .GT. 0.) THEN
3587
3588	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3589	wake_fip)
3590	ELSE
3591	wake_fip(:) = wake_fip_0(:)
3592	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3593
3594	ENDIF ! (iflag_wake>=1)
3595	!
3596	!===================================================================
3597	! Convection seche (thermiques ou ajustement)
3598	!===================================================================
3599	!
3600	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3601	,seuil_inversion,weak_inversion,dthmin)
3602
3603
3604
3605	d_t_ajsb(:,:)=0.
3606	d_q_ajsb(:,:)=0.
3607	d_t_ajs(:,:)=0.
3608	d_u_ajs(:,:)=0.
3609	d_v_ajs(:,:)=0.
3610	d_q_ajs(:,:)=0.
3611	clwcon0th(:,:)=0.
3612	!
3613	! fm_therm(:,:)=0.
3614	! entr_therm(:,:)=0.
3615	! detr_therm(:,:)=0.
3616	!
3617	IF (prt_level>9) WRITE(lunout,*) &
3618	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3619	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3620	IF (iflag_thermals<0) THEN
3621	! Rien
3622	! ====
3623	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3624	WRITE(lunout,*) 'WARNING : running without dry convection. Somme intermediate variables are not properly defined in physiq_mod.F90'
3625	! Reprendre proprement les initialisation ci dessouds si on veut vraiment utiliser l'option (FH)
3626	fraca(:,:)=0.
3627	fm_therm(:,:)=0.
3628	ztv(:,:)=t_seri(:,:)
3629	zqasc(:,:)=q_seri(:,:)
3630	ztla(:,:)=0.
3631	zthl(:,:)=0.
3632	zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
3633
3634
3635
3636	ELSE
3637
3638	! Thermiques
3639	! ==========
3640	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3641	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3642
3643
3644	!cc nrlmd le 10/04/2012
3645	DO k=1,klev+1
3646	DO i=1,klon
3647	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3648	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3649	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3650	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3651	ENDDO
3652	ENDDO
3653	!cc fin nrlmd le 10/04/2012
3654
3655	IF (iflag_thermals>=1) THEN
3656
3657	! Tests Fredho, instensibilite au pas de temps -------------------------------
3658	! A detruire en 2024 une fois les tests documentes et les choix faits !
3659	if (iflag_thermals_tenv /10 == 0 ) then !
3660	do k=1,klev !
3661	do i=1,klon !
3662	t_env(i,k)=t_seri(i,k) !
3663	q_env(i,k)=q_seri(i,k) !
3664	enddo !
3665	enddo !
3666	else if (iflag_thermals_tenv / 10 == 2 ) then !
3667	do k=1,klev !
3668	do i=1,klon !
3669	q_env(i,k)=q_seri(i,k) !
3670	enddo !
3671	enddo !
3672	else if (iflag_thermals_tenv / 10 == 3 ) then !
3673	do k=1,klev !
3674	do i=1,klon !
3675	t_env(i,k)=t(i,k) !
3676	q_env(i,k)=qx(i,k,1) !
3677	enddo !
3678	enddo !
3679	endif !
3680	! Tests Fredho, instensibilite au pas de temps ------------------------------
3681
3682	!jyg<
3683	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3684	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3685	! Appel des thermiques avec les profils exterieurs aux poches
3686	DO k=1,klev
3687	DO i=1,klon
3688	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3689	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3690	t_env(i,k) = t_env(i,k) - wake_s(i)*wake_deltat(i,k)
3691	q_env(i,k) = q_env(i,k) - wake_s(i)*wake_deltaq(i,k)
3692	u_therm(i,k) = u_seri(i,k)
3693	v_therm(i,k) = v_seri(i,k)
3694	ENDDO
3695	ENDDO
3696	ELSE
3697	! Appel des thermiques avec les profils moyens
3698	DO k=1,klev
3699	DO i=1,klon
3700	t_therm(i,k) = t_seri(i,k)
3701	q_therm(i,k) = q_seri(i,k)
3702	u_therm(i,k) = u_seri(i,k)
3703	v_therm(i,k) = v_seri(i,k)
3704	ENDDO
3705	ENDDO
3706	ENDIF
3707	!>jyg
3708	CALL calltherm(pdtphys &
3709	,pplay,paprs,pphi,weak_inversion &
3710	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3711	,u_therm,v_therm,t_therm,q_therm,t_env,q_env,zqsat,debut & !jyg
3712	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3713	,fm_therm,entr_therm,detr_therm &
3714	,zqasc,clwcon0th,lmax_th,ratqscth &
3715	,ratqsdiff,zqsatth &
3716	!on rajoute ale et alp, et les
3717	!caracteristiques de la couche alim
3718	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3719	,ztv,zpspsk,ztla,zthl &
3720	!cc nrlmd le 10/04/2012
3721	,pbl_tke_input,pctsrf,omega,cell_area &
3722	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3723	,n2,s2,strig,zcong,ale_bl_stat &
3724	,therm_tke_max,env_tke_max &
3725	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3726	,alp_bl_conv,alp_bl_stat &
3727	!cc fin nrlmd le 10/04/2012
3728	,zqla,ztva )
3729	!
3730	!jyg<
3731	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3732	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3733	! Si les thermiques ne sont presents que hors des
3734	! poches, la tendance moyenne associ\'ee doit etre
3735	! multipliee par la fraction surfacique qu'ils couvrent.
3736	DO k=1,klev
3737	DO i=1,klon
3738	!
3739	d_deltat_the(i,k) = - d_t_ajs(i,k)
3740	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3741	!
3742	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3743	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3744	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3745	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3746	!
3747	ENDDO
3748	ENDDO
3749	!
3750	IF (ok_bug_split_th) THEN
3751	CALL add_wake_tend &
3752	(d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3753	ELSE
3754	CALL add_wake_tend &
3755	(d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3756	ENDIF
3757	CALL prt_enerbil('the',itap)
3758	!
3759	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3760	!
3761	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3762	dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3763	CALL prt_enerbil('thermals',itap)
3764	!
3765	!
3766	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3767	cin, s2, n2, strig, &
3768	ale_bl_trig, ale_bl_stat, ale_bl, &
3769	alp_bl, alp_bl_stat, &
3770	proba_notrig, random_notrig, cv_gen)
3771	!>jyg
3772
3773	! ------------------------------------------------------------------
3774	! Transport de la TKE par les panaches thermiques.
3775	! FH : 2010/02/01
3776	if (iflag_thermcell_tke==1) then
3777	call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,rg,paprs,pbl_tke)
3778	endif
3779	! -------------------------------------------------------------------
3780
3781	DO i=1,klon
3782	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3783	!CR:04/05/12:correction calcul zmax
3784	zmax_th(i)=zmax0(i)
3785	ENDDO
3786
3787	ENDIF
3788
3789	! Ajustement sec
3790	! ==============
3791
3792	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3793	! a partir du sommet des thermiques.
3794	! Dans le cas contraire, on demarre au niveau 1.
3795
3796	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3797
3798	IF (iflag_thermals.eq.0) THEN
3799	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3800	limbas(:)=1
3801	ELSE
3802	limbas(:)=lmax_th(:)
3803	ENDIF
3804
3805	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3806	! pour des test de convergence numerique.
3807	! Le nouveau ajsec est a priori mieux, meme pour le cas
3808	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3809	! non nulles numeriquement pour des mailles non concernees.
3810
3811	IF (iflag_thermals==0) THEN
3812	! Calling adjustment alone (but not the thermal plume model)
3813	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3814	, d_t_ajsb, d_q_ajsb)
3815	ELSE IF (iflag_thermals>0) THEN
3816	! Calling adjustment above the top of thermal plumes
3817	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3818	, d_t_ajsb, d_q_ajsb)
3819	ENDIF
3820
3821	!--------------------------------------------------------------------
3822	! ajout des tendances de l'ajustement sec ou des thermiques
3823	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,dqbs0,paprs, &
3824	'ajsb',abortphy,flag_inhib_tend,itap,0)
3825	CALL prt_enerbil('ajsb',itap)
3826	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3827	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3828
3829	!---------------------------------------------------------------------
3830
3831	ENDIF
3832
3833	ENDIF
3834	!
3835	!===================================================================
3836	! Computation of subrgid scale near-surface wind distribution
3837	call surf_wind(klon,nsrfwnd,u10m,v10m,wake_s,wake_Cstar,ustar,wstar,surf_wind_value,surf_wind_proba)
3838
3839	!===================================================================
3840	! Computation of ratqs, the width (normalized) of the subrid scale
3841	! water distribution
3842
3843	l_mix_ave(:,:)=0.
3844	wprime_ave(:,:)=0.
3845
3846	DO nsrf = 1, nbsrf
3847	DO i = 1, klon
3848	l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3849	wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3850	ENDDO
3851	ENDDO
3852
3853	CALL ratqs_main(klon,klev,nbsrf,prt_level,lunout, &
3854	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3855	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3856	pctsrf,s_pblh,zstd, &
3857	tau_ratqs,fact_cldcon,wake_s, wake_deltaq, &
3858	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3859	paprs,pplay,t_seri,q_seri, &
3860	qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
3861	omega,pbl_tke(:,:,is_ave),pbl_eps(:,:,is_ave),l_mix_ave,wprime_ave, &
3862	t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
3863	ratqs,ratqsc,ratqs_inter_,sigma_qtherm)
3864
3865	!
3866	! Appeler le processus de condensation a grande echelle
3867	! et le processus de precipitation
3868	!-------------------------------------------------------------------------
3869	IF (prt_level .GE.10) THEN
3870	print *,'itap, ->fisrtilp ',itap
3871	ENDIF
3872	!
3873
3874	picefra(:,:)=0.
3875
3876	IF (ok_new_lscp) THEN
3877
3878
3879	DO k = 1, klev
3880	DO i = 1, klon
3881	ql_seri_lscp(i,k) = ratio_ql_qtot(i,k) * q_seri(i,k)
3882	qi_seri_lscp(i,k) = ratio_qi_qtot(i,k) * q_seri(i,k)
3883	ENDDO
3884	ENDDO
3885
3886
3887	!--mise à jour de flight_m et flight_h2o dans leur module
3888	!IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3889	! CALL airplane(debut,pphis,pplay,paprs,t_seri)
3890	!ENDIF
3891
3892	CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
3893	t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ptconv, ratqs, sigma_qtherm, &
3894	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
3895	pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
3896	radocond, picefra, rain_lsc, snow_lsc, &
3897	frac_impa, frac_nucl, beta_prec_fisrt, &
3898	prfl, psfl, rhcl, &
3899	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3900	iflag_ice_thermo, distcltop, temp_cltop, &
3901	pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
3902	cell_area, &
3903	cf_seri, rvc_seri, u_seri, v_seri, &
3904	qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
3905	dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
3906	dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
3907	rcont_seri, flight_dist, flight_h2o, &
3908	Tcritcont, qcritcont, potcontfraP, potcontfraNP, dcf_avi, dqi_avi, dqvc_avi, &
3909	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
3910	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
3911	dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
3912
3913
3914	ELSE
3915
3916	CALL fisrtilp(klon,klev,phys_tstep,paprs,pplay, &
3917	t_seri, q_seri,ptconv,ratqs,sigma_qtherm, &
3918	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, radocond, &
3919	rain_lsc, snow_lsc, &
3920	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3921	frac_impa, frac_nucl, beta_prec_fisrt, &
3922	prfl, psfl, rhcl, &
3923	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3924	iflag_ice_thermo, &
3925	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
3926
3927	ENDIF
3928	!
3929	WHERE (rain_lsc < 0) rain_lsc = 0.
3930	WHERE (snow_lsc < 0) snow_lsc = 0.
3931
3932	!+JLD
3933	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3934	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3935	! & ,rain_lsc,snow_lsc
3936	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3937	!-JLD
3938	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
3939	'lsc',abortphy,flag_inhib_tend,itap,0)
3940	CALL prt_enerbil('lsc',itap)
3941	rain_num(:)=0.
3942	DO k = 1, klev
3943	DO i = 1, klon
3944	IF (ql_seri(i,k)>oliqmax) THEN
3945	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3946	ql_seri(i,k)=oliqmax
3947	ENDIF
3948	ENDDO
3949	ENDDO
3950	IF (nqo >= 3) THEN
3951	DO k = 1, klev
3952	DO i = 1, klon
3953	IF (qs_seri(i,k)>oicemax) THEN
3954	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3955	qs_seri(i,k)=oicemax
3956	ENDIF
3957	ENDDO
3958	ENDDO
3959	ENDIF
3960
3961
3962	!---------------------------------------------------------------------------
3963	DO k = 1, klev
3964	DO i = 1, klon
3965	cldfra(i,k) = rneb(i,k)
3966	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3967	!EV: en effet etrange, j'ajouterais aussi qs_seri
3968	! plus largement, je nettoierais (enleverrais) ces lignes
3969	IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
3970	ENDDO
3971	ENDDO
3972
3973
3974	! Option to activate the radiative effect of blowing snow (ok_rad_bs)
3975	! makes sense only if the new large scale condensation scheme is active
3976	! with the ok_icefra_lscp flag active as well
3977
3978	IF (ok_bs .AND. ok_rad_bs) THEN
3979	IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
3980	DO k=1,klev
3981	DO i=1,klon
3982	radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
3983	picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
3984	qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
3985	cldfra(i,k)=max(cldfra(i,k),qbsfra)
3986	ENDDO
3987	ENDDO
3988	ELSE
3989	WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
3990	WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
3991	abort_message='inconsistency in cloud phase for blowing snow'
3992	CALL abort_physic(modname,abort_message,1)
3993	ENDIF
3994
3995	ENDIF
3996
3997	IF (mydebug) THEN
3998	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3999	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4000	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4001	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4002	ENDIF
4003
4004	!
4005	!-------------------------------------------------------------------
4006	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
4007	!-------------------------------------------------------------------
4008
4009	! 1. NUAGES CONVECTIFS
4010	!
4011	!IM cf FH
4012	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
4013	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
4014	snow_tiedtke=0.
4015	! print*,'avant calcul de la pseudo precip '
4016	! print*,'iflag_cld_th',iflag_cld_th
4017	IF (iflag_cld_th.eq.-1) THEN
4018	rain_tiedtke=rain_con
4019	ELSE
4020	! print*,'calcul de la pseudo precip '
4021	rain_tiedtke=0.
4022	! print*,'calcul de la pseudo precip 0'
4023	DO k=1,klev
4024	DO i=1,klon
4025	IF (d_q_con(i,k).lt.0.) THEN
4026	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
4027	*(paprs(i,k)-paprs(i,k+1))/rg
4028	ENDIF
4029	ENDDO
4030	ENDDO
4031	ENDIF
4032	!
4033	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
4034	!
4035
4036	! Nuages diagnostiques pour Tiedtke
4037	CALL diagcld1(paprs,pplay, &
4038	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
4039	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
4040	diafra,dialiq)
4041	DO k = 1, klev
4042	DO i = 1, klon
4043	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4044	radocond(i,k) = dialiq(i,k)
4045	cldfra(i,k) = diafra(i,k)
4046	ENDIF
4047	ENDDO
4048	ENDDO
4049
4050	ELSE IF (iflag_cld_th.ge.3) THEN
4051	! On prend pour les nuages convectifs le max du calcul de la
4052	! convection et du calcul du pas de temps precedent diminue d'un facteur
4053	! facttemps
4054	facteur = pdtphys *facttemps
4055	DO k=1,klev
4056	DO i=1,klon
4057	rnebcon(i,k)=rnebcon(i,k)*facteur
4058	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
4059	rnebcon(i,k)=rnebcon0(i,k)
4060	clwcon(i,k)=clwcon0(i,k)
4061	ENDIF
4062	ENDDO
4063	ENDDO
4064
4065	! On prend la somme des fractions nuageuses et des contenus en eau
4066
4067	IF (iflag_cld_th>=5) THEN
4068
4069	DO k=1,klev
4070	ptconvth(:,k)=fm_therm(:,k+1)>0.
4071	ENDDO
4072
4073	IF (iflag_coupl==4) THEN
4074
4075	! Dans le cas iflag_coupl==4, on prend la somme des convertures
4076	! convectives et lsc dans la partie des thermiques
4077	! Le controle par iflag_coupl est peut etre provisoire.
4078	DO k=1,klev
4079	DO i=1,klon
4080	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
4081	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4082	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4083	ELSE IF (ptconv(i,k)) THEN
4084	cldfra(i,k)=rnebcon(i,k)
4085	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4086	ENDIF
4087	ENDDO
4088	ENDDO
4089
4090	ELSE IF (iflag_coupl==5) THEN
4091	DO k=1,klev
4092	DO i=1,klon
4093	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4094	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4095	ENDDO
4096	ENDDO
4097
4098	ELSE
4099
4100	! Si on est sur un point touche par la convection
4101	! profonde et pas par les thermiques, on prend la
4102	! couverture nuageuse et l'eau nuageuse de la convection
4103	! profonde.
4104
4105	!IM/FH: 2011/02/23
4106	! definition des points sur lesquels ls thermiques sont actifs
4107
4108	DO k=1,klev
4109	DO i=1,klon
4110	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
4111	cldfra(i,k)=rnebcon(i,k)
4112	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4113	ENDIF
4114	ENDDO
4115	ENDDO
4116
4117	ENDIF
4118
4119	ELSE
4120
4121	! Ancienne version
4122	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
4123	radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
4124	ENDIF
4125
4126	ENDIF
4127
4128	! plulsc(:)=0.
4129	! do k=1,klev,-1
4130	! do i=1,klon
4131	! zzz=prfl(:,k)+psfl(:,k)
4132	! if (.not.ptconvth.zzz.gt.0.)
4133	! enddo prfl, psfl,
4134	! enddo
4135	!
4136	! 2. NUAGES STARTIFORMES
4137	!
4138	IF (ok_stratus) THEN
4139	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
4140	DO k = 1, klev
4141	DO i = 1, klon
4142	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4143	radocond(i,k) = dialiq(i,k)
4144	cldfra(i,k) = diafra(i,k)
4145	ENDIF
4146	ENDDO
4147	ENDDO
4148	ENDIF
4149	!
4150	! Precipitation totale
4151	!
4152	DO i = 1, klon
4153	rain_fall(i) = rain_con(i) + rain_lsc(i)
4154	snow_fall(i) = snow_con(i) + snow_lsc(i)
4155	ENDDO
4156	!
4157	! Calculer l'humidite relative pour diagnostique
4158	!
4159	DO k = 1, klev
4160	DO i = 1, klon
4161	zx_t = t_seri(i,k)
4162	IF (thermcep) THEN
4163	!! if (iflag_ice_thermo.eq.0) then !jyg
4164	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
4165	!! else !jyg
4166	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
4167	!! endif !jyg
4168	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
4169	zx_qs = MIN(0.5,zx_qs)
4170	zcor = 1./(1.-retv*zx_qs)
4171	zx_qs = zx_qs*zcor
4172	ELSE
4173	!! IF (zx_t.LT.t_coup) THEN !jyg
4174	IF (zx_t.LT.rtt) THEN !jyg
4175	zx_qs = qsats(zx_t)/pplay(i,k)
4176	ELSE
4177	zx_qs = qsatl(zx_t)/pplay(i,k)
4178	ENDIF
4179	ENDIF
4180	zx_rh(i,k) = q_seri(i,k)/zx_qs
4181	IF (iflag_ice_thermo .GT. 0) THEN
4182	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
4183	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
4184	ENDIF
4185	zqsat(i,k)=zx_qs
4186	ENDDO
4187	ENDDO
4188
4189	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
4190	! equivalente a 2m (tpote) pour diagnostique
4191	!
4192	DO i = 1, klon
4193	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
4194	IF (thermcep) THEN
4195	IF(zt2m(i).LT.RTT) then
4196	Lheat=RLSTT
4197	ELSE
4198	Lheat=RLVTT
4199	ENDIF
4200	ELSE
4201	IF (zt2m(i).LT.RTT) THEN
4202	Lheat=RLSTT
4203	ELSE
4204	Lheat=RLVTT
4205	ENDIF
4206	ENDIF
4207	tpote(i) = tpot(i)* &
4208	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
4209	ENDDO
4210
4211	IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
4212	IF (CPPKEY_INCA) THEN
4213	CALL VTe(VTphysiq)
4214	CALL VTb(VTinca)
4215	calday = REAL(days_elapsed + 1) + jH_cur
4216
4217	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
4218	CALL AEROSOL_METEO_CALC( &
4219	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
4220	prfl,psfl,pctsrf,cell_area, &
4221	latitude_deg,longitude_deg,u10m,v10m)
4222
4223	zxsnow_dummy(:) = 0.0
4224
4225	CALL chemhook_begin (calday, &
4226	days_elapsed+1, &
4227	jH_cur, &
4228	pctsrf(1,1), &
4229	latitude_deg, &
4230	longitude_deg, &
4231	cell_area, &
4232	paprs, &
4233	pplay, &
4234	coefh(1:klon,1:klev,is_ave), &
4235	pphi, &
4236	t_seri, &
4237	u, &
4238	v, &
4239	rot, &
4240	wo(:, :, 1), &
4241	q_seri, &
4242	zxtsol, &
4243	zt2m, &
4244	zxsnow_dummy, &
4245	solsw, &
4246	albsol1, &
4247	rain_fall, &
4248	snow_fall, &
4249	itop_con, &
4250	ibas_con, &
4251	cldfra, &
4252	nbp_lon, &
4253	nbp_lat-1, &
4254	tr_seri(:,:,1+nqCO2:nbtr), &
4255	ftsol, &
4256	paprs, &
4257	cdragh, &
4258	cdragm, &
4259	pctsrf, &
4260	pdtphys, &
4261	itap)
4262
4263	CALL VTe(VTinca)
4264	CALL VTb(VTphysiq)
4265	END IF
4266	ENDIF !type_trac = inca or inco
4267	IF (type_trac == 'repr') THEN
4268	IF (CPPKEY_REPROBUS) THEN
4269	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
4270	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
4271	END IF
4272	ENDIF
4273
4274	!
4275	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
4276	!
4277	IF (MOD(itaprad,radpas).EQ.0) THEN
4278
4279	!
4280	!jq - introduce the aerosol direct and first indirect radiative forcings
4281	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
4282	IF (flag_aerosol .GT. 0) THEN
4283	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4284	IF (.NOT. aerosol_couple) THEN
4285	!
4286	CALL readaerosol_optic( &
4287	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4288	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4289	mass_solu_aero, mass_solu_aero_pi, &
4290	tau_aero, piz_aero, cg_aero, &
4291	tausum_aero, tau3d_aero)
4292	ENDIF
4293	ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation
4294	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
4295	abort_message='config_inca=aero et rrtm=1 impossible'
4296	CALL abort_physic(modname,abort_message,1)
4297	ELSE
4298	!
4299	#ifdef CPP_RRTM
4300	IF (NSW.EQ.6) THEN
4301	!--new aerosol properties SW and LW
4302	!
4303	IF (CPPKEY_DUST) THEN
4304	!--SPL aerosol model
4305	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
4306	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4307	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4308	tausum_aero, tau3d_aero)
4309	ELSE
4310	!--climatologies or INCA aerosols
4311	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
4312	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4313	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4314	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4315	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4316	tausum_aero, drytausum_aero, tau3d_aero)
4317	END IF
4318
4319	IF (flag_aerosol .EQ. 7) THEN
4320	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
4321	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
4322	ENDIF
4323
4324	!
4325	ELSE IF (NSW.EQ.2) THEN
4326	!--for now we use the old aerosol properties
4327	!
4328	CALL readaerosol_optic( &
4329	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4330	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4331	mass_solu_aero, mass_solu_aero_pi, &
4332	tau_aero, piz_aero, cg_aero, &
4333	tausum_aero, tau3d_aero)
4334	!
4335	!--natural aerosols
4336	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4337	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4338	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4339	!--all aerosols
4340	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4341	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4342	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
4343	!
4344	!--no LW optics
4345	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4346	!
4347	ELSE
4348	abort_message='Only NSW=2 or 6 are possible with ' &
4349	// 'aerosols and iflag_rrtm=1'
4350	CALL abort_physic(modname,abort_message,1)
4351	ENDIF
4352	#else
4353	abort_message='You should compile with -rrtm if running ' &
4354	// 'with iflag_rrtm=1'
4355	CALL abort_physic(modname,abort_message,1)
4356	#endif
4357	!
4358	ENDIF
4359	ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION
4360	#ifdef CPP_ECRAD
4361	!--climatologies or INCA aerosols
4362	CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
4363	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4364	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4365	tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer)
4366	#else
4367	abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
4368	CALL abort_physic(modname,abort_message,1)
4369	#endif
4370	ENDIF
4371
4372	ELSE !--flag_aerosol = 0
4373	tausum_aero(:,:,:) = 0.
4374	drytausum_aero(:,:) = 0.
4375	mass_solu_aero(:,:) = 0.
4376	mass_solu_aero_pi(:,:) = 0.
4377	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4378	tau_aero(:,:,:,:) = 1.e-15
4379	piz_aero(:,:,:,:) = 1.
4380	cg_aero(:,:,:,:) = 0.
4381	ELSE
4382	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4383	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4384	piz_aero_sw_rrtm(:,:,:,:) = 1.0
4385	cg_aero_sw_rrtm(:,:,:,:) = 0.0
4386	ENDIF
4387	ENDIF
4388	!
4389	!--WMO criterion to determine tropopause
4390	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4391	!
4392	!--STRAT AEROSOL
4393	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
4394	IF (flag_aerosol_strat.GT.0) THEN
4395	IF (prt_level .GE.10) THEN
4396	PRINT *,'appel a readaerosolstrat', mth_cur
4397	ENDIF
4398	IF (iflag_rrtm.EQ.0) THEN
4399	IF (flag_aerosol_strat.EQ.1) THEN
4400	CALL readaerosolstrato(debut)
4401	ELSE
4402	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4403	CALL abort_physic(modname,abort_message,1)
4404	ENDIF
4405	ELSE
4406	#ifdef CPP_RRTM
4407	IF (.NOT. CPPKEY_STRATAER) THEN
4408	!--prescribed strat aerosols
4409	!--only in the case of non-interactive strat aerosols
4410	IF (flag_aerosol_strat.EQ.1) THEN
4411	CALL readaerosolstrato1_rrtm(debut)
4412	ELSEIF (flag_aerosol_strat.EQ.2) THEN
4413	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
4414	ELSE
4415	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4416	CALL abort_physic(modname,abort_message,1)
4417	ENDIF
4418	END IF
4419	#else
4420	abort_message='You should compile with -rrtm if running ' &
4421	// 'with iflag_rrtm=1'
4422	CALL abort_physic(modname,abort_message,1)
4423	#endif
4424	ENDIF
4425	ELSE
4426	tausum_aero(:,:,id_STRAT_phy) = 0.
4427	ENDIF
4428	!
4429	#ifdef CPP_RRTM
4430	IF (CPPKEY_STRATAER) THEN
4431	!--compute stratospheric mask
4432	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4433	!--interactive strat aerosols
4434	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4435	END IF
4436	#endif
4437	!--fin STRAT AEROSOL
4438	!
4439
4440	! Calculer les parametres optiques des nuages et quelques
4441	! parametres pour diagnostiques:
4442	!
4443	IF (aerosol_couple.AND.config_inca=='aero') THEN
4444	mass_solu_aero(:,:) = ccm(:,:,1)
4445	mass_solu_aero_pi(:,:) = ccm(:,:,2)
4446	ENDIF
4447
4448	!Rajout appel a interface calcul proprietes optiques des nuages
4449	CALL call_cloud_optics_prop(klon, klev, ok_newmicro, &
4450	paprs, pplay, t_seri, radocond, picefra, cldfra, &
4451	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4452	flwp, fiwp, flwc, fiwc, ok_aie, &
4453	mass_solu_aero, mass_solu_aero_pi, &
4454	cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
4455	ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
4456	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
4457	zfice, dNovrN, ptconv, rnebcon, clwcon)
4458
4459	!
4460	!IM betaCRF
4461	!
4462	cldtaurad = cldtau
4463	cldtaupirad = cldtaupi
4464	cldemirad = cldemi
4465	cldfrarad = cldfra
4466
4467	!
4468	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4469	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4470	!
4471	! global
4472	!
4473	!IM 251017 begin
4474	! print*,'physiq betaCRF global zdtime=',zdtime
4475	!IM 251017 end
4476	DO k=1, klev
4477	DO i=1, klon
4478	IF (pplay(i,k).GE.pfree) THEN
4479	beta(i,k) = beta_pbl
4480	ELSE
4481	beta(i,k) = beta_free
4482	ENDIF
4483	IF (mskocean_beta) THEN
4484	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4485	ENDIF
4486	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4487	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4488	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4489	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4490	ENDDO
4491	ENDDO
4492	!
4493	ELSE
4494	!
4495	! regional
4496	!
4497	DO k=1, klev
4498	DO i=1,klon
4499	!
4500	IF (longitude_deg(i).ge.lon1_beta.AND. &
4501	longitude_deg(i).le.lon2_beta.AND. &
4502	latitude_deg(i).le.lat1_beta.AND. &
4503	latitude_deg(i).ge.lat2_beta) THEN
4504	IF (pplay(i,k).GE.pfree) THEN
4505	beta(i,k) = beta_pbl
4506	ELSE
4507	beta(i,k) = beta_free
4508	ENDIF
4509	IF (mskocean_beta) THEN
4510	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4511	ENDIF
4512	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4513	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4514	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4515	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4516	ENDIF
4517	!
4518	ENDDO
4519	ENDDO
4520	!
4521	ENDIF
4522
4523	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4524	IF (ok_chlorophyll) THEN
4525	print*,"-- reading chlorophyll"
4526	CALL readchlorophyll(debut)
4527	ENDIF
4528
4529	!--if ok_suntime_rrtm we use ancillay data for RSUN
4530	!--previous values are therefore overwritten
4531	!--this is needed for CMIP6 runs
4532	!--and only possible for new radiation scheme
4533	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4534	#ifdef CPP_RRTM
4535	CALL read_rsun_rrtm(debut)
4536	#endif
4537	ENDIF
4538
4539	IF (mydebug) THEN
4540	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4541	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4542	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4543	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4544	ENDIF
4545
4546	!
4547	!sonia : If Iflag_radia >=2, pertubation of some variables
4548	!input to radiation (DICE)
4549	!
4550	IF (iflag_radia .ge. 2) THEN
4551	zsav_tsol (:) = zxtsol(:)
4552	CALL perturb_radlwsw(zxtsol,iflag_radia)
4553	ENDIF
4554
4555	IF (aerosol_couple.AND.config_inca=='aero') THEN
4556	IF (CPPKEY_INCA) THEN
4557	CALL radlwsw_inca &
4558	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4559	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4560	size(wo,3), wo, &
4561	cldfrarad, cldemirad, cldtaurad, &
4562	heat,heat0,cool,cool0,albpla, &
4563	topsw,toplw,solsw,sollw, &
4564	sollwdown, &
4565	topsw0,toplw0,solsw0,sollw0, &
4566	lwdn0, lwdn, lwup0, lwup, &
4567	swdn0, swdn, swup0, swup, &
4568	ok_ade, ok_aie, &
4569	tau_aero, piz_aero, cg_aero, &
4570	topswad_aero, solswad_aero, &
4571	topswad0_aero, solswad0_aero, &
4572	topsw_aero, topsw0_aero, &
4573	solsw_aero, solsw0_aero, &
4574	cldtaupirad, &
4575	topswai_aero, solswai_aero)
4576	END IF
4577	ELSE
4578	!
4579	!IM calcul radiatif pour le cas actuel
4580	!
4581	RCO2 = RCO2_act
4582	RCH4 = RCH4_act
4583	RN2O = RN2O_act
4584	RCFC11 = RCFC11_act
4585	RCFC12 = RCFC12_act
4586	!
4587	!--interactive CO2 in ppm from carbon cycle
4588	IF (carbon_cycle_rad) RCO2=RCO2_glo
4589	!
4590	IF (prt_level .GE.10) THEN
4591	print *,' ->radlwsw, number 1 '
4592	ENDIF
4593	!
4594	! AI namelist utilise pour l appel principal de radlwsw (ecrad)
4595	namelist_ecrad_file='namelist_ecrad'
4596	!
4597	CALL radlwsw &
4598	(debut, dist, rmu0, fract, &
4599	!albedo SB >>>
4600	! paprs, pplay,zxtsol,albsol1, albsol2, &
4601	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4602	!albedo SB <<<
4603	t_seri,q_seri,wo, &
4604	cldfrarad, cldemirad, cldtaurad, &
4605	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4606	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4607	tau_aero, piz_aero, cg_aero, &
4608	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4609	! Rajoute par OB pour RRTM
4610	tau_aero_lw_rrtm, &
4611	cldtaupirad, m_allaer, &
4612	! zqsat, flwcrad, fiwcrad, &
4613	zqsat, flwc, fiwc, &
4614	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4615	namelist_ecrad_file, &
4616	heat,heat0,cool,cool0,albpla, &
4617	heat_volc,cool_volc, &
4618	topsw,toplw,solsw,solswfdiff,sollw, &
4619	sollwdown, &
4620	topsw0,toplw0,solsw0,sollw0, &
4621	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4622	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4623	topswad_aero, solswad_aero, &
4624	topswai_aero, solswai_aero, &
4625	topswad0_aero, solswad0_aero, &
4626	topsw_aero, topsw0_aero, &
4627	solsw_aero, solsw0_aero, &
4628	topswcf_aero, solswcf_aero, &
4629	!-C. Kleinschmitt for LW diagnostics
4630	toplwad_aero, sollwad_aero,&
4631	toplwai_aero, sollwai_aero, &
4632	toplwad0_aero, sollwad0_aero,&
4633	!-end
4634	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4635	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4636	cloud_cover_sw)
4637
4638	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4639	!schemes
4640	toplw = toplw + betalwoff * (toplw0 - toplw)
4641	sollw = sollw + betalwoff * (sollw0 - sollw)
4642	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4643	lwup = lwup + betalwoff * (lwup0 - lwup)
4644	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4645	sollwdown(:))
4646	cool = cool + betalwoff * (cool0 - cool)
4647
4648	IF (.NOT. using_xios) THEN
4649	!
4650	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4651	!IM des taux doit etre different du taux actuel
4652	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4653	!
4654	IF (RCO2_per.NE.RCO2_act.OR. &
4655	RCH4_per.NE.RCH4_act.OR. &
4656	RN2O_per.NE.RN2O_act.OR. &
4657	RCFC11_per.NE.RCFC11_act.OR. &
4658	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4659	ENDIF
4660	!
4661	IF (ok_4xCO2atm) THEN
4662	!
4663	RCO2 = RCO2_per
4664	RCH4 = RCH4_per
4665	RN2O = RN2O_per
4666	RCFC11 = RCFC11_per
4667	RCFC12 = RCFC12_per
4668	!
4669	IF (prt_level .GE.10) THEN
4670	print *,' ->radlwsw, number 2 '
4671	ENDIF
4672	!
4673	namelist_ecrad_file='namelist_ecrad'
4674	!
4675	CALL radlwsw &
4676	(debut, dist, rmu0, fract, &
4677	!albedo SB >>>
4678	! paprs, pplay,zxtsol,albsol1, albsol2, &
4679	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4680	!albedo SB <<<
4681	t_seri,q_seri,wo, &
4682	cldfrarad, cldemirad, cldtaurad, &
4683	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4684	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4685	tau_aero, piz_aero, cg_aero, &
4686	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4687	! Rajoute par OB pour RRTM
4688	tau_aero_lw_rrtm, &
4689	cldtaupi, m_allaer, &
4690	! zqsat, flwcrad, fiwcrad, &
4691	zqsat, flwc, fiwc, &
4692	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4693	namelist_ecrad_file, &
4694	heatp,heat0p,coolp,cool0p,albplap, &
4695	heat_volc,cool_volc, &
4696	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4697	sollwdownp, &
4698	topsw0p,toplw0p,solsw0p,sollw0p, &
4699	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4700	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4701	topswad_aerop, solswad_aerop, &
4702	topswai_aerop, solswai_aerop, &
4703	topswad0_aerop, solswad0_aerop, &
4704	topsw_aerop, topsw0_aerop, &
4705	solsw_aerop, solsw0_aerop, &
4706	topswcf_aerop, solswcf_aerop, &
4707	!-C. Kleinschmitt for LW diagnostics
4708	toplwad_aerop, sollwad_aerop,&
4709	toplwai_aerop, sollwai_aerop, &
4710	toplwad0_aerop, sollwad0_aerop,&
4711	!-end
4712	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4713	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4714	cloud_cover_sw)
4715	ENDIF !ok_4xCO2atm
4716
4717	! A.I aout 2023
4718	! Effet 3D des nuages Ecrad
4719	! a passer : nom du ficher namelist et cles ok_3Deffect
4720	! a declarer comme iflag_rrtm et a lire dans physiq.def
4721	#ifdef CPP_ECRAD
4722	IF (ok_3Deffect) then
4723	! print*,'ok_3Deffect = ',ok_3Deffect
4724	namelist_ecrad_file='namelist_ecrad_s2'
4725	CALL radlwsw &
4726	(debut, dist, rmu0, fract, &
4727	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4728	t_seri,q_seri,wo, &
4729	cldfrarad, cldemirad, cldtaurad, &
4730	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4731	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4732	tau_aero, piz_aero, cg_aero, &
4733	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4734	tau_aero_lw_rrtm, &
4735	cldtaupi, m_allaer, &
4736	zqsat, flwc, fiwc, &
4737	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4738	namelist_ecrad_file, &
4739	! A modifier
4740	heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
4741	heat_volc,cool_volc, &
4742	topsw_s2,toplw_s2,solsw_s2,solswfdiff_s2,sollw_s2, &
4743	sollwdown_s2, &
4744	topsw0_s2,toplw0_s2,solsw0_s2,sollw0_s2, &
4745	lwdnc0_s2, lwdn0_s2, lwdn_s2, lwupc0_s2, lwup0_s2, lwup_s2, &
4746	swdnc0_s2, swdn0_s2, swdn_s2, swupc0_s2, swup0_s2, swup_s2, &
4747	topswad_aero_s2, solswad_aero_s2, &
4748	topswai_aero_s2, solswai_aero_s2, &
4749	topswad0_aero_s2, solswad0_aero_s2, &
4750	topsw_aero_s2, topsw0_aero_s2, &
4751	solsw_aero_s2, solsw0_aero_s2, &
4752	topswcf_aero_s2, solswcf_aero_s2, &
4753	!-C. Kleinschmitt for LW diagnostics
4754	toplwad_aero_s2, sollwad_aero_s2,&
4755	toplwai_aero_s2, sollwai_aero_s2, &
4756	toplwad0_aero_s2, sollwad0_aero_s2,&
4757	!-end
4758	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4759	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4760	cloud_cover_sw_s2)
4761	ENDIF ! ok_3Deffect
4762	#endif
4763
4764	ENDIF ! aerosol_couple
4765	itaprad = 0
4766	!
4767	! If Iflag_radia >=2, reset pertubed variables
4768	!
4769	IF (iflag_radia .ge. 2) THEN
4770	zxtsol(:) = zsav_tsol (:)
4771	ENDIF
4772	ENDIF ! MOD(itaprad,radpas)
4773	itaprad = itaprad + 1
4774
4775	IF (iflag_radia.eq.0) THEN
4776	IF (prt_level.ge.9) THEN
4777	PRINT *,'--------------------------------------------------'
4778	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4779	PRINT *,'>>>> heat et cool mis a zero '
4780	PRINT *,'--------------------------------------------------'
4781	ENDIF
4782	heat=0.
4783	cool=0.
4784	sollw=0. ! MPL 01032011
4785	solsw=0.
4786	radsol=0.
4787	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4788	swup0=0.
4789	lwup=0.
4790	lwup0=0.
4791	lwdn=0.
4792	lwdn0=0.
4793	ENDIF
4794
4795	!
4796	! Calculer radsol a l'exterieur de radlwsw
4797	! pour prendre en compte le cycle diurne
4798	! recode par Olivier Boucher en sept 2015
4799	!
4800	radsol=solsw*swradcorr+sollw
4801
4802	IF (ok_4xCO2atm) THEN
4803	radsolp=solswp*swradcorr+sollwp
4804	ENDIF
4805
4806	!
4807	! Ajouter la tendance des rayonnements (tous les pas)
4808	! avec une correction pour le cycle diurne dans le SW
4809	!
4810
4811	DO k=1, klev
4812	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4813	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4814	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4815	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4816	ENDDO
4817
4818	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4819	CALL prt_enerbil('SW',itap)
4820	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4821	CALL prt_enerbil('LW',itap)
4822
4823	!
4824	IF (mydebug) THEN
4825	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4826	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4827	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4828	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4829	ENDIF
4830
4831	! Calculer l'hydrologie de la surface
4832	!
4833	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4834	! . agesno, ftsol,fqsurf,fsnow, ruis)
4835	!
4836
4837	!
4838	! Calculer le bilan du sol et la derive de temperature (couplage)
4839	!
4840	DO i = 1, klon
4841	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4842	! a la demande de JLD
4843	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4844	ENDDO
4845	!
4846	!moddeblott(jan95)
4847	! Appeler le programme de parametrisation de l'orographie
4848	! a l'echelle sous-maille:
4849	!
4850	IF (prt_level .GE.10) THEN
4851	print *,' call orography ? ', ok_orodr
4852	ENDIF
4853	!
4854	IF (ok_orodr) THEN
4855	!
4856	! selection des points pour lesquels le shema est actif:
4857	igwd=0
4858	DO i=1,klon
4859	itest(i)=0
4860	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4861	!zrel_oro: relative mountain height wrt relief explained by mean slope
4862	! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
4863	! such as ice sheets (work by V. Wiener)
4864	! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
4865	! earn computation time but they are not physical.
4866	IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
4867	itest(i)=1
4868	igwd=igwd+1
4869	idx(igwd)=i
4870	ENDIF
4871	ENDDO
4872	! igwdim=MAX(1,igwd)
4873	!
4874	IF (ok_strato) THEN
4875
4876	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4877	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4878	igwd,idx,itest, &
4879	t_seri, u_seri, v_seri, &
4880	zulow, zvlow, zustrdr, zvstrdr, &
4881	d_t_oro, d_u_oro, d_v_oro)
4882
4883	ELSE
4884	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4885	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4886	igwd,idx,itest, &
4887	t_seri, u_seri, v_seri, &
4888	zulow, zvlow, zustrdr, zvstrdr, &
4889	d_t_oro, d_u_oro, d_v_oro)
4890	ENDIF
4891	!
4892	! ajout des tendances
4893	!-----------------------------------------------------------------------
4894	! ajout des tendances de la trainee de l'orographie
4895	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
4896	abortphy,flag_inhib_tend,itap,0)
4897	CALL prt_enerbil('oro',itap)
4898	!----------------------------------------------------------------------
4899	!
4900	ENDIF ! fin de test sur ok_orodr
4901	!
4902	IF (mydebug) THEN
4903	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4904	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4905	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4906	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4907	ENDIF
4908
4909	IF (ok_orolf) THEN
4910	!
4911	! selection des points pour lesquels le shema est actif:
4912	igwd=0
4913	DO i=1,klon
4914	itest(i)=0
4915	!zrel_oro: relative mountain height wrt relief explained by mean slope
4916	! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
4917	! such as ice sheets (work by V. Wiener)
4918	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4919	IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
4920	itest(i)=1
4921	igwd=igwd+1
4922	idx(igwd)=i
4923	ENDIF
4924	ENDDO
4925	! igwdim=MAX(1,igwd)
4926	!
4927	IF (ok_strato) THEN
4928
4929	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4930	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4931	igwd,idx,itest, &
4932	t_seri, u_seri, v_seri, &
4933	zulow, zvlow, zustrli, zvstrli, &
4934	d_t_lif, d_u_lif, d_v_lif )
4935
4936	ELSE
4937	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4938	latitude_deg,zmea,zstd,zpic, &
4939	itest, &
4940	t_seri, u_seri, v_seri, &
4941	zulow, zvlow, zustrli, zvstrli, &
4942	d_t_lif, d_u_lif, d_v_lif)
4943	ENDIF
4944
4945	! ajout des tendances de la portance de l'orographie
4946	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
4947	'lif', abortphy,flag_inhib_tend,itap,0)
4948	CALL prt_enerbil('lif',itap)
4949	ENDIF ! fin de test sur ok_orolf
4950
4951	IF (ok_hines) then
4952	! HINES GWD PARAMETRIZATION
4953	east_gwstress=0.
4954	west_gwstress=0.
4955	du_gwd_hines=0.
4956	dv_gwd_hines=0.
4957	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4958	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4959	du_gwd_hines, dv_gwd_hines)
4960	zustr_gwd_hines=0.
4961	zvstr_gwd_hines=0.
4962	DO k = 1, klev
4963	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4964	* (paprs(:, k)-paprs(:, k+1))/rg
4965	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4966	* (paprs(:, k)-paprs(:, k+1))/rg
4967	ENDDO
4968
4969	d_t_hin(:, :)=0.
4970	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4971	dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4972	CALL prt_enerbil('hin',itap)
4973	ENDIF
4974
4975	IF (.not. ok_hines .and. ok_gwd_rando) then
4976	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4977	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4978	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4979	dv_gwd_front, east_gwstress, west_gwstress)
4980	zustr_gwd_front=0.
4981	zvstr_gwd_front=0.
4982	DO k = 1, klev
4983	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4984	* (paprs(:, k)-paprs(:, k+1))/rg
4985	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4986	* (paprs(:, k)-paprs(:, k+1))/rg
4987	ENDDO
4988
4989	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
4990	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4991	CALL prt_enerbil('front_gwd_rando',itap)
4992	ENDIF
4993
4994	IF (ok_gwd_rando) THEN
4995	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4996	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4997	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4998	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
4999	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
5000	CALL prt_enerbil('flott_gwd_rando',itap)
5001	zustr_gwd_rando=0.
5002	zvstr_gwd_rando=0.
5003	DO k = 1, klev
5004	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
5005	* (paprs(:, k)-paprs(:, k+1))/rg
5006	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
5007	* (paprs(:, k)-paprs(:, k+1))/rg
5008	ENDDO
5009	ENDIF
5010
5011	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
5012
5013	IF (mydebug) THEN
5014	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5015	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5016	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5017	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5018	ENDIF
5019
5020	DO i = 1, klon
5021	zustrph(i)=0.
5022	zvstrph(i)=0.
5023	ENDDO
5024	DO k = 1, klev
5025	DO i = 1, klon
5026	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
5027	(paprs(i,k)-paprs(i,k+1))/rg
5028	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
5029	(paprs(i,k)-paprs(i,k+1))/rg
5030	ENDDO
5031	ENDDO
5032	!
5033	!IM calcul composantes axiales du moment angulaire et couple des montagnes
5034	!
5035	IF (is_sequential .and. ok_orodr) THEN
5036	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
5037	ra,rg,romega, &
5038	latitude_deg,longitude_deg,pphis, &
5039	zustrdr,zustrli,zustrph, &
5040	zvstrdr,zvstrli,zvstrph, &
5041	paprs,u,v, &
5042	aam, torsfc)
5043	ENDIF
5044	!IM cf. FLott END
5045	!DC Calcul de la tendance due au methane
5046	IF (ok_qch4) THEN
5047	! d_q_ch4: H2O source from CH4 in MMR/s (mass mixing ratio/s or kg H2O/kg air/s)
5048	IF (CPPKEY_STRATAER) THEN
5049
5050	CALL stratH2O_methox(debut,paprs,d_q_ch4)
5051	ELSE
5052	! ECMWF routine METHOX
5053	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
5054	END IF
5055	! add humidity tendency due to methane
5056	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
5057	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
5058	'q_ch4', abortphy,flag_inhib_tend,itap,0)
5059	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep ! update with H2O conserv done in add_phys_tend
5060	ENDIF
5061	!
5062	!
5063	IF (CPPKEY_STRATAER) THEN
5064	IF (ok_qemiss) THEN
5065	flh2o=1
5066	IF(flag_verbose_strataer) THEN
5067	print *,'IN physiq_mod: ok_qemiss =yes (',ok_qemiss,'), flh2o=',flh2o
5068	print *,'IN physiq_mod: flag_emit=',flag_emit,', nErupt=',nErupt
5069	print *,'IN physiq_mod: nAerErupt=',nAerErupt
5070	ENDIF
5071
5072	SELECT CASE(flag_emit)
5073	CASE(1) ! emission volc H2O in LMDZ
5074	DO ieru=1, nErupt
5075	IF (year_cur==year_emit_vol(ieru).AND.&
5076	mth_cur==mth_emit_vol(ieru).AND.&
5077	day_cur>=day_emit_vol(ieru).AND.&
5078	day_cur<(day_emit_vol(ieru)+injdur)) THEN
5079
5080	IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
5081	! initialisation of q tendency emission
5082	d_q_emiss(:,:)=0.
5083	! daily injection mass emission - NL
5084	m_H2O_emiss_vol_daily = m_H2O_emiss_vol(ieru)/(REAL(injdur)&
5085	*REAL(ponde_lonlat_vol(ieru)))
5086	!
5087	CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
5088	pplay,paprs,tr_seri,&
5089	m_H2O_emiss_vol_daily,&
5090	xlat_min_vol(ieru),xlat_max_vol(ieru),&
5091	xlon_min_vol(ieru),xlon_max_vol(ieru),&
5092	altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
5093	nAerErupt+1,0)
5094
5095	IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
5096	minval(d_q_emiss),maxval(d_q_emiss)
5097
5098	CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
5099	'q_emiss',abortphy,flag_inhib_tend,itap,0)
5100	IF (abortphy==1) Print*,'ERROR ABORT TEND EMISS'
5101	ENDIF
5102	ENDDO
5103	flh2o=0
5104	END SELECT ! emission scenario (flag_emit)
5105	ENDIF
5106	END IF
5107
5108	!===============================================================
5109	! Additional tendency of TKE due to orography
5110	!===============================================================
5111	!
5112	! Inititialization
5113	!------------------
5114
5115	addtkeoro=0
5116	CALL getin_p('addtkeoro',addtkeoro)
5117
5118	IF (prt_level.ge.5) &
5119	print*,'addtkeoro', addtkeoro
5120
5121	alphatkeoro=1.
5122	CALL getin_p('alphatkeoro',alphatkeoro)
5123	alphatkeoro=min(max(0.,alphatkeoro),1.)
5124
5125	smallscales_tkeoro=.FALSE.
5126	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
5127
5128
5129	dtadd(:,:)=0.
5130	duadd(:,:)=0.
5131	dvadd(:,:)=0.
5132
5133	! Choices for addtkeoro:
5134	! ** 0 no TKE tendency from orography
5135	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
5136	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
5137	!
5138
5139	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
5140	! -------------------------------------------
5141
5142
5143	! selection des points pour lesquels le schema est actif:
5144
5145
5146	IF (addtkeoro .EQ. 1 ) THEN
5147
5148	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
5149	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
5150
5151	ELSE IF (addtkeoro .EQ. 2) THEN
5152
5153	IF (smallscales_tkeoro) THEN
5154	igwd=0
5155	DO i=1,klon
5156	itest(i)=0
5157	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
5158	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
5159	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
5160	IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
5161	itest(i)=1
5162	igwd=igwd+1
5163	idx(igwd)=i
5164	ENDIF
5165	ENDDO
5166
5167	ELSE
5168
5169	igwd=0
5170	DO i=1,klon
5171	itest(i)=0
5172	IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
5173	itest(i)=1
5174	igwd=igwd+1
5175	idx(igwd)=i
5176	ENDIF
5177	ENDDO
5178
5179	ENDIF
5180
5181	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
5182	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
5183	igwd,idx,itest, &
5184	t_seri, u_seri, v_seri, &
5185	zulow, zvlow, zustrdr, zvstrdr, &
5186	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
5187
5188	zustrdr(:)=0.
5189	zvstrdr(:)=0.
5190	zulow(:)=0.
5191	zvlow(:)=0.
5192
5193	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
5194	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
5195	ENDIF
5196
5197
5198	! TKE update from subgrid temperature and wind tendencies
5199	!----------------------------------------------------------
5200	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5201
5202
5203	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
5204	!
5205	! Prevent pbl_tke_w from becoming negative
5206	wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
5207	!
5208
5209	ENDIF
5210	! -----
5211	!===============================================================
5212
5213
5214	!====================================================================
5215	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
5216	!====================================================================
5217	! Abderrahmane 24.08.09
5218
5219	IF (ok_cosp) THEN
5220	! adeclarer
5221	IF (CPPKEY_COSP) THEN
5222	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5223
5224	IF (prt_level .GE.10) THEN
5225	print*,'freq_cosp',freq_cosp
5226	ENDIF
5227	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5228	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
5229	! s ref_liq,ref_ice
5230	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
5231	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5232	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5233	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5234	JrNt,ref_liq,ref_ice, &
5235	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5236	zu10m,zv10m,pphis, &
5237	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5238	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5239	prfl(:,1:klev),psfl(:,1:klev), &
5240	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5241	mr_ozone,cldtau, cldemi)
5242
5243	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
5244	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
5245	! M clMISR,
5246	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
5247	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
5248
5249	ENDIF
5250	END IF
5251
5252	IF (CPPKEY_COSP2) THEN
5253	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5254
5255	IF (prt_level .GE.10) THEN
5256	print*,'freq_cosp',freq_cosp
5257	ENDIF
5258	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5259	print*,'Dans physiq.F avant appel '
5260	! s ref_liq,ref_ice
5261	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
5262	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5263	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5264	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5265	JrNt,ref_liq,ref_ice, &
5266	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5267	zu10m,zv10m,pphis, &
5268	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5269	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5270	prfl(:,1:klev),psfl(:,1:klev), &
5271	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5272	mr_ozone,cldtau, cldemi)
5273	ENDIF
5274	END IF
5275
5276	IF (CPPKEY_COSPV2) THEN
5277	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5278	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5279
5280	IF (prt_level .GE.10) THEN
5281	print*,'freq_cosp',freq_cosp
5282	ENDIF
5283	DO k = 1, klev
5284	DO i = 1, klon
5285	phicosp(i,k) = pphi(i,k) + pphis(i)
5286	ENDDO
5287	ENDDO
5288	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5289	print*,'Dans physiq.F avant appel '
5290	! s ref_liq,ref_ice
5291	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
5292	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5293	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5294	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5295	JrNt,ref_liq,ref_ice, &
5296	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5297	zu10m,zv10m,pphis, &
5298	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5299	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5300	prfl(:,1:klev),psfl(:,1:klev), &
5301	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5302	mr_ozone,cldtau, cldemi)
5303	ENDIF
5304	END IF
5305
5306	ENDIF !ok_cosp
5307
5308
5309	! Marine
5310
5311	IF (ok_airs) then
5312
5313	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
5314	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
5315	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
5316	& map_prop_hc,map_prop_hist,&
5317	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
5318	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
5319	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
5320	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
5321	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
5322	& map_ntot,map_hc,map_hist,&
5323	& map_Cb,map_ThCi,map_Anv,&
5324	& alt_tropo )
5325	ENDIF
5326
5327	ENDIF ! ok_airs
5328
5329
5330	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5331	!AA
5332	!AA Installation de l'interface online-offline pour traceurs
5333	!AA
5334	!====================================================================
5335	! Calcul des tendances traceurs
5336	!====================================================================
5337	!
5338
5339	IF (type_trac == 'repr') THEN
5340	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
5341	!MM dans Reprobus
5342	sh_in(:,:) = q_seri(:,:)
5343	IF (CPPKEY_REPROBUS) THEN
5344	d_q_rep(:,:) = 0.
5345	d_ql_rep(:,:) = 0.
5346	d_qi_rep(:,:) = 0.
5347	END IF
5348	ELSE
5349	sh_in(:,:) = qx(:,:,ivap)
5350	IF (nqo >= 3) THEN
5351	ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
5352	ELSE
5353	ch_in(:,:) = qx(:,:,iliq)
5354	ENDIF
5355	ENDIF
5356
5357	! Merge wdtrainA and wdtrainS in the total source of precipitation due to
5358	! adiabatic updraughts.
5359	!
5360	wdtrainAS(:,:) = wdtrainA(:,:) + wdtrainS(:,:)
5361
5362	IF (CPPKEY_DUST) THEN
5363	! Avec SPLA, iflag_phytrac est forcé =1
5364	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
5365	pdtphys,ftsol, & ! I
5366	t,q_seri,paprs,pplay,RHcl, & ! I
5367	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
5368	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
5369	u_seri, v_seri, latitude_deg, longitude_deg, &
5370	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
5371	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
5372	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
5373	ev,wdtrainAS, wdtrainM,wght_cvfd, & ! I
5374	fm_therm, entr_therm, rneb, & ! I
5375	beta_prec_fisrt,beta_prec, & !I
5376	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
5377	d_tr_dyn,tr_seri)
5378
5379	ELSE
5380	IF (iflag_phytrac == 1 ) THEN
5381	CALL phytrac ( &
5382	itap, days_elapsed+1, jH_cur, debut, &
5383	lafin, phys_tstep, u, v, t, &
5384	paprs, pplay, pmfu, pmfd, &
5385	pen_u, pde_u, pen_d, pde_d, &
5386	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
5387	u1, v1, ftsol, pctsrf, &
5388	zustar, zu10m, zv10m, &
5389	wstar(:,is_ave), ale_bl, ale_wake, &
5390	latitude_deg, longitude_deg, &
5391	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
5392	presnivs, pphis, pphi, albsol1, &
5393	sh_in, ch_in, rhcl, cldfra, rneb, &
5394	diafra, radocond, itop_con, ibas_con, &
5395	pmflxr, pmflxs, prfl, psfl, &
5396	da, phi, mp, upwd, &
5397	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
5398	wdtrainAS, wdtrainM, sigd, clw,elij, & !<<RomP
5399	ev, ep, epmlmMm, eplaMm, & !<<RomP
5400	dnwd, aerosol_couple, flxmass_w, &
5401	tau_aero, piz_aero, cg_aero, ccm, &
5402	rfname, &
5403	d_tr_dyn, & !<<RomP
5404	tr_seri, init_source)
5405	IF (CPPKEY_REPROBUS) THEN
5406
5407
5408	print*,'avt add phys rep',abortphy
5409
5410	CALL add_phys_tend &
5411	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
5412	'rep',abortphy,flag_inhib_tend,itap,0)
5413	IF (abortphy==1) Print*,'ERROR ABORT REP'
5414
5415	print*,'apr add phys rep',abortphy
5416
5417	END IF
5418	ENDIF ! (iflag_phytrac=1)
5419
5420	END IF
5421	!ENDIF ! (iflag_phytrac=1)
5422
5423	IF (offline) THEN
5424
5425	IF (prt_level.ge.9) &
5426	print*,'Attention on met a 0 les thermiques pour phystoke'
5427	CALL phystokenc ( &
5428	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5429	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5430	fm_therm,entr_therm, &
5431	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5432	frac_impa, frac_nucl, &
5433	pphis,cell_area,phys_tstep,itap, &
5434	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
5435
5436
5437	ENDIF
5438
5439	!
5440	! Calculer le transport de l'eau et de l'energie (diagnostique)
5441	!
5442	CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5443	ue, ve, uq, vq, uwat, vwat)
5444	!
5445	!IM global posePB BEG
5446	IF(1.EQ.0) THEN
5447	!
5448	CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
5449	ve_lay, vq_lay, ue_lay, uq_lay)
5450	!
5451	ENDIF !(1.EQ.0) THEN
5452	!IM global posePB END
5453	!
5454	! Accumuler les variables a stocker dans les fichiers histoire:
5455	!
5456
5457	!================================================================
5458	! Conversion of kinetic and potential energy into heat, for
5459	! parameterisation of subgrid-scale motions
5460	!================================================================
5461
5462	d_t_ec(:,:)=0.
5463	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5464	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
5465	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
5466	zmasse,exner,d_t_ec)
5467	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
5468
5469	!==================================================================
5470	!--OB water mass fixer for the physics
5471	!--water profiles are corrected to force mass conservation of water
5472	!--currently flag is turned off
5473	!==================================================================
5474	IF (ok_water_mass_fixer) THEN
5475	qql2(:)=0.0
5476	DO k = 1, klev
5477	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
5478	IF (nqo >= 3) THEN
5479	qql2(:)=qql2(:)+qs_seri(:,k)*zmasse(:,k)
5480	ENDIF
5481	IF (ok_bs) THEN
5482	qql2(:)=qql2(:)+qbs_seri(:,k)*zmasse(:,k)
5483	ENDIF
5484	ENDDO
5485
5486	IF (CPPKEY_STRATAER) THEN
5487	IF (ok_qemiss) THEN
5488	DO k = 1, klev
5489	qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
5490	ENDDO
5491	ENDIF
5492	END IF
5493	IF (ok_qch4) THEN
5494	DO k = 1, klev
5495	qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
5496	ENDDO
5497	ENDIF
5498
5499	DO i = 1, klon
5500	!--compute ratio of what q+ql should be with conservation to what it is
5501	IF (ok_bs) THEN
5502	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i)-bs_fall(i))*pdtphys)/qql2(i)
5503	ELSE
5504	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5505	ENDIF
5506	DO k = 1, klev
5507	q_seri(i,k) =q_seri(i,k)*corrqql
5508	ql_seri(i,k)=ql_seri(i,k)*corrqql
5509	IF (nqo >= 3) THEN
5510	qs_seri(i,k)=qs_seri(i,k)*corrqql
5511	ENDIF
5512	IF (ok_bs) THEN
5513	qbs_seri(i,k)=qbs_seri(i,k)*corrqql
5514	ENDIF
5515	ENDDO
5516	ENDDO
5517	ENDIF
5518	!--fin mass fixer
5519
5520	!cc prw = eau precipitable
5521	! prlw = colonne eau liquide
5522	! prlw = colonne eau solide
5523	! prbsw = colonne neige soufflee
5524	! water_budget = non-conservation residual from the LMDZ physics
5525	! (should be equal to machine precision if mass fixer is activated)
5526	prw(:) = 0.
5527	prlw(:) = 0.
5528	prsw(:) = 0.
5529	prbsw(:) = 0.
5530	water_budget(:) = 0.0
5531	DO k = 1, klev
5532	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
5533	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5534	water_budget(:) = water_budget(:) + (q_seri(:,k)-qx(:,k,ivap)+ql_seri(:,k)-qx(:,k,iliq))*zmasse(:,k)
5535	IF (nqo >= 3) THEN
5536	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
5537	water_budget(:) = water_budget(:) + (qs_seri(:,k)-qx(:,k,isol))*zmasse(:,k)
5538	ENDIF
5539	IF (nqo >= 4 .AND. ok_bs) THEN
5540	prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
5541	water_budget(:) = water_budget(:) + (qbs_seri(:,k)-qx(:,k,ibs))*zmasse(:,k)
5542	ENDIF
5543	ENDDO
5544	water_budget(:)=water_budget(:)+(rain_fall(:)+snow_fall(:)-evap(:))*pdtphys
5545	IF (ok_bs) THEN
5546	water_budget(:)=water_budget(:)+bs_fall(:)*pdtphys
5547	ENDIF
5548
5549	!=======================================================================
5550	! SORTIES
5551	!=======================================================================
5552	!
5553	!IM initialisation + calculs divers diag AMIP2
5554	CALL calcul_divers(itap, itapm1, un_jour)
5555	!
5556	!IM Interpolation sur les niveaux de pression du NMC
5557	! -------------------------------------------------
5558	!
5559	include "calcul_STDlev.h"
5560	!
5561	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5562	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5563	!
5564	!
5565	IF (ANY(type_trac == ['inca','inco'])) THEN
5566	IF (CPPKEY_INCA) THEN
5567	CALL VTe(VTphysiq)
5568	CALL VTb(VTinca)
5569
5570	CALL chemhook_end ( &
5571	phys_tstep, &
5572	pplay, &
5573	t_seri, &
5574	tr_seri(:,:,1+nqCO2:nbtr), &
5575	nbtr, &
5576	paprs, &
5577	q_seri, &
5578	cell_area, &
5579	pphi, &
5580	pphis, &
5581	zx_rh, &
5582	aps, bps, ap, bp, lafin)
5583
5584	CALL VTe(VTinca)
5585	CALL VTb(VTphysiq)
5586	END IF
5587	ENDIF
5588
5589	IF (type_trac == 'repr') THEN
5590	IF (CPPKEY_REPROBUS) THEN
5591	CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
5592	END IF
5593	ENDIF
5594
5595	!
5596	! Convertir les incrementations en tendances
5597	!
5598	IF (prt_level .GE.10) THEN
5599	print *,'Convertir les incrementations en tendances '
5600	ENDIF
5601	!
5602	IF (mydebug) THEN
5603	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5604	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5605	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5606	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5607	ENDIF
5608
5609	DO k = 1, klev
5610	DO i = 1, klon
5611	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5612	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5613	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5614	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5615	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
5616	!CR: on ajoute le contenu en glace
5617	IF (nqo >= 3) THEN
5618	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5619	ENDIF
5620	!--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
5621	IF (nqo.ge.5 .and. ok_ice_supersat) THEN
5622	d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
5623	d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
5624	IF (nqo.ge.6 .and. ok_plane_contrail) THEN
5625	d_qx(i,k,ircont) = ( rcont_seri(i,k) - qx(i,k,ircont) ) / phys_tstep
5626	ENDIF
5627	ENDIF
5628
5629	IF (nqo.ge.4 .and. ok_bs) THEN
5630	d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
5631	ENDIF
5632
5633	ENDDO
5634	ENDDO
5635	!
5636	! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
5637	itr = 0
5638	DO iq = 1, nqtot
5639	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
5640	itr = itr+1
5641	DO k = 1, klev
5642	DO i = 1, klon
5643	d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
5644	ENDDO
5645	ENDDO
5646	ENDDO
5647	!
5648	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5649	!IM global posePB include "write_bilKP_ins.h"
5650	!IM global posePB include "write_bilKP_ave.h"
5651	!
5652	! Sauvegarder les valeurs de t et q a la fin de la physique:
5653	!
5654	u_ancien(:,:) = u_seri(:,:)
5655	v_ancien(:,:) = v_seri(:,:)
5656	t_ancien(:,:) = t_seri(:,:)
5657	q_ancien(:,:) = q_seri(:,:)
5658	ql_ancien(:,:) = ql_seri(:,:)
5659	qs_ancien(:,:) = qs_seri(:,:)
5660	qbs_ancien(:,:)= qbs_seri(:,:)
5661	cf_ancien(:,:) = cf_seri(:,:)
5662	rvc_ancien(:,:)= rvc_seri(:,:)
5663	rcont_ancien(:,:)= rcont_seri(:,:)
5664	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5665	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5666	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5667	CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
5668	! !! RomP >>>
5669	IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
5670	! !! RomP <<<
5671	!==========================================================================
5672	! Sorties des tendances pour un point particulier
5673	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5674	! pour le debug
5675	! La valeur de igout est attribuee plus haut dans le programme
5676	!==========================================================================
5677
5678	IF (prt_level.ge.1) THEN
5679	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5680	write(lunout,*) &
5681	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5682	write(lunout,*) &
5683	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5684	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5685	pctsrf(igout,is_sic)
5686	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5687	DO k=1,klev
5688	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5689	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5690	d_t_eva(igout,k)
5691	ENDDO
5692	write(lunout,*) 'cool,heat'
5693	DO k=1,klev
5694	write(lunout,*) cool(igout,k),heat(igout,k)
5695	ENDDO
5696
5697	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5698	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5699	!jyg! do k=1,klev
5700	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5701	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5702	!jyg! enddo
5703	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5704	DO k=1,klev
5705	write(lunout,*) d_t_vdf(igout,k), &
5706	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5707	ENDDO
5708	!>jyg
5709
5710	write(lunout,*) 'd_ps ',d_ps(igout)
5711	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5712	DO k=1,klev
5713	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5714	d_qx(igout,k,1),d_qx(igout,k,2)
5715	ENDDO
5716	ENDIF
5717
5718	!============================================================
5719	! Calcul de la temperature potentielle
5720	!============================================================
5721	DO k = 1, klev
5722	DO i = 1, klon
5723	!JYG/IM theta en debut du pas de temps
5724	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5725	!JYG/IM theta en fin de pas de temps de physique
5726	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5727	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5728	! MPL 20130625
5729	! fth_fonctions.F90 et parkind1.F90
5730	! sinon thetal=theta
5731	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5732	! : ql_seri(i,k))
5733	thetal(i,k)=theta(i,k)
5734	ENDDO
5735	ENDDO
5736	!
5737
5738	! 22.03.04 BEG
5739	!=============================================================
5740	! Ecriture des sorties
5741	!=============================================================
5742
5743	! Recupere des varibles calcule dans differents modules
5744	! pour ecriture dans histxxx.nc
5745
5746	! Get some variables from module fonte_neige_mod
5747	CALL fonte_neige_get_vars(pctsrf, &
5748	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5749
5750
5751	!=============================================================
5752	! Separation entre thermiques et non thermiques dans les sorties
5753	! de fisrtilp
5754	!=============================================================
5755
5756	IF (iflag_thermals>=1) THEN
5757	d_t_lscth=0.
5758	d_t_lscst=0.
5759	d_q_lscth=0.
5760	d_q_lscst=0.
5761	DO k=1,klev
5762	DO i=1,klon
5763	IF (ptconvth(i,k)) THEN
5764	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5765	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5766	ELSE
5767	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5768	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5769	ENDIF
5770	ENDDO
5771	ENDDO
5772
5773	DO i=1,klon
5774	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5775	plul_th(i)=prfl(i,1)+psfl(i,1)
5776	ENDDO
5777	ENDIF
5778
5779	!On effectue les sorties:
5780
5781	IF (CPPKEY_DUST) THEN
5782	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5783	pplay, lmax_th, aerosol_couple, &
5784	ok_ade, ok_aie, ivap, ok_sync, &
5785	ptconv, read_climoz, clevSTD, &
5786	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5787	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5788	ELSE
5789	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5790	pplay, lmax_th, aerosol_couple, &
5791	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs, &
5792	ok_sync, ptconv, read_climoz, clevSTD, &
5793	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5794	flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
5795	END IF
5796
5797	#ifndef CPP_XIOS
5798	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5799	#endif
5800
5801	! Petit appelle de sorties pour accompagner le travail sur phyex
5802	if ( iflag_physiq == 1 ) then
5803	call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.t,0.t,0.*t,pbl_tke,theta)
5804	endif
5805
5806	!====================================================================
5807	! Arret du modele apres hgardfou en cas de detection d'un
5808	! plantage par hgardfou
5809	!====================================================================
5810
5811	IF (abortphy==1) THEN
5812	abort_message ='Plantage hgardfou'
5813	CALL abort_physic (modname,abort_message,1)
5814	ENDIF
5815
5816	! 22.03.04 END
5817	!
5818	!====================================================================
5819	! Si c'est la fin, il faut conserver l'etat de redemarrage
5820	!====================================================================
5821	!
5822
5823	! Disabling calls to the prt_alerte function
5824	alert_first_call = .FALSE.
5825
5826
5827	IF (lafin) THEN
5828	itau_phy = itau_phy + itap
5829	CALL phyredem ("restartphy.nc")
5830	! open(97,form="unformatted",file="finbin")
5831	! write(97) u_seri,v_seri,t_seri,q_seri
5832	! close(97)
5833
5834	IF (is_omp_master) THEN
5835
5836	IF (read_climoz >= 1) THEN
5837	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5838	DEALLOCATE(press_edg_climoz)
5839	DEALLOCATE(press_cen_climoz)
5840	ENDIF
5841
5842	ENDIF
5843
5844	IF (using_xios) THEN
5845
5846	IF (CPPKEY_INCA) THEN
5847	IF (type_trac == 'inca') THEN
5848	IF (is_omp_master .AND. grid_type==unstructured) THEN
5849	CALL finalize_inca
5850	ENDIF
5851	ENDIF
5852	END IF
5853
5854	! close xios physiq context (call LMDZ)
5855	IF (is_omp_master) CALL xios_context_finalize
5856	ENDIF
5857
5858	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5859
5860	ENDIF
5861
5862	! first=.false.
5863
5864	END SUBROUTINE physiq
5865
5866	END MODULE physiq_mod

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

Download in other formats:

Original Format