Context Navigation

physiq_mod.F90 @ 5435

Last change on this file since 5435 was 5203, checked in by Laurent Fairhead, 3 months ago
Merge with trunk revision 5202 before reintegration to trunk
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: 214.2 KB

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

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

Download in other formats:

Original Format