Context Navigation

physiq_mod.F90 @ 5445

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

Line
1
2	! $Id: physiq_mod.F90 5431 2024-12-19 14:46:10Z 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, wdtrainAS, &
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	wdtrainAS(:,:) = 0.
3087	upwd(:,:) = 0.
3088	dnwd(:,:) = 0.
3089	ep(:,:) = 0.
3090	da(:,:)=0.
3091	mp(:,:)=0.
3092	wght_cvfd(:,:)=0.
3093	phi(:,:,:)=0.
3094	phi2(:,:,:)=0.
3095	epmlmMm(:,:,:)=0.
3096	eplaMm(:,:)=0.
3097	d1a(:,:)=0.
3098	dam(:,:)=0.
3099	elij(:,:,:)=0.
3100	ev(:,:)=0.
3101	qtaa(:,:)=0.
3102	clw(:,:)=0.
3103	sij(:,:,:)=0.
3104	!
3105	IF (iflag_con.EQ.1) THEN
3106	abort_message ='reactiver le call conlmd dans physiq.F'
3107	CALL abort_physic (modname,abort_message,1)
3108	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
3109	! . d_t_con, d_q_con,
3110	! . rain_con, snow_con, ibas_con, itop_con)
3111	ELSE IF (iflag_con.EQ.2) THEN
3112	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
3113	conv_t, conv_q, -evap, omega, &
3114	d_t_con, d_q_con, rain_con, snow_con, &
3115	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3116	kcbot, kctop, kdtop, pmflxr, pmflxs)
3117	d_u_con = 0.
3118	d_v_con = 0.
3119
3120	WHERE (rain_con < 0.) rain_con = 0.
3121	WHERE (snow_con < 0.) snow_con = 0.
3122	DO i = 1, klon
3123	ibas_con(i) = klev+1 - kcbot(i)
3124	itop_con(i) = klev+1 - kctop(i)
3125	ENDDO
3126	ELSE IF (iflag_con.GE.3) THEN
3127	! nb of tracers for the KE convection:
3128	! MAF la partie traceurs est faite dans phytrac
3129	! on met ntra=1 pour limiter les appels mais on peut
3130	! supprimer les calculs / ftra.
3131	ntra = 1
3132
3133	!=======================================================================
3134	!ajout pour la parametrisation des poches froides: calcul de
3135	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
3136	IF (iflag_wake>=1) THEN
3137	DO k=1,klev
3138	DO i=1,klon
3139	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
3140	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
3141	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3142	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3143	ENDDO
3144	ENDDO
3145	ELSE
3146	t_w(:,:) = t_seri(:,:)
3147	q_w(:,:) = q_seri(:,:)
3148	t_x(:,:) = t_seri(:,:)
3149	q_x(:,:) = q_seri(:,:)
3150	ENDIF
3151	!
3152	!jyg<
3153	! Perform dry adiabatic adjustment on wake profile
3154	! The corresponding tendencies are added to the convective tendencies
3155	! after the call to the convective scheme.
3156	IF (iflag_wake>=1) then
3157	IF (iflag_adjwk >= 1) THEN
3158	limbas(:) = 1
3159	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
3160	d_t_adjwk, d_q_adjwk)
3161	!
3162	DO k=1,klev
3163	DO i=1,klon
3164	IF (wake_s(i) .GT. 1.e-3) THEN
3165	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
3166	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
3167	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
3168	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
3169	ELSE
3170	d_deltat_ajs_cv(i,k) = 0.
3171	d_deltaq_ajs_cv(i,k) = 0.
3172	ENDIF
3173	ENDDO
3174	ENDDO
3175	IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
3176	CALL add_wake_tend &
3177	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
3178	ENDIF ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
3179	ENDIF ! (iflag_adjwk >= 1)
3180	ENDIF ! (iflag_wake>=1)
3181	!>jyg
3182	!
3183
3184	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
3185	!! (k, q_w(1,k), q_x(1,k),k=1,25)
3186
3187	!jyg<
3188	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
3189	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
3190	ale_bl_prescr, alp_bl_prescr, &
3191	wake_pe, wake_fip, &
3192	Ale_bl, Ale_bl_trig, Alp_bl, &
3193	Ale, Alp , Ale_wake, Alp_wake)
3194	!>jyg
3195	!
3196	! sb, oct02:
3197	! Schema de convection modularise et vectorise:
3198	! (driver commun aux versions 3 et 4)
3199	!
3200	IF (ok_cvl) THEN ! new driver for convectL
3201	!
3202	!jyg<
3203	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3204	! Calculate the upmost level of deep convection loops: k_upper_cv
3205	! (near 22 km)
3206	k_upper_cv = klev
3207	!izero = klon/2+1/klon
3208	!DO k = klev,1,-1
3209	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
3210	!ENDDO
3211	! FH : nouveau calcul base sur un profil global sans quoi
3212	! le modele etait sensible au decoupage de domaines
3213	DO k = klev,1,-1
3214	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
3215	ENDDO
3216	IF (prt_level .ge. 5) THEN
3217	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
3218	k_upper_cv
3219	ENDIF
3220	!
3221	!>jyg
3222	IF (type_trac == 'repr') THEN
3223	nbtr_tmp=ntra
3224	ELSE
3225	nbtr_tmp=nbtr
3226	ENDIF
3227	!jyg iflag_con est dans clesphys
3228	!c CALL concvl (iflag_con,iflag_clos,
3229	CALL concvl (iflag_clos, &
3230	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
3231	t_w,q_w,wake_s, &
3232	u_seri,v_seri,tr_seri,nbtr_tmp, &
3233	ALE,ALP, &
3234	sig1,w01, &
3235	d_t_con,d_q_con,fqcomp,d_u_con,d_v_con,d_tr, &
3236	rain_con, snow_con, ibas_con, itop_con, sigd, &
3237	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
3238	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
3239	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
3240	! RomP >>>
3241	!! . pmflxr,pmflxs,da,phi,mp,
3242	!! . ftd,fqd,lalim_conv,wght_th)
3243	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
3244	ftd,fqd,lalim_conv,wght_th, &
3245	ev, ep,epmlmMm,eplaMm, &
3246	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv,detrain_cv, &
3247	tau_cld_cv,coefw_cld_cv,epmax_diag)
3248
3249	! RomP <<<
3250
3251	!IM begin
3252	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
3253	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
3254	!IM end
3255	!IM cf. FH
3256	clwcon0=qcondc
3257	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
3258	fm_cv(:,:)=upwd(:,:)+dnwd(:,:)+dnwd0(:,:)
3259	!
3260	!jyg<
3261	! If convective tendencies are too large, then call convection
3262	! every time step
3263	cvpas = cvpas_0
3264	DO k=1,k_upper_cv
3265	DO i=1,klon
3266	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3267	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3268	dtcon_multistep_max = 3.
3269	dqcon_multistep_max = 0.02
3270	ENDIF
3271	ENDDO
3272	ENDDO
3273	!
3274	DO k=1,k_upper_cv
3275	DO i=1,klon
3276	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3277	!! abs(d_q_con(i,k)) > 2.e-2) THEN
3278	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3279	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3280	cvpas = 1
3281	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3282	!! i,k,d_t_con(i,k),d_q_con(i,k)
3283	ENDIF
3284	ENDDO
3285	ENDDO
3286	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3287	!!! call bcast(cvpas)
3288	!!! ------------------------------------------------------------
3289	!>jyg
3290	!
3291	DO i = 1, klon
3292	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
3293	ENDDO
3294	!
3295	!jyg<
3296	! Add the tendency due to the dry adjustment of the wake profile
3297	IF (iflag_wake>=1) THEN
3298	IF (iflag_adjwk == 2) THEN
3299	DO k=1,klev
3300	DO i=1,klon
3301	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3302	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
3303	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3304	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3305	ENDDO
3306	ENDDO
3307	ENDIF ! (iflag_adjwk = 2)
3308	ENDIF ! (iflag_wake>=1)
3309	!>jyg
3310	!
3311	ELSE ! ok_cvl
3312
3313	! MAF conema3 ne contient pas les traceurs
3314	CALL conema3 (phys_tstep, &
3315	paprs,pplay,t_seri,q_seri, &
3316	u_seri,v_seri,tr_seri,ntra, &
3317	sig1,w01, &
3318	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3319	rain_con, snow_con, ibas_con, itop_con, &
3320	upwd,dnwd,dnwd0,bas,top, &
3321	Ma,cape,tvp,rflag, &
3322	pbase &
3323	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3324	,clwcon0)
3325
3326	ENDIF ! ok_cvl
3327
3328	!
3329	! Correction precip
3330	rain_con = rain_con * cvl_corr
3331	snow_con = snow_con * cvl_corr
3332	!
3333
3334	IF (.NOT. ok_gust) THEN
3335	do i = 1, klon
3336	wd(i)=0.0
3337	enddo
3338	ENDIF
3339
3340	! =================================================================== c
3341	! Calcul des proprietes des nuages convectifs
3342	!
3343
3344	! calcul des proprietes des nuages convectifs
3345	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3346	IF (iflag_cld_cv == 0) THEN
3347	CALL clouds_gno &
3348	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3349	ELSE
3350	CALL clouds_bigauss &
3351	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3352	ENDIF
3353
3354
3355	! =================================================================== c
3356
3357	DO i = 1, klon
3358	itop_con(i) = min(max(itop_con(i),1),klev)
3359	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3360	ENDDO
3361
3362	DO i = 1, klon
3363	! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3364	! où i n'est pas un point convectif et donc ibas_con(i)=0
3365	! c'est un pb indépendant des isotopes
3366	if (ibas_con(i) > 0) then
3367	ema_pcb(i) = paprs(i,ibas_con(i))
3368	else
3369	ema_pcb(i) = 0.0
3370	endif
3371	ENDDO
3372	DO i = 1, klon
3373	! L'idicage de itop_con peut cacher un pb potentiel
3374	! FH sous la dictee de JYG, CR
3375	ema_pct(i) = paprs(i,itop_con(i)+1)
3376
3377	IF (itop_con(i).gt.klev-3) THEN
3378	IF (prt_level >= 9) THEN
3379	write(lunout,*)'La convection monte trop haut '
3380	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
3381	ENDIF
3382	ENDIF
3383	ENDDO
3384	ELSE IF (iflag_con.eq.0) THEN
3385	write(lunout,*) 'On n appelle pas la convection'
3386	clwcon0=0.
3387	rnebcon0=0.
3388	d_t_con=0.
3389	d_q_con=0.
3390	d_u_con=0.
3391	d_v_con=0.
3392	rain_con=0.
3393	snow_con=0.
3394	bas=1
3395	top=1
3396	ELSE
3397	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
3398	CALL abort_physic("physiq", "", 1)
3399	ENDIF
3400
3401	!--saving d_q_con * zmass for next timestep if convection is not called every timestep
3402	IF (ok_conserv_d_q_con) THEN
3403	d_q_con_zmasse(:,:) = d_q_con(:,:) * zmasse(:,:)
3404	ENDIF
3405
3406	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3407	! . d_u_con, d_v_con)
3408
3409	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3410	proba_notrig(:) = 1.
3411	itapcv = 0
3412	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3413	!
3414	itapcv = itapcv+1
3415	!
3416	! Compter les steps ou cvpas=1
3417	IF (cvpas == 1) THEN
3418	Ncvpaseq1 = Ncvpaseq1+1
3419	ENDIF
3420	IF (mod(itap,1000) == 0) THEN
3421	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3422	ENDIF
3423
3424	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3425	!!! l'energie dans les courants satures.
3426	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3427	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3428	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3429	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3430	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3431	!! itap, 1)
3432	!! call prt_enerbil('convection_sat',itap)
3433	!!
3434	!!
3435
3436	!--recompute d_q_con with zmasse from new timestep
3437	IF (ok_conserv_d_q_con) THEN
3438	d_q_con(:,:)=d_q_con_zmasse(:,:)/zmasse(:,:)
3439	ENDIF
3440
3441	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, dqbs0, paprs, &
3442	'convection',abortphy,flag_inhib_tend,itap,0)
3443	CALL prt_enerbil('convection',itap)
3444
3445	!-------------------------------------------------------------------------
3446
3447	IF (mydebug) THEN
3448	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3449	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3450	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3451	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3452	ENDIF
3453
3454	!
3455	!==========================================================================
3456	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3457	!pour la couche limite diffuse pour l instant
3458	!
3459	!
3460	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3461	! il faut rajouter cette tendance calcul\'ee hors des poches
3462	! froides
3463	!
3464	IF (iflag_wake>=1) THEN
3465	!
3466	!
3467	! Call wakes every "wkpas" step
3468	!
3469	IF (MOD(itapwk,wkpas).EQ.0) THEN
3470	!
3471	DO k=1,klev
3472	DO i=1,klon
3473	dt_dwn(i,k) = ftd(i,k)
3474	dq_dwn(i,k) = fqd(i,k)
3475	M_dwn(i,k) = dnwd0(i,k)
3476	M_up(i,k) = upwd(i,k)
3477	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3478	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3479	ENDDO
3480	ENDDO
3481
3482	IF (iflag_wake==2) THEN
3483	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3484	DO k = 1,klev
3485	dt_dwn(:,k)= dt_dwn(:,k)+ &
3486	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3487	dq_dwn(:,k)= dq_dwn(:,k)+ &
3488	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3489	ENDDO
3490	ELSEIF (iflag_wake==3) THEN
3491	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3492	DO k = 1,klev
3493	DO i=1,klon
3494	IF (rneb(i,k)==0.) THEN
3495	! On ne tient compte des tendances qu'en dehors des
3496	! nuages (c'est-\`a-dire a priri dans une region ou
3497	! l'eau se reevapore).
3498	dt_dwn(i,k)= dt_dwn(i,k)+ &
3499	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3500	dq_dwn(i,k)= dq_dwn(i,k)+ &
3501	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3502	ENDIF
3503	ENDDO
3504	ENDDO
3505	ENDIF
3506
3507	!
3508	!calcul caracteristiques de la poche froide
3509	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3510	t_seri, q_seri, omega, &
3511	dt_dwn, dq_dwn, M_dwn, M_up, &
3512	dt_a, dq_a, cv_gen, &
3513	sigd, cin, &
3514	wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, &
3515	wake_dth, wake_h, &
3516	!! wake_pe, wake_fip, wake_gfl, &
3517	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3518	d_t_wake, d_q_wake, &
3519	wake_k, t_x, q_x, &
3520	wake_omgbdth, wake_dp_omgb, &
3521	wake_dtKE, wake_dqKE, &
3522	wake_omg, wake_dp_deltomg, &
3523	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3524	d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk)
3525	!
3526	!jyg Reinitialize itapwk when wakes have been called
3527	itapwk = 0
3528	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3529	!
3530	itapwk = itapwk+1
3531	!
3532	!-----------------------------------------------------------------------
3533	! ajout des tendances des poches froides
3534	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,dqbs0,paprs,'wake', &
3535	abortphy,flag_inhib_tend,itap,0)
3536	CALL prt_enerbil('wake',itap)
3537	!------------------------------------------------------------------------
3538
3539	! Increment Wake state variables
3540	IF (iflag_wake_tend .GT. 0.) THEN
3541
3542	CALL add_wake_tend &
3543	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &
3544	'wake', abortphy)
3545	CALL prt_enerbil('wake',itap)
3546	ENDIF ! (iflag_wake_tend .GT. 0.)
3547	!
3548	IF (prt_level .GE. 10) THEN
3549	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3550	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3551	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3552	ENDIF
3553
3554	IF (iflag_alp_wk_cond .GT. 0.) THEN
3555
3556	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3557	wake_fip)
3558	ELSE
3559	wake_fip(:) = wake_fip_0(:)
3560	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3561
3562	ENDIF ! (iflag_wake>=1)
3563	!
3564	!===================================================================
3565	! Convection seche (thermiques ou ajustement)
3566	!===================================================================
3567	!
3568	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3569	,seuil_inversion,weak_inversion,dthmin)
3570
3571
3572
3573	d_t_ajsb(:,:)=0.
3574	d_q_ajsb(:,:)=0.
3575	d_t_ajs(:,:)=0.
3576	d_u_ajs(:,:)=0.
3577	d_v_ajs(:,:)=0.
3578	d_q_ajs(:,:)=0.
3579	clwcon0th(:,:)=0.
3580	!
3581	! fm_therm(:,:)=0.
3582	! entr_therm(:,:)=0.
3583	! detr_therm(:,:)=0.
3584	!
3585	IF (prt_level>9) WRITE(lunout,*) &
3586	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3587	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3588	IF (iflag_thermals<0) THEN
3589	! Rien
3590	! ====
3591	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3592	WRITE(lunout,*) 'WARNING : running without dry convection. Somme intermediate variables are not properly defined in physiq_mod.F90'
3593	! Reprendre proprement les initialisation ci dessouds si on veut vraiment utiliser l'option (FH)
3594	fraca(:,:)=0.
3595	fm_therm(:,:)=0.
3596	ztv(:,:)=t_seri(:,:)
3597	zqasc(:,:)=q_seri(:,:)
3598	ztla(:,:)=0.
3599	zthl(:,:)=0.
3600	zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
3601
3602
3603
3604	ELSE
3605
3606	! Thermiques
3607	! ==========
3608	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3609	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3610
3611
3612	!cc nrlmd le 10/04/2012
3613	DO k=1,klev+1
3614	DO i=1,klon
3615	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3616	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3617	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3618	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3619	ENDDO
3620	ENDDO
3621	!cc fin nrlmd le 10/04/2012
3622
3623	IF (iflag_thermals>=1) THEN
3624
3625	! Tests Fredho, instensibilite au pas de temps -------------------------------
3626	! A detruire en 2024 une fois les tests documentes et les choix faits !
3627	if (iflag_thermals_tenv /10 == 0 ) then !
3628	do k=1,klev !
3629	do i=1,klon !
3630	t_env(i,k)=t_seri(i,k) !
3631	q_env(i,k)=q_seri(i,k) !
3632	enddo !
3633	enddo !
3634	else if (iflag_thermals_tenv / 10 == 2 ) then !
3635	do k=1,klev !
3636	do i=1,klon !
3637	q_env(i,k)=q_seri(i,k) !
3638	enddo !
3639	enddo !
3640	else if (iflag_thermals_tenv / 10 == 3 ) then !
3641	do k=1,klev !
3642	do i=1,klon !
3643	t_env(i,k)=t(i,k) !
3644	q_env(i,k)=qx(i,k,1) !
3645	enddo !
3646	enddo !
3647	endif !
3648	! Tests Fredho, instensibilite au pas de temps ------------------------------
3649
3650	!jyg<
3651	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3652	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3653	! Appel des thermiques avec les profils exterieurs aux poches
3654	DO k=1,klev
3655	DO i=1,klon
3656	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3657	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3658	t_env(i,k) = t_env(i,k) - wake_s(i)*wake_deltat(i,k)
3659	q_env(i,k) = q_env(i,k) - wake_s(i)*wake_deltaq(i,k)
3660	u_therm(i,k) = u_seri(i,k)
3661	v_therm(i,k) = v_seri(i,k)
3662	ENDDO
3663	ENDDO
3664	ELSE
3665	! Appel des thermiques avec les profils moyens
3666	DO k=1,klev
3667	DO i=1,klon
3668	t_therm(i,k) = t_seri(i,k)
3669	q_therm(i,k) = q_seri(i,k)
3670	u_therm(i,k) = u_seri(i,k)
3671	v_therm(i,k) = v_seri(i,k)
3672	ENDDO
3673	ENDDO
3674	ENDIF
3675	!>jyg
3676	CALL calltherm(pdtphys &
3677	,pplay,paprs,pphi,weak_inversion &
3678	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3679	,u_therm,v_therm,t_therm,q_therm,t_env,q_env,zqsat,debut & !jyg
3680	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3681	,fm_therm,entr_therm,detr_therm &
3682	,zqasc,clwcon0th,lmax_th,ratqscth &
3683	,ratqsdiff,zqsatth &
3684	!on rajoute ale et alp, et les
3685	!caracteristiques de la couche alim
3686	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3687	,ztv,zpspsk,ztla,zthl &
3688	!cc nrlmd le 10/04/2012
3689	,pbl_tke_input,pctsrf,omega,cell_area &
3690	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3691	,n2,s2,strig,zcong,ale_bl_stat &
3692	,therm_tke_max,env_tke_max &
3693	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3694	,alp_bl_conv,alp_bl_stat &
3695	!cc fin nrlmd le 10/04/2012
3696	,zqla,ztva )
3697	!
3698	!jyg<
3699	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3700	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3701	! Si les thermiques ne sont presents que hors des
3702	! poches, la tendance moyenne associ\'ee doit etre
3703	! multipliee par la fraction surfacique qu'ils couvrent.
3704	DO k=1,klev
3705	DO i=1,klon
3706	!
3707	d_deltat_the(i,k) = - d_t_ajs(i,k)
3708	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3709	!
3710	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3711	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3712	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3713	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3714	!
3715	ENDDO
3716	ENDDO
3717	!
3718	IF (ok_bug_split_th) THEN
3719	CALL add_wake_tend &
3720	(d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3721	ELSE
3722	CALL add_wake_tend &
3723	(d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3724	ENDIF
3725	CALL prt_enerbil('the',itap)
3726	!
3727	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3728	!
3729	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3730	dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3731	CALL prt_enerbil('thermals',itap)
3732	!
3733	!
3734	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3735	cin, s2, n2, strig, &
3736	ale_bl_trig, ale_bl_stat, ale_bl, &
3737	alp_bl, alp_bl_stat, &
3738	proba_notrig, random_notrig, cv_gen)
3739	!>jyg
3740
3741	! ------------------------------------------------------------------
3742	! Transport de la TKE par les panaches thermiques.
3743	! FH : 2010/02/01
3744	if (iflag_thermcell_tke==1) then
3745	call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,rg,paprs,pbl_tke)
3746	endif
3747	! -------------------------------------------------------------------
3748
3749	DO i=1,klon
3750	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3751	!CR:04/05/12:correction calcul zmax
3752	zmax_th(i)=zmax0(i)
3753	ENDDO
3754
3755	ENDIF
3756
3757	! Ajustement sec
3758	! ==============
3759
3760	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3761	! a partir du sommet des thermiques.
3762	! Dans le cas contraire, on demarre au niveau 1.
3763
3764	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3765
3766	IF (iflag_thermals.eq.0) THEN
3767	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3768	limbas(:)=1
3769	ELSE
3770	limbas(:)=lmax_th(:)
3771	ENDIF
3772
3773	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3774	! pour des test de convergence numerique.
3775	! Le nouveau ajsec est a priori mieux, meme pour le cas
3776	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3777	! non nulles numeriquement pour des mailles non concernees.
3778
3779	IF (iflag_thermals==0) THEN
3780	! Calling adjustment alone (but not the thermal plume model)
3781	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3782	, d_t_ajsb, d_q_ajsb)
3783	ELSE IF (iflag_thermals>0) THEN
3784	! Calling adjustment above the top of thermal plumes
3785	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3786	, d_t_ajsb, d_q_ajsb)
3787	ENDIF
3788
3789	!--------------------------------------------------------------------
3790	! ajout des tendances de l'ajustement sec ou des thermiques
3791	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,dqbs0,paprs, &
3792	'ajsb',abortphy,flag_inhib_tend,itap,0)
3793	CALL prt_enerbil('ajsb',itap)
3794	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3795	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3796
3797	!---------------------------------------------------------------------
3798
3799	ENDIF
3800
3801	ENDIF
3802	!
3803	!===================================================================
3804	! Computation of ratqs, the width (normalized) of the subrid scale
3805	! water distribution
3806
3807	l_mix_ave(:,:)=0.
3808	wprime_ave(:,:)=0.
3809
3810	DO nsrf = 1, nbsrf
3811	DO i = 1, klon
3812	l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3813	wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3814	ENDDO
3815	ENDDO
3816
3817	CALL ratqs_main(klon,klev,nbsrf,prt_level,lunout, &
3818	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3819	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3820	pctsrf,s_pblh,zstd, &
3821	tau_ratqs,fact_cldcon,wake_s, wake_deltaq, &
3822	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3823	paprs,pplay,t_seri,q_seri, &
3824	qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
3825	omega,pbl_tke(:,:,is_ave),pbl_eps(:,:,is_ave),l_mix_ave,wprime_ave, &
3826	t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
3827	ratqs,ratqsc,ratqs_inter_,sigma_qtherm)
3828
3829	!
3830	! Appeler le processus de condensation a grande echelle
3831	! et le processus de precipitation
3832	!-------------------------------------------------------------------------
3833	IF (prt_level .GE.10) THEN
3834	print *,'itap, ->fisrtilp ',itap
3835	ENDIF
3836	!
3837
3838	picefra(:,:)=0.
3839
3840	IF (ok_new_lscp) THEN
3841
3842
3843	DO k = 1, klev
3844	DO i = 1, klon
3845	ql_seri_lscp(i,k) = ratio_ql_qtot(i,k) * q_seri(i,k)
3846	qi_seri_lscp(i,k) = ratio_qi_qtot(i,k) * q_seri(i,k)
3847	ENDDO
3848	ENDDO
3849
3850
3851	!--mise à jour de flight_m et flight_h2o dans leur module
3852	!IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3853	! CALL airplane(debut,pphis,pplay,paprs,t_seri)
3854	!ENDIF
3855
3856	CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
3857	t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ptconv, ratqs, sigma_qtherm, &
3858	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
3859	pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
3860	radocond, picefra, rain_lsc, snow_lsc, &
3861	frac_impa, frac_nucl, beta_prec_fisrt, &
3862	prfl, psfl, rhcl, &
3863	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3864	iflag_ice_thermo, distcltop, temp_cltop, &
3865	pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
3866	cell_area, &
3867	cf_seri, rvc_seri, u_seri, v_seri, &
3868	qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
3869	dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
3870	dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
3871	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
3872	dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
3873	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
3874	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
3875	dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
3876
3877
3878	ELSE
3879
3880	CALL fisrtilp(klon,klev,phys_tstep,paprs,pplay, &
3881	t_seri, q_seri,ptconv,ratqs,sigma_qtherm, &
3882	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, radocond, &
3883	rain_lsc, snow_lsc, &
3884	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3885	frac_impa, frac_nucl, beta_prec_fisrt, &
3886	prfl, psfl, rhcl, &
3887	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3888	iflag_ice_thermo, &
3889	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
3890
3891	ENDIF
3892	!
3893	WHERE (rain_lsc < 0) rain_lsc = 0.
3894	WHERE (snow_lsc < 0) snow_lsc = 0.
3895
3896	!+JLD
3897	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3898	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3899	! & ,rain_lsc,snow_lsc
3900	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3901	!-JLD
3902	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
3903	'lsc',abortphy,flag_inhib_tend,itap,0)
3904	CALL prt_enerbil('lsc',itap)
3905	rain_num(:)=0.
3906	DO k = 1, klev
3907	DO i = 1, klon
3908	IF (ql_seri(i,k)>oliqmax) THEN
3909	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3910	ql_seri(i,k)=oliqmax
3911	ENDIF
3912	ENDDO
3913	ENDDO
3914	IF (nqo >= 3) THEN
3915	DO k = 1, klev
3916	DO i = 1, klon
3917	IF (qs_seri(i,k)>oicemax) THEN
3918	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3919	qs_seri(i,k)=oicemax
3920	ENDIF
3921	ENDDO
3922	ENDDO
3923	ENDIF
3924
3925
3926	!---------------------------------------------------------------------------
3927	DO k = 1, klev
3928	DO i = 1, klon
3929	cldfra(i,k) = rneb(i,k)
3930	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3931	!EV: en effet etrange, j'ajouterais aussi qs_seri
3932	! plus largement, je nettoierais (enleverrais) ces lignes
3933	IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
3934	ENDDO
3935	ENDDO
3936
3937
3938	! Option to activate the radiative effect of blowing snow (ok_rad_bs)
3939	! makes sense only if the new large scale condensation scheme is active
3940	! with the ok_icefra_lscp flag active as well
3941
3942	IF (ok_bs .AND. ok_rad_bs) THEN
3943	IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
3944	DO k=1,klev
3945	DO i=1,klon
3946	radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
3947	picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
3948	qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
3949	cldfra(i,k)=max(cldfra(i,k),qbsfra)
3950	ENDDO
3951	ENDDO
3952	ELSE
3953	WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
3954	WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
3955	abort_message='inconsistency in cloud phase for blowing snow'
3956	CALL abort_physic(modname,abort_message,1)
3957	ENDIF
3958
3959	ENDIF
3960
3961	IF (mydebug) THEN
3962	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3963	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3964	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3965	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3966	ENDIF
3967
3968	!
3969	!-------------------------------------------------------------------
3970	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3971	!-------------------------------------------------------------------
3972
3973	! 1. NUAGES CONVECTIFS
3974	!
3975	!IM cf FH
3976	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3977	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3978	snow_tiedtke=0.
3979	! print*,'avant calcul de la pseudo precip '
3980	! print*,'iflag_cld_th',iflag_cld_th
3981	IF (iflag_cld_th.eq.-1) THEN
3982	rain_tiedtke=rain_con
3983	ELSE
3984	! print*,'calcul de la pseudo precip '
3985	rain_tiedtke=0.
3986	! print*,'calcul de la pseudo precip 0'
3987	DO k=1,klev
3988	DO i=1,klon
3989	IF (d_q_con(i,k).lt.0.) THEN
3990	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3991	*(paprs(i,k)-paprs(i,k+1))/rg
3992	ENDIF
3993	ENDDO
3994	ENDDO
3995	ENDIF
3996	!
3997	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3998	!
3999
4000	! Nuages diagnostiques pour Tiedtke
4001	CALL diagcld1(paprs,pplay, &
4002	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
4003	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
4004	diafra,dialiq)
4005	DO k = 1, klev
4006	DO i = 1, klon
4007	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4008	radocond(i,k) = dialiq(i,k)
4009	cldfra(i,k) = diafra(i,k)
4010	ENDIF
4011	ENDDO
4012	ENDDO
4013
4014	ELSE IF (iflag_cld_th.ge.3) THEN
4015	! On prend pour les nuages convectifs le max du calcul de la
4016	! convection et du calcul du pas de temps precedent diminue d'un facteur
4017	! facttemps
4018	facteur = pdtphys *facttemps
4019	DO k=1,klev
4020	DO i=1,klon
4021	rnebcon(i,k)=rnebcon(i,k)*facteur
4022	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
4023	rnebcon(i,k)=rnebcon0(i,k)
4024	clwcon(i,k)=clwcon0(i,k)
4025	ENDIF
4026	ENDDO
4027	ENDDO
4028
4029	! On prend la somme des fractions nuageuses et des contenus en eau
4030
4031	IF (iflag_cld_th>=5) THEN
4032
4033	DO k=1,klev
4034	ptconvth(:,k)=fm_therm(:,k+1)>0.
4035	ENDDO
4036
4037	IF (iflag_coupl==4) THEN
4038
4039	! Dans le cas iflag_coupl==4, on prend la somme des convertures
4040	! convectives et lsc dans la partie des thermiques
4041	! Le controle par iflag_coupl est peut etre provisoire.
4042	DO k=1,klev
4043	DO i=1,klon
4044	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
4045	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4046	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4047	ELSE IF (ptconv(i,k)) THEN
4048	cldfra(i,k)=rnebcon(i,k)
4049	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4050	ENDIF
4051	ENDDO
4052	ENDDO
4053
4054	ELSE IF (iflag_coupl==5) THEN
4055	DO k=1,klev
4056	DO i=1,klon
4057	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4058	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4059	ENDDO
4060	ENDDO
4061
4062	ELSE
4063
4064	! Si on est sur un point touche par la convection
4065	! profonde et pas par les thermiques, on prend la
4066	! couverture nuageuse et l'eau nuageuse de la convection
4067	! profonde.
4068
4069	!IM/FH: 2011/02/23
4070	! definition des points sur lesquels ls thermiques sont actifs
4071
4072	DO k=1,klev
4073	DO i=1,klon
4074	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
4075	cldfra(i,k)=rnebcon(i,k)
4076	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4077	ENDIF
4078	ENDDO
4079	ENDDO
4080
4081	ENDIF
4082
4083	ELSE
4084
4085	! Ancienne version
4086	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
4087	radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
4088	ENDIF
4089
4090	ENDIF
4091
4092	! plulsc(:)=0.
4093	! do k=1,klev,-1
4094	! do i=1,klon
4095	! zzz=prfl(:,k)+psfl(:,k)
4096	! if (.not.ptconvth.zzz.gt.0.)
4097	! enddo prfl, psfl,
4098	! enddo
4099	!
4100	! 2. NUAGES STARTIFORMES
4101	!
4102	IF (ok_stratus) THEN
4103	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
4104	DO k = 1, klev
4105	DO i = 1, klon
4106	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4107	radocond(i,k) = dialiq(i,k)
4108	cldfra(i,k) = diafra(i,k)
4109	ENDIF
4110	ENDDO
4111	ENDDO
4112	ENDIF
4113	!
4114	! Precipitation totale
4115	!
4116	DO i = 1, klon
4117	rain_fall(i) = rain_con(i) + rain_lsc(i)
4118	snow_fall(i) = snow_con(i) + snow_lsc(i)
4119	ENDDO
4120	!
4121	! Calculer l'humidite relative pour diagnostique
4122	!
4123	DO k = 1, klev
4124	DO i = 1, klon
4125	zx_t = t_seri(i,k)
4126	IF (thermcep) THEN
4127	!! if (iflag_ice_thermo.eq.0) then !jyg
4128	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
4129	!! else !jyg
4130	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
4131	!! endif !jyg
4132	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
4133	zx_qs = MIN(0.5,zx_qs)
4134	zcor = 1./(1.-retv*zx_qs)
4135	zx_qs = zx_qs*zcor
4136	ELSE
4137	!! IF (zx_t.LT.t_coup) THEN !jyg
4138	IF (zx_t.LT.rtt) THEN !jyg
4139	zx_qs = qsats(zx_t)/pplay(i,k)
4140	ELSE
4141	zx_qs = qsatl(zx_t)/pplay(i,k)
4142	ENDIF
4143	ENDIF
4144	zx_rh(i,k) = q_seri(i,k)/zx_qs
4145	IF (iflag_ice_thermo .GT. 0) THEN
4146	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
4147	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
4148	ENDIF
4149	zqsat(i,k)=zx_qs
4150	ENDDO
4151	ENDDO
4152
4153	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
4154	! equivalente a 2m (tpote) pour diagnostique
4155	!
4156	DO i = 1, klon
4157	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
4158	IF (thermcep) THEN
4159	IF(zt2m(i).LT.RTT) then
4160	Lheat=RLSTT
4161	ELSE
4162	Lheat=RLVTT
4163	ENDIF
4164	ELSE
4165	IF (zt2m(i).LT.RTT) THEN
4166	Lheat=RLSTT
4167	ELSE
4168	Lheat=RLVTT
4169	ENDIF
4170	ENDIF
4171	tpote(i) = tpot(i)* &
4172	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
4173	ENDDO
4174
4175	IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
4176	IF (CPPKEY_INCA) THEN
4177	CALL VTe(VTphysiq)
4178	CALL VTb(VTinca)
4179	calday = REAL(days_elapsed + 1) + jH_cur
4180
4181	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
4182	CALL AEROSOL_METEO_CALC( &
4183	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
4184	prfl,psfl,pctsrf,cell_area, &
4185	latitude_deg,longitude_deg,u10m,v10m)
4186
4187	zxsnow_dummy(:) = 0.0
4188
4189	CALL chemhook_begin (calday, &
4190	days_elapsed+1, &
4191	jH_cur, &
4192	pctsrf(1,1), &
4193	latitude_deg, &
4194	longitude_deg, &
4195	cell_area, &
4196	paprs, &
4197	pplay, &
4198	coefh(1:klon,1:klev,is_ave), &
4199	pphi, &
4200	t_seri, &
4201	u, &
4202	v, &
4203	rot, &
4204	wo(:, :, 1), &
4205	q_seri, &
4206	zxtsol, &
4207	zt2m, &
4208	zxsnow_dummy, &
4209	solsw, &
4210	albsol1, &
4211	rain_fall, &
4212	snow_fall, &
4213	itop_con, &
4214	ibas_con, &
4215	cldfra, &
4216	nbp_lon, &
4217	nbp_lat-1, &
4218	tr_seri(:,:,1+nqCO2:nbtr), &
4219	ftsol, &
4220	paprs, &
4221	cdragh, &
4222	cdragm, &
4223	pctsrf, &
4224	pdtphys, &
4225	itap)
4226
4227	CALL VTe(VTinca)
4228	CALL VTb(VTphysiq)
4229	END IF
4230	ENDIF !type_trac = inca or inco
4231	IF (type_trac == 'repr') THEN
4232	IF (CPPKEY_REPROBUS) THEN
4233	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
4234	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
4235	END IF
4236	ENDIF
4237
4238	!
4239	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
4240	!
4241	IF (MOD(itaprad,radpas).EQ.0) THEN
4242
4243	!
4244	!jq - introduce the aerosol direct and first indirect radiative forcings
4245	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
4246	IF (flag_aerosol .GT. 0) THEN
4247	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4248	IF (.NOT. aerosol_couple) THEN
4249	!
4250	CALL readaerosol_optic( &
4251	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4252	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4253	mass_solu_aero, mass_solu_aero_pi, &
4254	tau_aero, piz_aero, cg_aero, &
4255	tausum_aero, tau3d_aero)
4256	ENDIF
4257	ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation
4258	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
4259	abort_message='config_inca=aero et rrtm=1 impossible'
4260	CALL abort_physic(modname,abort_message,1)
4261	ELSE
4262	!
4263	#ifdef CPP_RRTM
4264	IF (NSW.EQ.6) THEN
4265	!--new aerosol properties SW and LW
4266	!
4267	IF (CPPKEY_DUST) THEN
4268	!--SPL aerosol model
4269	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
4270	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4271	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4272	tausum_aero, tau3d_aero)
4273	ELSE
4274	!--climatologies or INCA aerosols
4275	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
4276	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4277	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4278	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4279	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4280	tausum_aero, drytausum_aero, tau3d_aero)
4281	END IF
4282
4283	IF (flag_aerosol .EQ. 7) THEN
4284	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
4285	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
4286	ENDIF
4287
4288	!
4289	ELSE IF (NSW.EQ.2) THEN
4290	!--for now we use the old aerosol properties
4291	!
4292	CALL readaerosol_optic( &
4293	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4294	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4295	mass_solu_aero, mass_solu_aero_pi, &
4296	tau_aero, piz_aero, cg_aero, &
4297	tausum_aero, tau3d_aero)
4298	!
4299	!--natural aerosols
4300	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4301	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4302	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4303	!--all aerosols
4304	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4305	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4306	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
4307	!
4308	!--no LW optics
4309	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4310	!
4311	ELSE
4312	abort_message='Only NSW=2 or 6 are possible with ' &
4313	// 'aerosols and iflag_rrtm=1'
4314	CALL abort_physic(modname,abort_message,1)
4315	ENDIF
4316	#else
4317	abort_message='You should compile with -rrtm if running ' &
4318	// 'with iflag_rrtm=1'
4319	CALL abort_physic(modname,abort_message,1)
4320	#endif
4321	!
4322	ENDIF
4323	ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION
4324	#ifdef CPP_ECRAD
4325	!--climatologies or INCA aerosols
4326	CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
4327	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4328	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4329	tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer)
4330	#else
4331	abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
4332	CALL abort_physic(modname,abort_message,1)
4333	#endif
4334	ENDIF
4335
4336	ELSE !--flag_aerosol = 0
4337	tausum_aero(:,:,:) = 0.
4338	drytausum_aero(:,:) = 0.
4339	mass_solu_aero(:,:) = 0.
4340	mass_solu_aero_pi(:,:) = 0.
4341	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4342	tau_aero(:,:,:,:) = 1.e-15
4343	piz_aero(:,:,:,:) = 1.
4344	cg_aero(:,:,:,:) = 0.
4345	ELSE
4346	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4347	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4348	piz_aero_sw_rrtm(:,:,:,:) = 1.0
4349	cg_aero_sw_rrtm(:,:,:,:) = 0.0
4350	ENDIF
4351	ENDIF
4352	!
4353	!--WMO criterion to determine tropopause
4354	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4355	!
4356	!--STRAT AEROSOL
4357	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
4358	IF (flag_aerosol_strat.GT.0) THEN
4359	IF (prt_level .GE.10) THEN
4360	PRINT *,'appel a readaerosolstrat', mth_cur
4361	ENDIF
4362	IF (iflag_rrtm.EQ.0) THEN
4363	IF (flag_aerosol_strat.EQ.1) THEN
4364	CALL readaerosolstrato(debut)
4365	ELSE
4366	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4367	CALL abort_physic(modname,abort_message,1)
4368	ENDIF
4369	ELSE
4370	#ifdef CPP_RRTM
4371	IF (.NOT. CPPKEY_STRATAER) THEN
4372	!--prescribed strat aerosols
4373	!--only in the case of non-interactive strat aerosols
4374	IF (flag_aerosol_strat.EQ.1) THEN
4375	CALL readaerosolstrato1_rrtm(debut)
4376	ELSEIF (flag_aerosol_strat.EQ.2) THEN
4377	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
4378	ELSE
4379	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4380	CALL abort_physic(modname,abort_message,1)
4381	ENDIF
4382	END IF
4383	#else
4384	abort_message='You should compile with -rrtm if running ' &
4385	// 'with iflag_rrtm=1'
4386	CALL abort_physic(modname,abort_message,1)
4387	#endif
4388	ENDIF
4389	ELSE
4390	tausum_aero(:,:,id_STRAT_phy) = 0.
4391	ENDIF
4392	!
4393	#ifdef CPP_RRTM
4394	IF (CPPKEY_STRATAER) THEN
4395	!--compute stratospheric mask
4396	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4397	!--interactive strat aerosols
4398	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4399	END IF
4400	#endif
4401	!--fin STRAT AEROSOL
4402	!
4403
4404	! Calculer les parametres optiques des nuages et quelques
4405	! parametres pour diagnostiques:
4406	!
4407	IF (aerosol_couple.AND.config_inca=='aero') THEN
4408	mass_solu_aero(:,:) = ccm(:,:,1)
4409	mass_solu_aero_pi(:,:) = ccm(:,:,2)
4410	ENDIF
4411
4412	!Rajout appel a interface calcul proprietes optiques des nuages
4413	CALL call_cloud_optics_prop(klon, klev, ok_newmicro, &
4414	paprs, pplay, t_seri, radocond, picefra, cldfra, &
4415	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4416	flwp, fiwp, flwc, fiwc, ok_aie, &
4417	mass_solu_aero, mass_solu_aero_pi, &
4418	cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
4419	ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
4420	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
4421	zfice, dNovrN, ptconv, rnebcon, clwcon)
4422
4423	!
4424	!IM betaCRF
4425	!
4426	cldtaurad = cldtau
4427	cldtaupirad = cldtaupi
4428	cldemirad = cldemi
4429	cldfrarad = cldfra
4430
4431	!
4432	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4433	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4434	!
4435	! global
4436	!
4437	!IM 251017 begin
4438	! print*,'physiq betaCRF global zdtime=',zdtime
4439	!IM 251017 end
4440	DO k=1, klev
4441	DO i=1, klon
4442	IF (pplay(i,k).GE.pfree) THEN
4443	beta(i,k) = beta_pbl
4444	ELSE
4445	beta(i,k) = beta_free
4446	ENDIF
4447	IF (mskocean_beta) THEN
4448	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4449	ENDIF
4450	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4451	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4452	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4453	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4454	ENDDO
4455	ENDDO
4456	!
4457	ELSE
4458	!
4459	! regional
4460	!
4461	DO k=1, klev
4462	DO i=1,klon
4463	!
4464	IF (longitude_deg(i).ge.lon1_beta.AND. &
4465	longitude_deg(i).le.lon2_beta.AND. &
4466	latitude_deg(i).le.lat1_beta.AND. &
4467	latitude_deg(i).ge.lat2_beta) THEN
4468	IF (pplay(i,k).GE.pfree) THEN
4469	beta(i,k) = beta_pbl
4470	ELSE
4471	beta(i,k) = beta_free
4472	ENDIF
4473	IF (mskocean_beta) THEN
4474	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4475	ENDIF
4476	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4477	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4478	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4479	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4480	ENDIF
4481	!
4482	ENDDO
4483	ENDDO
4484	!
4485	ENDIF
4486
4487	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4488	IF (ok_chlorophyll) THEN
4489	print*,"-- reading chlorophyll"
4490	CALL readchlorophyll(debut)
4491	ENDIF
4492
4493	!--if ok_suntime_rrtm we use ancillay data for RSUN
4494	!--previous values are therefore overwritten
4495	!--this is needed for CMIP6 runs
4496	!--and only possible for new radiation scheme
4497	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4498	#ifdef CPP_RRTM
4499	CALL read_rsun_rrtm(debut)
4500	#endif
4501	ENDIF
4502
4503	IF (mydebug) THEN
4504	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4505	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4506	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4507	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4508	ENDIF
4509
4510	!
4511	!sonia : If Iflag_radia >=2, pertubation of some variables
4512	!input to radiation (DICE)
4513	!
4514	IF (iflag_radia .ge. 2) THEN
4515	zsav_tsol (:) = zxtsol(:)
4516	CALL perturb_radlwsw(zxtsol,iflag_radia)
4517	ENDIF
4518
4519	IF (aerosol_couple.AND.config_inca=='aero') THEN
4520	IF (CPPKEY_INCA) THEN
4521	CALL radlwsw_inca &
4522	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4523	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4524	size(wo,3), wo, &
4525	cldfrarad, cldemirad, cldtaurad, &
4526	heat,heat0,cool,cool0,albpla, &
4527	topsw,toplw,solsw,sollw, &
4528	sollwdown, &
4529	topsw0,toplw0,solsw0,sollw0, &
4530	lwdn0, lwdn, lwup0, lwup, &
4531	swdn0, swdn, swup0, swup, &
4532	ok_ade, ok_aie, &
4533	tau_aero, piz_aero, cg_aero, &
4534	topswad_aero, solswad_aero, &
4535	topswad0_aero, solswad0_aero, &
4536	topsw_aero, topsw0_aero, &
4537	solsw_aero, solsw0_aero, &
4538	cldtaupirad, &
4539	topswai_aero, solswai_aero)
4540	END IF
4541	ELSE
4542	!
4543	!IM calcul radiatif pour le cas actuel
4544	!
4545	RCO2 = RCO2_act
4546	RCH4 = RCH4_act
4547	RN2O = RN2O_act
4548	RCFC11 = RCFC11_act
4549	RCFC12 = RCFC12_act
4550	!
4551	!--interactive CO2 in ppm from carbon cycle
4552	IF (carbon_cycle_rad) RCO2=RCO2_glo
4553	!
4554	IF (prt_level .GE.10) THEN
4555	print *,' ->radlwsw, number 1 '
4556	ENDIF
4557	!
4558	! AI namelist utilise pour l appel principal de radlwsw (ecrad)
4559	namelist_ecrad_file='namelist_ecrad'
4560	!
4561	CALL radlwsw &
4562	(debut, dist, rmu0, fract, &
4563	!albedo SB >>>
4564	! paprs, pplay,zxtsol,albsol1, albsol2, &
4565	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4566	!albedo SB <<<
4567	t_seri,q_seri,wo, &
4568	cldfrarad, cldemirad, cldtaurad, &
4569	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4570	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4571	tau_aero, piz_aero, cg_aero, &
4572	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4573	! Rajoute par OB pour RRTM
4574	tau_aero_lw_rrtm, &
4575	cldtaupirad, m_allaer, &
4576	! zqsat, flwcrad, fiwcrad, &
4577	zqsat, flwc, fiwc, &
4578	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4579	namelist_ecrad_file, &
4580	heat,heat0,cool,cool0,albpla, &
4581	heat_volc,cool_volc, &
4582	topsw,toplw,solsw,solswfdiff,sollw, &
4583	sollwdown, &
4584	topsw0,toplw0,solsw0,sollw0, &
4585	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4586	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4587	topswad_aero, solswad_aero, &
4588	topswai_aero, solswai_aero, &
4589	topswad0_aero, solswad0_aero, &
4590	topsw_aero, topsw0_aero, &
4591	solsw_aero, solsw0_aero, &
4592	topswcf_aero, solswcf_aero, &
4593	!-C. Kleinschmitt for LW diagnostics
4594	toplwad_aero, sollwad_aero,&
4595	toplwai_aero, sollwai_aero, &
4596	toplwad0_aero, sollwad0_aero,&
4597	!-end
4598	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4599	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4600	cloud_cover_sw)
4601
4602	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4603	!schemes
4604	toplw = toplw + betalwoff * (toplw0 - toplw)
4605	sollw = sollw + betalwoff * (sollw0 - sollw)
4606	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4607	lwup = lwup + betalwoff * (lwup0 - lwup)
4608	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4609	sollwdown(:))
4610	cool = cool + betalwoff * (cool0 - cool)
4611
4612	IF (.NOT. using_xios) THEN
4613	!
4614	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4615	!IM des taux doit etre different du taux actuel
4616	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4617	!
4618	IF (RCO2_per.NE.RCO2_act.OR. &
4619	RCH4_per.NE.RCH4_act.OR. &
4620	RN2O_per.NE.RN2O_act.OR. &
4621	RCFC11_per.NE.RCFC11_act.OR. &
4622	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4623	ENDIF
4624	!
4625	IF (ok_4xCO2atm) THEN
4626	!
4627	RCO2 = RCO2_per
4628	RCH4 = RCH4_per
4629	RN2O = RN2O_per
4630	RCFC11 = RCFC11_per
4631	RCFC12 = RCFC12_per
4632	!
4633	IF (prt_level .GE.10) THEN
4634	print *,' ->radlwsw, number 2 '
4635	ENDIF
4636	!
4637	namelist_ecrad_file='namelist_ecrad'
4638	!
4639	CALL radlwsw &
4640	(debut, dist, rmu0, fract, &
4641	!albedo SB >>>
4642	! paprs, pplay,zxtsol,albsol1, albsol2, &
4643	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4644	!albedo SB <<<
4645	t_seri,q_seri,wo, &
4646	cldfrarad, cldemirad, cldtaurad, &
4647	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4648	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4649	tau_aero, piz_aero, cg_aero, &
4650	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4651	! Rajoute par OB pour RRTM
4652	tau_aero_lw_rrtm, &
4653	cldtaupi, m_allaer, &
4654	! zqsat, flwcrad, fiwcrad, &
4655	zqsat, flwc, fiwc, &
4656	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4657	namelist_ecrad_file, &
4658	heatp,heat0p,coolp,cool0p,albplap, &
4659	heat_volc,cool_volc, &
4660	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4661	sollwdownp, &
4662	topsw0p,toplw0p,solsw0p,sollw0p, &
4663	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4664	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4665	topswad_aerop, solswad_aerop, &
4666	topswai_aerop, solswai_aerop, &
4667	topswad0_aerop, solswad0_aerop, &
4668	topsw_aerop, topsw0_aerop, &
4669	solsw_aerop, solsw0_aerop, &
4670	topswcf_aerop, solswcf_aerop, &
4671	!-C. Kleinschmitt for LW diagnostics
4672	toplwad_aerop, sollwad_aerop,&
4673	toplwai_aerop, sollwai_aerop, &
4674	toplwad0_aerop, sollwad0_aerop,&
4675	!-end
4676	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4677	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4678	cloud_cover_sw)
4679	ENDIF !ok_4xCO2atm
4680
4681	! A.I aout 2023
4682	! Effet 3D des nuages Ecrad
4683	! a passer : nom du ficher namelist et cles ok_3Deffect
4684	! a declarer comme iflag_rrtm et a lire dans physiq.def
4685	#ifdef CPP_ECRAD
4686	IF (ok_3Deffect) then
4687	! print*,'ok_3Deffect = ',ok_3Deffect
4688	namelist_ecrad_file='namelist_ecrad_s2'
4689	CALL radlwsw &
4690	(debut, dist, rmu0, fract, &
4691	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4692	t_seri,q_seri,wo, &
4693	cldfrarad, cldemirad, cldtaurad, &
4694	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4695	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4696	tau_aero, piz_aero, cg_aero, &
4697	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4698	tau_aero_lw_rrtm, &
4699	cldtaupi, m_allaer, &
4700	zqsat, flwc, fiwc, &
4701	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4702	namelist_ecrad_file, &
4703	! A modifier
4704	heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
4705	heat_volc,cool_volc, &
4706	topsw_s2,toplw_s2,solsw_s2,solswfdiff_s2,sollw_s2, &
4707	sollwdown_s2, &
4708	topsw0_s2,toplw0_s2,solsw0_s2,sollw0_s2, &
4709	lwdnc0_s2, lwdn0_s2, lwdn_s2, lwupc0_s2, lwup0_s2, lwup_s2, &
4710	swdnc0_s2, swdn0_s2, swdn_s2, swupc0_s2, swup0_s2, swup_s2, &
4711	topswad_aero_s2, solswad_aero_s2, &
4712	topswai_aero_s2, solswai_aero_s2, &
4713	topswad0_aero_s2, solswad0_aero_s2, &
4714	topsw_aero_s2, topsw0_aero_s2, &
4715	solsw_aero_s2, solsw0_aero_s2, &
4716	topswcf_aero_s2, solswcf_aero_s2, &
4717	!-C. Kleinschmitt for LW diagnostics
4718	toplwad_aero_s2, sollwad_aero_s2,&
4719	toplwai_aero_s2, sollwai_aero_s2, &
4720	toplwad0_aero_s2, sollwad0_aero_s2,&
4721	!-end
4722	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4723	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4724	cloud_cover_sw_s2)
4725	ENDIF ! ok_3Deffect
4726	#endif
4727
4728	ENDIF ! aerosol_couple
4729	itaprad = 0
4730	!
4731	! If Iflag_radia >=2, reset pertubed variables
4732	!
4733	IF (iflag_radia .ge. 2) THEN
4734	zxtsol(:) = zsav_tsol (:)
4735	ENDIF
4736	ENDIF ! MOD(itaprad,radpas)
4737	itaprad = itaprad + 1
4738
4739	IF (iflag_radia.eq.0) THEN
4740	IF (prt_level.ge.9) THEN
4741	PRINT *,'--------------------------------------------------'
4742	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4743	PRINT *,'>>>> heat et cool mis a zero '
4744	PRINT *,'--------------------------------------------------'
4745	ENDIF
4746	heat=0.
4747	cool=0.
4748	sollw=0. ! MPL 01032011
4749	solsw=0.
4750	radsol=0.
4751	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4752	swup0=0.
4753	lwup=0.
4754	lwup0=0.
4755	lwdn=0.
4756	lwdn0=0.
4757	ENDIF
4758
4759	!
4760	! Calculer radsol a l'exterieur de radlwsw
4761	! pour prendre en compte le cycle diurne
4762	! recode par Olivier Boucher en sept 2015
4763	!
4764	radsol=solsw*swradcorr+sollw
4765
4766	IF (ok_4xCO2atm) THEN
4767	radsolp=solswp*swradcorr+sollwp
4768	ENDIF
4769
4770	!
4771	! Ajouter la tendance des rayonnements (tous les pas)
4772	! avec une correction pour le cycle diurne dans le SW
4773	!
4774
4775	DO k=1, klev
4776	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4777	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4778	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4779	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4780	ENDDO
4781
4782	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4783	CALL prt_enerbil('SW',itap)
4784	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4785	CALL prt_enerbil('LW',itap)
4786
4787	!
4788	IF (mydebug) THEN
4789	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4790	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4791	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4792	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4793	ENDIF
4794
4795	! Calculer l'hydrologie de la surface
4796	!
4797	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4798	! . agesno, ftsol,fqsurf,fsnow, ruis)
4799	!
4800
4801	!
4802	! Calculer le bilan du sol et la derive de temperature (couplage)
4803	!
4804	DO i = 1, klon
4805	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4806	! a la demande de JLD
4807	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4808	ENDDO
4809	!
4810	!moddeblott(jan95)
4811	! Appeler le programme de parametrisation de l'orographie
4812	! a l'echelle sous-maille:
4813	!
4814	IF (prt_level .GE.10) THEN
4815	print *,' call orography ? ', ok_orodr
4816	ENDIF
4817	!
4818	IF (ok_orodr) THEN
4819	!
4820	! selection des points pour lesquels le shema est actif:
4821	igwd=0
4822	DO i=1,klon
4823	itest(i)=0
4824	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4825	!zrel_oro: relative mountain height wrt relief explained by mean slope
4826	! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
4827	! such as ice sheets (work by V. Wiener)
4828	! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
4829	! earn computation time but they are not physical.
4830	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
4831	itest(i)=1
4832	igwd=igwd+1
4833	idx(igwd)=i
4834	ENDIF
4835	ENDDO
4836	! igwdim=MAX(1,igwd)
4837	!
4838	IF (ok_strato) THEN
4839
4840	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4841	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4842	igwd,idx,itest, &
4843	t_seri, u_seri, v_seri, &
4844	zulow, zvlow, zustrdr, zvstrdr, &
4845	d_t_oro, d_u_oro, d_v_oro)
4846
4847	ELSE
4848	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4849	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4850	igwd,idx,itest, &
4851	t_seri, u_seri, v_seri, &
4852	zulow, zvlow, zustrdr, zvstrdr, &
4853	d_t_oro, d_u_oro, d_v_oro)
4854	ENDIF
4855	!
4856	! ajout des tendances
4857	!-----------------------------------------------------------------------
4858	! ajout des tendances de la trainee de l'orographie
4859	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
4860	abortphy,flag_inhib_tend,itap,0)
4861	CALL prt_enerbil('oro',itap)
4862	!----------------------------------------------------------------------
4863	!
4864	ENDIF ! fin de test sur ok_orodr
4865	!
4866	IF (mydebug) THEN
4867	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4868	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4869	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4870	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4871	ENDIF
4872
4873	IF (ok_orolf) THEN
4874	!
4875	! selection des points pour lesquels le shema est actif:
4876	igwd=0
4877	DO i=1,klon
4878	itest(i)=0
4879	!zrel_oro: relative mountain height wrt relief explained by mean slope
4880	! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
4881	! such as ice sheets (work by V. Wiener)
4882	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4883	IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
4884	itest(i)=1
4885	igwd=igwd+1
4886	idx(igwd)=i
4887	ENDIF
4888	ENDDO
4889	! igwdim=MAX(1,igwd)
4890	!
4891	IF (ok_strato) THEN
4892
4893	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4894	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4895	igwd,idx,itest, &
4896	t_seri, u_seri, v_seri, &
4897	zulow, zvlow, zustrli, zvstrli, &
4898	d_t_lif, d_u_lif, d_v_lif )
4899
4900	ELSE
4901	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4902	latitude_deg,zmea,zstd,zpic, &
4903	itest, &
4904	t_seri, u_seri, v_seri, &
4905	zulow, zvlow, zustrli, zvstrli, &
4906	d_t_lif, d_u_lif, d_v_lif)
4907	ENDIF
4908
4909	! ajout des tendances de la portance de l'orographie
4910	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
4911	'lif', abortphy,flag_inhib_tend,itap,0)
4912	CALL prt_enerbil('lif',itap)
4913	ENDIF ! fin de test sur ok_orolf
4914
4915	IF (ok_hines) then
4916	! HINES GWD PARAMETRIZATION
4917	east_gwstress=0.
4918	west_gwstress=0.
4919	du_gwd_hines=0.
4920	dv_gwd_hines=0.
4921	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4922	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4923	du_gwd_hines, dv_gwd_hines)
4924	zustr_gwd_hines=0.
4925	zvstr_gwd_hines=0.
4926	DO k = 1, klev
4927	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4928	* (paprs(:, k)-paprs(:, k+1))/rg
4929	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4930	* (paprs(:, k)-paprs(:, k+1))/rg
4931	ENDDO
4932
4933	d_t_hin(:, :)=0.
4934	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4935	dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4936	CALL prt_enerbil('hin',itap)
4937	ENDIF
4938
4939	IF (.not. ok_hines .and. ok_gwd_rando) then
4940	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4941	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4942	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4943	dv_gwd_front, east_gwstress, west_gwstress)
4944	zustr_gwd_front=0.
4945	zvstr_gwd_front=0.
4946	DO k = 1, klev
4947	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4948	* (paprs(:, k)-paprs(:, k+1))/rg
4949	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4950	* (paprs(:, k)-paprs(:, k+1))/rg
4951	ENDDO
4952
4953	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
4954	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4955	CALL prt_enerbil('front_gwd_rando',itap)
4956	ENDIF
4957
4958	IF (ok_gwd_rando) THEN
4959	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4960	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4961	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4962	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
4963	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4964	CALL prt_enerbil('flott_gwd_rando',itap)
4965	zustr_gwd_rando=0.
4966	zvstr_gwd_rando=0.
4967	DO k = 1, klev
4968	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4969	* (paprs(:, k)-paprs(:, k+1))/rg
4970	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4971	* (paprs(:, k)-paprs(:, k+1))/rg
4972	ENDDO
4973	ENDIF
4974
4975	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4976
4977	IF (mydebug) THEN
4978	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4979	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4980	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4981	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4982	ENDIF
4983
4984	DO i = 1, klon
4985	zustrph(i)=0.
4986	zvstrph(i)=0.
4987	ENDDO
4988	DO k = 1, klev
4989	DO i = 1, klon
4990	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4991	(paprs(i,k)-paprs(i,k+1))/rg
4992	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4993	(paprs(i,k)-paprs(i,k+1))/rg
4994	ENDDO
4995	ENDDO
4996	!
4997	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4998	!
4999	IF (is_sequential .and. ok_orodr) THEN
5000	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
5001	ra,rg,romega, &
5002	latitude_deg,longitude_deg,pphis, &
5003	zustrdr,zustrli,zustrph, &
5004	zvstrdr,zvstrli,zvstrph, &
5005	paprs,u,v, &
5006	aam, torsfc)
5007	ENDIF
5008	!IM cf. FLott END
5009	!DC Calcul de la tendance due au methane
5010	IF (ok_qch4) THEN
5011	! d_q_ch4: H2O source from CH4 in MMR/s (mass mixing ratio/s or kg H2O/kg air/s)
5012	IF (CPPKEY_STRATAER) THEN
5013
5014	CALL stratH2O_methox(debut,paprs,d_q_ch4)
5015	ELSE
5016	! ECMWF routine METHOX
5017	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
5018	END IF
5019	! add humidity tendency due to methane
5020	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
5021	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
5022	'q_ch4', abortphy,flag_inhib_tend,itap,0)
5023	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep ! update with H2O conserv done in add_phys_tend
5024	ENDIF
5025	!
5026	!
5027	IF (CPPKEY_STRATAER) THEN
5028	IF (ok_qemiss) THEN
5029	flh2o=1
5030	IF(flag_verbose_strataer) THEN
5031	print *,'IN physiq_mod: ok_qemiss =yes (',ok_qemiss,'), flh2o=',flh2o
5032	print *,'IN physiq_mod: flag_emit=',flag_emit,', nErupt=',nErupt
5033	print *,'IN physiq_mod: nAerErupt=',nAerErupt
5034	ENDIF
5035
5036	SELECT CASE(flag_emit)
5037	CASE(1) ! emission volc H2O in LMDZ
5038	DO ieru=1, nErupt
5039	IF (year_cur==year_emit_vol(ieru).AND.&
5040	mth_cur==mth_emit_vol(ieru).AND.&
5041	day_cur>=day_emit_vol(ieru).AND.&
5042	day_cur<(day_emit_vol(ieru)+injdur)) THEN
5043
5044	IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
5045	! initialisation of q tendency emission
5046	d_q_emiss(:,:)=0.
5047	! daily injection mass emission - NL
5048	m_H2O_emiss_vol_daily = m_H2O_emiss_vol(ieru)/(REAL(injdur)&
5049	*REAL(ponde_lonlat_vol(ieru)))
5050	!
5051	CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
5052	pplay,paprs,tr_seri,&
5053	m_H2O_emiss_vol_daily,&
5054	xlat_min_vol(ieru),xlat_max_vol(ieru),&
5055	xlon_min_vol(ieru),xlon_max_vol(ieru),&
5056	altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
5057	nAerErupt+1,0)
5058
5059	IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
5060	minval(d_q_emiss),maxval(d_q_emiss)
5061
5062	CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
5063	'q_emiss',abortphy,flag_inhib_tend,itap,0)
5064	IF (abortphy==1) Print*,'ERROR ABORT TEND EMISS'
5065	ENDIF
5066	ENDDO
5067	flh2o=0
5068	END SELECT ! emission scenario (flag_emit)
5069	ENDIF
5070	END IF
5071
5072	!===============================================================
5073	! Additional tendency of TKE due to orography
5074	!===============================================================
5075	!
5076	! Inititialization
5077	!------------------
5078
5079	addtkeoro=0
5080	CALL getin_p('addtkeoro',addtkeoro)
5081
5082	IF (prt_level.ge.5) &
5083	print*,'addtkeoro', addtkeoro
5084
5085	alphatkeoro=1.
5086	CALL getin_p('alphatkeoro',alphatkeoro)
5087	alphatkeoro=min(max(0.,alphatkeoro),1.)
5088
5089	smallscales_tkeoro=.FALSE.
5090	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
5091
5092
5093	dtadd(:,:)=0.
5094	duadd(:,:)=0.
5095	dvadd(:,:)=0.
5096
5097	! Choices for addtkeoro:
5098	! ** 0 no TKE tendency from orography
5099	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
5100	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
5101	!
5102
5103	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
5104	! -------------------------------------------
5105
5106
5107	! selection des points pour lesquels le schema est actif:
5108
5109
5110	IF (addtkeoro .EQ. 1 ) THEN
5111
5112	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
5113	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
5114
5115	ELSE IF (addtkeoro .EQ. 2) THEN
5116
5117	IF (smallscales_tkeoro) THEN
5118	igwd=0
5119	DO i=1,klon
5120	itest(i)=0
5121	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
5122	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
5123	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
5124	IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
5125	itest(i)=1
5126	igwd=igwd+1
5127	idx(igwd)=i
5128	ENDIF
5129	ENDDO
5130
5131	ELSE
5132
5133	igwd=0
5134	DO i=1,klon
5135	itest(i)=0
5136	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
5137	itest(i)=1
5138	igwd=igwd+1
5139	idx(igwd)=i
5140	ENDIF
5141	ENDDO
5142
5143	ENDIF
5144
5145	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
5146	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
5147	igwd,idx,itest, &
5148	t_seri, u_seri, v_seri, &
5149	zulow, zvlow, zustrdr, zvstrdr, &
5150	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
5151
5152	zustrdr(:)=0.
5153	zvstrdr(:)=0.
5154	zulow(:)=0.
5155	zvlow(:)=0.
5156
5157	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
5158	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
5159	ENDIF
5160
5161
5162	! TKE update from subgrid temperature and wind tendencies
5163	!----------------------------------------------------------
5164	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5165
5166
5167	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
5168	!
5169	! Prevent pbl_tke_w from becoming negative
5170	wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
5171	!
5172
5173	ENDIF
5174	! -----
5175	!===============================================================
5176
5177
5178	!====================================================================
5179	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
5180	!====================================================================
5181	! Abderrahmane 24.08.09
5182
5183	IF (ok_cosp) THEN
5184	! adeclarer
5185	IF (CPPKEY_COSP) THEN
5186	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5187
5188	IF (prt_level .GE.10) THEN
5189	print*,'freq_cosp',freq_cosp
5190	ENDIF
5191	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5192	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
5193	! s ref_liq,ref_ice
5194	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
5195	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5196	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5197	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5198	JrNt,ref_liq,ref_ice, &
5199	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5200	zu10m,zv10m,pphis, &
5201	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5202	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5203	prfl(:,1:klev),psfl(:,1:klev), &
5204	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5205	mr_ozone,cldtau, cldemi)
5206
5207	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
5208	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
5209	! M clMISR,
5210	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
5211	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
5212
5213	ENDIF
5214	END IF
5215
5216	IF (CPPKEY_COSP2) THEN
5217	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5218
5219	IF (prt_level .GE.10) THEN
5220	print*,'freq_cosp',freq_cosp
5221	ENDIF
5222	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5223	print*,'Dans physiq.F avant appel '
5224	! s ref_liq,ref_ice
5225	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
5226	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5227	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5228	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5229	JrNt,ref_liq,ref_ice, &
5230	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5231	zu10m,zv10m,pphis, &
5232	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5233	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5234	prfl(:,1:klev),psfl(:,1:klev), &
5235	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5236	mr_ozone,cldtau, cldemi)
5237	ENDIF
5238	END IF
5239
5240	IF (CPPKEY_COSPV2) THEN
5241	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5242	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5243
5244	IF (prt_level .GE.10) THEN
5245	print*,'freq_cosp',freq_cosp
5246	ENDIF
5247	DO k = 1, klev
5248	DO i = 1, klon
5249	phicosp(i,k) = pphi(i,k) + pphis(i)
5250	ENDDO
5251	ENDDO
5252	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5253	print*,'Dans physiq.F avant appel '
5254	! s ref_liq,ref_ice
5255	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
5256	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5257	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5258	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5259	JrNt,ref_liq,ref_ice, &
5260	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5261	zu10m,zv10m,pphis, &
5262	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5263	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5264	prfl(:,1:klev),psfl(:,1:klev), &
5265	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5266	mr_ozone,cldtau, cldemi)
5267	ENDIF
5268	END IF
5269
5270	ENDIF !ok_cosp
5271
5272
5273	! Marine
5274
5275	IF (ok_airs) then
5276
5277	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
5278	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
5279	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
5280	& map_prop_hc,map_prop_hist,&
5281	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
5282	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
5283	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
5284	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
5285	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
5286	& map_ntot,map_hc,map_hist,&
5287	& map_Cb,map_ThCi,map_Anv,&
5288	& alt_tropo )
5289	ENDIF
5290
5291	ENDIF ! ok_airs
5292
5293
5294	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5295	!AA
5296	!AA Installation de l'interface online-offline pour traceurs
5297	!AA
5298	!====================================================================
5299	! Calcul des tendances traceurs
5300	!====================================================================
5301	!
5302
5303	IF (type_trac == 'repr') THEN
5304	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
5305	!MM dans Reprobus
5306	sh_in(:,:) = q_seri(:,:)
5307	IF (CPPKEY_REPROBUS) THEN
5308	d_q_rep(:,:) = 0.
5309	d_ql_rep(:,:) = 0.
5310	d_qi_rep(:,:) = 0.
5311	END IF
5312	ELSE
5313	sh_in(:,:) = qx(:,:,ivap)
5314	IF (nqo >= 3) THEN
5315	ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
5316	ELSE
5317	ch_in(:,:) = qx(:,:,iliq)
5318	ENDIF
5319	ENDIF
5320
5321	! Merge wdtrainA and wdtrainS in the total source of precipitation due to
5322	! adiabatic updraughts.
5323	!
5324	wdtrainAS(:,:) = wdtrainA(:,:) + wdtrainS(:,:)
5325
5326	IF (CPPKEY_DUST) THEN
5327	! Avec SPLA, iflag_phytrac est forcé =1
5328	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
5329	pdtphys,ftsol, & ! I
5330	t,q_seri,paprs,pplay,RHcl, & ! I
5331	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
5332	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
5333	u_seri, v_seri, latitude_deg, longitude_deg, &
5334	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
5335	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
5336	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
5337	ev,wdtrainAS, wdtrainM,wght_cvfd, & ! I
5338	fm_therm, entr_therm, rneb, & ! I
5339	beta_prec_fisrt,beta_prec, & !I
5340	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
5341	d_tr_dyn,tr_seri)
5342
5343	ELSE
5344	IF (iflag_phytrac == 1 ) THEN
5345	CALL phytrac ( &
5346	itap, days_elapsed+1, jH_cur, debut, &
5347	lafin, phys_tstep, u, v, t, &
5348	paprs, pplay, pmfu, pmfd, &
5349	pen_u, pde_u, pen_d, pde_d, &
5350	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
5351	u1, v1, ftsol, pctsrf, &
5352	zustar, zu10m, zv10m, &
5353	wstar(:,is_ave), ale_bl, ale_wake, &
5354	latitude_deg, longitude_deg, &
5355	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
5356	presnivs, pphis, pphi, albsol1, &
5357	sh_in, ch_in, rhcl, cldfra, rneb, &
5358	diafra, radocond, itop_con, ibas_con, &
5359	pmflxr, pmflxs, prfl, psfl, &
5360	da, phi, mp, upwd, &
5361	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
5362	wdtrainAS, wdtrainM, sigd, clw,elij, & !<<RomP
5363	ev, ep, epmlmMm, eplaMm, & !<<RomP
5364	dnwd, aerosol_couple, flxmass_w, &
5365	tau_aero, piz_aero, cg_aero, ccm, &
5366	rfname, &
5367	d_tr_dyn, & !<<RomP
5368	tr_seri, init_source)
5369	IF (CPPKEY_REPROBUS) THEN
5370
5371
5372	print*,'avt add phys rep',abortphy
5373
5374	CALL add_phys_tend &
5375	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
5376	'rep',abortphy,flag_inhib_tend,itap,0)
5377	IF (abortphy==1) Print*,'ERROR ABORT REP'
5378
5379	print*,'apr add phys rep',abortphy
5380
5381	END IF
5382	ENDIF ! (iflag_phytrac=1)
5383
5384	END IF
5385	!ENDIF ! (iflag_phytrac=1)
5386
5387	IF (offline) THEN
5388
5389	IF (prt_level.ge.9) &
5390	print*,'Attention on met a 0 les thermiques pour phystoke'
5391	CALL phystokenc ( &
5392	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5393	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5394	fm_therm,entr_therm, &
5395	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5396	frac_impa, frac_nucl, &
5397	pphis,cell_area,phys_tstep,itap, &
5398	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
5399
5400
5401	ENDIF
5402
5403	!
5404	! Calculer le transport de l'eau et de l'energie (diagnostique)
5405	!
5406	CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5407	ue, ve, uq, vq, uwat, vwat)
5408	!
5409	!IM global posePB BEG
5410	IF(1.EQ.0) THEN
5411	!
5412	CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
5413	ve_lay, vq_lay, ue_lay, uq_lay)
5414	!
5415	ENDIF !(1.EQ.0) THEN
5416	!IM global posePB END
5417	!
5418	! Accumuler les variables a stocker dans les fichiers histoire:
5419	!
5420
5421	!================================================================
5422	! Conversion of kinetic and potential energy into heat, for
5423	! parameterisation of subgrid-scale motions
5424	!================================================================
5425
5426	d_t_ec(:,:)=0.
5427	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5428	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
5429	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
5430	zmasse,exner,d_t_ec)
5431	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
5432
5433	!==================================================================
5434	!--OB water mass fixer for the physics
5435	!--water profiles are corrected to force mass conservation of water
5436	!--currently flag is turned off
5437	!==================================================================
5438	IF (ok_water_mass_fixer) THEN
5439	qql2(:)=0.0
5440	DO k = 1, klev
5441	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
5442	IF (nqo >= 3) THEN
5443	qql2(:)=qql2(:)+qs_seri(:,k)*zmasse(:,k)
5444	ENDIF
5445	IF (ok_bs) THEN
5446	qql2(:)=qql2(:)+qbs_seri(:,k)*zmasse(:,k)
5447	ENDIF
5448	ENDDO
5449
5450	IF (CPPKEY_STRATAER) THEN
5451	IF (ok_qemiss) THEN
5452	DO k = 1, klev
5453	qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
5454	ENDDO
5455	ENDIF
5456	END IF
5457	IF (ok_qch4) THEN
5458	DO k = 1, klev
5459	qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
5460	ENDDO
5461	ENDIF
5462
5463	DO i = 1, klon
5464	!--compute ratio of what q+ql should be with conservation to what it is
5465	IF (ok_bs) THEN
5466	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i)-bs_fall(i))*pdtphys)/qql2(i)
5467	ELSE
5468	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5469	ENDIF
5470	DO k = 1, klev
5471	q_seri(i,k) =q_seri(i,k)*corrqql
5472	ql_seri(i,k)=ql_seri(i,k)*corrqql
5473	IF (nqo >= 3) THEN
5474	qs_seri(i,k)=qs_seri(i,k)*corrqql
5475	ENDIF
5476	IF (ok_bs) THEN
5477	qbs_seri(i,k)=qbs_seri(i,k)*corrqql
5478	ENDIF
5479	ENDDO
5480	ENDDO
5481	ENDIF
5482	!--fin mass fixer
5483
5484	!cc prw = eau precipitable
5485	! prlw = colonne eau liquide
5486	! prlw = colonne eau solide
5487	! prbsw = colonne neige soufflee
5488	! water_budget = non-conservation residual from the LMDZ physics
5489	! (should be equal to machine precision if mass fixer is activated)
5490	prw(:) = 0.
5491	prlw(:) = 0.
5492	prsw(:) = 0.
5493	prbsw(:) = 0.
5494	water_budget(:) = 0.0
5495	DO k = 1, klev
5496	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
5497	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5498	water_budget(:) = water_budget(:) + (q_seri(:,k)-qx(:,k,ivap)+ql_seri(:,k)-qx(:,k,iliq))*zmasse(:,k)
5499	IF (nqo >= 3) THEN
5500	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
5501	water_budget(:) = water_budget(:) + (qs_seri(:,k)-qx(:,k,isol))*zmasse(:,k)
5502	ENDIF
5503	IF (nqo >= 4 .AND. ok_bs) THEN
5504	prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
5505	water_budget(:) = water_budget(:) + (qbs_seri(:,k)-qx(:,k,ibs))*zmasse(:,k)
5506	ENDIF
5507	ENDDO
5508	water_budget(:)=water_budget(:)+(rain_fall(:)+snow_fall(:)-evap(:))*pdtphys
5509	IF (ok_bs) THEN
5510	water_budget(:)=water_budget(:)+bs_fall(:)*pdtphys
5511	ENDIF
5512
5513	!=======================================================================
5514	! SORTIES
5515	!=======================================================================
5516	!
5517	!IM initialisation + calculs divers diag AMIP2
5518	CALL calcul_divers(itap, itapm1, un_jour)
5519	!
5520	!IM Interpolation sur les niveaux de pression du NMC
5521	! -------------------------------------------------
5522	!
5523	include "calcul_STDlev.h"
5524	!
5525	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5526	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5527	!
5528	!
5529	IF (ANY(type_trac == ['inca','inco'])) THEN
5530	IF (CPPKEY_INCA) THEN
5531	CALL VTe(VTphysiq)
5532	CALL VTb(VTinca)
5533
5534	CALL chemhook_end ( &
5535	phys_tstep, &
5536	pplay, &
5537	t_seri, &
5538	tr_seri(:,:,1+nqCO2:nbtr), &
5539	nbtr, &
5540	paprs, &
5541	q_seri, &
5542	cell_area, &
5543	pphi, &
5544	pphis, &
5545	zx_rh, &
5546	aps, bps, ap, bp, lafin)
5547
5548	CALL VTe(VTinca)
5549	CALL VTb(VTphysiq)
5550	END IF
5551	ENDIF
5552
5553	IF (type_trac == 'repr') THEN
5554	IF (CPPKEY_REPROBUS) THEN
5555	CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
5556	END IF
5557	ENDIF
5558
5559	!
5560	! Convertir les incrementations en tendances
5561	!
5562	IF (prt_level .GE.10) THEN
5563	print *,'Convertir les incrementations en tendances '
5564	ENDIF
5565	!
5566	IF (mydebug) THEN
5567	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5568	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5569	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5570	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5571	ENDIF
5572
5573	DO k = 1, klev
5574	DO i = 1, klon
5575	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5576	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5577	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5578	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5579	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
5580	!CR: on ajoute le contenu en glace
5581	IF (nqo >= 3) THEN
5582	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5583	ENDIF
5584	!--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
5585	IF (nqo.ge.5 .and. ok_ice_supersat) THEN
5586	d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
5587	d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
5588	ENDIF
5589
5590	IF (nqo.ge.4 .and. ok_bs) THEN
5591	d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
5592	ENDIF
5593
5594	ENDDO
5595	ENDDO
5596	!
5597	! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
5598	itr = 0
5599	DO iq = 1, nqtot
5600	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
5601	itr = itr+1
5602	DO k = 1, klev
5603	DO i = 1, klon
5604	d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
5605	ENDDO
5606	ENDDO
5607	ENDDO
5608	!
5609	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5610	!IM global posePB include "write_bilKP_ins.h"
5611	!IM global posePB include "write_bilKP_ave.h"
5612	!
5613	! Sauvegarder les valeurs de t et q a la fin de la physique:
5614	!
5615	u_ancien(:,:) = u_seri(:,:)
5616	v_ancien(:,:) = v_seri(:,:)
5617	t_ancien(:,:) = t_seri(:,:)
5618	q_ancien(:,:) = q_seri(:,:)
5619	ql_ancien(:,:) = ql_seri(:,:)
5620	qs_ancien(:,:) = qs_seri(:,:)
5621	qbs_ancien(:,:)= qbs_seri(:,:)
5622	cf_ancien(:,:) = cf_seri(:,:)
5623	rvc_ancien(:,:)= rvc_seri(:,:)
5624	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5625	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5626	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5627	CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
5628	! !! RomP >>>
5629	IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
5630	! !! RomP <<<
5631	!==========================================================================
5632	! Sorties des tendances pour un point particulier
5633	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5634	! pour le debug
5635	! La valeur de igout est attribuee plus haut dans le programme
5636	!==========================================================================
5637
5638	IF (prt_level.ge.1) THEN
5639	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5640	write(lunout,*) &
5641	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5642	write(lunout,*) &
5643	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5644	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5645	pctsrf(igout,is_sic)
5646	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5647	DO k=1,klev
5648	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5649	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5650	d_t_eva(igout,k)
5651	ENDDO
5652	write(lunout,*) 'cool,heat'
5653	DO k=1,klev
5654	write(lunout,*) cool(igout,k),heat(igout,k)
5655	ENDDO
5656
5657	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5658	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5659	!jyg! do k=1,klev
5660	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5661	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5662	!jyg! enddo
5663	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5664	DO k=1,klev
5665	write(lunout,*) d_t_vdf(igout,k), &
5666	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5667	ENDDO
5668	!>jyg
5669
5670	write(lunout,*) 'd_ps ',d_ps(igout)
5671	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5672	DO k=1,klev
5673	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5674	d_qx(igout,k,1),d_qx(igout,k,2)
5675	ENDDO
5676	ENDIF
5677
5678	!============================================================
5679	! Calcul de la temperature potentielle
5680	!============================================================
5681	DO k = 1, klev
5682	DO i = 1, klon
5683	!JYG/IM theta en debut du pas de temps
5684	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5685	!JYG/IM theta en fin de pas de temps de physique
5686	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5687	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5688	! MPL 20130625
5689	! fth_fonctions.F90 et parkind1.F90
5690	! sinon thetal=theta
5691	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5692	! : ql_seri(i,k))
5693	thetal(i,k)=theta(i,k)
5694	ENDDO
5695	ENDDO
5696	!
5697
5698	! 22.03.04 BEG
5699	!=============================================================
5700	! Ecriture des sorties
5701	!=============================================================
5702
5703	! Recupere des varibles calcule dans differents modules
5704	! pour ecriture dans histxxx.nc
5705
5706	! Get some variables from module fonte_neige_mod
5707	CALL fonte_neige_get_vars(pctsrf, &
5708	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5709
5710
5711	!=============================================================
5712	! Separation entre thermiques et non thermiques dans les sorties
5713	! de fisrtilp
5714	!=============================================================
5715
5716	IF (iflag_thermals>=1) THEN
5717	d_t_lscth=0.
5718	d_t_lscst=0.
5719	d_q_lscth=0.
5720	d_q_lscst=0.
5721	DO k=1,klev
5722	DO i=1,klon
5723	IF (ptconvth(i,k)) THEN
5724	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5725	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5726	ELSE
5727	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5728	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5729	ENDIF
5730	ENDDO
5731	ENDDO
5732
5733	DO i=1,klon
5734	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5735	plul_th(i)=prfl(i,1)+psfl(i,1)
5736	ENDDO
5737	ENDIF
5738
5739	!On effectue les sorties:
5740
5741	IF (CPPKEY_DUST) THEN
5742	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5743	pplay, lmax_th, aerosol_couple, &
5744	ok_ade, ok_aie, ivap, ok_sync, &
5745	ptconv, read_climoz, clevSTD, &
5746	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5747	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5748	ELSE
5749	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5750	pplay, lmax_th, aerosol_couple, &
5751	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs, &
5752	ok_sync, ptconv, read_climoz, clevSTD, &
5753	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5754	flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
5755	END IF
5756
5757	#ifndef CPP_XIOS
5758	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5759	#endif
5760
5761	! Petit appelle de sorties pour accompagner le travail sur phyex
5762	if ( iflag_physiq == 1 ) then
5763	call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.t,0.t,0.*t,pbl_tke,theta)
5764	endif
5765
5766	!====================================================================
5767	! Arret du modele apres hgardfou en cas de detection d'un
5768	! plantage par hgardfou
5769	!====================================================================
5770
5771	IF (abortphy==1) THEN
5772	abort_message ='Plantage hgardfou'
5773	CALL abort_physic (modname,abort_message,1)
5774	ENDIF
5775
5776	! 22.03.04 END
5777	!
5778	!====================================================================
5779	! Si c'est la fin, il faut conserver l'etat de redemarrage
5780	!====================================================================
5781	!
5782
5783	! Disabling calls to the prt_alerte function
5784	alert_first_call = .FALSE.
5785
5786
5787	IF (lafin) THEN
5788	itau_phy = itau_phy + itap
5789	CALL phyredem ("restartphy.nc")
5790	! open(97,form="unformatted",file="finbin")
5791	! write(97) u_seri,v_seri,t_seri,q_seri
5792	! close(97)
5793
5794	IF (is_omp_master) THEN
5795
5796	IF (read_climoz >= 1) THEN
5797	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5798	DEALLOCATE(press_edg_climoz)
5799	DEALLOCATE(press_cen_climoz)
5800	ENDIF
5801
5802	ENDIF
5803
5804	IF (using_xios) THEN
5805
5806	IF (CPPKEY_INCA) THEN
5807	IF (type_trac == 'inca') THEN
5808	IF (is_omp_master .AND. grid_type==unstructured) THEN
5809	CALL finalize_inca
5810	ENDIF
5811	ENDIF
5812	END IF
5813
5814	! close xios physiq context (call LMDZ)
5815	IF (is_omp_master) CALL xios_context_finalize
5816	ENDIF
5817
5818	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5819
5820	ENDIF
5821
5822	! first=.false.
5823
5824	END SUBROUTINE physiq
5825
5826	END MODULE physiq_mod

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

Download in other formats:

Original Format