Context Navigation

physiq_mod.F90 @ 5523

Last change on this file since 5523 was 5518, checked in by Laurent Fairhead, 5 days ago
Remains of replay mode. Will not work with DYNAMICO as is LF/AB
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 215.7 KB

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

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

Download in other formats:

Original Format