Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq_mod.F90

Last change on this file was 4887, checked in by oboucher, 5 weeks ago
Correction of another small bug on the water mass fixer (not activated) Introduction of a water_budget diagnostic that quantifies non-conservation from the LMDZ physics
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.4 KB

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

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

Download in other formats:

Original Format