Context Navigation

physiq_mod.F90 @ 4283

Last change on this file since 4283 was 4236, checked in by Laurent Fairhead, 21 months ago
Need to switch to the LMDZ context each time the physics is called if xios output is activated in the dynamics (forgot to commit) LF
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.2 KB

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

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

Download in other formats:

Original Format