Context Navigation

physiq_mod.F90 @ 4449

Last change on this file since 4449 was 4449, checked in by evignon, 15 months ago
commission du nouveau schema de turbulence developpe dans le cadre de l'atelier tke
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 `Id`
File size: 192.4 KB

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

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

Download in other formats:

Original Format