Context Navigation

physiq_mod.F90 @ 5506

Last change on this file since 5506 was 5506, checked in by evignon, 12 hours ago
petite correction pour le bilan d'eau en cas de neige soufflee du au changement de convention des flux
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.7 KB

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

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

Download in other formats:

Original Format