Context Navigation

physiq_mod.F90 @ 5477

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

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

Download in other formats:

Original Format