Context Navigation

physiq_mod.F90 @ 5450

Last change on this file since 5450 was 5445, checked in by fhourdin, 5 days ago
Adding wind subrid scale distribution
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: 215.6 KB

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

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

Download in other formats:

Original Format