Context Navigation

physiq_mod.F90 @ 4358

Last change on this file since 4358 was 4358, checked in by dcugnet, 19 months ago
remove "config_inca" variable from "control_mod" and "infotrac_phy" (read in infotrac) only kept version of "type_trac" is in tracinca ; few tests are moved from infotrac to this module. simplify and generalize a bit the routines "phyetat0_get" and "phyetat0_srf" from phyetat0, converted to a module. fix the isotopic version: few "USE … » were misplaced between ISOVERIF CPP keys fix the old water and derived isotopes names in the ISOTRAC case
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.8 KB

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

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

Download in other formats:

Original Format