Context Navigation

physiq_mod.F90 @ 4629

Last change on this file since 4629 was 4629, checked in by dcugnet, 15 months ago
Fix for r4625
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 202.2 KB

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

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

Download in other formats:

Original Format