Context Navigation

physiq_mod.F90 @ 3522

Last change on this file since 3522 was 3522, checked in by oboucher, 5 years ago
Call to initialisation of S3A model
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.6 KB

Line
1	!
2	! $Id: physiq_mod.F90 3522 2019-05-22 21:41:07Z oboucher $
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	STOP
1273	ENDIF
1274
1275	Ncvpaseq1 = 0
1276	dnwd0=0.0
1277	ftd=0.0
1278	fqd=0.0
1279	cin=0.
1280	!ym Attention pbase pas initialise dans concvl !!!!
1281	pbase=0
1282	!IM 180608
1283
1284	itau_con=0
1285	first=.FALSE.
1286
1287	ENDIF ! first
1288
1289	!ym => necessaire pour iflag_con != 2
1290	pmfd(:,:) = 0.
1291	pen_u(:,:) = 0.
1292	pen_d(:,:) = 0.
1293	pde_d(:,:) = 0.
1294	pde_u(:,:) = 0.
1295	aam=0.
1296	d_t_adjwk(:,:)=0
1297	d_q_adjwk(:,:)=0
1298
1299	alp_bl_conv(:)=0.
1300
1301	torsfc=0.
1302	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1303
1304	modname = 'physiq'
1305
1306	IF (debut) THEN
1307	CALL suphel ! initialiser constantes et parametres phys.
1308	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1309	tau_gl=5.
1310	CALL getin_p('tau_gl', tau_gl)
1311	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1312	! secondes
1313	tau_gl=86400.*tau_gl
1314	WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1315
1316	iflag_bug_t2m_ipslcm61 = 1
1317	CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
1318	iflag_bug_t2m_stab_ipslcm61 = -1
1319	CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
1320
1321	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1322	CALL getin_p('random_notrig_max',random_notrig_max)
1323	CALL getin_p('ok_adjwk',ok_adjwk)
1324	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1325	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1326	! 1 => convective adjustment but state variables are unchanged
1327	! 2 => convective adjustment and state variables are changed
1328	CALL getin_p('iflag_adjwk',iflag_adjwk)
1329	CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1330	CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
1331	CALL getin_p('oliqmax',oliqmax)
1332	CALL getin_p('oicemax',oicemax)
1333	CALL getin_p('ratqsp0',ratqsp0)
1334	CALL getin_p('ratqsdp',ratqsdp)
1335	iflag_wake_tend = 0
1336	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1337	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1338	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1339	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1340	CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
1341	CALL getin_p('ok_bug_split_th',ok_bug_split_th)
1342	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1343	CALL getin_p('fl_ebil',fl_ebil)
1344	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1345	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1346	iflag_phytrac = 1 ! by default we do want to call phytrac
1347	CALL getin_p('iflag_phytrac',iflag_phytrac)
1348	nvm_lmdz = 13
1349	CALL getin_p('NVM',nvm_lmdz)
1350
1351	WRITE(lunout,*) 'iflag_alp_wk_cond=', iflag_alp_wk_cond
1352	WRITE(lunout,*) 'random_ntrig_max=', random_notrig_max
1353	WRITE(lunout,*) 'ok_adjwk=', ok_adjwk
1354	WRITE(lunout,*) 'iflag_adjwk=', iflag_adjwk
1355	WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1356	WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1357	WRITE(lunout,*) 'ratqsp0=', ratqsp0
1358	WRITE(lunout,*) 'ratqsdp=', ratqsdp
1359	WRITE(lunout,*) 'iflag_wake_tend=', iflag_wake_tend
1360	WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1361	WRITE(lunout,*) 'ok_bug_cv_trac=', ok_bug_cv_trac
1362	WRITE(lunout,*) 'ok_bug_split_th=', ok_bug_split_th
1363	WRITE(lunout,*) 'fl_ebil=', fl_ebil
1364	WRITE(lunout,*) 'fl_cor_ebil=', fl_cor_ebil
1365	WRITE(lunout,*) 'iflag_phytrac=', iflag_phytrac
1366	WRITE(lunout,*) 'NVM=', nvm_lmdz
1367
1368	!--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1369	WRITE(lunout,*) 'Call to infocfields from physiq'
1370	CALL infocfields_init
1371
1372	ENDIF
1373
1374	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1375
1376	!======================================================================
1377	! Gestion calendrier : mise a jour du module phys_cal_mod
1378	!
1379	! CALL phys_cal_update(jD_cur,jH_cur)
1380
1381	!
1382	! Si c'est le debut, il faut initialiser plusieurs choses
1383	! ********
1384	!
1385	IF (debut) THEN
1386	!rv CRinitialisation de wght_th et lalim_conv pour la
1387	!definition de la couche alimentation de la convection a partir
1388	!des caracteristiques du thermique
1389	wght_th(:,:)=1.
1390	lalim_conv(:)=1
1391	!RC
1392	ustar(:,:)=0.
1393	! u10m(:,:)=0.
1394	! v10m(:,:)=0.
1395	rain_con(:)=0.
1396	snow_con(:)=0.
1397	topswai(:)=0.
1398	topswad(:)=0.
1399	solswai(:)=0.
1400	solswad(:)=0.
1401
1402	wmax_th(:)=0.
1403	tau_overturning_th(:)=0.
1404
1405	IF (type_trac == 'inca') THEN
1406	! jg : initialisation jusqu'au ces variables sont dans restart
1407	ccm(:,:,:) = 0.
1408	tau_aero(:,:,:,:) = 0.
1409	piz_aero(:,:,:,:) = 0.
1410	cg_aero(:,:,:,:) = 0.
1411
1412	config_inca='none' ! default
1413	CALL getin_p('config_inca',config_inca)
1414
1415	ELSE
1416	config_inca='none' ! default
1417	ENDIF
1418
1419	tau_aero(:,:,:,:) = 1.e-15
1420	piz_aero(:,:,:,:) = 1.
1421	cg_aero(:,:,:,:) = 0.
1422
1423	IF (aerosol_couple .AND. (config_inca /= "aero" &
1424	.AND. config_inca /= "aeNP ")) THEN
1425	abort_message &
1426	= 'if aerosol_couple is activated, config_inca need to be ' &
1427	// 'aero or aeNP'
1428	CALL abort_physic (modname,abort_message,1)
1429	ENDIF
1430
1431	rnebcon0(:,:) = 0.0
1432	clwcon0(:,:) = 0.0
1433	rnebcon(:,:) = 0.0
1434	clwcon(:,:) = 0.0
1435
1436	!
1437	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1438	iflag_coupl,iflag_clos,iflag_wake
1439	print*,'iflag_cycle_diurne', iflag_cycle_diurne
1440	!
1441	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1442	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1443	CALL abort_physic (modname,abort_message,1)
1444	ENDIF
1445	!
1446	!
1447	! Initialiser les compteurs:
1448	!
1449	itap = 0
1450	itaprad = 0
1451	itapcv = 0
1452	itapwk = 0
1453
1454	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1455	!! Un petit travail \`a faire ici.
1456	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1457
1458	IF (iflag_pbl>1) THEN
1459	PRINT*, "Using method MELLOR&YAMADA"
1460	ENDIF
1461
1462	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1464	! phylmd plutot que dyn3d
1465	! Attention : la version precedente n'etait pas tres propre.
1466	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1467	! pour obtenir le meme resultat.
1468	!jyg for fh<
1469	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1470	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
1471	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1472	CALL abort_physic(modname,abort_message,1)
1473	ENDIF
1474	!>jyg
1475	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1476	radpas = NINT( 86400./phys_tstep)/nbapp_rad
1477	ELSE
1478	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1479	'multiple de nbapp_rad'
1480	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1481	'mais 1+1<>2'
1482	abort_message='nbre de pas de temps physique n est pas multiple ' &
1483	// 'de nbapp_rad'
1484	CALL abort_physic(modname,abort_message,1)
1485	ENDIF
1486	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1487	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
1488	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1489	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1490	cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
1491	cvpas = cvpas_0
1492	print *,'physiq, cvpas ',cvpas
1493	ELSE
1494	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1495	'multiple de nbapp_cv'
1496	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1497	'mais 1+1<>2'
1498	abort_message='nbre de pas de temps physique n est pas multiple ' &
1499	// 'de nbapp_cv'
1500	CALL abort_physic(modname,abort_message,1)
1501	ENDIF
1502	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1503	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
1504	! print *,'physiq, wkpas ',wkpas
1505	ELSE
1506	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1507	'multiple de nbapp_wk'
1508	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1509	'mais 1+1<>2'
1510	abort_message='nbre de pas de temps physique n est pas multiple ' &
1511	// 'de nbapp_wk'
1512	CALL abort_physic(modname,abort_message,1)
1513	ENDIF
1514	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1515	CALL init_iophy_new(latitude_deg,longitude_deg)
1516
1517	!===================================================================
1518	!IM stations CFMIP
1519	nCFMIP=npCFMIP
1520	OPEN(98,file='npCFMIP_param.data',status='old', &
1521	form='formatted',iostat=iostat)
1522	IF (iostat == 0) THEN
1523	READ(98,*,end=998) nCFMIP
1524	998 CONTINUE
1525	CLOSE(98)
1526	CONTINUE
1527	IF(nCFMIP.GT.npCFMIP) THEN
1528	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1529	CALL abort_physic("physiq", "", 1)
1530	ELSE
1531	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1532	ENDIF
1533
1534	!
1535	ALLOCATE(tabCFMIP(nCFMIP))
1536	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1537	ALLOCATE(tabijGCM(nCFMIP))
1538	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1539	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1540	!
1541	! lecture des nCFMIP stations CFMIP, de leur numero
1542	! et des coordonnees geographiques lonCFMIP, latCFMIP
1543	!
1544	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1545	lonCFMIP, latCFMIP)
1546	!
1547	! identification des
1548	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1549	! grille de LMDZ
1550	! 2) indices points tabijGCM de la grille physique 1d sur
1551	! klon points
1552	! 3) indices iGCM, jGCM de la grille physique 2d
1553	!
1554	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1555	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1556	!
1557	ELSE
1558	ALLOCATE(tabijGCM(0))
1559	ALLOCATE(lonGCM(0), latGCM(0))
1560	ALLOCATE(iGCM(0), jGCM(0))
1561	ENDIF
1562
1563	#ifdef CPP_IOIPSL
1564
1565	!$OMP MASTER
1566	! FH : if ok_sync=.true. , the time axis is written at each time step
1567	! in the output files. Only at the end in the opposite case
1568	ok_sync_omp=.FALSE.
1569	CALL getin('ok_sync',ok_sync_omp)
1570	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1571	iGCM,jGCM,lonGCM,latGCM, &
1572	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1573	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1574	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1575	read_climoz, phys_out_filestations, &
1576	new_aod, aerosol_couple, &
1577	flag_aerosol_strat, pdtphys, paprs, pphis, &
1578	pplay, lmax_th, ptconv, ptconvth, ivap, &
1579	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1580	!$OMP END MASTER
1581	!$OMP BARRIER
1582	ok_sync=ok_sync_omp
1583
1584	freq_outNMC(1) = ecrit_files(7)
1585	freq_outNMC(2) = ecrit_files(8)
1586	freq_outNMC(3) = ecrit_files(9)
1587	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1588	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1589	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1590
1591	#ifndef CPP_XIOS
1592	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1593	#endif
1594
1595	#endif
1596	ecrit_reg = ecrit_reg * un_jour
1597	ecrit_tra = ecrit_tra * un_jour
1598
1599	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1600	date0 = jD_ref
1601	WRITE(,) 'physiq date0 : ',date0
1602	!
1603
1604	! CALL create_climoz(read_climoz)
1605	IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol) !! initialise aero from file for XIOS interpolation (unstructured_grid)
1606	IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat) !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
1607
1608	#ifdef CPP_COSP
1609	IF (ok_cosp) THEN
1610	DO k = 1, klev
1611	DO i = 1, klon
1612	phicosp(i,k) = pphi(i,k) + pphis(i)
1613	ENDDO
1614	ENDDO
1615	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
1616	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1617	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1618	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1619	JrNt,ref_liq,ref_ice, &
1620	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1621	zu10m,zv10m,pphis, &
1622	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1623	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1624	prfl(:,1:klev),psfl(:,1:klev), &
1625	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1626	mr_ozone,cldtau, cldemi)
1627	ENDIF
1628	#endif
1629
1630	#ifdef CPP_COSP2
1631	IF (ok_cosp) THEN
1632	DO k = 1, klev
1633	DO i = 1, klon
1634	phicosp(i,k) = pphi(i,k) + pphis(i)
1635	ENDDO
1636	ENDDO
1637	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
1638	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1639	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1640	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1641	JrNt,ref_liq,ref_ice, &
1642	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1643	zu10m,zv10m,pphis, &
1644	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1645	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1646	prfl(:,1:klev),psfl(:,1:klev), &
1647	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1648	mr_ozone,cldtau, cldemi)
1649	ENDIF
1650	#endif
1651
1652	#ifdef CPP_COSPV2
1653	IF (ok_cosp) THEN
1654	DO k = 1, klev
1655	DO i = 1, klon
1656	phicosp(i,k) = pphi(i,k) + pphis(i)
1657	ENDDO
1658	ENDDO
1659	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1660	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1661	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1662	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1663	JrNt,ref_liq,ref_ice, &
1664	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1665	zu10m,zv10m,pphis, &
1666	phicosp,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1667	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1668	prfl(:,1:klev),psfl(:,1:klev), &
1669	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1670	mr_ozone,cldtau, cldemi)
1671	ENDIF
1672	#endif
1673
1674	!
1675	!
1676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1677	! Nouvelle initialisation pour le rayonnement RRTM
1678	!
1679	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1680
1681	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1682	! Initialisation des champs dans phytrac qui sont utilisés par phys_output_write
1683	IF (iflag_phytrac == 1 ) THEN
1684	CALL phytrac_init()
1685	ENDIF
1686
1687	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
1688	pplay, lmax_th, aerosol_couple, &
1689	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, ok_sync,&
1690	ptconv, read_climoz, clevSTD, &
1691	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
1692	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1693
1694	#ifdef CPP_XIOS
1695	IF (is_omp_master) CALL xios_update_calendar(1)
1696	#endif
1697	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1698	CALL create_etat0_limit_unstruct
1699	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1700
1701	!jyg<
1702	IF (klon_glo==1) THEN
1703	pbl_tke(:,:,is_ave) = 0.
1704	DO nsrf=1,nbsrf
1705	DO k = 1,klev+1
1706	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1707	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1708	ENDDO
1709	ENDDO
1710	ELSE
1711	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1712	!>jyg
1713	ENDIF
1714	!IM begin
1715	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1716	,ratqs(1,1)
1717	!IM end
1718
1719
1720	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1721	!
1722	! on remet le calendrier a zero
1723	!
1724	IF (raz_date .eq. 1) THEN
1725	itau_phy = 0
1726	ENDIF
1727
1728	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1729	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1730	! pdtphys
1731	! abort_message='Pas physique n est pas correct '
1732	! ! call abort_physic(modname,abort_message,1)
1733	! phys_tstep=pdtphys
1734	! ENDIF
1735	IF (nlon .NE. klon) THEN
1736	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1737	klon
1738	abort_message='nlon et klon ne sont pas coherents'
1739	CALL abort_physic(modname,abort_message,1)
1740	ENDIF
1741	IF (nlev .NE. klev) THEN
1742	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1743	klev
1744	abort_message='nlev et klev ne sont pas coherents'
1745	CALL abort_physic(modname,abort_message,1)
1746	ENDIF
1747	!
1748	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1749	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1750	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1751	abort_message='Nbre d appels au rayonnement insuffisant'
1752	CALL abort_physic(modname,abort_message,1)
1753	ENDIF
1754	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1755	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1756	ok_cvl
1757	!
1758	!KE43
1759	! Initialisation pour la convection de K.E. (sb):
1760	IF (iflag_con.GE.3) THEN
1761
1762	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1763	WRITE(lunout,*) &
1764	"On va utiliser le melange convectif des traceurs qui"
1765	WRITE(lunout,*)"est calcule dans convect4.3"
1766	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1767
1768	DO i = 1, klon
1769	ema_cbmf(i) = 0.
1770	ema_pcb(i) = 0.
1771	ema_pct(i) = 0.
1772	! ema_workcbmf(i) = 0.
1773	ENDDO
1774	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1775	DO i = 1, klon
1776	ibas_con(i) = 1
1777	itop_con(i) = 1
1778	ENDDO
1779	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1780	!================================================================
1781	!CR:04.12.07: initialisations poches froides
1782	! Controle de ALE et ALP pour la fermeture convective (jyg)
1783	IF (iflag_wake>=1) THEN
1784	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1785	,alp_bl_prescr, ale_bl_prescr)
1786	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1787	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1788	!
1789	! Initialize tendencies of wake state variables (for some flag values
1790	! they are not computed).
1791	d_deltat_wk(:,:) = 0.
1792	d_deltaq_wk(:,:) = 0.
1793	d_deltat_wk_gw(:,:) = 0.
1794	d_deltaq_wk_gw(:,:) = 0.
1795	d_deltat_vdf(:,:) = 0.
1796	d_deltaq_vdf(:,:) = 0.
1797	d_deltat_the(:,:) = 0.
1798	d_deltaq_the(:,:) = 0.
1799	d_deltat_ajs_cv(:,:) = 0.
1800	d_deltaq_ajs_cv(:,:) = 0.
1801	d_s_wk(:) = 0.
1802	d_dens_wk(:) = 0.
1803	ENDIF
1804
1805	! do i = 1,klon
1806	! Ale_bl(i)=0.
1807	! Alp_bl(i)=0.
1808	! enddo
1809
1810	!ELSE
1811	! ALLOCATE(tabijGCM(0))
1812	! ALLOCATE(lonGCM(0), latGCM(0))
1813	! ALLOCATE(iGCM(0), jGCM(0))
1814	ENDIF
1815
1816	DO i=1,klon
1817	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1818	ENDDO
1819
1820	!34EK
1821	IF (ok_orodr) THEN
1822
1823	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1824	! FH sans doute a enlever de finitivement ou, si on le
1825	! garde, l'activer justement quand ok_orodr = false.
1826	! ce rugoro est utilise par la couche limite et fait double emploi
1827	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1828	! DO i=1,klon
1829	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1830	! ENDDO
1831	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1832	IF (ok_strato) THEN
1833	CALL SUGWD_strato(klon,klev,paprs,pplay)
1834	ELSE
1835	CALL SUGWD(klon,klev,paprs,pplay)
1836	ENDIF
1837
1838	DO i=1,klon
1839	zuthe(i)=0.
1840	zvthe(i)=0.
1841	IF (zstd(i).gt.10.) THEN
1842	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1843	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1844	ENDIF
1845	ENDDO
1846	ENDIF
1847	!
1848	!
1849	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
1850	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1851	lmt_pas
1852	!
1853	capemaxcels = 't_max(X)'
1854	t2mincels = 't_min(X)'
1855	t2maxcels = 't_max(X)'
1856	tinst = 'inst(X)'
1857	tave = 'ave(X)'
1858	!IM cf. AM 081204 BEG
1859	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1860	!IM cf. AM 081204 END
1861	!
1862	!=============================================================
1863	! Initialisation des sorties
1864	!=============================================================
1865
1866	#ifdef CPP_XIOS
1867	! Get "missing_val" value from XML files (from temperature variable)
1868	!$OMP MASTER
1869	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1870	!$OMP END MASTER
1871	!$OMP BARRIER
1872	missing_val=missing_val_omp
1873	#endif
1874
1875	#ifdef CPP_XIOS
1876	! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1877	! initialised at that moment
1878	! Get "missing_val" value from XML files (from temperature variable)
1879	!$OMP MASTER
1880	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1881	!$OMP END MASTER
1882	!$OMP BARRIER
1883	missing_val=missing_val_omp
1884	#endif
1885
1886
1887	CALL printflag( tabcntr0,radpas,ok_journe, &
1888	ok_instan, ok_region )
1889	!
1890	!
1891	!
1892	! Prescrire l'ozone dans l'atmosphere
1893	!
1894	!
1895	!c DO i = 1, klon
1896	!c DO k = 1, klev
1897	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1898	!c ENDDO
1899	!c ENDDO
1900	!
1901	IF (type_trac == 'inca') THEN
1902	#ifdef INCA
1903	CALL VTe(VTphysiq)
1904	CALL VTb(VTinca)
1905	calday = REAL(days_elapsed) + jH_cur
1906	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1907
1908	CALL chemini( &
1909	rg, &
1910	ra, &
1911	cell_area, &
1912	latitude_deg, &
1913	longitude_deg, &
1914	presnivs, &
1915	calday, &
1916	klon, &
1917	nqtot, &
1918	nqo, &
1919	pdtphys, &
1920	annee_ref, &
1921	year_cur, &
1922	day_ref, &
1923	day_ini, &
1924	start_time, &
1925	itau_phy, &
1926	date0, &
1927	io_lon, &
1928	io_lat, &
1929	chemistry_couple, &
1930	init_source, &
1931	init_tauinca, &
1932	init_pizinca, &
1933	init_cginca, &
1934	init_ccminca)
1935
1936
1937	! initialisation des variables depuis le restart de inca
1938	ccm(:,:,:) = init_ccminca
1939	tau_aero(:,:,:,:) = init_tauinca
1940	piz_aero(:,:,:,:) = init_pizinca
1941	cg_aero(:,:,:,:) = init_cginca
1942	!
1943
1944
1945	CALL VTe(VTinca)
1946	CALL VTb(VTphysiq)
1947	#endif
1948	ENDIF
1949
1950	!$omp single
1951	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1952	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1953	!$omp end single
1954	!
1955	!IM betaCRF
1956	pfree=70000. !Pa
1957	beta_pbl=1.
1958	beta_free=1.
1959	lon1_beta=-180.
1960	lon2_beta=+180.
1961	lat1_beta=90.
1962	lat2_beta=-90.
1963	mskocean_beta=.FALSE.
1964
1965	!albedo SB >>>
1966	SELECT CASE(nsw)
1967	CASE(2)
1968	SFRWL(1)=0.45538747
1969	SFRWL(2)=0.54461211
1970	CASE(4)
1971	SFRWL(1)=0.45538747
1972	SFRWL(2)=0.32870591
1973	SFRWL(3)=0.18568763
1974	SFRWL(4)=3.02191470E-02
1975	CASE(6)
1976	SFRWL(1)=1.28432794E-03
1977	SFRWL(2)=0.12304168
1978	SFRWL(3)=0.33106142
1979	SFRWL(4)=0.32870591
1980	SFRWL(5)=0.18568763
1981	SFRWL(6)=3.02191470E-02
1982	END SELECT
1983
1984
1985	!albedo SB <<<
1986
1987	OPEN(99,file='beta_crf.data',status='old', &
1988	form='formatted',err=9999)
1989	READ(99,*,end=9998) pfree
1990	READ(99,*,end=9998) beta_pbl
1991	READ(99,*,end=9998) beta_free
1992	READ(99,*,end=9998) lon1_beta
1993	READ(99,*,end=9998) lon2_beta
1994	READ(99,*,end=9998) lat1_beta
1995	READ(99,*,end=9998) lat2_beta
1996	READ(99,*,end=9998) mskocean_beta
1997	9998 Continue
1998	CLOSE(99)
1999	9999 Continue
2000	WRITE(,)'pfree=',pfree
2001	WRITE(,)'beta_pbl=',beta_pbl
2002	WRITE(,)'beta_free=',beta_free
2003	WRITE(,)'lon1_beta=',lon1_beta
2004	WRITE(,)'lon2_beta=',lon2_beta
2005	WRITE(,)'lat1_beta=',lat1_beta
2006	WRITE(,)'lat2_beta=',lat2_beta
2007	WRITE(,)'mskocean_beta=',mskocean_beta
2008
2009	!lwoff=y : offset LW CRE for radiation code and other schemes
2010	!lwoff=y : betalwoff=1.
2011	betalwoff=0.
2012	IF (ok_lwoff) THEN
2013	betalwoff=1.
2014	ENDIF
2015	WRITE(,)'ok_lwoff=',ok_lwoff
2016	!
2017	!lwoff=y to begin only sollw and sollwdown are set up to CS values
2018	sollw = sollw + betalwoff * (sollw0 - sollw)
2019	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
2020	sollwdown(:))
2021	ENDIF
2022	!
2023	! ************** Fin de IF ( debut ) *************
2024	!
2025	!
2026	! Incrementer le compteur de la physique
2027	!
2028	itap = itap + 1
2029	IF (is_master .OR. prt_level > 9) THEN
2030	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
2031	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2032	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2033	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
2034	ENDIF
2035	ENDIF
2036	!
2037	!
2038	! Update fraction of the sub-surfaces (pctsrf) and
2039	! initialize, where a new fraction has appeared, all variables depending
2040	! on the surface fraction.
2041	!
2042	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
2043	pctsrf, fevap, z0m, z0h, agesno, &
2044	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
2045
2046	! Update time and other variables in Reprobus
2047	IF (type_trac == 'repr') THEN
2048	#ifdef REPROBUS
2049	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2050	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2051	CALL Rtime(debut)
2052	#endif
2053	ENDIF
2054
2055	! Tendances bidons pour les processus qui n'affectent pas certaines
2056	! variables.
2057	du0(:,:)=0.
2058	dv0(:,:)=0.
2059	dt0 = 0.
2060	dq0(:,:)=0.
2061	dql0(:,:)=0.
2062	dqi0(:,:)=0.
2063	dsig0(:) = 0.
2064	ddens0(:) = 0.
2065	wkoccur1(:)=1
2066	!
2067	! Mettre a zero des variables de sortie (pour securite)
2068	!
2069	DO i = 1, klon
2070	d_ps(i) = 0.0
2071	ENDDO
2072	DO k = 1, klev
2073	DO i = 1, klon
2074	d_t(i,k) = 0.0
2075	d_u(i,k) = 0.0
2076	d_v(i,k) = 0.0
2077	ENDDO
2078	ENDDO
2079	DO iq = 1, nqtot
2080	DO k = 1, klev
2081	DO i = 1, klon
2082	d_qx(i,k,iq) = 0.0
2083	ENDDO
2084	ENDDO
2085	ENDDO
2086	beta_prec_fisrt(:,:)=0.
2087	beta_prec(:,:)=0.
2088	!
2089	! Output variables from the convective scheme should not be set to 0
2090	! since convection is not always called at every time step.
2091	IF (ok_bug_cv_trac) THEN
2092	da(:,:)=0.
2093	mp(:,:)=0.
2094	phi(:,:,:)=0.
2095	! RomP >>>
2096	phi2(:,:,:)=0.
2097	epmlmMm(:,:,:)=0.
2098	eplaMm(:,:)=0.
2099	d1a(:,:)=0.
2100	dam(:,:)=0.
2101	pmflxr(:,:)=0.
2102	pmflxs(:,:)=0.
2103	! RomP <<<
2104	ENDIF
2105
2106	!
2107	! Ne pas affecter les valeurs entrees de u, v, h, et q
2108	!
2109	DO k = 1, klev
2110	DO i = 1, klon
2111	t_seri(i,k) = t(i,k)
2112	u_seri(i,k) = u(i,k)
2113	v_seri(i,k) = v(i,k)
2114	q_seri(i,k) = qx(i,k,ivap)
2115	ql_seri(i,k) = qx(i,k,iliq)
2116	!CR: ATTENTION, on rajoute la variable glace
2117	IF (nqo.eq.2) THEN
2118	qs_seri(i,k) = 0.
2119	ELSE IF (nqo.eq.3) THEN
2120	qs_seri(i,k) = qx(i,k,isol)
2121	ENDIF
2122	ENDDO
2123	ENDDO
2124	!
2125	!--OB mass fixer
2126	IF (mass_fixer) THEN
2127	!--store initial water burden
2128	qql1(:)=0.0
2129	DO k = 1, klev
2130	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2131	ENDDO
2132	ENDIF
2133	!--fin mass fixer
2134
2135	tke0(:,:)=pbl_tke(:,:,is_ave)
2136	!CR:Nombre de traceurs de l'eau: nqo
2137	! IF (nqtot.GE.3) THEN
2138	IF (nqtot.GE.(nqo+1)) THEN
2139	! DO iq = 3, nqtot
2140	DO iq = nqo+1, nqtot
2141	DO k = 1, klev
2142	DO i = 1, klon
2143	! tr_seri(i,k,iq-2) = qx(i,k,iq)
2144	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
2145	ENDDO
2146	ENDDO
2147	ENDDO
2148	ELSE
2149	DO k = 1, klev
2150	DO i = 1, klon
2151	tr_seri(i,k,1) = 0.0
2152	ENDDO
2153	ENDDO
2154	ENDIF
2155	!
2156	DO i = 1, klon
2157	ztsol(i) = 0.
2158	ENDDO
2159	DO nsrf = 1, nbsrf
2160	DO i = 1, klon
2161	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2162	ENDDO
2163	ENDDO
2164	! Initialize variables used for diagnostic purpose
2165	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2166
2167	! Diagnostiquer la tendance dynamique
2168	!
2169	IF (ancien_ok) THEN
2170	!
2171	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2172	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2173	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2174	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2175	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2176	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2177	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2178	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2179	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2180	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2181	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2182	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2183	! !! RomP >>> td dyn traceur
2184	IF (nqtot.GT.nqo) THEN ! jyg
2185	DO iq = nqo+1, nqtot ! jyg
2186	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
2187	ENDDO
2188	ENDIF
2189	! !! RomP <<<
2190	ELSE
2191	d_u_dyn(:,:) = 0.0
2192	d_v_dyn(:,:) = 0.0
2193	d_t_dyn(:,:) = 0.0
2194	d_q_dyn(:,:) = 0.0
2195	d_ql_dyn(:,:) = 0.0
2196	d_qs_dyn(:,:) = 0.0
2197	d_q_dyn2d(:) = 0.0
2198	d_ql_dyn2d(:) = 0.0
2199	d_qs_dyn2d(:) = 0.0
2200	! !! RomP >>> td dyn traceur
2201	IF (nqtot.GT.nqo) THEN ! jyg
2202	DO iq = nqo+1, nqtot ! jyg
2203	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
2204	ENDDO
2205	ENDIF
2206	! !! RomP <<<
2207	ancien_ok = .TRUE.
2208	ENDIF
2209	!
2210	! Ajouter le geopotentiel du sol:
2211	!
2212	DO k = 1, klev
2213	DO i = 1, klon
2214	zphi(i,k) = pphi(i,k) + pphis(i)
2215	ENDDO
2216	ENDDO
2217	!
2218	! Verifier les temperatures
2219	!
2220	!IM BEG
2221	IF (check) THEN
2222	amn=MIN(ftsol(1,is_ter),1000.)
2223	amx=MAX(ftsol(1,is_ter),-1000.)
2224	DO i=2, klon
2225	amn=MIN(ftsol(i,is_ter),amn)
2226	amx=MAX(ftsol(i,is_ter),amx)
2227	ENDDO
2228	!
2229	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2230	ENDIF !(check) THEN
2231	!IM END
2232	!
2233	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2234	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2235
2236	!
2237	!IM BEG
2238	IF (check) THEN
2239	amn=MIN(ftsol(1,is_ter),1000.)
2240	amx=MAX(ftsol(1,is_ter),-1000.)
2241	DO i=2, klon
2242	amn=MIN(ftsol(i,is_ter),amn)
2243	amx=MAX(ftsol(i,is_ter),amx)
2244	ENDDO
2245	!
2246	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2247	ENDIF !(check) THEN
2248	!IM END
2249	!
2250	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2251	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2252	!
2253	! Update ozone if day change
2254	IF (MOD(itap-1,lmt_pas) == 0) THEN
2255	IF (read_climoz <= 0) THEN
2256	! Once per day, update ozone from Royer:
2257	IF (solarlong0<-999.) then
2258	! Generic case with evolvoing season
2259	zzz=real(days_elapsed+1)
2260	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2261	! Particular case with annual mean insolation
2262	zzz=real(90) ! could be revisited
2263	IF (read_climoz/=-1) THEN
2264	abort_message ='read_climoz=-1 is recommended when ' &
2265	// 'solarlong0=1000.'
2266	CALL abort_physic (modname,abort_message,1)
2267	ENDIF
2268	ELSE
2269	! Case where the season is imposed with solarlong0
2270	zzz=real(90) ! could be revisited
2271	ENDIF
2272
2273	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2274	ELSE
2275	!--- ro3i = elapsed days number since current year 1st january, 0h
2276	ro3i=days_elapsed+jh_cur-jh_1jan
2277	!--- scaling for old style files (360 records)
2278	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2279	IF(adjust_tropopause) THEN
2280	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2281	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2282	time_climoz , longitude_deg, latitude_deg, &
2283	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2284	ELSE
2285	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2286	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2287	time_climoz )
2288	ENDIF
2289	! Convert from mole fraction of ozone to column density of ozone in a
2290	! cell, in kDU:
2291	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2292	* zmasse / dobson_u / 1e3
2293	! (By regridding ozone values for LMDZ only once a day, we
2294	! have already neglected the variation of pressure in one
2295	! day. So do not recompute "wo" at each time step even if
2296	! "zmasse" changes a little.)
2297	ENDIF
2298	ENDIF
2299	!
2300	! Re-evaporer l'eau liquide nuageuse
2301	!
2302	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2303	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2304
2305	CALL add_phys_tend &
2306	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2307	'eva',abortphy,flag_inhib_tend,itap,0)
2308	CALL prt_enerbil('eva',itap)
2309
2310	!=========================================================================
2311	! Calculs de l'orbite.
2312	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2313	! doit donc etre plac\'e avant radlwsw et pbl_surface
2314
2315	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2316	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2317	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2318	!
2319	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2320	! solarlong0
2321	IF (solarlong0<-999.) THEN
2322	IF (new_orbit) THEN
2323	! calcul selon la routine utilisee pour les planetes
2324	CALL solarlong(day_since_equinox, zlongi, dist)
2325	ELSE
2326	! calcul selon la routine utilisee pour l'AR4
2327	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2328	ENDIF
2329	ELSE
2330	zlongi=solarlong0 ! longitude solaire vraie
2331	dist=1. ! distance au soleil / moyenne
2332	ENDIF
2333
2334	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2335
2336
2337	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2338	! Calcul de l'ensoleillement :
2339	! ============================
2340	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2341	! l'annee a partir d'une formule analytique.
2342	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2343	! non nul aux poles.
2344	IF (abs(solarlong0-1000.)<1.e-4) THEN
2345	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2346	latitude_deg,longitude_deg,rmu0,fract)
2347	swradcorr(:) = 1.0
2348	JrNt(:) = 1.0
2349	zrmu0(:) = rmu0(:)
2350	ELSE
2351	! recode par Olivier Boucher en sept 2015
2352	SELECT CASE (iflag_cycle_diurne)
2353	CASE(0)
2354	! Sans cycle diurne
2355	CALL angle(zlongi, latitude_deg, fract, rmu0)
2356	swradcorr = 1.0
2357	JrNt = 1.0
2358	zrmu0 = rmu0
2359	CASE(1)
2360	! Avec cycle diurne sans application des poids
2361	! bit comparable a l ancienne formulation cycle_diurne=true
2362	! on integre entre gmtime et gmtime+radpas
2363	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2364	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2365	latitude_deg,longitude_deg,rmu0,fract)
2366	zrmu0 = rmu0
2367	swradcorr = 1.0
2368	! Calcul du flag jour-nuit
2369	JrNt = 0.0
2370	WHERE (fract.GT.0.0) JrNt = 1.0
2371	CASE(2)
2372	! Avec cycle diurne sans application des poids
2373	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2374	! Comme cette routine est appele a tous les pas de temps de
2375	! la physique meme si le rayonnement n'est pas appele je
2376	! remonte en arriere les radpas-1 pas de temps
2377	! suivant. Petite ruse avec MOD pour prendre en compte le
2378	! premier pas de temps de la physique pendant lequel
2379	! itaprad=0
2380	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2381	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2382	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2383	latitude_deg,longitude_deg,rmu0,fract)
2384	!
2385	! Calcul des poids
2386	!
2387	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2388	zdtime2=0.0 !--pas de temps de la physique qui se termine
2389	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2390	latitude_deg,longitude_deg,zrmu0,zfract)
2391	swradcorr = 0.0
2392	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2393	swradcorr=zfract/fract*zrmu0/rmu0
2394	! Calcul du flag jour-nuit
2395	JrNt = 0.0
2396	WHERE (zfract.GT.0.0) JrNt = 1.0
2397	END SELECT
2398	ENDIF
2399	sza_o = ACOS (rmu0) *180./pi
2400
2401	IF (mydebug) THEN
2402	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2403	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2404	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2405	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2406	ENDIF
2407
2408	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2409	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2410	! Cela implique tous les interactions des sous-surfaces et la
2411	! partie diffusion turbulent du couche limit.
2412	!
2413	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2414	! ecriture des fihiers hist_XXXX.nc, ces sont :
2415	! qsol, zq2m, s_pblh, s_lcl,
2416	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2417	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2418	! zu10m, zv10m, fder,
2419	! zxqsurf, rh2m, zxfluxu, zxfluxv,
2420	! frugs, agesno, fsollw, fsolsw,
2421	! d_ts, fevap, fluxlat, t2m,
2422	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2423	!
2424	! Certains ne sont pas utiliser du tout :
2425	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2426	!
2427
2428	! Calcul de l'humidite de saturation au niveau du sol
2429
2430
2431
2432	IF (iflag_pbl/=0) THEN
2433
2434	!jyg+nrlmd<
2435	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2436	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2437	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2438	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2439	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2440	ENDIF
2441	! !!
2442	!>jyg+nrlmd
2443	!
2444	!-------gustiness calculation-------!
2445	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2446	gustiness=0 !ym missing init
2447
2448	IF (iflag_gusts==0) THEN
2449	gustiness(1:klon)=0
2450	ELSE IF (iflag_gusts==1) THEN
2451	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2452	ELSE IF (iflag_gusts==2) THEN
2453	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2454	! ELSE IF (iflag_gusts==2) THEN
2455	! do i = 1, klon
2456	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2457	! *ale_wake(i) !! need to make sigma_wk accessible here
2458	! enddo
2459	! ELSE IF (iflag_gusts==3) THEN
2460	! do i = 1, klon
2461	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2462	! enddo
2463	ENDIF
2464
2465	CALL pbl_surface( &
2466	phys_tstep, date0, itap, days_elapsed+1, &
2467	debut, lafin, &
2468	longitude_deg, latitude_deg, rugoro, zrmu0, &
2469	zsig, sollwdown, pphi, cldt, &
2470	rain_fall, snow_fall, solsw, sollw, &
2471	gustiness, &
2472	t_seri, q_seri, u_seri, v_seri, &
2473	!nrlmd+jyg<
2474	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2475	!>nrlmd+jyg
2476	pplay, paprs, pctsrf, &
2477	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2478	!albedo SB <<<
2479	cdragh, cdragm, u1, v1, &
2480	!albedo SB >>>
2481	! albsol1, albsol2, sens, evap, &
2482	albsol_dir, albsol_dif, sens, evap, &
2483	!albedo SB <<<
2484	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2485	zxtsol, zxfluxlat, zt2m, qsat2m, &
2486	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2487	!nrlmd<
2488	!jyg<
2489	d_t_vdf_w, d_q_vdf_w, &
2490	d_t_vdf_x, d_q_vdf_x, &
2491	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2492	!>jyg
2493	delta_tsurf,wake_dens, &
2494	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2495	kh,kh_x,kh_w, &
2496	!>nrlmd
2497	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2498	slab_wfbils, &
2499	qsol, zq2m, s_pblh, s_lcl, &
2500	!jyg<
2501	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2502	!>jyg
2503	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2504	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2505	zustar, zu10m, zv10m, fder, &
2506	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2507	z0m, z0h, agesno, fsollw, fsolsw, &
2508	d_ts, fevap, fluxlat, t2m, &
2509	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2510	fluxt, fluxu, fluxv, &
2511	dsens, devap, zxsnow, &
2512	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2513	!nrlmd+jyg<
2514	wake_delta_pbl_TKE, &
2515	!>nrlmd+jyg
2516	treedrg )
2517	!FC
2518	!
2519	! Add turbulent diffusion tendency to the wake difference variables
2520	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2521	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2522	!jyg<
2523	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2524	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2525	CALL add_wake_tend &
2526	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2527	ELSE
2528	d_deltat_vdf(:,:) = 0.
2529	d_deltaq_vdf(:,:) = 0.
2530	!>jyg
2531	ENDIF
2532
2533	!add limitation for t,q at and wind at 10m
2534	if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
2535	CALL borne_var_surf( klon,klev,nbsrf, &
2536	iflag_bug_t2m_stab_ipslcm61, &
2537	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1), &
2538	ftsol,zxqsurf,pctsrf,paprs, &
2539	t2m, q2m, u10m, v10m, &
2540	zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor, &
2541	zrh2m_cor, zqsat2m_cor)
2542	ELSE
2543	zt2m_cor(:)=zt2m(:)
2544	zq2m_cor(:)=zq2m(:)
2545	zu10m_cor(:)=zu10m(:)
2546	zv10m_cor(:)=zv10m(:)
2547	ENDIF
2548
2549	!---------------------------------------------------------------------
2550	! ajout des tendances de la diffusion turbulente
2551	IF (klon_glo==1) THEN
2552	CALL add_pbl_tend &
2553	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2554	'vdf',abortphy,flag_inhib_tend,itap)
2555	ELSE
2556	CALL add_phys_tend &
2557	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2558	'vdf',abortphy,flag_inhib_tend,itap,0)
2559	ENDIF
2560	CALL prt_enerbil('vdf',itap)
2561	!--------------------------------------------------------------------
2562
2563	IF (mydebug) THEN
2564	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2565	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2566	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2567	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2568	ENDIF
2569
2570	!albedo SB >>>
2571	albsol1=0.
2572	albsol2=0.
2573	falb1=0.
2574	falb2=0.
2575	SELECT CASE(nsw)
2576	CASE(2)
2577	albsol1=albsol_dir(:,1)
2578	albsol2=albsol_dir(:,2)
2579	falb1=falb_dir(:,1,:)
2580	falb2=falb_dir(:,2,:)
2581	CASE(4)
2582	albsol1=albsol_dir(:,1)
2583	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2584	+albsol_dir(:,4)*SFRWL(4)
2585	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2586	falb1=falb_dir(:,1,:)
2587	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2588	+falb_dir(:,4,:)*SFRWL(4)
2589	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2590	CASE(6)
2591	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2592	+albsol_dir(:,3)*SFRWL(3)
2593	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2594	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2595	+albsol_dir(:,6)*SFRWL(6)
2596	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2597	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2598	+falb_dir(:,3,:)*SFRWL(3)
2599	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2600	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2601	+falb_dir(:,6,:)*SFRWL(6)
2602	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2603	END SELECt
2604	!albedo SB <<<
2605
2606
2607	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2608	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2609
2610	ENDIF
2611	! =================================================================== c
2612	! Calcul de Qsat
2613
2614	DO k = 1, klev
2615	DO i = 1, klon
2616	zx_t = t_seri(i,k)
2617	IF (thermcep) THEN
2618	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2619	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2620	zx_qs = MIN(0.5,zx_qs)
2621	zcor = 1./(1.-retv*zx_qs)
2622	zx_qs = zx_qs*zcor
2623	ELSE
2624	!! IF (zx_t.LT.t_coup) THEN !jyg
2625	IF (zx_t.LT.rtt) THEN !jyg
2626	zx_qs = qsats(zx_t)/pplay(i,k)
2627	ELSE
2628	zx_qs = qsatl(zx_t)/pplay(i,k)
2629	ENDIF
2630	ENDIF
2631	zqsat(i,k)=zx_qs
2632	ENDDO
2633	ENDDO
2634
2635	IF (prt_level.ge.1) THEN
2636	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2637	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2638	ENDIF
2639	!
2640	! Appeler la convection (au choix)
2641	!
2642	DO k = 1, klev
2643	DO i = 1, klon
2644	conv_q(i,k) = d_q_dyn(i,k) &
2645	+ d_q_vdf(i,k)/phys_tstep
2646	conv_t(i,k) = d_t_dyn(i,k) &
2647	+ d_t_vdf(i,k)/phys_tstep
2648	ENDDO
2649	ENDDO
2650	IF (check) THEN
2651	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2652	WRITE(lunout,*) "avantcon=", za
2653	ENDIF
2654	zx_ajustq = .FALSE.
2655	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2656	IF (zx_ajustq) THEN
2657	DO i = 1, klon
2658	z_avant(i) = 0.0
2659	ENDDO
2660	DO k = 1, klev
2661	DO i = 1, klon
2662	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2663	*(paprs(i,k)-paprs(i,k+1))/RG
2664	ENDDO
2665	ENDDO
2666	ENDIF
2667
2668	! Calcule de vitesse verticale a partir de flux de masse verticale
2669	DO k = 1, klev
2670	DO i = 1, klon
2671	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2672	ENDDO
2673	ENDDO
2674
2675	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2676	omega(igout, :)
2677	!
2678	! Appel de la convection tous les "cvpas"
2679	!
2680	!!jyg IF (MOD(itapcv,cvpas).EQ.0) THEN
2681	!! print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2682	!! itapcv, cvpas, itap-1, cvpas_0
2683	IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
2684
2685	!
2686	! Mettre a zero des variables de sortie (pour securite)
2687	!
2688	pmflxr(:,:) = 0.
2689	pmflxs(:,:) = 0.
2690	wdtrainA(:,:) = 0.
2691	wdtrainS(:,:) = 0.
2692	wdtrainM(:,:) = 0.
2693	upwd(:,:) = 0.
2694	dnwd(:,:) = 0.
2695	ep(:,:) = 0.
2696	da(:,:)=0.
2697	mp(:,:)=0.
2698	wght_cvfd(:,:)=0.
2699	phi(:,:,:)=0.
2700	phi2(:,:,:)=0.
2701	epmlmMm(:,:,:)=0.
2702	eplaMm(:,:)=0.
2703	d1a(:,:)=0.
2704	dam(:,:)=0.
2705	elij(:,:,:)=0.
2706	ev(:,:)=0.
2707	qtaa(:,:)=0.
2708	clw(:,:)=0.
2709	sij(:,:,:)=0.
2710	!
2711	IF (iflag_con.EQ.1) THEN
2712	abort_message ='reactiver le call conlmd dans physiq.F'
2713	CALL abort_physic (modname,abort_message,1)
2714	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
2715	! . d_t_con, d_q_con,
2716	! . rain_con, snow_con, ibas_con, itop_con)
2717	ELSE IF (iflag_con.EQ.2) THEN
2718	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
2719	conv_t, conv_q, -evap, omega, &
2720	d_t_con, d_q_con, rain_con, snow_con, &
2721	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2722	kcbot, kctop, kdtop, pmflxr, pmflxs)
2723	d_u_con = 0.
2724	d_v_con = 0.
2725
2726	WHERE (rain_con < 0.) rain_con = 0.
2727	WHERE (snow_con < 0.) snow_con = 0.
2728	DO i = 1, klon
2729	ibas_con(i) = klev+1 - kcbot(i)
2730	itop_con(i) = klev+1 - kctop(i)
2731	ENDDO
2732	ELSE IF (iflag_con.GE.3) THEN
2733	! nb of tracers for the KE convection:
2734	! MAF la partie traceurs est faite dans phytrac
2735	! on met ntra=1 pour limiter les appels mais on peut
2736	! supprimer les calculs / ftra.
2737	ntra = 1
2738
2739	!=======================================================================
2740	!ajout pour la parametrisation des poches froides: calcul de
2741	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2742	IF (iflag_wake>=1) THEN
2743	DO k=1,klev
2744	DO i=1,klon
2745	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2746	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2747	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2748	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2749	ENDDO
2750	ENDDO
2751	ELSE
2752	t_w(:,:) = t_seri(:,:)
2753	q_w(:,:) = q_seri(:,:)
2754	t_x(:,:) = t_seri(:,:)
2755	q_x(:,:) = q_seri(:,:)
2756	ENDIF
2757	!
2758	!jyg<
2759	! Perform dry adiabatic adjustment on wake profile
2760	! The corresponding tendencies are added to the convective tendencies
2761	! after the call to the convective scheme.
2762	IF (iflag_wake>=1) then
2763	IF (iflag_adjwk >= 1) THEN
2764	limbas(:) = 1
2765	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2766	d_t_adjwk, d_q_adjwk)
2767	!
2768	DO k=1,klev
2769	DO i=1,klon
2770	IF (wake_s(i) .GT. 1.e-3) THEN
2771	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2772	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2773	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2774	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2775	ELSE
2776	d_deltat_ajs_cv(i,k) = 0.
2777	d_deltaq_ajs_cv(i,k) = 0.
2778	ENDIF
2779	ENDDO
2780	ENDDO
2781	IF (iflag_adjwk == 2) THEN
2782	CALL add_wake_tend &
2783	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2784	ENDIF ! (iflag_adjwk == 2)
2785	ENDIF ! (iflag_adjwk >= 1)
2786	ENDIF ! (iflag_wake>=1)
2787	!>jyg
2788	!
2789
2790	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
2791	!! (k, q_w(1,k), q_x(1,k),k=1,25)
2792
2793	!jyg<
2794	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
2795	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
2796	ale_bl_prescr, alp_bl_prescr, &
2797	wake_pe, wake_fip, &
2798	Ale_bl, Ale_bl_trig, Alp_bl, &
2799	Ale, Alp , Ale_wake, Alp_wake)
2800	!>jyg
2801	!
2802	! sb, oct02:
2803	! Schema de convection modularise et vectorise:
2804	! (driver commun aux versions 3 et 4)
2805	!
2806	IF (ok_cvl) THEN ! new driver for convectL
2807	!
2808	!jyg<
2809	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2810	! Calculate the upmost level of deep convection loops: k_upper_cv
2811	! (near 22 km)
2812	k_upper_cv = klev
2813	!izero = klon/2+1/klon
2814	!DO k = klev,1,-1
2815	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2816	!ENDDO
2817	! FH : nouveau calcul base sur un profil global sans quoi
2818	! le modele etait sensible au decoupage de domaines
2819	DO k = klev,1,-1
2820	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
2821	ENDDO
2822	IF (prt_level .ge. 5) THEN
2823	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2824	k_upper_cv
2825	ENDIF
2826	!
2827	!>jyg
2828	IF (type_trac == 'repr') THEN
2829	nbtr_tmp=ntra
2830	ELSE
2831	nbtr_tmp=nbtr
2832	ENDIF
2833	!jyg iflag_con est dans clesphys
2834	!c CALL concvl (iflag_con,iflag_clos,
2835	CALL concvl (iflag_clos, &
2836	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
2837	t_w,q_w,wake_s, &
2838	u_seri,v_seri,tr_seri,nbtr_tmp, &
2839	ALE,ALP, &
2840	sig1,w01, &
2841	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2842	rain_con, snow_con, ibas_con, itop_con, sigd, &
2843	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2844	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2845	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2846	! RomP >>>
2847	!! . pmflxr,pmflxs,da,phi,mp,
2848	!! . ftd,fqd,lalim_conv,wght_th)
2849	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
2850	ftd,fqd,lalim_conv,wght_th, &
2851	ev, ep,epmlmMm,eplaMm, &
2852	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2853	tau_cld_cv,coefw_cld_cv,epmax_diag)
2854
2855	! RomP <<<
2856
2857	!IM begin
2858	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2859	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2860	!IM end
2861	!IM cf. FH
2862	clwcon0=qcondc
2863	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2864	!
2865	!jyg<
2866	! If convective tendencies are too large, then call convection
2867	! every time step
2868	cvpas = cvpas_0
2869	DO k=1,k_upper_cv
2870	DO i=1,klon
2871	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
2872	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
2873	dtcon_multistep_max = 3.
2874	dqcon_multistep_max = 0.02
2875	ENDIF
2876	ENDDO
2877	ENDDO
2878	!
2879	DO k=1,k_upper_cv
2880	DO i=1,klon
2881	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
2882	!! abs(d_q_con(i,k)) > 2.e-2) THEN
2883	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
2884	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
2885	cvpas = 1
2886	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
2887	!! i,k,d_t_con(i,k),d_q_con(i,k)
2888	ENDIF
2889	ENDDO
2890	ENDDO
2891	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
2892	!!! call bcast(cvpas)
2893	!!! ------------------------------------------------------------
2894	!>jyg
2895	!
2896	DO i = 1, klon
2897	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
2898	ENDDO
2899	!
2900	!jyg<
2901	! Add the tendency due to the dry adjustment of the wake profile
2902	IF (iflag_wake>=1) THEN
2903	IF (iflag_adjwk == 2) THEN
2904	DO k=1,klev
2905	DO i=1,klon
2906	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
2907	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
2908	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2909	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2910	ENDDO
2911	ENDDO
2912	ENDIF ! (iflag_adjwk = 2)
2913	ENDIF ! (iflag_wake>=1)
2914	!>jyg
2915	!
2916	ELSE ! ok_cvl
2917
2918	! MAF conema3 ne contient pas les traceurs
2919	CALL conema3 (phys_tstep, &
2920	paprs,pplay,t_seri,q_seri, &
2921	u_seri,v_seri,tr_seri,ntra, &
2922	sig1,w01, &
2923	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2924	rain_con, snow_con, ibas_con, itop_con, &
2925	upwd,dnwd,dnwd0,bas,top, &
2926	Ma,cape,tvp,rflag, &
2927	pbase &
2928	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2929	,clwcon0)
2930
2931	ENDIF ! ok_cvl
2932
2933	!
2934	! Correction precip
2935	rain_con = rain_con * cvl_corr
2936	snow_con = snow_con * cvl_corr
2937	!
2938
2939	IF (.NOT. ok_gust) THEN
2940	do i = 1, klon
2941	wd(i)=0.0
2942	enddo
2943	ENDIF
2944
2945	! =================================================================== c
2946	! Calcul des proprietes des nuages convectifs
2947	!
2948
2949	! calcul des proprietes des nuages convectifs
2950	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2951	IF (iflag_cld_cv == 0) THEN
2952	CALL clouds_gno &
2953	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2954	ELSE
2955	CALL clouds_bigauss &
2956	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2957	ENDIF
2958
2959
2960	! =================================================================== c
2961
2962	DO i = 1, klon
2963	itop_con(i) = min(max(itop_con(i),1),klev)
2964	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2965	ENDDO
2966
2967	DO i = 1, klon
2968	ema_pcb(i) = paprs(i,ibas_con(i))
2969	ENDDO
2970	DO i = 1, klon
2971	! L'idicage de itop_con peut cacher un pb potentiel
2972	! FH sous la dictee de JYG, CR
2973	ema_pct(i) = paprs(i,itop_con(i)+1)
2974
2975	IF (itop_con(i).gt.klev-3) THEN
2976	IF (prt_level >= 9) THEN
2977	write(lunout,*)'La convection monte trop haut '
2978	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2979	ENDIF
2980	ENDIF
2981	ENDDO
2982	ELSE IF (iflag_con.eq.0) THEN
2983	write(lunout,*) 'On n appelle pas la convection'
2984	clwcon0=0.
2985	rnebcon0=0.
2986	d_t_con=0.
2987	d_q_con=0.
2988	d_u_con=0.
2989	d_v_con=0.
2990	rain_con=0.
2991	snow_con=0.
2992	bas=1
2993	top=1
2994	ELSE
2995	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2996	CALL abort_physic("physiq", "", 1)
2997	ENDIF
2998
2999	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3000	! . d_u_con, d_v_con)
3001
3002	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3003	proba_notrig(:) = 1.
3004	itapcv = 0
3005	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3006	!
3007	itapcv = itapcv+1
3008	!
3009	! Compter les steps ou cvpas=1
3010	IF (cvpas == 1) THEN
3011	Ncvpaseq1 = Ncvpaseq1+1
3012	ENDIF
3013	IF (mod(itap,1000) == 0) THEN
3014	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3015	ENDIF
3016
3017	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3018	!!! l'energie dans les courants satures.
3019	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3020	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3021	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3022	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3023	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3024	!! itap, 1)
3025	!! call prt_enerbil('convection_sat',itap)
3026	!!
3027	!!
3028	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
3029	'convection',abortphy,flag_inhib_tend,itap,0)
3030	CALL prt_enerbil('convection',itap)
3031
3032	!-------------------------------------------------------------------------
3033
3034	IF (mydebug) THEN
3035	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3036	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3037	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3038	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3039	ENDIF
3040
3041	IF (check) THEN
3042	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3043	WRITE(lunout,*)"aprescon=", za
3044	zx_t = 0.0
3045	za = 0.0
3046	DO i = 1, klon
3047	za = za + cell_area(i)/REAL(klon)
3048	zx_t = zx_t + (rain_con(i)+ &
3049	snow_con(i))*cell_area(i)/REAL(klon)
3050	ENDDO
3051	zx_t = zx_t/za*phys_tstep
3052	WRITE(lunout,*)"Precip=", zx_t
3053	ENDIF
3054	IF (zx_ajustq) THEN
3055	DO i = 1, klon
3056	z_apres(i) = 0.0
3057	ENDDO
3058	DO k = 1, klev
3059	DO i = 1, klon
3060	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3061	*(paprs(i,k)-paprs(i,k+1))/RG
3062	ENDDO
3063	ENDDO
3064	DO i = 1, klon
3065	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
3066	/z_apres(i)
3067	ENDDO
3068	DO k = 1, klev
3069	DO i = 1, klon
3070	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3071	z_factor(i).LT.(1.0-1.0E-08)) THEN
3072	q_seri(i,k) = q_seri(i,k) * z_factor(i)
3073	ENDIF
3074	ENDDO
3075	ENDDO
3076	ENDIF
3077	zx_ajustq=.FALSE.
3078
3079	!
3080	!==========================================================================
3081	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3082	!pour la couche limite diffuse pour l instant
3083	!
3084	!
3085	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3086	! il faut rajouter cette tendance calcul\'ee hors des poches
3087	! froides
3088	!
3089	IF (iflag_wake>=1) THEN
3090	!
3091	!
3092	! Call wakes every "wkpas" step
3093	!
3094	IF (MOD(itapwk,wkpas).EQ.0) THEN
3095	!
3096	DO k=1,klev
3097	DO i=1,klon
3098	dt_dwn(i,k) = ftd(i,k)
3099	dq_dwn(i,k) = fqd(i,k)
3100	M_dwn(i,k) = dnwd0(i,k)
3101	M_up(i,k) = upwd(i,k)
3102	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3103	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3104	ENDDO
3105	ENDDO
3106
3107	IF (iflag_wake==2) THEN
3108	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3109	DO k = 1,klev
3110	dt_dwn(:,k)= dt_dwn(:,k)+ &
3111	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3112	dq_dwn(:,k)= dq_dwn(:,k)+ &
3113	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3114	ENDDO
3115	ELSEIF (iflag_wake==3) THEN
3116	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3117	DO k = 1,klev
3118	DO i=1,klon
3119	IF (rneb(i,k)==0.) THEN
3120	! On ne tient compte des tendances qu'en dehors des
3121	! nuages (c'est-\`a-dire a priri dans une region ou
3122	! l'eau se reevapore).
3123	dt_dwn(i,k)= dt_dwn(i,k)+ &
3124	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3125	dq_dwn(i,k)= dq_dwn(i,k)+ &
3126	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3127	ENDIF
3128	ENDDO
3129	ENDDO
3130	ENDIF
3131
3132	!
3133	!calcul caracteristiques de la poche froide
3134	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3135	t_seri, q_seri, omega, &
3136	dt_dwn, dq_dwn, M_dwn, M_up, &
3137	dt_a, dq_a, cv_gen, &
3138	sigd, cin, &
3139	wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
3140	wake_dth, wake_h, &
3141	!! wake_pe, wake_fip, wake_gfl, &
3142	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3143	d_t_wake, d_q_wake, &
3144	wake_k, t_x, q_x, &
3145	wake_omgbdth, wake_dp_omgb, &
3146	wake_dtKE, wake_dqKE, &
3147	wake_omg, wake_dp_deltomg, &
3148	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3149	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
3150	!
3151	!jyg Reinitialize itapwk when wakes have been called
3152	itapwk = 0
3153	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3154	!
3155	itapwk = itapwk+1
3156	!
3157	!-----------------------------------------------------------------------
3158	! ajout des tendances des poches froides
3159	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
3160	abortphy,flag_inhib_tend,itap,0)
3161	CALL prt_enerbil('wake',itap)
3162	!------------------------------------------------------------------------
3163
3164	! Increment Wake state variables
3165	IF (iflag_wake_tend .GT. 0.) THEN
3166
3167	CALL add_wake_tend &
3168	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
3169	'wake', abortphy)
3170	CALL prt_enerbil('wake',itap)
3171	ENDIF ! (iflag_wake_tend .GT. 0.)
3172	!
3173	IF (prt_level .GE. 10) THEN
3174	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3175	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3176	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3177	ENDIF
3178
3179	IF (iflag_alp_wk_cond .GT. 0.) THEN
3180
3181	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3182	wake_fip)
3183	ELSE
3184	wake_fip(:) = wake_fip_0(:)
3185	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3186
3187	ENDIF ! (iflag_wake>=1)
3188	!
3189	!===================================================================
3190	! Convection seche (thermiques ou ajustement)
3191	!===================================================================
3192	!
3193	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3194	,seuil_inversion,weak_inversion,dthmin)
3195
3196
3197
3198	d_t_ajsb(:,:)=0.
3199	d_q_ajsb(:,:)=0.
3200	d_t_ajs(:,:)=0.
3201	d_u_ajs(:,:)=0.
3202	d_v_ajs(:,:)=0.
3203	d_q_ajs(:,:)=0.
3204	clwcon0th(:,:)=0.
3205	!
3206	! fm_therm(:,:)=0.
3207	! entr_therm(:,:)=0.
3208	! detr_therm(:,:)=0.
3209	!
3210	IF (prt_level>9) WRITE(lunout,*) &
3211	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3212	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3213	IF (iflag_thermals<0) THEN
3214	! Rien
3215	! ====
3216	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3217
3218
3219	ELSE
3220
3221	! Thermiques
3222	! ==========
3223	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3224	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3225
3226
3227	!cc nrlmd le 10/04/2012
3228	DO k=1,klev+1
3229	DO i=1,klon
3230	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3231	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3232	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3233	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3234	ENDDO
3235	ENDDO
3236	!cc fin nrlmd le 10/04/2012
3237
3238	IF (iflag_thermals>=1) THEN
3239	!jyg<
3240	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3241	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3242	! Appel des thermiques avec les profils exterieurs aux poches
3243	DO k=1,klev
3244	DO i=1,klon
3245	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3246	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3247	u_therm(i,k) = u_seri(i,k)
3248	v_therm(i,k) = v_seri(i,k)
3249	ENDDO
3250	ENDDO
3251	ELSE
3252	! Appel des thermiques avec les profils moyens
3253	DO k=1,klev
3254	DO i=1,klon
3255	t_therm(i,k) = t_seri(i,k)
3256	q_therm(i,k) = q_seri(i,k)
3257	u_therm(i,k) = u_seri(i,k)
3258	v_therm(i,k) = v_seri(i,k)
3259	ENDDO
3260	ENDDO
3261	ENDIF
3262	!>jyg
3263	CALL calltherm(pdtphys &
3264	,pplay,paprs,pphi,weak_inversion &
3265	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3266	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
3267	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3268	,fm_therm,entr_therm,detr_therm &
3269	,zqasc,clwcon0th,lmax_th,ratqscth &
3270	,ratqsdiff,zqsatth &
3271	!on rajoute ale et alp, et les
3272	!caracteristiques de la couche alim
3273	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3274	,ztv,zpspsk,ztla,zthl &
3275	!cc nrlmd le 10/04/2012
3276	,pbl_tke_input,pctsrf,omega,cell_area &
3277	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3278	,n2,s2,ale_bl_stat &
3279	,therm_tke_max,env_tke_max &
3280	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3281	,alp_bl_conv,alp_bl_stat &
3282	!cc fin nrlmd le 10/04/2012
3283	,zqla,ztva )
3284	!
3285	!jyg<
3286	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3287	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3288	! Si les thermiques ne sont presents que hors des
3289	! poches, la tendance moyenne associ\'ee doit etre
3290	! multipliee par la fraction surfacique qu'ils couvrent.
3291	DO k=1,klev
3292	DO i=1,klon
3293	!
3294	d_deltat_the(i,k) = - d_t_ajs(i,k)
3295	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3296	!
3297	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3298	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3299	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3300	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3301	!
3302	ENDDO
3303	ENDDO
3304	!
3305	IF (ok_bug_split_th) THEN
3306	CALL add_wake_tend &
3307	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3308	ELSE
3309	CALL add_wake_tend &
3310	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3311	ENDIF
3312	CALL prt_enerbil('the',itap)
3313	!
3314	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3315	!
3316	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3317	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3318	CALL prt_enerbil('thermals',itap)
3319	!
3320	!
3321	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3322	cin, s2, n2, &
3323	ale_bl_trig, ale_bl_stat, ale_bl, &
3324	alp_bl, alp_bl_stat, &
3325	proba_notrig, random_notrig, cv_gen)
3326	!>jyg
3327
3328	! ------------------------------------------------------------------
3329	! Transport de la TKE par les panaches thermiques.
3330	! FH : 2010/02/01
3331	! if (iflag_pbl.eq.10) then
3332	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3333	! s rg,paprs,pbl_tke)
3334	! endif
3335	! -------------------------------------------------------------------
3336
3337	DO i=1,klon
3338	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3339	!CR:04/05/12:correction calcul zmax
3340	zmax_th(i)=zmax0(i)
3341	ENDDO
3342
3343	ENDIF
3344
3345	! Ajustement sec
3346	! ==============
3347
3348	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3349	! a partir du sommet des thermiques.
3350	! Dans le cas contraire, on demarre au niveau 1.
3351
3352	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3353
3354	IF (iflag_thermals.eq.0) THEN
3355	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3356	limbas(:)=1
3357	ELSE
3358	limbas(:)=lmax_th(:)
3359	ENDIF
3360
3361	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3362	! pour des test de convergence numerique.
3363	! Le nouveau ajsec est a priori mieux, meme pour le cas
3364	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3365	! non nulles numeriquement pour des mailles non concernees.
3366
3367	IF (iflag_thermals==0) THEN
3368	! Calling adjustment alone (but not the thermal plume model)
3369	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3370	, d_t_ajsb, d_q_ajsb)
3371	ELSE IF (iflag_thermals>0) THEN
3372	! Calling adjustment above the top of thermal plumes
3373	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3374	, d_t_ajsb, d_q_ajsb)
3375	ENDIF
3376
3377	!--------------------------------------------------------------------
3378	! ajout des tendances de l'ajustement sec ou des thermiques
3379	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3380	'ajsb',abortphy,flag_inhib_tend,itap,0)
3381	CALL prt_enerbil('ajsb',itap)
3382	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3383	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3384
3385	!---------------------------------------------------------------------
3386
3387	ENDIF
3388
3389	ENDIF
3390	!
3391	!===================================================================
3392	! Computation of ratqs, the width (normalized) of the subrid scale
3393	! water distribution
3394	CALL calcratqs(klon,klev,prt_level,lunout, &
3395	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3396	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3397	tau_ratqs,fact_cldcon, &
3398	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3399	paprs,pplay,q_seri,zqsat,fm_therm, &
3400	ratqs,ratqsc)
3401
3402
3403	!
3404	! Appeler le processus de condensation a grande echelle
3405	! et le processus de precipitation
3406	!-------------------------------------------------------------------------
3407	IF (prt_level .GE.10) THEN
3408	print *,'itap, ->fisrtilp ',itap
3409	ENDIF
3410	!
3411	CALL fisrtilp(phys_tstep,paprs,pplay, &
3412	t_seri, q_seri,ptconv,ratqs, &
3413	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3414	rain_lsc, snow_lsc, &
3415	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3416	frac_impa, frac_nucl, beta_prec_fisrt, &
3417	prfl, psfl, rhcl, &
3418	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3419	iflag_ice_thermo)
3420	!
3421	WHERE (rain_lsc < 0) rain_lsc = 0.
3422	WHERE (snow_lsc < 0) snow_lsc = 0.
3423
3424	!+JLD
3425	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3426	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3427	! & ,rain_lsc,snow_lsc
3428	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3429	!-JLD
3430	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3431	'lsc',abortphy,flag_inhib_tend,itap,0)
3432	CALL prt_enerbil('lsc',itap)
3433	rain_num(:)=0.
3434	DO k = 1, klev
3435	DO i = 1, klon
3436	IF (ql_seri(i,k)>oliqmax) THEN
3437	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3438	ql_seri(i,k)=oliqmax
3439	ENDIF
3440	ENDDO
3441	ENDDO
3442	IF (nqo==3) THEN
3443	DO k = 1, klev
3444	DO i = 1, klon
3445	IF (qs_seri(i,k)>oicemax) THEN
3446	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3447	qs_seri(i,k)=oicemax
3448	ENDIF
3449	ENDDO
3450	ENDDO
3451	ENDIF
3452
3453	!---------------------------------------------------------------------------
3454	DO k = 1, klev
3455	DO i = 1, klon
3456	cldfra(i,k) = rneb(i,k)
3457	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3458	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3459	ENDDO
3460	ENDDO
3461	IF (check) THEN
3462	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3463	WRITE(lunout,*)"apresilp=", za
3464	zx_t = 0.0
3465	za = 0.0
3466	DO i = 1, klon
3467	za = za + cell_area(i)/REAL(klon)
3468	zx_t = zx_t + (rain_lsc(i) &
3469	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3470	ENDDO
3471	zx_t = zx_t/za*phys_tstep
3472	WRITE(lunout,*)"Precip=", zx_t
3473	ENDIF
3474
3475	IF (mydebug) THEN
3476	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3477	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3478	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3479	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3480	ENDIF
3481
3482	!
3483	!-------------------------------------------------------------------
3484	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3485	!-------------------------------------------------------------------
3486
3487	! 1. NUAGES CONVECTIFS
3488	!
3489	!IM cf FH
3490	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3491	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3492	snow_tiedtke=0.
3493	! print*,'avant calcul de la pseudo precip '
3494	! print*,'iflag_cld_th',iflag_cld_th
3495	IF (iflag_cld_th.eq.-1) THEN
3496	rain_tiedtke=rain_con
3497	ELSE
3498	! print*,'calcul de la pseudo precip '
3499	rain_tiedtke=0.
3500	! print*,'calcul de la pseudo precip 0'
3501	DO k=1,klev
3502	DO i=1,klon
3503	IF (d_q_con(i,k).lt.0.) THEN
3504	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3505	*(paprs(i,k)-paprs(i,k+1))/rg
3506	ENDIF
3507	ENDDO
3508	ENDDO
3509	ENDIF
3510	!
3511	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3512	!
3513
3514	! Nuages diagnostiques pour Tiedtke
3515	CALL diagcld1(paprs,pplay, &
3516	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3517	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3518	diafra,dialiq)
3519	DO k = 1, klev
3520	DO i = 1, klon
3521	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3522	cldliq(i,k) = dialiq(i,k)
3523	cldfra(i,k) = diafra(i,k)
3524	ENDIF
3525	ENDDO
3526	ENDDO
3527
3528	ELSE IF (iflag_cld_th.ge.3) THEN
3529	! On prend pour les nuages convectifs le max du calcul de la
3530	! convection et du calcul du pas de temps precedent diminue d'un facteur
3531	! facttemps
3532	facteur = pdtphys *facttemps
3533	DO k=1,klev
3534	DO i=1,klon
3535	rnebcon(i,k)=rnebcon(i,k)*facteur
3536	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3537	rnebcon(i,k)=rnebcon0(i,k)
3538	clwcon(i,k)=clwcon0(i,k)
3539	ENDIF
3540	ENDDO
3541	ENDDO
3542
3543	! On prend la somme des fractions nuageuses et des contenus en eau
3544
3545	IF (iflag_cld_th>=5) THEN
3546
3547	DO k=1,klev
3548	ptconvth(:,k)=fm_therm(:,k+1)>0.
3549	ENDDO
3550
3551	IF (iflag_coupl==4) THEN
3552
3553	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3554	! convectives et lsc dans la partie des thermiques
3555	! Le controle par iflag_coupl est peut etre provisoire.
3556	DO k=1,klev
3557	DO i=1,klon
3558	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3559	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3560	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3561	ELSE IF (ptconv(i,k)) THEN
3562	cldfra(i,k)=rnebcon(i,k)
3563	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3564	ENDIF
3565	ENDDO
3566	ENDDO
3567
3568	ELSE IF (iflag_coupl==5) THEN
3569	DO k=1,klev
3570	DO i=1,klon
3571	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3572	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3573	ENDDO
3574	ENDDO
3575
3576	ELSE
3577
3578	! Si on est sur un point touche par la convection
3579	! profonde et pas par les thermiques, on prend la
3580	! couverture nuageuse et l'eau nuageuse de la convection
3581	! profonde.
3582
3583	!IM/FH: 2011/02/23
3584	! definition des points sur lesquels ls thermiques sont actifs
3585
3586	DO k=1,klev
3587	DO i=1,klon
3588	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3589	cldfra(i,k)=rnebcon(i,k)
3590	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3591	ENDIF
3592	ENDDO
3593	ENDDO
3594
3595	ENDIF
3596
3597	ELSE
3598
3599	! Ancienne version
3600	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3601	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3602	ENDIF
3603
3604	ENDIF
3605
3606	! plulsc(:)=0.
3607	! do k=1,klev,-1
3608	! do i=1,klon
3609	! zzz=prfl(:,k)+psfl(:,k)
3610	! if (.not.ptconvth.zzz.gt.0.)
3611	! enddo prfl, psfl,
3612	! enddo
3613	!
3614	! 2. NUAGES STARTIFORMES
3615	!
3616	IF (ok_stratus) THEN
3617	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3618	DO k = 1, klev
3619	DO i = 1, klon
3620	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3621	cldliq(i,k) = dialiq(i,k)
3622	cldfra(i,k) = diafra(i,k)
3623	ENDIF
3624	ENDDO
3625	ENDDO
3626	ENDIF
3627	!
3628	! Precipitation totale
3629	!
3630	DO i = 1, klon
3631	rain_fall(i) = rain_con(i) + rain_lsc(i)
3632	snow_fall(i) = snow_con(i) + snow_lsc(i)
3633	ENDDO
3634	!
3635	! Calculer l'humidite relative pour diagnostique
3636	!
3637	DO k = 1, klev
3638	DO i = 1, klon
3639	zx_t = t_seri(i,k)
3640	IF (thermcep) THEN
3641	!! if (iflag_ice_thermo.eq.0) then !jyg
3642	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3643	!! else !jyg
3644	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3645	!! endif !jyg
3646	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3647	zx_qs = MIN(0.5,zx_qs)
3648	zcor = 1./(1.-retv*zx_qs)
3649	zx_qs = zx_qs*zcor
3650	ELSE
3651	!! IF (zx_t.LT.t_coup) THEN !jyg
3652	IF (zx_t.LT.rtt) THEN !jyg
3653	zx_qs = qsats(zx_t)/pplay(i,k)
3654	ELSE
3655	zx_qs = qsatl(zx_t)/pplay(i,k)
3656	ENDIF
3657	ENDIF
3658	zx_rh(i,k) = q_seri(i,k)/zx_qs
3659	zqsat(i,k)=zx_qs
3660	ENDDO
3661	ENDDO
3662
3663	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3664	! equivalente a 2m (tpote) pour diagnostique
3665	!
3666	DO i = 1, klon
3667	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3668	IF (thermcep) THEN
3669	IF(zt2m(i).LT.RTT) then
3670	Lheat=RLSTT
3671	ELSE
3672	Lheat=RLVTT
3673	ENDIF
3674	ELSE
3675	IF (zt2m(i).LT.RTT) THEN
3676	Lheat=RLSTT
3677	ELSE
3678	Lheat=RLVTT
3679	ENDIF
3680	ENDIF
3681	tpote(i) = tpot(i)* &
3682	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3683	ENDDO
3684
3685	IF (type_trac == 'inca') THEN
3686	#ifdef INCA
3687	CALL VTe(VTphysiq)
3688	CALL VTb(VTinca)
3689	calday = REAL(days_elapsed + 1) + jH_cur
3690
3691	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
3692	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3693	CALL AEROSOL_METEO_CALC( &
3694	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3695	prfl,psfl,pctsrf,cell_area, &
3696	latitude_deg,longitude_deg,u10m,v10m)
3697	ENDIF
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	ENDIF
3883	!
3884	#ifdef CPP_RRTM
3885	#ifdef CPP_StratAer
3886	!--compute stratospheric mask
3887	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3888	!--interactive strat aerosols
3889	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3890	#endif
3891	#endif
3892	!--fin STRAT AEROSOL
3893	!
3894
3895	! Calculer les parametres optiques des nuages et quelques
3896	! parametres pour diagnostiques:
3897	!
3898	IF (aerosol_couple.AND.config_inca=='aero') THEN
3899	mass_solu_aero(:,:) = ccm(:,:,1)
3900	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3901	ENDIF
3902
3903	IF (ok_newmicro) then
3904	IF (iflag_rrtm.NE.0) THEN
3905	#ifdef CPP_RRTM
3906	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3907	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3908	// 'pour ok_cdnc'
3909	CALL abort_physic(modname,abort_message,1)
3910	ENDIF
3911	#else
3912
3913	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3914	CALL abort_physic(modname,abort_message,1)
3915	#endif
3916	ENDIF
3917	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
3918	paprs, pplay, t_seri, cldliq, cldfra, &
3919	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3920	flwp, fiwp, flwc, fiwc, &
3921	mass_solu_aero, mass_solu_aero_pi, &
3922	cldtaupi, re, fl, ref_liq, ref_ice, &
3923	ref_liq_pi, ref_ice_pi)
3924	ELSE
3925	CALL nuage (paprs, pplay, &
3926	t_seri, cldliq, cldfra, cldtau, cldemi, &
3927	cldh, cldl, cldm, cldt, cldq, &
3928	ok_aie, &
3929	mass_solu_aero, mass_solu_aero_pi, &
3930	bl95_b0, bl95_b1, &
3931	cldtaupi, re, fl)
3932	ENDIF
3933	!
3934	!IM betaCRF
3935	!
3936	cldtaurad = cldtau
3937	cldtaupirad = cldtaupi
3938	cldemirad = cldemi
3939	cldfrarad = cldfra
3940
3941	!
3942	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3943	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3944	!
3945	! global
3946	!
3947	!IM 251017 begin
3948	! print*,'physiq betaCRF global zdtime=',zdtime
3949	!IM 251017 end
3950	DO k=1, klev
3951	DO i=1, klon
3952	IF (pplay(i,k).GE.pfree) THEN
3953	beta(i,k) = beta_pbl
3954	ELSE
3955	beta(i,k) = beta_free
3956	ENDIF
3957	IF (mskocean_beta) THEN
3958	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3959	ENDIF
3960	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3961	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3962	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3963	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3964	ENDDO
3965	ENDDO
3966	!
3967	ELSE
3968	!
3969	! regional
3970	!
3971	DO k=1, klev
3972	DO i=1,klon
3973	!
3974	IF (longitude_deg(i).ge.lon1_beta.AND. &
3975	longitude_deg(i).le.lon2_beta.AND. &
3976	latitude_deg(i).le.lat1_beta.AND. &
3977	latitude_deg(i).ge.lat2_beta) THEN
3978	IF (pplay(i,k).GE.pfree) THEN
3979	beta(i,k) = beta_pbl
3980	ELSE
3981	beta(i,k) = beta_free
3982	ENDIF
3983	IF (mskocean_beta) THEN
3984	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3985	ENDIF
3986	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3987	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3988	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3989	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3990	ENDIF
3991	!
3992	ENDDO
3993	ENDDO
3994	!
3995	ENDIF
3996
3997	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3998	IF (ok_chlorophyll) THEN
3999	print*,"-- reading chlorophyll"
4000	CALL readchlorophyll(debut)
4001	ENDIF
4002
4003	!--if ok_suntime_rrtm we use ancillay data for RSUN
4004	!--previous values are therefore overwritten
4005	!--this is needed for CMIP6 runs
4006	!--and only possible for new radiation scheme
4007	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4008	#ifdef CPP_RRTM
4009	CALL read_rsun_rrtm(debut)
4010	#endif
4011	ENDIF
4012
4013	IF (mydebug) THEN
4014	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4015	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4016	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4017	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4018	ENDIF
4019
4020	!
4021	!sonia : If Iflag_radia >=2, pertubation of some variables
4022	!input to radiation (DICE)
4023	!
4024	IF (iflag_radia .ge. 2) THEN
4025	zsav_tsol (:) = zxtsol(:)
4026	CALL perturb_radlwsw(zxtsol,iflag_radia)
4027	ENDIF
4028
4029	IF (aerosol_couple.AND.config_inca=='aero') THEN
4030	#ifdef INCA
4031	CALL radlwsw_inca &
4032	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4033	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4034	size(wo,3), wo, &
4035	cldfrarad, cldemirad, cldtaurad, &
4036	heat,heat0,cool,cool0,albpla, &
4037	topsw,toplw,solsw,sollw, &
4038	sollwdown, &
4039	topsw0,toplw0,solsw0,sollw0, &
4040	lwdn0, lwdn, lwup0, lwup, &
4041	swdn0, swdn, swup0, swup, &
4042	ok_ade, ok_aie, &
4043	tau_aero, piz_aero, cg_aero, &
4044	topswad_aero, solswad_aero, &
4045	topswad0_aero, solswad0_aero, &
4046	topsw_aero, topsw0_aero, &
4047	solsw_aero, solsw0_aero, &
4048	cldtaupirad, &
4049	topswai_aero, solswai_aero)
4050	#endif
4051	ELSE
4052	!
4053	!IM calcul radiatif pour le cas actuel
4054	!
4055	RCO2 = RCO2_act
4056	RCH4 = RCH4_act
4057	RN2O = RN2O_act
4058	RCFC11 = RCFC11_act
4059	RCFC12 = RCFC12_act
4060	!
4061	!--interactive CO2 in ppm from carbon cycle
4062	IF (carbon_cycle_rad.AND..NOT.debut) THEN
4063	RCO2=RCO2_glo
4064	ENDIF
4065	!
4066	IF (prt_level .GE.10) THEN
4067	print *,' ->radlwsw, number 1 '
4068	ENDIF
4069	!
4070	CALL radlwsw &
4071	(dist, rmu0, fract, &
4072	!albedo SB >>>
4073	! paprs, pplay,zxtsol,albsol1, albsol2, &
4074	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4075	!albedo SB <<<
4076	t_seri,q_seri,wo, &
4077	cldfrarad, cldemirad, cldtaurad, &
4078	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4079	flag_aerosol, &
4080	flag_aerosol_strat, flag_aer_feedback, &
4081	tau_aero, piz_aero, cg_aero, &
4082	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4083	! Rajoute par OB pour RRTM
4084	tau_aero_lw_rrtm, &
4085	cldtaupirad,new_aod, &
4086	! zqsat, flwcrad, fiwcrad, &
4087	zqsat, flwc, fiwc, &
4088	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4089	heat,heat0,cool,cool0,albpla, &
4090	heat_volc,cool_volc, &
4091	topsw,toplw,solsw,sollw, &
4092	sollwdown, &
4093	topsw0,toplw0,solsw0,sollw0, &
4094	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4095	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4096	topswad_aero, solswad_aero, &
4097	topswai_aero, solswai_aero, &
4098	topswad0_aero, solswad0_aero, &
4099	topsw_aero, topsw0_aero, &
4100	solsw_aero, solsw0_aero, &
4101	topswcf_aero, solswcf_aero, &
4102	!-C. Kleinschmitt for LW diagnostics
4103	toplwad_aero, sollwad_aero,&
4104	toplwai_aero, sollwai_aero, &
4105	toplwad0_aero, sollwad0_aero,&
4106	!-end
4107	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4108	ZSWFT0_i, ZFSDN0, ZFSUP0)
4109
4110	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4111	!schemes
4112	toplw = toplw + betalwoff * (toplw0 - toplw)
4113	sollw = sollw + betalwoff * (sollw0 - sollw)
4114	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4115	lwup = lwup + betalwoff * (lwup0 - lwup)
4116	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4117	sollwdown(:))
4118	cool = cool + betalwoff * (cool0 - cool)
4119
4120	#ifndef CPP_XIOS
4121	!--OB 30/05/2016 modified 21/10/2016
4122	!--here we return swaero_diag and dryaod_diag to FALSE
4123	!--and histdef will switch it back to TRUE if necessary
4124	!--this is necessary to get the right swaero at first step
4125	!--but only in the case of no XIOS as XIOS is covered elsewhere
4126	IF (debut) swaerofree_diag = .FALSE.
4127	IF (debut) swaero_diag = .FALSE.
4128	IF (debut) dryaod_diag = .FALSE.
4129	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4130	!--as for swaero_diag, see above
4131	IF (debut) ok_4xCO2atm = .FALSE.
4132
4133	!
4134	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4135	!IM des taux doit etre different du taux actuel
4136	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4137	!
4138	IF (RCO2_per.NE.RCO2_act.OR. &
4139	RCH4_per.NE.RCH4_act.OR. &
4140	RN2O_per.NE.RN2O_act.OR. &
4141	RCFC11_per.NE.RCFC11_act.OR. &
4142	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4143	#endif
4144	!
4145	IF (ok_4xCO2atm) THEN
4146	!
4147	RCO2 = RCO2_per
4148	RCH4 = RCH4_per
4149	RN2O = RN2O_per
4150	RCFC11 = RCFC11_per
4151	RCFC12 = RCFC12_per
4152	!
4153	IF (prt_level .GE.10) THEN
4154	print *,' ->radlwsw, number 2 '
4155	ENDIF
4156	!
4157	CALL radlwsw &
4158	(dist, rmu0, fract, &
4159	!albedo SB >>>
4160	! paprs, pplay,zxtsol,albsol1, albsol2, &
4161	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4162	!albedo SB <<<
4163	t_seri,q_seri,wo, &
4164	cldfrarad, cldemirad, cldtaurad, &
4165	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4166	flag_aerosol, &
4167	flag_aerosol_strat, flag_aer_feedback, &
4168	tau_aero, piz_aero, cg_aero, &
4169	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4170	! Rajoute par OB pour RRTM
4171	tau_aero_lw_rrtm, &
4172	cldtaupi,new_aod, &
4173	! zqsat, flwcrad, fiwcrad, &
4174	zqsat, flwc, fiwc, &
4175	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4176	heatp,heat0p,coolp,cool0p,albplap, &
4177	heat_volc,cool_volc, &
4178	topswp,toplwp,solswp,sollwp, &
4179	sollwdownp, &
4180	topsw0p,toplw0p,solsw0p,sollw0p, &
4181	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4182	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4183	topswad_aerop, solswad_aerop, &
4184	topswai_aerop, solswai_aerop, &
4185	topswad0_aerop, solswad0_aerop, &
4186	topsw_aerop, topsw0_aerop, &
4187	solsw_aerop, solsw0_aerop, &
4188	topswcf_aerop, solswcf_aerop, &
4189	!-C. Kleinschmitt for LW diagnostics
4190	toplwad_aerop, sollwad_aerop,&
4191	toplwai_aerop, sollwai_aerop, &
4192	toplwad0_aerop, sollwad0_aerop,&
4193	!-end
4194	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4195	ZSWFT0_i, ZFSDN0, ZFSUP0)
4196	endif !ok_4xCO2atm
4197	ENDIF ! aerosol_couple
4198	itaprad = 0
4199	!
4200	! If Iflag_radia >=2, reset pertubed variables
4201	!
4202	IF (iflag_radia .ge. 2) THEN
4203	zxtsol(:) = zsav_tsol (:)
4204	ENDIF
4205	ENDIF ! MOD(itaprad,radpas)
4206	itaprad = itaprad + 1
4207
4208	IF (iflag_radia.eq.0) THEN
4209	IF (prt_level.ge.9) THEN
4210	PRINT *,'--------------------------------------------------'
4211	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4212	PRINT *,'>>>> heat et cool mis a zero '
4213	PRINT *,'--------------------------------------------------'
4214	ENDIF
4215	heat=0.
4216	cool=0.
4217	sollw=0. ! MPL 01032011
4218	solsw=0.
4219	radsol=0.
4220	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4221	swup0=0.
4222	lwup=0.
4223	lwup0=0.
4224	lwdn=0.
4225	lwdn0=0.
4226	ENDIF
4227
4228	!
4229	! Calculer radsol a l'exterieur de radlwsw
4230	! pour prendre en compte le cycle diurne
4231	! recode par Olivier Boucher en sept 2015
4232	!
4233	radsol=solsw*swradcorr+sollw
4234
4235	IF (ok_4xCO2atm) THEN
4236	radsolp=solswp*swradcorr+sollwp
4237	ENDIF
4238
4239	!
4240	! Ajouter la tendance des rayonnements (tous les pas)
4241	! avec une correction pour le cycle diurne dans le SW
4242	!
4243
4244	DO k=1, klev
4245	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4246	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4247	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4248	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4249	ENDDO
4250
4251	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4252	CALL prt_enerbil('SW',itap)
4253	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4254	CALL prt_enerbil('LW',itap)
4255
4256	!
4257	IF (mydebug) THEN
4258	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4259	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4260	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4261	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4262	ENDIF
4263
4264	! Calculer l'hydrologie de la surface
4265	!
4266	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4267	! . agesno, ftsol,fqsurf,fsnow, ruis)
4268	!
4269
4270	!
4271	! Calculer le bilan du sol et la derive de temperature (couplage)
4272	!
4273	DO i = 1, klon
4274	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4275	! a la demande de JLD
4276	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4277	ENDDO
4278	!
4279	!moddeblott(jan95)
4280	! Appeler le programme de parametrisation de l'orographie
4281	! a l'echelle sous-maille:
4282	!
4283	IF (prt_level .GE.10) THEN
4284	print *,' call orography ? ', ok_orodr
4285	ENDIF
4286	!
4287	IF (ok_orodr) THEN
4288	!
4289	! selection des points pour lesquels le shema est actif:
4290	igwd=0
4291	DO i=1,klon
4292	itest(i)=0
4293	! IF ((zstd(i).gt.10.0)) THEN
4294	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4295	itest(i)=1
4296	igwd=igwd+1
4297	idx(igwd)=i
4298	ENDIF
4299	ENDDO
4300	! igwdim=MAX(1,igwd)
4301	!
4302	IF (ok_strato) THEN
4303
4304	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4305	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4306	igwd,idx,itest, &
4307	t_seri, u_seri, v_seri, &
4308	zulow, zvlow, zustrdr, zvstrdr, &
4309	d_t_oro, d_u_oro, d_v_oro)
4310
4311	ELSE
4312	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4313	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4314	igwd,idx,itest, &
4315	t_seri, u_seri, v_seri, &
4316	zulow, zvlow, zustrdr, zvstrdr, &
4317	d_t_oro, d_u_oro, d_v_oro)
4318	ENDIF
4319	!
4320	! ajout des tendances
4321	!-----------------------------------------------------------------------
4322	! ajout des tendances de la trainee de l'orographie
4323	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4324	abortphy,flag_inhib_tend,itap,0)
4325	CALL prt_enerbil('oro',itap)
4326	!----------------------------------------------------------------------
4327	!
4328	ENDIF ! fin de test sur ok_orodr
4329	!
4330	IF (mydebug) THEN
4331	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4332	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4333	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4334	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4335	ENDIF
4336
4337	IF (ok_orolf) THEN
4338	!
4339	! selection des points pour lesquels le shema est actif:
4340	igwd=0
4341	DO i=1,klon
4342	itest(i)=0
4343	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4344	itest(i)=1
4345	igwd=igwd+1
4346	idx(igwd)=i
4347	ENDIF
4348	ENDDO
4349	! igwdim=MAX(1,igwd)
4350	!
4351	IF (ok_strato) THEN
4352
4353	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4354	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4355	igwd,idx,itest, &
4356	t_seri, u_seri, v_seri, &
4357	zulow, zvlow, zustrli, zvstrli, &
4358	d_t_lif, d_u_lif, d_v_lif )
4359
4360	ELSE
4361	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4362	latitude_deg,zmea,zstd,zpic, &
4363	itest, &
4364	t_seri, u_seri, v_seri, &
4365	zulow, zvlow, zustrli, zvstrli, &
4366	d_t_lif, d_u_lif, d_v_lif)
4367	ENDIF
4368
4369	! ajout des tendances de la portance de l'orographie
4370	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4371	'lif', abortphy,flag_inhib_tend,itap,0)
4372	CALL prt_enerbil('lif',itap)
4373	ENDIF ! fin de test sur ok_orolf
4374
4375	IF (ok_hines) then
4376	! HINES GWD PARAMETRIZATION
4377	east_gwstress=0.
4378	west_gwstress=0.
4379	du_gwd_hines=0.
4380	dv_gwd_hines=0.
4381	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4382	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4383	du_gwd_hines, dv_gwd_hines)
4384	zustr_gwd_hines=0.
4385	zvstr_gwd_hines=0.
4386	DO k = 1, klev
4387	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4388	* (paprs(:, k)-paprs(:, k+1))/rg
4389	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4390	* (paprs(:, k)-paprs(:, k+1))/rg
4391	ENDDO
4392
4393	d_t_hin(:, :)=0.
4394	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4395	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4396	CALL prt_enerbil('hin',itap)
4397	ENDIF
4398
4399	IF (.not. ok_hines .and. ok_gwd_rando) then
4400	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4401	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4402	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4403	dv_gwd_front, east_gwstress, west_gwstress)
4404	zustr_gwd_front=0.
4405	zvstr_gwd_front=0.
4406	DO k = 1, klev
4407	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4408	* (paprs(:, k)-paprs(:, k+1))/rg
4409	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4410	* (paprs(:, k)-paprs(:, k+1))/rg
4411	ENDDO
4412
4413	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4414	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4415	CALL prt_enerbil('front_gwd_rando',itap)
4416	ENDIF
4417
4418	IF (ok_gwd_rando) THEN
4419	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4420	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4421	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4422	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4423	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4424	CALL prt_enerbil('flott_gwd_rando',itap)
4425	zustr_gwd_rando=0.
4426	zvstr_gwd_rando=0.
4427	DO k = 1, klev
4428	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4429	* (paprs(:, k)-paprs(:, k+1))/rg
4430	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4431	* (paprs(:, k)-paprs(:, k+1))/rg
4432	ENDDO
4433	ENDIF
4434
4435	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4436
4437	IF (mydebug) THEN
4438	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4439	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4440	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4441	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4442	ENDIF
4443
4444	DO i = 1, klon
4445	zustrph(i)=0.
4446	zvstrph(i)=0.
4447	ENDDO
4448	DO k = 1, klev
4449	DO i = 1, klon
4450	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4451	(paprs(i,k)-paprs(i,k+1))/rg
4452	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4453	(paprs(i,k)-paprs(i,k+1))/rg
4454	ENDDO
4455	ENDDO
4456	!
4457	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4458	!
4459	IF (is_sequential .and. ok_orodr) THEN
4460	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4461	ra,rg,romega, &
4462	latitude_deg,longitude_deg,pphis, &
4463	zustrdr,zustrli,zustrph, &
4464	zvstrdr,zvstrli,zvstrph, &
4465	paprs,u,v, &
4466	aam, torsfc)
4467	ENDIF
4468	!IM cf. FLott END
4469	!DC Calcul de la tendance due au methane
4470	IF (ok_qch4) THEN
4471	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4472	! ajout de la tendance d'humidite due au methane
4473	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4474	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4475	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4476	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4477	ENDIF
4478	!
4479	!
4480
4481	!===============================================================
4482	! Additional tendency of TKE due to orography
4483	!===============================================================
4484	!
4485	! Inititialization
4486	!------------------
4487
4488	addtkeoro=0
4489	CALL getin_p('addtkeoro',addtkeoro)
4490
4491	IF (prt_level.ge.5) &
4492	print*,'addtkeoro', addtkeoro
4493
4494	alphatkeoro=1.
4495	CALL getin_p('alphatkeoro',alphatkeoro)
4496	alphatkeoro=min(max(0.,alphatkeoro),1.)
4497
4498	smallscales_tkeoro=.FALSE.
4499	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4500
4501
4502	dtadd(:,:)=0.
4503	duadd(:,:)=0.
4504	dvadd(:,:)=0.
4505
4506	! Choices for addtkeoro:
4507	! ** 0 no TKE tendency from orography
4508	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4509	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4510	!
4511
4512	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4513	! -------------------------------------------
4514
4515
4516	! selection des points pour lesquels le schema est actif:
4517
4518
4519	IF (addtkeoro .EQ. 1 ) THEN
4520
4521	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4522	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4523
4524	ELSE IF (addtkeoro .EQ. 2) THEN
4525
4526	IF (smallscales_tkeoro) THEN
4527	igwd=0
4528	DO i=1,klon
4529	itest(i)=0
4530	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4531	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4532	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4533	IF (zstd(i).GT.1.0) THEN
4534	itest(i)=1
4535	igwd=igwd+1
4536	idx(igwd)=i
4537	ENDIF
4538	ENDDO
4539
4540	ELSE
4541
4542	igwd=0
4543	DO i=1,klon
4544	itest(i)=0
4545	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4546	itest(i)=1
4547	igwd=igwd+1
4548	idx(igwd)=i
4549	ENDIF
4550	ENDDO
4551
4552	ENDIF
4553
4554	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4555	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4556	igwd,idx,itest, &
4557	t_seri, u_seri, v_seri, &
4558	zulow, zvlow, zustrdr, zvstrdr, &
4559	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4560
4561	zustrdr(:)=0.
4562	zvstrdr(:)=0.
4563	zulow(:)=0.
4564	zvlow(:)=0.
4565
4566	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4567	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4568	ENDIF
4569
4570
4571	! TKE update from subgrid temperature and wind tendencies
4572	!----------------------------------------------------------
4573	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4574
4575
4576	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4577
4578
4579	ENDIF
4580	! -----
4581	!===============================================================
4582
4583
4584	!====================================================================
4585	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4586	!====================================================================
4587	! Abderrahmane 24.08.09
4588
4589	IF (ok_cosp) THEN
4590	! adeclarer
4591	#ifdef CPP_COSP
4592	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4593
4594	IF (prt_level .GE.10) THEN
4595	print*,'freq_cosp',freq_cosp
4596	ENDIF
4597	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4598	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4599	! s ref_liq,ref_ice
4600	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4601	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4602	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4603	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4604	JrNt,ref_liq,ref_ice, &
4605	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4606	zu10m,zv10m,pphis, &
4607	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4608	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4609	prfl(:,1:klev),psfl(:,1:klev), &
4610	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4611	mr_ozone,cldtau, cldemi)
4612
4613	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4614	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4615	! M clMISR,
4616	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4617	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4618
4619	ENDIF
4620	#endif
4621
4622	#ifdef CPP_COSP2
4623	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4624
4625	IF (prt_level .GE.10) THEN
4626	print*,'freq_cosp',freq_cosp
4627	ENDIF
4628	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4629	print*,'Dans physiq.F avant appel '
4630	! s ref_liq,ref_ice
4631	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4632	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4633	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4634	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4635	JrNt,ref_liq,ref_ice, &
4636	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4637	zu10m,zv10m,pphis, &
4638	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4639	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4640	prfl(:,1:klev),psfl(:,1:klev), &
4641	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4642	mr_ozone,cldtau, cldemi)
4643	ENDIF
4644	#endif
4645
4646	#ifdef CPP_COSPV2
4647	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4648
4649	IF (prt_level .GE.10) THEN
4650	print*,'freq_cosp',freq_cosp
4651	ENDIF
4652	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4653	print*,'Dans physiq.F avant appel '
4654	! s ref_liq,ref_ice
4655	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4656	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4657	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4658	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4659	JrNt,ref_liq,ref_ice, &
4660	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4661	zu10m,zv10m,pphis, &
4662	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4663	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4664	prfl(:,1:klev),psfl(:,1:klev), &
4665	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4666	mr_ozone,cldtau, cldemi)
4667	ENDIF
4668	#endif
4669
4670	ENDIF !ok_cosp
4671
4672
4673	! Marine
4674
4675	IF (ok_airs) then
4676
4677	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4678	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4679	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4680	& map_prop_hc,map_prop_hist,&
4681	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4682	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4683	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4684	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4685	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4686	& map_ntot,map_hc,map_hist,&
4687	& map_Cb,map_ThCi,map_Anv,&
4688	& alt_tropo )
4689	ENDIF
4690
4691	ENDIF ! ok_airs
4692
4693
4694	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4695	!AA
4696	!AA Installation de l'interface online-offline pour traceurs
4697	!AA
4698	!====================================================================
4699	! Calcul des tendances traceurs
4700	!====================================================================
4701	!
4702
4703	IF (type_trac=='repr') THEN
4704	sh_in(:,:) = q_seri(:,:)
4705	ELSE
4706	sh_in(:,:) = qx(:,:,ivap)
4707	ch_in(:,:) = qx(:,:,iliq)
4708	ENDIF
4709
4710	IF (iflag_phytrac == 1 ) THEN
4711
4712	#ifdef CPP_Dust
4713	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4714	pdtphys,ftsol, & ! I
4715	t,q_seri,paprs,pplay,RHcl, & ! I
4716	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4717	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4718	u_seri, v_seri, latitude_deg, longitude_deg, &
4719	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4720	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4721	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4722	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4723	fm_therm, entr_therm, rneb, & ! I
4724	beta_prec_fisrt,beta_prec, & !I
4725	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4726	d_tr_dyn,tr_seri)
4727
4728	#else
4729
4730	CALL phytrac ( &
4731	itap, days_elapsed+1, jH_cur, debut, &
4732	lafin, phys_tstep, u, v, t, &
4733	paprs, pplay, pmfu, pmfd, &
4734	pen_u, pde_u, pen_d, pde_d, &
4735	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4736	u1, v1, ftsol, pctsrf, &
4737	zustar, zu10m, zv10m, &
4738	wstar(:,is_ave), ale_bl, ale_wake, &
4739	latitude_deg, longitude_deg, &
4740	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4741	presnivs, pphis, pphi, albsol1, &
4742	sh_in, ch_in, rhcl, cldfra, rneb, &
4743	diafra, cldliq, itop_con, ibas_con, &
4744	pmflxr, pmflxs, prfl, psfl, &
4745	da, phi, mp, upwd, &
4746	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4747	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4748	ev, ep, epmlmMm, eplaMm, & !<<RomP
4749	dnwd, aerosol_couple, flxmass_w, &
4750	tau_aero, piz_aero, cg_aero, ccm, &
4751	rfname, &
4752	d_tr_dyn, & !<<RomP
4753	tr_seri, init_source)
4754	#endif
4755	ENDIF ! (iflag_phytrac=1)
4756
4757	IF (offline) THEN
4758
4759	IF (prt_level.ge.9) &
4760	print*,'Attention on met a 0 les thermiques pour phystoke'
4761	CALL phystokenc ( &
4762	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4763	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4764	fm_therm,entr_therm, &
4765	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4766	frac_impa, frac_nucl, &
4767	pphis,cell_area,phys_tstep,itap, &
4768	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4769
4770
4771	ENDIF
4772
4773	!
4774	! Calculer le transport de l'eau et de l'energie (diagnostique)
4775	!
4776	CALL transp (paprs,zxtsol, &
4777	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4778	ve, vq, ue, uq, vwat, uwat)
4779	!
4780	!IM global posePB BEG
4781	IF(1.EQ.0) THEN
4782	!
4783	CALL transp_lay (paprs,zxtsol, &
4784	t_seri, q_seri, u_seri, v_seri, zphi, &
4785	ve_lay, vq_lay, ue_lay, uq_lay)
4786	!
4787	ENDIF !(1.EQ.0) THEN
4788	!IM global posePB END
4789	! Accumuler les variables a stocker dans les fichiers histoire:
4790	!
4791
4792	!================================================================
4793	! Conversion of kinetic and potential energy into heat, for
4794	! parameterisation of subgrid-scale motions
4795	!================================================================
4796
4797	d_t_ec(:,:)=0.
4798	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4799	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4800	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4801	zmasse,exner,d_t_ec)
4802	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4803
4804	!=======================================================================
4805	! SORTIES
4806	!=======================================================================
4807	!
4808	!IM initialisation + calculs divers diag AMIP2
4809	!
4810	include "calcul_divers.h"
4811	!
4812	!IM Interpolation sur les niveaux de pression du NMC
4813	! -------------------------------------------------
4814	!
4815	include "calcul_STDlev.h"
4816	!
4817	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4818	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4819	!
4820	!cc prw = eau precipitable
4821	! prlw = colonne eau liquide
4822	! prlw = colonne eau solide
4823	prw(:) = 0.
4824	prlw(:) = 0.
4825	prsw(:) = 0.
4826	DO k = 1, klev
4827	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4828	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4829	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4830	ENDDO
4831	!
4832	IF (type_trac == 'inca') THEN
4833	#ifdef INCA
4834	CALL VTe(VTphysiq)
4835	CALL VTb(VTinca)
4836
4837	CALL chemhook_end ( &
4838	phys_tstep, &
4839	pplay, &
4840	t_seri, &
4841	tr_seri, &
4842	nbtr, &
4843	paprs, &
4844	q_seri, &
4845	cell_area, &
4846	pphi, &
4847	pphis, &
4848	zx_rh, &
4849	aps, bps)
4850
4851	CALL VTe(VTinca)
4852	CALL VTb(VTphysiq)
4853	#endif
4854	ENDIF
4855
4856
4857	!
4858	! Convertir les incrementations en tendances
4859	!
4860	IF (prt_level .GE.10) THEN
4861	print *,'Convertir les incrementations en tendances '
4862	ENDIF
4863	!
4864	IF (mydebug) THEN
4865	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4866	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4867	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4868	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4869	ENDIF
4870
4871	DO k = 1, klev
4872	DO i = 1, klon
4873	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
4874	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
4875	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
4876	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
4877	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
4878	!CR: on ajoute le contenu en glace
4879	IF (nqo.eq.3) THEN
4880	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
4881	ENDIF
4882	ENDDO
4883	ENDDO
4884	!
4885	!CR: nb de traceurs eau: nqo
4886	! IF (nqtot.GE.3) THEN
4887	IF (nqtot.GE.(nqo+1)) THEN
4888	! DO iq = 3, nqtot
4889	DO iq = nqo+1, nqtot
4890	DO k = 1, klev
4891	DO i = 1, klon
4892	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
4893	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
4894	ENDDO
4895	ENDDO
4896	ENDDO
4897	ENDIF
4898	!
4899	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4900	!IM global posePB include "write_bilKP_ins.h"
4901	!IM global posePB include "write_bilKP_ave.h"
4902	!
4903
4904	!--OB mass fixer
4905	!--profile is corrected to force mass conservation of water
4906	IF (mass_fixer) THEN
4907	qql2(:)=0.0
4908	DO k = 1, klev
4909	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4910	ENDDO
4911	DO i = 1, klon
4912	!--compute ratio of what q+ql should be with conservation to what it is
4913	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4914	DO k = 1, klev
4915	q_seri(i,k) =q_seri(i,k)*corrqql
4916	ql_seri(i,k)=ql_seri(i,k)*corrqql
4917	ENDDO
4918	ENDDO
4919	ENDIF
4920	!--fin mass fixer
4921
4922	! Sauvegarder les valeurs de t et q a la fin de la physique:
4923	!
4924	u_ancien(:,:) = u_seri(:,:)
4925	v_ancien(:,:) = v_seri(:,:)
4926	t_ancien(:,:) = t_seri(:,:)
4927	q_ancien(:,:) = q_seri(:,:)
4928	ql_ancien(:,:) = ql_seri(:,:)
4929	qs_ancien(:,:) = qs_seri(:,:)
4930	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4931	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4932	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4933	! !! RomP >>>
4934	!CR: nb de traceurs eau: nqo
4935	IF (nqtot.GT.nqo) THEN
4936	DO iq = nqo+1, nqtot
4937	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4938	ENDDO
4939	ENDIF
4940	! !! RomP <<<
4941	!==========================================================================
4942	! Sorties des tendances pour un point particulier
4943	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4944	! pour le debug
4945	! La valeur de igout est attribuee plus haut dans le programme
4946	!==========================================================================
4947
4948	IF (prt_level.ge.1) THEN
4949	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4950	write(lunout,*) &
4951	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4952	write(lunout,*) &
4953	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4954	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4955	pctsrf(igout,is_sic)
4956	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4957	DO k=1,klev
4958	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4959	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4960	d_t_eva(igout,k)
4961	ENDDO
4962	write(lunout,*) 'cool,heat'
4963	DO k=1,klev
4964	write(lunout,*) cool(igout,k),heat(igout,k)
4965	ENDDO
4966
4967	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4968	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4969	!jyg! do k=1,klev
4970	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4971	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4972	!jyg! enddo
4973	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4974	DO k=1,klev
4975	write(lunout,*) d_t_vdf(igout,k), &
4976	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4977	ENDDO
4978	!>jyg
4979
4980	write(lunout,*) 'd_ps ',d_ps(igout)
4981	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4982	DO k=1,klev
4983	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4984	d_qx(igout,k,1),d_qx(igout,k,2)
4985	ENDDO
4986	ENDIF
4987
4988	!============================================================
4989	! Calcul de la temperature potentielle
4990	!============================================================
4991	DO k = 1, klev
4992	DO i = 1, klon
4993	!JYG/IM theta en debut du pas de temps
4994	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4995	!JYG/IM theta en fin de pas de temps de physique
4996	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4997	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4998	! MPL 20130625
4999	! fth_fonctions.F90 et parkind1.F90
5000	! sinon thetal=theta
5001	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5002	! : ql_seri(i,k))
5003	thetal(i,k)=theta(i,k)
5004	ENDDO
5005	ENDDO
5006	!
5007
5008	! 22.03.04 BEG
5009	!=============================================================
5010	! Ecriture des sorties
5011	!=============================================================
5012	#ifdef CPP_IOIPSL
5013
5014	! Recupere des varibles calcule dans differents modules
5015	! pour ecriture dans histxxx.nc
5016
5017	! Get some variables from module fonte_neige_mod
5018	CALL fonte_neige_get_vars(pctsrf, &
5019	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5020
5021
5022	!=============================================================
5023	! Separation entre thermiques et non thermiques dans les sorties
5024	! de fisrtilp
5025	!=============================================================
5026
5027	IF (iflag_thermals>=1) THEN
5028	d_t_lscth=0.
5029	d_t_lscst=0.
5030	d_q_lscth=0.
5031	d_q_lscst=0.
5032	DO k=1,klev
5033	DO i=1,klon
5034	IF (ptconvth(i,k)) THEN
5035	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5036	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5037	ELSE
5038	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5039	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5040	ENDIF
5041	ENDDO
5042	ENDDO
5043
5044	DO i=1,klon
5045	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5046	plul_th(i)=prfl(i,1)+psfl(i,1)
5047	ENDDO
5048	ENDIF
5049
5050	!On effectue les sorties:
5051
5052	#ifdef CPP_Dust
5053	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5054	pplay, lmax_th, aerosol_couple, &
5055	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
5056	ptconv, read_climoz, clevSTD, &
5057	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5058	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5059	#else
5060	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5061	pplay, lmax_th, aerosol_couple, &
5062	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, &
5063	ok_sync, ptconv, read_climoz, clevSTD, &
5064	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5065	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5066	#endif
5067
5068	#ifndef CPP_XIOS
5069	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5070	#endif
5071
5072	#endif
5073
5074	! On remet des variables a .false. apres un premier appel
5075	IF (debut) THEN
5076	#ifdef CPP_XIOS
5077	swaero_diag=.FALSE.
5078	swaerofree_diag=.FALSE.
5079	dryaod_diag=.FALSE.
5080	ok_4xCO2atm= .FALSE.
5081	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5082
5083	IF (is_master) THEN
5084	!--setting up swaero_diag to TRUE in XIOS case
5085	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
5086	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
5087	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
5088	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
5089	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
5090	!!!--for now these fields are not in the XML files so they are omitted
5091	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
5092	swaero_diag=.TRUE.
5093
5094	!--setting up swaerofree_diag to TRUE in XIOS case
5095	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
5096	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
5097	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
5098	xios_field_is_active("LWupTOAcleanclr")) &
5099	swaerofree_diag=.TRUE.
5100
5101	!--setting up dryaod_diag to TRUE in XIOS case
5102	DO naero = 1, naero_tot-1
5103	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
5104	ENDDO
5105	!
5106	!--setting up ok_4xCO2atm to TRUE in XIOS case
5107	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
5108	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
5109	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
5110	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
5111	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
5112	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
5113	ok_4xCO2atm=.TRUE.
5114	ENDIF
5115	!$OMP BARRIER
5116	CALL bcast(swaero_diag)
5117	CALL bcast(swaerofree_diag)
5118	CALL bcast(dryaod_diag)
5119	CALL bcast(ok_4xCO2atm)
5120	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5121	#endif
5122	ENDIF
5123
5124	!====================================================================
5125	! Arret du modele apres hgardfou en cas de detection d'un
5126	! plantage par hgardfou
5127	!====================================================================
5128
5129	IF (abortphy==1) THEN
5130	abort_message ='Plantage hgardfou'
5131	CALL abort_physic (modname,abort_message,1)
5132	ENDIF
5133
5134	! 22.03.04 END
5135	!
5136	!====================================================================
5137	! Si c'est la fin, il faut conserver l'etat de redemarrage
5138	!====================================================================
5139	!
5140
5141	IF (lafin) THEN
5142	itau_phy = itau_phy + itap
5143	CALL phyredem ("restartphy.nc")
5144	! open(97,form="unformatted",file="finbin")
5145	! write(97) u_seri,v_seri,t_seri,q_seri
5146	! close(97)
5147
5148	IF (is_omp_master) THEN
5149
5150	IF (read_climoz >= 1) THEN
5151	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5152	DEALLOCATE(press_edg_climoz) ! pointer
5153	DEALLOCATE(press_cen_climoz) ! pointer
5154	ENDIF
5155
5156	ENDIF
5157	#ifdef CPP_XIOS
5158	IF (is_omp_master) CALL xios_context_finalize
5159	#endif
5160	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5161	ENDIF
5162
5163	! first=.false.
5164
5165	END SUBROUTINE physiq
5166
5167	END MODULE physiq_mod

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

Download in other formats:

Original Format