Context Navigation

physiq_mod.F90 @ 5489

Last change on this file since 5489 was 5489, checked in by aborella, 13 days ago
Merge with trunk
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 216.9 KB

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

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

Download in other formats:

Original Format