Context Navigation

physiq_mod.F90 @ 3575

Last change on this file since 3575 was 3575, checked in by fhourdin, 5 years ago
Correction suite a la detections de problemes en 1D avec le mode debug sur jean-zay.
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: 182.7 KB

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

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

Download in other formats:

Original Format