Context Navigation

physiq_mod.F90 @ 5151

Last change on this file since 5151 was 5151, checked in by abarral, 3 months ago
Move CPP_COSP* in lmdz_cppkeys_wrapper.F90 Change calcul_divers.h into lmdz_calcul_divers.f90
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: 211.4 KB

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

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

Download in other formats:

Original Format