Context Navigation

physiq_mod.F90 @ 4660

Last change on this file since 4660 was 4657, checked in by fhourdin, 16 months ago
Poursuite nettoyage replauisation
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: 203.1 KB

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

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

Download in other formats:

Original Format