Context Navigation

physiq_mod.F90 @ 5112

Last change on this file since 5112 was 5112, checked in by abarral, 4 months ago
Rename modules in phy_common from _mod > lmdz_
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 210.9 KB

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

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

Download in other formats:

Original Format