Context Navigation

physiq_mod.F90 @ 4229

Last change on this file since 4229 was 4229, checked in by oboucher, 2 years ago
Cleaning up the energy and humidity diag transport routine
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: 191.7 KB

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

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

Download in other formats:

Original Format