Context Navigation

physiq_mod.F90 @ 5322

Last change on this file since 5322 was 5315, checked in by abarral, 10 days ago
Turn regdim.h into modules
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: 213.5 KB

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

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

Download in other formats:

Original Format