Context Navigation

physiq_mod.F90 @ 3577

Last change on this file since 3577 was 3577, checked in by fhourdin, 5 years ago
Corrections d'erreurs identifiees en debug ancienne physique sur Jean-Zay Vive le poihl interactif.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 182.8 KB

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

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

Download in other formats:

Original Format