Context Navigation

physiq_mod.F90 @ 5171

Last change on this file since 5171 was 5169, checked in by idelkadi, 10 months ago

Update of the COSP part in LMDZ:

initialization of input fields for the first call to Cospv2

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: 212.5 KB

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

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

Download in other formats:

Original Format