Context Navigation

physiq_mod.F90 @ 4601

Last change on this file since 4601 was 4601, checked in by dcugnet, 17 months ago
StratAer? commit, N. Lebas
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: 201.6 KB

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

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

Download in other formats:

Original Format