Context Navigation

physiq_mod.F90 @ 4448

Last change on this file since 4448 was 4448, checked in by fhourdin, 15 months ago
Debut de replay_isation de yamada4
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: 192.3 KB

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

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

Download in other formats:

Original Format