Context Navigation

physiq_mod.F90 @ 4707

Last change on this file since 4707 was 4707, checked in by Laurent Fairhead, 12 months ago
Continuing on from the morning's poihl workshop: getting rid of the includes in the routine
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.4 KB

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

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

Download in other formats:

Original Format