Context Navigation

physiq_mod.F90 @ 5486

Last change on this file since 5486 was 5486, checked in by evignon, 2 days ago
inclusion d'un diagnostique de la sublimation de la glace sur les landice pour la conservation de l'eau
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 215.4 KB

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

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

Download in other formats:

Original Format