Context Navigation

physiq_mod.F90 @ 3627

Last change on this file since 3627 was 3617, checked in by lguez, 4 years ago
Sync latest trunk changes to branch `Ocean_skin`
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: 183.1 KB

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

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

Download in other formats:

Original Format