Context Navigation

physiq_mod.F90 @ 3436

Last change on this file since 3436 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, 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 Author Date Id Revision

File size: 177.0 KB

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

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

Download in other formats:

Original Format