Context Navigation

physiq_mod.F90 @ 5221

Last change on this file since 5221 was 5221, checked in by abarral, 4 weeks ago
Merge r5196, r5199 Nb: skipping r5197 r5198 as they revert one-another
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 210.8 KB

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

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

Download in other formats:

Original Format