Context Navigation

physiq_mod.F90 @ 4062

Last change on this file since 4062 was 4062, checked in by oboucher, 2 years ago
more keys for plane effects to avoid pb on outputs
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: 189.5 KB

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

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

Download in other formats:

Original Format