Context Navigation

physiq_mod.F90 @ 5480

Last change on this file since 5480 was 5480, checked in by yann meurdesoif, 14 hours ago

For GPU porting :

add wrapper for source to source tools
Separate initialization phase of lscp_old (firstilp), SAVE variable in compute subroutine cannot be managed properly.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 215.6 KB

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

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

Download in other formats:

Original Format