Context Navigation

physiq_mod.F90 @ 4666

Last change on this file since 4666 was 4666, checked in by fhourdin, 9 months ago
Replayisation lmdz_lscp_old
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: 202.9 KB

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

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

Download in other formats:

Original Format