Context Navigation

physiq_mod.F90 @ 4744

Last change on this file since 4744 was 4744, checked in by jyg, 11 months ago
Implementation of a two radii model for wake population dynamics. It is activated with : iflag_wk_pop_dyn=3
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 211.4 KB

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

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

Download in other formats:

Original Format