Context Navigation

physiq_mod.F90 @ 5396

Last change on this file since 5396 was 5396, checked in by evignon, 5 weeks ago
ajout de ql_seri et qi_seri dans lmdz_lscp
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: 214.4 KB

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

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

Download in other formats:

Original Format