Context Navigation

physiq_mod.F90 @ 5649

Last change on this file since 5649 was 5645, checked in by evignon, 2 months ago
Le nom du parametre new_oliq est tres trompeur. Je le change en liqice_in_radocond
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.4 KB

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

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

Download in other formats:

Original Format