Context Navigation

physiq_mod.F90 @ 4651

Last change on this file since 4651 was 4651, checked in by fhourdin, 9 months ago

Pre-replayisation pour cloudth

Cration de lmdz_cloudth_ini
cloudth_mod.F90 -> lmdz_cloudth.F90
Supression des constantes dans lmdz_cloudth
En fait, cloudth lui même n'est pas vraiment replayisable car appeler ligne par ligne.

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: 202.6 KB

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

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

Download in other formats:

Original Format