Context Navigation

physiq_mod.F90 @ 4537

Last change on this file since 4537 was 4537, checked in by fhourdin, 13 months ago

Travail preparatoire au couplage avec mesoNH

Travail preparatoire aux test d'integration de la physique de MesoNH
et plus generalement a la reecriture du moniteur de la physique.
phylmd appelle phylmdex si une cle iflag_physiq=2 est ajoutee dans run.def
Il a ete necessaire en plus d'eliminer un certain nombre d'appels dans
phyetat0 et phyredem.

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: 198.7 KB

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

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

Download in other formats:

Original Format