Context Navigation

physiq_mod.F90 @ 4636

Last change on this file since 4636 was 4636, checked in by lfalletti, 12 months ago
REPROBUS update: to use a copy of StratAer? emission (volcanic eruption)
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.4 KB

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

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

Download in other formats:

Original Format