Context Navigation

physiq_mod.F90 @ 5428

Last change on this file since 5428 was 5402, checked in by evignon, 2 weeks ago
petite correction liee a la specificite des variables ozone
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 214.5 KB

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

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

Download in other formats:

Original Format