Context Navigation

physiq_mod.F90 @ 5327

Last change on this file since 5327 was 5326, checked in by abarral, 22 hours ago
Turn condsurfs_new.f90 into module Add own handle_err subroutine to condsurfs_new_mod.f90 phytracr_spl_mod.F90, as they were using incorrectly the one from grilles_gcm_netcdf_sub
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 213.7 KB

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

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

Download in other formats:

Original Format