Context Navigation

physiq_mod.F90 @ 4678

Last change on this file since 4678 was 4678, checked in by fhourdin, 8 months ago
Petits nettoyages sur les thermiques
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: 206.2 KB

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

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

Download in other formats:

Original Format