Context Navigation

physiq_mod.F90 @ 5287

Last change on this file since 5287 was 5285, checked in by abarral, 9 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h 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: 214.0 KB

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

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

Download in other formats:

Original Format