Context Navigation

physiq_mod.F90 @ 5006

Last change on this file since 5006 was 4976, checked in by idelkadi, 7 months ago
Error correction introduced in previous version
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.7 KB

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

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

Download in other formats:

Original Format