Context Navigation

physiq_mod.F90 @ 4886

Last change on this file since 4886 was 4886, checked in by oboucher, 7 weeks ago
Correcting a couple of bugs in the mass fixer for water vapour in the physics (the mass_fixer is switched to FALSE by default)
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: 213.0 KB

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

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

Download in other formats:

Original Format