Context Navigation

physiq_mod.F90 @ 5143

Last change on this file since 5143 was 5143, checked in by abarral, 3 months ago
Put YOEGWD.h, FCTTRE.h into modules
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.0 KB

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

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

Download in other formats:

Original Format