Context Navigation

physiq_mod.F90 @ 4738

Last change on this file since 4738 was 4727, checked in by idelkadi, 10 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
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: 210.9 KB

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

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

Download in other formats:

Original Format