Context Navigation

physiq_mod.F90 @ 5185

Last change on this file since 5185 was 4951, checked in by aborella, 7 months ago
New version of condensation and ice supersaturation in LSCP. Multiple changes troughout the code (in particular, two new water phase tracers).
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: 216.0 KB

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

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

Download in other formats:

Original Format