Context Navigation

physiq_mod.F90 @ 5491

Last change on this file since 5491 was 5491, checked in by jyg, 5 hours ago

New outputs :

+ coef_clos = [conv mass flux given by Alp closure]/[conv mass flux given by Emanuel scheme closure]
+ coef_clos_eff = effective coefficient used in the convective scheme.

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.5 KB

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

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

Download in other formats:

Original Format