Context Navigation

physiq_mod.F90 @ 5110

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

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

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

Download in other formats:

Original Format