Context Navigation

physiq_mod.F90 @ 5555

Last change on this file since 5555 was 5554, checked in by asima, 3 days ago
Encapsulating in modules 3 Dust (SPLA) subroutines with identical names in INCA (see also r5505); Related changes in other files, especially some cleaning of chem_mod_f.90. Everything cf the new coding conventions bien sûr !
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.9 KB

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

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

Download in other formats:

Original Format