Context Navigation

physiq_mod.F90 @ 5499

Last change on this file since 5499 was 5499, checked in by yann meurdesoif, 18 hours ago
Adapt calwake and wake for automatic GPU port. YM
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: 215.6 KB

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

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

Download in other formats:

Original Format