Context Navigation

physiq_mod.F90 @ 4790

Last change on this file since 4790 was 4790, checked in by idelkadi, 4 months ago
Continued work on implementing the Ecrad code in LMDZ (Integration of aerosols):
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.1 KB

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

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

Download in other formats:

Original Format