Context Navigation

physiq_mod.F90 @ 4069

Last change on this file since 4069 was 4069, checked in by dcugnet, 2 years ago
Fix in strings_mod ; parser routines taken from version 9 of https://svn.lmd.jussieu.fr/tracers-parser
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.6 KB

Line
1	!
2	! $Id: physiq_mod.F90 4069 2022-01-30 17:58:23Z dcugnet $
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: addPhase
43	USE strings_mod, ONLY: strIdx
44	USE iophy
45	USE limit_read_mod, ONLY : init_limit_read
46	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
47	USE mod_phys_lmdz_mpi_data, only: is_mpi_root
48	USE mod_phys_lmdz_para
49	USE netcdf95, only: nf95_close
50	USE netcdf, only: nf90_fill_real ! IM for NMC files
51	USE open_climoz_m, only: open_climoz ! ozone climatology from a file
52	USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
53	USE pbl_surface_mod, ONLY : pbl_surface
54	USE phyaqua_mod, only: zenang_an
55	USE 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	rneb_seri(i,k) = 0.
2247	ELSE IF (nqo.eq.3) THEN
2248	qs_seri(i,k) = qx(i,k,isol)
2249	rneb_seri(i,k) = 0.
2250	ELSE IF (nqo.eq.4) THEN
2251	qs_seri(i,k) = qx(i,k,isol)
2252	rneb_seri(i,k) = qx(i,k,irneb)
2253	ENDIF
2254	ENDDO
2255	ENDDO
2256	!
2257	!--OB mass fixer
2258	IF (mass_fixer) THEN
2259	!--store initial water burden
2260	qql1(:)=0.0
2261	DO k = 1, klev
2262	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2263	ENDDO
2264	ENDIF
2265	!--fin mass fixer
2266
2267	tke0(:,:)=pbl_tke(:,:,is_ave)
2268	IF (nqtot > nqo) THEN
2269	! water isotopes are not included in tr_seri
2270	itr = 0
2271	DO iq = 1, nqtot
2272	IF(tracers(iq)%isH2Ofamily) CYCLE
2273	itr = itr+1
2274	DO k = 1, klev
2275	DO i = 1, klon
2276	tr_seri(i,k,itr) = qx(i,k,iq)
2277	ENDDO
2278	ENDDO
2279	ENDDO
2280	ELSE
2281	! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
2282	! tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
2283	tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','v',''))) = 0.0
2284	ENDIF
2285	!
2286	! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2287	! LF
2288	IF (debut) THEN
2289	WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
2290	itr = 0
2291	do iq = 1, nqtot
2292	IF(tracers(iq)%isH2Ofamily) CYCLE
2293	itr = itr+1
2294	tr_ancien(:,:,itr)=tr_seri(:,:,itr)
2295	enddo
2296	ENDIF
2297	!
2298	DO i = 1, klon
2299	ztsol(i) = 0.
2300	ENDDO
2301	DO nsrf = 1, nbsrf
2302	DO i = 1, klon
2303	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2304	ENDDO
2305	ENDDO
2306	! Initialize variables used for diagnostic purpose
2307	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2308
2309	! Diagnostiquer la tendance dynamique
2310	!
2311	IF (ancien_ok) THEN
2312	!
2313	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2314	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2315	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2316	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2317	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2318	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2319	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2320	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2321	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2322	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2323	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2324	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2325	! !! RomP >>> td dyn traceur
2326	IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
2327	! !! RomP <<<
2328	!!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
2329	d_rneb_dyn(:,:)=0.0
2330	ELSE
2331	d_u_dyn(:,:) = 0.0
2332	d_v_dyn(:,:) = 0.0
2333	d_t_dyn(:,:) = 0.0
2334	d_q_dyn(:,:) = 0.0
2335	d_ql_dyn(:,:) = 0.0
2336	d_qs_dyn(:,:) = 0.0
2337	d_q_dyn2d(:) = 0.0
2338	d_ql_dyn2d(:) = 0.0
2339	d_qs_dyn2d(:) = 0.0
2340	! !! RomP >>> td dyn traceur
2341	IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
2342	! !! RomP <<<
2343	d_rneb_dyn(:,:)=0.0
2344	ancien_ok = .TRUE.
2345	ENDIF
2346	!
2347	! Ajouter le geopotentiel du sol:
2348	!
2349	DO k = 1, klev
2350	DO i = 1, klon
2351	zphi(i,k) = pphi(i,k) + pphis(i)
2352	ENDDO
2353	ENDDO
2354	!
2355	! Verifier les temperatures
2356	!
2357	!IM BEG
2358	IF (check) THEN
2359	amn=MIN(ftsol(1,is_ter),1000.)
2360	amx=MAX(ftsol(1,is_ter),-1000.)
2361	DO i=2, klon
2362	amn=MIN(ftsol(i,is_ter),amn)
2363	amx=MAX(ftsol(i,is_ter),amx)
2364	ENDDO
2365	!
2366	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2367	ENDIF !(check) THEN
2368	!IM END
2369	!
2370	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2371	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2372
2373	!
2374	!IM BEG
2375	IF (check) THEN
2376	amn=MIN(ftsol(1,is_ter),1000.)
2377	amx=MAX(ftsol(1,is_ter),-1000.)
2378	DO i=2, klon
2379	amn=MIN(ftsol(i,is_ter),amn)
2380	amx=MAX(ftsol(i,is_ter),amx)
2381	ENDDO
2382	!
2383	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2384	ENDIF !(check) THEN
2385	!IM END
2386	!
2387	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2388	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2389	!
2390	! Update ozone if day change
2391	IF (MOD(itap-1,lmt_pas) == 0) THEN
2392	IF (read_climoz <= 0) THEN
2393	! Once per day, update ozone from Royer:
2394	IF (solarlong0<-999.) then
2395	! Generic case with evolvoing season
2396	zzz=real(days_elapsed+1)
2397	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2398	! Particular case with annual mean insolation
2399	zzz=real(90) ! could be revisited
2400	IF (read_climoz/=-1) THEN
2401	abort_message ='read_climoz=-1 is recommended when ' &
2402	// 'solarlong0=1000.'
2403	CALL abort_physic (modname,abort_message,1)
2404	ENDIF
2405	ELSE
2406	! Case where the season is imposed with solarlong0
2407	zzz=real(90) ! could be revisited
2408	ENDIF
2409
2410	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2411	#ifdef REPROBUS
2412	ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2413	DO i = 1, klon
2414	Z1=t_seri(i,itroprep(i)+1)
2415	Z2=t_seri(i,itroprep(i))
2416	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2417	B=Z2-fac*alog(pplay(i,itroprep(i)))
2418	ttrop(i)= fac*alog(ptrop(i))+B
2419	!
2420	Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2421	Z2= 1.e-3 * ( pphi(i,itroprep(i)) +pphis(i) ) / gravit
2422	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2423	B=Z2-fac*alog(pplay(i,itroprep(i)))
2424	ztrop(i)=fac*alog(ptrop(i))+B
2425	ENDDO
2426	#endif
2427	ELSE
2428	!--- ro3i = elapsed days number since current year 1st january, 0h
2429	ro3i=days_elapsed+jh_cur-jh_1jan
2430	!--- scaling for old style files (360 records)
2431	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2432	IF(adjust_tropopause) THEN
2433	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2434	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2435	time_climoz , longitude_deg, latitude_deg, &
2436	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2437	ELSE
2438	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2439	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2440	time_climoz )
2441	ENDIF
2442	! Convert from mole fraction of ozone to column density of ozone in a
2443	! cell, in kDU:
2444	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2445	* zmasse / dobson_u / 1e3
2446	! (By regridding ozone values for LMDZ only once a day, we
2447	! have already neglected the variation of pressure in one
2448	! day. So do not recompute "wo" at each time step even if
2449	! "zmasse" changes a little.)
2450	ENDIF
2451	ENDIF
2452	!
2453	! Re-evaporer l'eau liquide nuageuse
2454	!
2455	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2456	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2457
2458	CALL add_phys_tend &
2459	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2460	'eva',abortphy,flag_inhib_tend,itap,0)
2461	CALL prt_enerbil('eva',itap)
2462
2463	!=========================================================================
2464	! Calculs de l'orbite.
2465	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2466	! doit donc etre plac\'e avant radlwsw et pbl_surface
2467
2468	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2469	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2470	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2471	!
2472	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2473	! solarlong0
2474	IF (solarlong0<-999.) THEN
2475	IF (new_orbit) THEN
2476	! calcul selon la routine utilisee pour les planetes
2477	CALL solarlong(day_since_equinox, zlongi, dist)
2478	ELSE
2479	! calcul selon la routine utilisee pour l'AR4
2480	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2481	ENDIF
2482	ELSE
2483	zlongi=solarlong0 ! longitude solaire vraie
2484	dist=1. ! distance au soleil / moyenne
2485	ENDIF
2486
2487	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2488
2489
2490	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2491	! Calcul de l'ensoleillement :
2492	! ============================
2493	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2494	! l'annee a partir d'une formule analytique.
2495	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2496	! non nul aux poles.
2497	IF (abs(solarlong0-1000.)<1.e-4) THEN
2498	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2499	latitude_deg,longitude_deg,rmu0,fract)
2500	swradcorr(:) = 1.0
2501	JrNt(:) = 1.0
2502	zrmu0(:) = rmu0(:)
2503	ELSE
2504	! recode par Olivier Boucher en sept 2015
2505	SELECT CASE (iflag_cycle_diurne)
2506	CASE(0)
2507	! Sans cycle diurne
2508	CALL angle(zlongi, latitude_deg, fract, rmu0)
2509	swradcorr = 1.0
2510	JrNt = 1.0
2511	zrmu0 = rmu0
2512	CASE(1)
2513	! Avec cycle diurne sans application des poids
2514	! bit comparable a l ancienne formulation cycle_diurne=true
2515	! on integre entre gmtime et gmtime+radpas
2516	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2517	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2518	latitude_deg,longitude_deg,rmu0,fract)
2519	zrmu0 = rmu0
2520	swradcorr = 1.0
2521	! Calcul du flag jour-nuit
2522	JrNt = 0.0
2523	WHERE (fract.GT.0.0) JrNt = 1.0
2524	CASE(2)
2525	! Avec cycle diurne sans application des poids
2526	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2527	! Comme cette routine est appele a tous les pas de temps de
2528	! la physique meme si le rayonnement n'est pas appele je
2529	! remonte en arriere les radpas-1 pas de temps
2530	! suivant. Petite ruse avec MOD pour prendre en compte le
2531	! premier pas de temps de la physique pendant lequel
2532	! itaprad=0
2533	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2534	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2535	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2536	latitude_deg,longitude_deg,rmu0,fract)
2537	!
2538	! Calcul des poids
2539	!
2540	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2541	zdtime2=0.0 !--pas de temps de la physique qui se termine
2542	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2543	latitude_deg,longitude_deg,zrmu0,zfract)
2544	swradcorr = 0.0
2545	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2546	swradcorr=zfract/fract*zrmu0/rmu0
2547	! Calcul du flag jour-nuit
2548	JrNt = 0.0
2549	WHERE (zfract.GT.0.0) JrNt = 1.0
2550	END SELECT
2551	ENDIF
2552	sza_o = ACOS (rmu0) *180./pi
2553
2554	IF (mydebug) THEN
2555	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2556	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2557	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2558	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2559	ENDIF
2560
2561	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2562	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2563	! Cela implique tous les interactions des sous-surfaces et la
2564	! partie diffusion turbulent du couche limit.
2565	!
2566	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2567	! ecriture des fihiers hist_XXXX.nc, ces sont :
2568	! qsol, zq2m, s_pblh, s_lcl,
2569	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2570	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2571	! zu10m, zv10m, fder,
2572	! zxqsurf, delta_qsurf,
2573	! rh2m, zxfluxu, zxfluxv,
2574	! frugs, agesno, fsollw, fsolsw,
2575	! d_ts, fevap, fluxlat, t2m,
2576	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2577	!
2578	! Certains ne sont pas utiliser du tout :
2579	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2580	!
2581
2582	! Calcul de l'humidite de saturation au niveau du sol
2583
2584
2585
2586	IF (iflag_pbl/=0) THEN
2587
2588	!jyg+nrlmd<
2589	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2590	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2591	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2592	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2593	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2594	ENDIF
2595	! !!
2596	!>jyg+nrlmd
2597	!
2598	!-------gustiness calculation-------!
2599	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2600	gustiness=0 !ym missing init
2601
2602	IF (iflag_gusts==0) THEN
2603	gustiness(1:klon)=0
2604	ELSE IF (iflag_gusts==1) THEN
2605	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2606	ELSE IF (iflag_gusts==2) THEN
2607	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2608	! ELSE IF (iflag_gusts==2) THEN
2609	! do i = 1, klon
2610	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2611	! *ale_wake(i) !! need to make sigma_wk accessible here
2612	! enddo
2613	! ELSE IF (iflag_gusts==3) THEN
2614	! do i = 1, klon
2615	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2616	! enddo
2617	ENDIF
2618
2619	CALL pbl_surface( &
2620	phys_tstep, date0, itap, days_elapsed+1, &
2621	debut, lafin, &
2622	longitude_deg, latitude_deg, rugoro, zrmu0, &
2623	sollwdown, cldt, &
2624	rain_fall, snow_fall, solsw, solswfdiff, sollw, &
2625	gustiness, &
2626	t_seri, q_seri, u_seri, v_seri, &
2627	!nrlmd+jyg<
2628	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2629	!>nrlmd+jyg
2630	pplay, paprs, pctsrf, &
2631	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2632	!albedo SB <<<
2633	cdragh, cdragm, u1, v1, &
2634	beta_aridity, &
2635	!albedo SB >>>
2636	! albsol1, albsol2, sens, evap, &
2637	albsol_dir, albsol_dif, sens, evap, &
2638	!albedo SB <<<
2639	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2640	zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, &
2641	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2642	!nrlmd<
2643	!jyg<
2644	d_t_vdf_w, d_q_vdf_w, &
2645	d_t_vdf_x, d_q_vdf_x, &
2646	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2647	!>jyg
2648	delta_tsurf,wake_dens, &
2649	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2650	kh,kh_x,kh_w, &
2651	!>nrlmd
2652	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2653	slab_wfbils, &
2654	qsol, zq2m, s_pblh, s_lcl, &
2655	!jyg<
2656	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2657	!>jyg
2658	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2659	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2660	zustar, zu10m, zv10m, fder, &
2661	zxqsurf, delta_qsurf, rh2m, zxfluxu, zxfluxv, &
2662	z0m, z0h, agesno, fsollw, fsolsw, &
2663	d_ts, fevap, fluxlat, t2m, &
2664	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2665	fluxt, fluxu, fluxv, &
2666	dsens, devap, zxsnow, &
2667	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2668	!nrlmd+jyg<
2669	wake_delta_pbl_TKE, &
2670	!>nrlmd+jyg
2671	treedrg )
2672	!FC
2673	!
2674	! Add turbulent diffusion tendency to the wake difference variables
2675	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2676	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2677	!jyg<
2678	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2679	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2680	CALL add_wake_tend &
2681	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2682	ELSE
2683	d_deltat_vdf(:,:) = 0.
2684	d_deltaq_vdf(:,:) = 0.
2685	!>jyg
2686	ENDIF
2687
2688	!---------------------------------------------------------------------
2689	! ajout des tendances de la diffusion turbulente
2690	IF (klon_glo==1) THEN
2691	CALL add_pbl_tend &
2692	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2693	'vdf',abortphy,flag_inhib_tend,itap)
2694	ELSE
2695	CALL add_phys_tend &
2696	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2697	'vdf',abortphy,flag_inhib_tend,itap,0)
2698	ENDIF
2699	CALL prt_enerbil('vdf',itap)
2700	!--------------------------------------------------------------------
2701
2702	IF (mydebug) THEN
2703	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2704	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2705	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2706	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2707	ENDIF
2708
2709	!albedo SB >>>
2710	albsol1=0.
2711	albsol2=0.
2712	falb1=0.
2713	falb2=0.
2714	SELECT CASE(nsw)
2715	CASE(2)
2716	albsol1=albsol_dir(:,1)
2717	albsol2=albsol_dir(:,2)
2718	falb1=falb_dir(:,1,:)
2719	falb2=falb_dir(:,2,:)
2720	CASE(4)
2721	albsol1=albsol_dir(:,1)
2722	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2723	+albsol_dir(:,4)*SFRWL(4)
2724	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2725	falb1=falb_dir(:,1,:)
2726	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2727	+falb_dir(:,4,:)*SFRWL(4)
2728	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2729	CASE(6)
2730	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2731	+albsol_dir(:,3)*SFRWL(3)
2732	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2733	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2734	+albsol_dir(:,6)*SFRWL(6)
2735	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2736	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2737	+falb_dir(:,3,:)*SFRWL(3)
2738	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2739	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2740	+falb_dir(:,6,:)*SFRWL(6)
2741	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2742	END SELECt
2743	!albedo SB <<<
2744
2745
2746	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2747	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2748
2749	ENDIF
2750	! =================================================================== c
2751	! Calcul de Qsat
2752
2753	DO k = 1, klev
2754	DO i = 1, klon
2755	zx_t = t_seri(i,k)
2756	IF (thermcep) THEN
2757	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2758	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2759	zx_qs = MIN(0.5,zx_qs)
2760	zcor = 1./(1.-retv*zx_qs)
2761	zx_qs = zx_qs*zcor
2762	ELSE
2763	!! IF (zx_t.LT.t_coup) THEN !jyg
2764	IF (zx_t.LT.rtt) THEN !jyg
2765	zx_qs = qsats(zx_t)/pplay(i,k)
2766	ELSE
2767	zx_qs = qsatl(zx_t)/pplay(i,k)
2768	ENDIF
2769	ENDIF
2770	zqsat(i,k)=zx_qs
2771	ENDDO
2772	ENDDO
2773
2774	IF (prt_level.ge.1) THEN
2775	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2776	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2777	ENDIF
2778	!
2779	! Appeler la convection (au choix)
2780	!
2781	DO k = 1, klev
2782	DO i = 1, klon
2783	conv_q(i,k) = d_q_dyn(i,k) &
2784	+ d_q_vdf(i,k)/phys_tstep
2785	conv_t(i,k) = d_t_dyn(i,k) &
2786	+ d_t_vdf(i,k)/phys_tstep
2787	ENDDO
2788	ENDDO
2789	IF (check) THEN
2790	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2791	WRITE(lunout,*) "avantcon=", za
2792	ENDIF
2793	zx_ajustq = .FALSE.
2794	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2795	IF (zx_ajustq) THEN
2796	DO i = 1, klon
2797	z_avant(i) = 0.0
2798	ENDDO
2799	DO k = 1, klev
2800	DO i = 1, klon
2801	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2802	*(paprs(i,k)-paprs(i,k+1))/RG
2803	ENDDO
2804	ENDDO
2805	ENDIF
2806
2807	! Calcule de vitesse verticale a partir de flux de masse verticale
2808	DO k = 1, klev
2809	DO i = 1, klon
2810	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2811	ENDDO
2812	ENDDO
2813
2814	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2815	omega(igout, :)
2816	!
2817	! Appel de la convection tous les "cvpas"
2818	!
2819	!!jyg IF (MOD(itapcv,cvpas).EQ.0) THEN
2820	!! print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2821	!! itapcv, cvpas, itap-1, cvpas_0
2822	IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
2823
2824	!
2825	! Mettre a zero des variables de sortie (pour securite)
2826	!
2827	pmflxr(:,:) = 0.
2828	pmflxs(:,:) = 0.
2829	wdtrainA(:,:) = 0.
2830	wdtrainS(:,:) = 0.
2831	wdtrainM(:,:) = 0.
2832	upwd(:,:) = 0.
2833	dnwd(:,:) = 0.
2834	ep(:,:) = 0.
2835	da(:,:)=0.
2836	mp(:,:)=0.
2837	wght_cvfd(:,:)=0.
2838	phi(:,:,:)=0.
2839	phi2(:,:,:)=0.
2840	epmlmMm(:,:,:)=0.
2841	eplaMm(:,:)=0.
2842	d1a(:,:)=0.
2843	dam(:,:)=0.
2844	elij(:,:,:)=0.
2845	ev(:,:)=0.
2846	qtaa(:,:)=0.
2847	clw(:,:)=0.
2848	sij(:,:,:)=0.
2849	!
2850	IF (iflag_con.EQ.1) THEN
2851	abort_message ='reactiver le call conlmd dans physiq.F'
2852	CALL abort_physic (modname,abort_message,1)
2853	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
2854	! . d_t_con, d_q_con,
2855	! . rain_con, snow_con, ibas_con, itop_con)
2856	ELSE IF (iflag_con.EQ.2) THEN
2857	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
2858	conv_t, conv_q, -evap, omega, &
2859	d_t_con, d_q_con, rain_con, snow_con, &
2860	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2861	kcbot, kctop, kdtop, pmflxr, pmflxs)
2862	d_u_con = 0.
2863	d_v_con = 0.
2864
2865	WHERE (rain_con < 0.) rain_con = 0.
2866	WHERE (snow_con < 0.) snow_con = 0.
2867	DO i = 1, klon
2868	ibas_con(i) = klev+1 - kcbot(i)
2869	itop_con(i) = klev+1 - kctop(i)
2870	ENDDO
2871	ELSE IF (iflag_con.GE.3) THEN
2872	! nb of tracers for the KE convection:
2873	! MAF la partie traceurs est faite dans phytrac
2874	! on met ntra=1 pour limiter les appels mais on peut
2875	! supprimer les calculs / ftra.
2876	ntra = 1
2877
2878	!=======================================================================
2879	!ajout pour la parametrisation des poches froides: calcul de
2880	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2881	IF (iflag_wake>=1) THEN
2882	DO k=1,klev
2883	DO i=1,klon
2884	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2885	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2886	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2887	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2888	ENDDO
2889	ENDDO
2890	ELSE
2891	t_w(:,:) = t_seri(:,:)
2892	q_w(:,:) = q_seri(:,:)
2893	t_x(:,:) = t_seri(:,:)
2894	q_x(:,:) = q_seri(:,:)
2895	ENDIF
2896	!
2897	!jyg<
2898	! Perform dry adiabatic adjustment on wake profile
2899	! The corresponding tendencies are added to the convective tendencies
2900	! after the call to the convective scheme.
2901	IF (iflag_wake>=1) then
2902	IF (iflag_adjwk >= 1) THEN
2903	limbas(:) = 1
2904	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2905	d_t_adjwk, d_q_adjwk)
2906	!
2907	DO k=1,klev
2908	DO i=1,klon
2909	IF (wake_s(i) .GT. 1.e-3) THEN
2910	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2911	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2912	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2913	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2914	ELSE
2915	d_deltat_ajs_cv(i,k) = 0.
2916	d_deltaq_ajs_cv(i,k) = 0.
2917	ENDIF
2918	ENDDO
2919	ENDDO
2920	IF (iflag_adjwk == 2) THEN
2921	CALL add_wake_tend &
2922	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2923	ENDIF ! (iflag_adjwk == 2)
2924	ENDIF ! (iflag_adjwk >= 1)
2925	ENDIF ! (iflag_wake>=1)
2926	!>jyg
2927	!
2928
2929	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
2930	!! (k, q_w(1,k), q_x(1,k),k=1,25)
2931
2932	!jyg<
2933	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
2934	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
2935	ale_bl_prescr, alp_bl_prescr, &
2936	wake_pe, wake_fip, &
2937	Ale_bl, Ale_bl_trig, Alp_bl, &
2938	Ale, Alp , Ale_wake, Alp_wake)
2939	!>jyg
2940	!
2941	! sb, oct02:
2942	! Schema de convection modularise et vectorise:
2943	! (driver commun aux versions 3 et 4)
2944	!
2945	IF (ok_cvl) THEN ! new driver for convectL
2946	!
2947	!jyg<
2948	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2949	! Calculate the upmost level of deep convection loops: k_upper_cv
2950	! (near 22 km)
2951	k_upper_cv = klev
2952	!izero = klon/2+1/klon
2953	!DO k = klev,1,-1
2954	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2955	!ENDDO
2956	! FH : nouveau calcul base sur un profil global sans quoi
2957	! le modele etait sensible au decoupage de domaines
2958	DO k = klev,1,-1
2959	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
2960	ENDDO
2961	IF (prt_level .ge. 5) THEN
2962	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2963	k_upper_cv
2964	ENDIF
2965	!
2966	!>jyg
2967	IF (type_trac == 'repr') THEN
2968	nbtr_tmp=ntra
2969	ELSE
2970	nbtr_tmp=nbtr
2971	ENDIF
2972	!jyg iflag_con est dans clesphys
2973	!c CALL concvl (iflag_con,iflag_clos,
2974	CALL concvl (iflag_clos, &
2975	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
2976	t_w,q_w,wake_s, &
2977	u_seri,v_seri,tr_seri,nbtr_tmp, &
2978	ALE,ALP, &
2979	sig1,w01, &
2980	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2981	rain_con, snow_con, ibas_con, itop_con, sigd, &
2982	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2983	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2984	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2985	! RomP >>>
2986	!! . pmflxr,pmflxs,da,phi,mp,
2987	!! . ftd,fqd,lalim_conv,wght_th)
2988	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
2989	ftd,fqd,lalim_conv,wght_th, &
2990	ev, ep,epmlmMm,eplaMm, &
2991	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2992	tau_cld_cv,coefw_cld_cv,epmax_diag)
2993
2994	! RomP <<<
2995
2996	!IM begin
2997	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2998	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2999	!IM end
3000	!IM cf. FH
3001	clwcon0=qcondc
3002	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
3003	!
3004	!jyg<
3005	! If convective tendencies are too large, then call convection
3006	! every time step
3007	cvpas = cvpas_0
3008	DO k=1,k_upper_cv
3009	DO i=1,klon
3010	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3011	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3012	dtcon_multistep_max = 3.
3013	dqcon_multistep_max = 0.02
3014	ENDIF
3015	ENDDO
3016	ENDDO
3017	!
3018	DO k=1,k_upper_cv
3019	DO i=1,klon
3020	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3021	!! abs(d_q_con(i,k)) > 2.e-2) THEN
3022	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3023	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3024	cvpas = 1
3025	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3026	!! i,k,d_t_con(i,k),d_q_con(i,k)
3027	ENDIF
3028	ENDDO
3029	ENDDO
3030	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3031	!!! call bcast(cvpas)
3032	!!! ------------------------------------------------------------
3033	!>jyg
3034	!
3035	DO i = 1, klon
3036	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
3037	ENDDO
3038	!
3039	!jyg<
3040	! Add the tendency due to the dry adjustment of the wake profile
3041	IF (iflag_wake>=1) THEN
3042	IF (iflag_adjwk == 2) THEN
3043	DO k=1,klev
3044	DO i=1,klon
3045	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3046	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
3047	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3048	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3049	ENDDO
3050	ENDDO
3051	ENDIF ! (iflag_adjwk = 2)
3052	ENDIF ! (iflag_wake>=1)
3053	!>jyg
3054	!
3055	ELSE ! ok_cvl
3056
3057	! MAF conema3 ne contient pas les traceurs
3058	CALL conema3 (phys_tstep, &
3059	paprs,pplay,t_seri,q_seri, &
3060	u_seri,v_seri,tr_seri,ntra, &
3061	sig1,w01, &
3062	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3063	rain_con, snow_con, ibas_con, itop_con, &
3064	upwd,dnwd,dnwd0,bas,top, &
3065	Ma,cape,tvp,rflag, &
3066	pbase &
3067	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3068	,clwcon0)
3069
3070	ENDIF ! ok_cvl
3071
3072	!
3073	! Correction precip
3074	rain_con = rain_con * cvl_corr
3075	snow_con = snow_con * cvl_corr
3076	!
3077
3078	IF (.NOT. ok_gust) THEN
3079	do i = 1, klon
3080	wd(i)=0.0
3081	enddo
3082	ENDIF
3083
3084	! =================================================================== c
3085	! Calcul des proprietes des nuages convectifs
3086	!
3087
3088	! calcul des proprietes des nuages convectifs
3089	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3090	IF (iflag_cld_cv == 0) THEN
3091	CALL clouds_gno &
3092	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3093	ELSE
3094	CALL clouds_bigauss &
3095	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3096	ENDIF
3097
3098
3099	! =================================================================== c
3100
3101	DO i = 1, klon
3102	itop_con(i) = min(max(itop_con(i),1),klev)
3103	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3104	ENDDO
3105
3106	DO i = 1, klon
3107	! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3108	! où i n'est pas un point convectif et donc ibas_con(i)=0
3109	! c'est un pb indépendant des isotopes
3110	if (ibas_con(i) > 0) then
3111	ema_pcb(i) = paprs(i,ibas_con(i))
3112	else
3113	ema_pcb(i) = 0.0
3114	endif
3115	ENDDO
3116	DO i = 1, klon
3117	! L'idicage de itop_con peut cacher un pb potentiel
3118	! FH sous la dictee de JYG, CR
3119	ema_pct(i) = paprs(i,itop_con(i)+1)
3120
3121	IF (itop_con(i).gt.klev-3) THEN
3122	IF (prt_level >= 9) THEN
3123	write(lunout,*)'La convection monte trop haut '
3124	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
3125	ENDIF
3126	ENDIF
3127	ENDDO
3128	ELSE IF (iflag_con.eq.0) THEN
3129	write(lunout,*) 'On n appelle pas la convection'
3130	clwcon0=0.
3131	rnebcon0=0.
3132	d_t_con=0.
3133	d_q_con=0.
3134	d_u_con=0.
3135	d_v_con=0.
3136	rain_con=0.
3137	snow_con=0.
3138	bas=1
3139	top=1
3140	ELSE
3141	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
3142	CALL abort_physic("physiq", "", 1)
3143	ENDIF
3144
3145	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3146	! . d_u_con, d_v_con)
3147
3148	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3149	proba_notrig(:) = 1.
3150	itapcv = 0
3151	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3152	!
3153	itapcv = itapcv+1
3154	!
3155	! Compter les steps ou cvpas=1
3156	IF (cvpas == 1) THEN
3157	Ncvpaseq1 = Ncvpaseq1+1
3158	ENDIF
3159	IF (mod(itap,1000) == 0) THEN
3160	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3161	ENDIF
3162
3163	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3164	!!! l'energie dans les courants satures.
3165	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3166	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3167	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3168	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3169	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3170	!! itap, 1)
3171	!! call prt_enerbil('convection_sat',itap)
3172	!!
3173	!!
3174	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
3175	'convection',abortphy,flag_inhib_tend,itap,0)
3176	CALL prt_enerbil('convection',itap)
3177
3178	!-------------------------------------------------------------------------
3179
3180	IF (mydebug) THEN
3181	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3182	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3183	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3184	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3185	ENDIF
3186
3187	IF (check) THEN
3188	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3189	WRITE(lunout,*)"aprescon=", za
3190	zx_t = 0.0
3191	za = 0.0
3192	DO i = 1, klon
3193	za = za + cell_area(i)/REAL(klon)
3194	zx_t = zx_t + (rain_con(i)+ &
3195	snow_con(i))*cell_area(i)/REAL(klon)
3196	ENDDO
3197	zx_t = zx_t/za*phys_tstep
3198	WRITE(lunout,*)"Precip=", zx_t
3199	ENDIF
3200	IF (zx_ajustq) THEN
3201	DO i = 1, klon
3202	z_apres(i) = 0.0
3203	ENDDO
3204	DO k = 1, klev
3205	DO i = 1, klon
3206	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3207	*(paprs(i,k)-paprs(i,k+1))/RG
3208	ENDDO
3209	ENDDO
3210	DO i = 1, klon
3211	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
3212	/z_apres(i)
3213	ENDDO
3214	DO k = 1, klev
3215	DO i = 1, klon
3216	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3217	z_factor(i).LT.(1.0-1.0E-08)) THEN
3218	q_seri(i,k) = q_seri(i,k) * z_factor(i)
3219	ENDIF
3220	ENDDO
3221	ENDDO
3222	ENDIF
3223	zx_ajustq=.FALSE.
3224
3225	!
3226	!==========================================================================
3227	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3228	!pour la couche limite diffuse pour l instant
3229	!
3230	!
3231	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3232	! il faut rajouter cette tendance calcul\'ee hors des poches
3233	! froides
3234	!
3235	IF (iflag_wake>=1) THEN
3236	!
3237	!
3238	! Call wakes every "wkpas" step
3239	!
3240	IF (MOD(itapwk,wkpas).EQ.0) THEN
3241	!
3242	DO k=1,klev
3243	DO i=1,klon
3244	dt_dwn(i,k) = ftd(i,k)
3245	dq_dwn(i,k) = fqd(i,k)
3246	M_dwn(i,k) = dnwd0(i,k)
3247	M_up(i,k) = upwd(i,k)
3248	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3249	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3250	ENDDO
3251	ENDDO
3252
3253	IF (iflag_wake==2) THEN
3254	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3255	DO k = 1,klev
3256	dt_dwn(:,k)= dt_dwn(:,k)+ &
3257	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3258	dq_dwn(:,k)= dq_dwn(:,k)+ &
3259	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3260	ENDDO
3261	ELSEIF (iflag_wake==3) THEN
3262	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3263	DO k = 1,klev
3264	DO i=1,klon
3265	IF (rneb(i,k)==0.) THEN
3266	! On ne tient compte des tendances qu'en dehors des
3267	! nuages (c'est-\`a-dire a priri dans une region ou
3268	! l'eau se reevapore).
3269	dt_dwn(i,k)= dt_dwn(i,k)+ &
3270	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3271	dq_dwn(i,k)= dq_dwn(i,k)+ &
3272	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3273	ENDIF
3274	ENDDO
3275	ENDDO
3276	ENDIF
3277
3278	!
3279	!calcul caracteristiques de la poche froide
3280	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3281	t_seri, q_seri, omega, &
3282	dt_dwn, dq_dwn, M_dwn, M_up, &
3283	dt_a, dq_a, cv_gen, &
3284	sigd, cin, &
3285	wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
3286	wake_dth, wake_h, &
3287	!! wake_pe, wake_fip, wake_gfl, &
3288	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3289	d_t_wake, d_q_wake, &
3290	wake_k, t_x, q_x, &
3291	wake_omgbdth, wake_dp_omgb, &
3292	wake_dtKE, wake_dqKE, &
3293	wake_omg, wake_dp_deltomg, &
3294	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3295	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
3296	!
3297	!jyg Reinitialize itapwk when wakes have been called
3298	itapwk = 0
3299	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3300	!
3301	itapwk = itapwk+1
3302	!
3303	!-----------------------------------------------------------------------
3304	! ajout des tendances des poches froides
3305	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
3306	abortphy,flag_inhib_tend,itap,0)
3307	CALL prt_enerbil('wake',itap)
3308	!------------------------------------------------------------------------
3309
3310	! Increment Wake state variables
3311	IF (iflag_wake_tend .GT. 0.) THEN
3312
3313	CALL add_wake_tend &
3314	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
3315	'wake', abortphy)
3316	CALL prt_enerbil('wake',itap)
3317	ENDIF ! (iflag_wake_tend .GT. 0.)
3318	!
3319	IF (prt_level .GE. 10) THEN
3320	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3321	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3322	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3323	ENDIF
3324
3325	IF (iflag_alp_wk_cond .GT. 0.) THEN
3326
3327	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3328	wake_fip)
3329	ELSE
3330	wake_fip(:) = wake_fip_0(:)
3331	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3332
3333	ENDIF ! (iflag_wake>=1)
3334	!
3335	!===================================================================
3336	! Convection seche (thermiques ou ajustement)
3337	!===================================================================
3338	!
3339	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3340	,seuil_inversion,weak_inversion,dthmin)
3341
3342
3343
3344	d_t_ajsb(:,:)=0.
3345	d_q_ajsb(:,:)=0.
3346	d_t_ajs(:,:)=0.
3347	d_u_ajs(:,:)=0.
3348	d_v_ajs(:,:)=0.
3349	d_q_ajs(:,:)=0.
3350	clwcon0th(:,:)=0.
3351	!
3352	! fm_therm(:,:)=0.
3353	! entr_therm(:,:)=0.
3354	! detr_therm(:,:)=0.
3355	!
3356	IF (prt_level>9) WRITE(lunout,*) &
3357	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3358	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3359	IF (iflag_thermals<0) THEN
3360	! Rien
3361	! ====
3362	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3363
3364
3365	ELSE
3366
3367	! Thermiques
3368	! ==========
3369	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3370	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3371
3372
3373	!cc nrlmd le 10/04/2012
3374	DO k=1,klev+1
3375	DO i=1,klon
3376	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3377	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3378	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3379	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3380	ENDDO
3381	ENDDO
3382	!cc fin nrlmd le 10/04/2012
3383
3384	IF (iflag_thermals>=1) THEN
3385	!jyg<
3386	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3387	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3388	! Appel des thermiques avec les profils exterieurs aux poches
3389	DO k=1,klev
3390	DO i=1,klon
3391	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3392	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3393	u_therm(i,k) = u_seri(i,k)
3394	v_therm(i,k) = v_seri(i,k)
3395	ENDDO
3396	ENDDO
3397	ELSE
3398	! Appel des thermiques avec les profils moyens
3399	DO k=1,klev
3400	DO i=1,klon
3401	t_therm(i,k) = t_seri(i,k)
3402	q_therm(i,k) = q_seri(i,k)
3403	u_therm(i,k) = u_seri(i,k)
3404	v_therm(i,k) = v_seri(i,k)
3405	ENDDO
3406	ENDDO
3407	ENDIF
3408	!>jyg
3409	CALL calltherm(pdtphys &
3410	,pplay,paprs,pphi,weak_inversion &
3411	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3412	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
3413	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3414	,fm_therm,entr_therm,detr_therm &
3415	,zqasc,clwcon0th,lmax_th,ratqscth &
3416	,ratqsdiff,zqsatth &
3417	!on rajoute ale et alp, et les
3418	!caracteristiques de la couche alim
3419	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3420	,ztv,zpspsk,ztla,zthl &
3421	!cc nrlmd le 10/04/2012
3422	,pbl_tke_input,pctsrf,omega,cell_area &
3423	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3424	,n2,s2,ale_bl_stat &
3425	,therm_tke_max,env_tke_max &
3426	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3427	,alp_bl_conv,alp_bl_stat &
3428	!cc fin nrlmd le 10/04/2012
3429	,zqla,ztva )
3430	!
3431	!jyg<
3432	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3433	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3434	! Si les thermiques ne sont presents que hors des
3435	! poches, la tendance moyenne associ\'ee doit etre
3436	! multipliee par la fraction surfacique qu'ils couvrent.
3437	DO k=1,klev
3438	DO i=1,klon
3439	!
3440	d_deltat_the(i,k) = - d_t_ajs(i,k)
3441	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3442	!
3443	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3444	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3445	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3446	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3447	!
3448	ENDDO
3449	ENDDO
3450	!
3451	IF (ok_bug_split_th) THEN
3452	CALL add_wake_tend &
3453	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3454	ELSE
3455	CALL add_wake_tend &
3456	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3457	ENDIF
3458	CALL prt_enerbil('the',itap)
3459	!
3460	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3461	!
3462	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3463	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3464	CALL prt_enerbil('thermals',itap)
3465	!
3466	!
3467	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3468	cin, s2, n2, &
3469	ale_bl_trig, ale_bl_stat, ale_bl, &
3470	alp_bl, alp_bl_stat, &
3471	proba_notrig, random_notrig, cv_gen)
3472	!>jyg
3473
3474	! ------------------------------------------------------------------
3475	! Transport de la TKE par les panaches thermiques.
3476	! FH : 2010/02/01
3477	! if (iflag_pbl.eq.10) then
3478	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3479	! s rg,paprs,pbl_tke)
3480	! endif
3481	! -------------------------------------------------------------------
3482
3483	DO i=1,klon
3484	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3485	!CR:04/05/12:correction calcul zmax
3486	zmax_th(i)=zmax0(i)
3487	ENDDO
3488
3489	ENDIF
3490
3491	! Ajustement sec
3492	! ==============
3493
3494	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3495	! a partir du sommet des thermiques.
3496	! Dans le cas contraire, on demarre au niveau 1.
3497
3498	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3499
3500	IF (iflag_thermals.eq.0) THEN
3501	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3502	limbas(:)=1
3503	ELSE
3504	limbas(:)=lmax_th(:)
3505	ENDIF
3506
3507	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3508	! pour des test de convergence numerique.
3509	! Le nouveau ajsec est a priori mieux, meme pour le cas
3510	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3511	! non nulles numeriquement pour des mailles non concernees.
3512
3513	IF (iflag_thermals==0) THEN
3514	! Calling adjustment alone (but not the thermal plume model)
3515	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3516	, d_t_ajsb, d_q_ajsb)
3517	ELSE IF (iflag_thermals>0) THEN
3518	! Calling adjustment above the top of thermal plumes
3519	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3520	, d_t_ajsb, d_q_ajsb)
3521	ENDIF
3522
3523	!--------------------------------------------------------------------
3524	! ajout des tendances de l'ajustement sec ou des thermiques
3525	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3526	'ajsb',abortphy,flag_inhib_tend,itap,0)
3527	CALL prt_enerbil('ajsb',itap)
3528	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3529	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3530
3531	!---------------------------------------------------------------------
3532
3533	ENDIF
3534
3535	ENDIF
3536	!
3537	!===================================================================
3538	! Computation of ratqs, the width (normalized) of the subrid scale
3539	! water distribution
3540
3541	tke_dissip_ave(:,:)=0.
3542	l_mix_ave(:,:)=0.
3543	wprime_ave(:,:)=0.
3544
3545	DO nsrf = 1, nbsrf
3546	DO i = 1, klon
3547	tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
3548	l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3549	wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3550	ENDDO
3551	ENDDO
3552
3553	CALL calcratqs(klon,klev,prt_level,lunout, &
3554	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3555	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3556	tau_ratqs,fact_cldcon,wake_s, wake_deltaq, &
3557	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3558	paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
3559	pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
3560	ratqs,ratqsc,ratqs_inter)
3561
3562	!
3563	! Appeler le processus de condensation a grande echelle
3564	! et le processus de precipitation
3565	!-------------------------------------------------------------------------
3566	IF (prt_level .GE.10) THEN
3567	print *,'itap, ->fisrtilp ',itap
3568	ENDIF
3569	!
3570
3571	picefra(:,:)=0.
3572
3573	IF (ok_new_lscp) THEN
3574
3575	!--mise à jour de flight_m et flight_h2o dans leur module
3576	IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3577	CALL airplane(debut,pphis,pplay,paprs,t_seri)
3578	ENDIF
3579
3580	CALL lscp(phys_tstep,missing_val,paprs,pplay, &
3581	t_seri, q_seri,ptconv,ratqs, &
3582	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, &
3583	cldliq, picefra, rain_lsc, snow_lsc, &
3584	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3585	frac_impa, frac_nucl, beta_prec_fisrt, &
3586	prfl, psfl, rhcl, &
3587	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3588	iflag_ice_thermo, ok_ice_sursat)
3589
3590	ELSE
3591
3592	CALL fisrtilp(phys_tstep,paprs,pplay, &
3593	t_seri, q_seri,ptconv,ratqs, &
3594	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3595	rain_lsc, snow_lsc, &
3596	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3597	frac_impa, frac_nucl, beta_prec_fisrt, &
3598	prfl, psfl, rhcl, &
3599	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3600	iflag_ice_thermo)
3601
3602	ENDIF
3603	!
3604	WHERE (rain_lsc < 0) rain_lsc = 0.
3605	WHERE (snow_lsc < 0) snow_lsc = 0.
3606
3607	!+JLD
3608	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3609	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3610	! & ,rain_lsc,snow_lsc
3611	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3612	!-JLD
3613	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3614	'lsc',abortphy,flag_inhib_tend,itap,0)
3615	CALL prt_enerbil('lsc',itap)
3616	rain_num(:)=0.
3617	DO k = 1, klev
3618	DO i = 1, klon
3619	IF (ql_seri(i,k)>oliqmax) THEN
3620	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3621	ql_seri(i,k)=oliqmax
3622	ENDIF
3623	ENDDO
3624	ENDDO
3625	IF (nqo==3) THEN
3626	DO k = 1, klev
3627	DO i = 1, klon
3628	IF (qs_seri(i,k)>oicemax) THEN
3629	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3630	qs_seri(i,k)=oicemax
3631	ENDIF
3632	ENDDO
3633	ENDDO
3634	ENDIF
3635
3636	!---------------------------------------------------------------------------
3637	DO k = 1, klev
3638	DO i = 1, klon
3639	cldfra(i,k) = rneb(i,k)
3640	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3641	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3642	ENDDO
3643	ENDDO
3644	IF (check) THEN
3645	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3646	WRITE(lunout,*)"apresilp=", za
3647	zx_t = 0.0
3648	za = 0.0
3649	DO i = 1, klon
3650	za = za + cell_area(i)/REAL(klon)
3651	zx_t = zx_t + (rain_lsc(i) &
3652	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3653	ENDDO
3654	zx_t = zx_t/za*phys_tstep
3655	WRITE(lunout,*)"Precip=", zx_t
3656	ENDIF
3657
3658	IF (mydebug) THEN
3659	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3660	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3661	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3662	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3663	ENDIF
3664
3665	!
3666	!-------------------------------------------------------------------
3667	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3668	!-------------------------------------------------------------------
3669
3670	! 1. NUAGES CONVECTIFS
3671	!
3672	!IM cf FH
3673	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3674	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3675	snow_tiedtke=0.
3676	! print*,'avant calcul de la pseudo precip '
3677	! print*,'iflag_cld_th',iflag_cld_th
3678	IF (iflag_cld_th.eq.-1) THEN
3679	rain_tiedtke=rain_con
3680	ELSE
3681	! print*,'calcul de la pseudo precip '
3682	rain_tiedtke=0.
3683	! print*,'calcul de la pseudo precip 0'
3684	DO k=1,klev
3685	DO i=1,klon
3686	IF (d_q_con(i,k).lt.0.) THEN
3687	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3688	*(paprs(i,k)-paprs(i,k+1))/rg
3689	ENDIF
3690	ENDDO
3691	ENDDO
3692	ENDIF
3693	!
3694	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3695	!
3696
3697	! Nuages diagnostiques pour Tiedtke
3698	CALL diagcld1(paprs,pplay, &
3699	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3700	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3701	diafra,dialiq)
3702	DO k = 1, klev
3703	DO i = 1, klon
3704	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3705	cldliq(i,k) = dialiq(i,k)
3706	cldfra(i,k) = diafra(i,k)
3707	ENDIF
3708	ENDDO
3709	ENDDO
3710
3711	ELSE IF (iflag_cld_th.ge.3) THEN
3712	! On prend pour les nuages convectifs le max du calcul de la
3713	! convection et du calcul du pas de temps precedent diminue d'un facteur
3714	! facttemps
3715	facteur = pdtphys *facttemps
3716	DO k=1,klev
3717	DO i=1,klon
3718	rnebcon(i,k)=rnebcon(i,k)*facteur
3719	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3720	rnebcon(i,k)=rnebcon0(i,k)
3721	clwcon(i,k)=clwcon0(i,k)
3722	ENDIF
3723	ENDDO
3724	ENDDO
3725
3726	! On prend la somme des fractions nuageuses et des contenus en eau
3727
3728	IF (iflag_cld_th>=5) THEN
3729
3730	DO k=1,klev
3731	ptconvth(:,k)=fm_therm(:,k+1)>0.
3732	ENDDO
3733
3734	IF (iflag_coupl==4) THEN
3735
3736	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3737	! convectives et lsc dans la partie des thermiques
3738	! Le controle par iflag_coupl est peut etre provisoire.
3739	DO k=1,klev
3740	DO i=1,klon
3741	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3742	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3743	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3744	ELSE IF (ptconv(i,k)) THEN
3745	cldfra(i,k)=rnebcon(i,k)
3746	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3747	ENDIF
3748	ENDDO
3749	ENDDO
3750
3751	ELSE IF (iflag_coupl==5) THEN
3752	DO k=1,klev
3753	DO i=1,klon
3754	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3755	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3756	ENDDO
3757	ENDDO
3758
3759	ELSE
3760
3761	! Si on est sur un point touche par la convection
3762	! profonde et pas par les thermiques, on prend la
3763	! couverture nuageuse et l'eau nuageuse de la convection
3764	! profonde.
3765
3766	!IM/FH: 2011/02/23
3767	! definition des points sur lesquels ls thermiques sont actifs
3768
3769	DO k=1,klev
3770	DO i=1,klon
3771	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3772	cldfra(i,k)=rnebcon(i,k)
3773	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3774	ENDIF
3775	ENDDO
3776	ENDDO
3777
3778	ENDIF
3779
3780	ELSE
3781
3782	! Ancienne version
3783	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3784	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3785	ENDIF
3786
3787	ENDIF
3788
3789	! plulsc(:)=0.
3790	! do k=1,klev,-1
3791	! do i=1,klon
3792	! zzz=prfl(:,k)+psfl(:,k)
3793	! if (.not.ptconvth.zzz.gt.0.)
3794	! enddo prfl, psfl,
3795	! enddo
3796	!
3797	! 2. NUAGES STARTIFORMES
3798	!
3799	IF (ok_stratus) THEN
3800	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3801	DO k = 1, klev
3802	DO i = 1, klon
3803	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3804	cldliq(i,k) = dialiq(i,k)
3805	cldfra(i,k) = diafra(i,k)
3806	ENDIF
3807	ENDDO
3808	ENDDO
3809	ENDIF
3810	!
3811	! Precipitation totale
3812	!
3813	DO i = 1, klon
3814	rain_fall(i) = rain_con(i) + rain_lsc(i)
3815	snow_fall(i) = snow_con(i) + snow_lsc(i)
3816	ENDDO
3817	!
3818	! Calculer l'humidite relative pour diagnostique
3819	!
3820	DO k = 1, klev
3821	DO i = 1, klon
3822	zx_t = t_seri(i,k)
3823	IF (thermcep) THEN
3824	!! if (iflag_ice_thermo.eq.0) then !jyg
3825	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3826	!! else !jyg
3827	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3828	!! endif !jyg
3829	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3830	zx_qs = MIN(0.5,zx_qs)
3831	zcor = 1./(1.-retv*zx_qs)
3832	zx_qs = zx_qs*zcor
3833	ELSE
3834	!! IF (zx_t.LT.t_coup) THEN !jyg
3835	IF (zx_t.LT.rtt) THEN !jyg
3836	zx_qs = qsats(zx_t)/pplay(i,k)
3837	ELSE
3838	zx_qs = qsatl(zx_t)/pplay(i,k)
3839	ENDIF
3840	ENDIF
3841	zx_rh(i,k) = q_seri(i,k)/zx_qs
3842	IF (iflag_ice_thermo .GT. 0) THEN
3843	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
3844	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
3845	ENDIF
3846	zqsat(i,k)=zx_qs
3847	ENDDO
3848	ENDDO
3849
3850	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3851	! equivalente a 2m (tpote) pour diagnostique
3852	!
3853	DO i = 1, klon
3854	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3855	IF (thermcep) THEN
3856	IF(zt2m(i).LT.RTT) then
3857	Lheat=RLSTT
3858	ELSE
3859	Lheat=RLVTT
3860	ENDIF
3861	ELSE
3862	IF (zt2m(i).LT.RTT) THEN
3863	Lheat=RLSTT
3864	ELSE
3865	Lheat=RLVTT
3866	ENDIF
3867	ENDIF
3868	tpote(i) = tpot(i)* &
3869	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3870	ENDDO
3871
3872	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN ! ModThL
3873	#ifdef INCA
3874	CALL VTe(VTphysiq)
3875	CALL VTb(VTinca)
3876	calday = REAL(days_elapsed + 1) + jH_cur
3877
3878	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
3879	CALL AEROSOL_METEO_CALC( &
3880	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3881	prfl,psfl,pctsrf,cell_area, &
3882	latitude_deg,longitude_deg,u10m,v10m)
3883
3884	zxsnow_dummy(:) = 0.0
3885
3886	CALL chemhook_begin (calday, &
3887	days_elapsed+1, &
3888	jH_cur, &
3889	pctsrf(1,1), &
3890	latitude_deg, &
3891	longitude_deg, &
3892	cell_area, &
3893	paprs, &
3894	pplay, &
3895	coefh(1:klon,1:klev,is_ave), &
3896	pphi, &
3897	t_seri, &
3898	u, &
3899	v, &
3900	rot, &
3901	wo(:, :, 1), &
3902	q_seri, &
3903	zxtsol, &
3904	zt2m, &
3905	zxsnow_dummy, &
3906	solsw, &
3907	albsol1, &
3908	rain_fall, &
3909	snow_fall, &
3910	itop_con, &
3911	ibas_con, &
3912	cldfra, &
3913	nbp_lon, &
3914	nbp_lat-1, &
3915	tr_seri(:,:,1+nqCO2:nbtr), &
3916	ftsol, &
3917	paprs, &
3918	cdragh, &
3919	cdragm, &
3920	pctsrf, &
3921	pdtphys, &
3922	itap)
3923
3924	CALL VTe(VTinca)
3925	CALL VTb(VTphysiq)
3926	#endif
3927	ENDIF !type_trac = inca or inco
3928	IF (type_trac == 'repr') THEN
3929	#ifdef REPROBUS
3930	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
3931	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
3932	#endif
3933	ENDIF
3934
3935	!
3936	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3937	!
3938	IF (MOD(itaprad,radpas).EQ.0) THEN
3939
3940	!
3941	!jq - introduce the aerosol direct and first indirect radiative forcings
3942	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3943	IF (flag_aerosol .GT. 0) THEN
3944	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3945	IF (.NOT. aerosol_couple) THEN
3946	!
3947	CALL readaerosol_optic( &
3948	debut, flag_aerosol, itap, jD_cur-jD_ref, &
3949	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3950	mass_solu_aero, mass_solu_aero_pi, &
3951	tau_aero, piz_aero, cg_aero, &
3952	tausum_aero, tau3d_aero)
3953	ENDIF
3954	ELSE ! RRTM radiation
3955	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3956	abort_message='config_inca=aero et rrtm=1 impossible'
3957	CALL abort_physic(modname,abort_message,1)
3958	ELSE
3959	!
3960	#ifdef CPP_RRTM
3961	IF (NSW.EQ.6) THEN
3962	!--new aerosol properties SW and LW
3963	!
3964	#ifdef CPP_Dust
3965	!--SPL aerosol model
3966	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3967	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3968	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3969	tausum_aero, tau3d_aero)
3970	#else
3971	!--climatologies or INCA aerosols
3972	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
3973	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3974	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3975	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3976	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3977	tausum_aero, drytausum_aero, tau3d_aero)
3978	#endif
3979
3980	IF (flag_aerosol .EQ. 7) THEN
3981	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
3982	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3983	ENDIF
3984
3985	!
3986	ELSE IF (NSW.EQ.2) THEN
3987	!--for now we use the old aerosol properties
3988	!
3989	CALL readaerosol_optic( &
3990	debut, flag_aerosol, itap, jD_cur-jD_ref, &
3991	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3992	mass_solu_aero, mass_solu_aero_pi, &
3993	tau_aero, piz_aero, cg_aero, &
3994	tausum_aero, tau3d_aero)
3995	!
3996	!--natural aerosols
3997	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3998	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3999	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4000	!--all aerosols
4001	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4002	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4003	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
4004	!
4005	!--no LW optics
4006	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4007	!
4008	ELSE
4009	abort_message='Only NSW=2 or 6 are possible with ' &
4010	// 'aerosols and iflag_rrtm=1'
4011	CALL abort_physic(modname,abort_message,1)
4012	ENDIF
4013	#else
4014	abort_message='You should compile with -rrtm if running ' &
4015	// 'with iflag_rrtm=1'
4016	CALL abort_physic(modname,abort_message,1)
4017	#endif
4018	!
4019	ENDIF
4020	ENDIF
4021	ELSE !--flag_aerosol = 0
4022	tausum_aero(:,:,:) = 0.
4023	drytausum_aero(:,:) = 0.
4024	mass_solu_aero(:,:) = 0.
4025	mass_solu_aero_pi(:,:) = 0.
4026	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4027	tau_aero(:,:,:,:) = 1.e-15
4028	piz_aero(:,:,:,:) = 1.
4029	cg_aero(:,:,:,:) = 0.
4030	ELSE
4031	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4032	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4033	piz_aero_sw_rrtm(:,:,:,:) = 1.0
4034	cg_aero_sw_rrtm(:,:,:,:) = 0.0
4035	ENDIF
4036	ENDIF
4037	!
4038	!--WMO criterion to determine tropopause
4039	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4040	!
4041	!--STRAT AEROSOL
4042	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
4043	IF (flag_aerosol_strat.GT.0) THEN
4044	IF (prt_level .GE.10) THEN
4045	PRINT *,'appel a readaerosolstrat', mth_cur
4046	ENDIF
4047	IF (iflag_rrtm.EQ.0) THEN
4048	IF (flag_aerosol_strat.EQ.1) THEN
4049	CALL readaerosolstrato(debut)
4050	ELSE
4051	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4052	CALL abort_physic(modname,abort_message,1)
4053	ENDIF
4054	ELSE
4055	#ifdef CPP_RRTM
4056	#ifndef CPP_StratAer
4057	!--prescribed strat aerosols
4058	!--only in the case of non-interactive strat aerosols
4059	IF (flag_aerosol_strat.EQ.1) THEN
4060	CALL readaerosolstrato1_rrtm(debut)
4061	ELSEIF (flag_aerosol_strat.EQ.2) THEN
4062	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
4063	ELSE
4064	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4065	CALL abort_physic(modname,abort_message,1)
4066	ENDIF
4067	#endif
4068	#else
4069	abort_message='You should compile with -rrtm if running ' &
4070	// 'with iflag_rrtm=1'
4071	CALL abort_physic(modname,abort_message,1)
4072	#endif
4073	ENDIF
4074	ELSE
4075	tausum_aero(:,:,id_STRAT_phy) = 0.
4076	ENDIF
4077	!
4078	#ifdef CPP_RRTM
4079	#ifdef CPP_StratAer
4080	!--compute stratospheric mask
4081	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4082	!--interactive strat aerosols
4083	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4084	#endif
4085	#endif
4086	!--fin STRAT AEROSOL
4087	!
4088
4089	! Calculer les parametres optiques des nuages et quelques
4090	! parametres pour diagnostiques:
4091	!
4092	IF (aerosol_couple.AND.config_inca=='aero') THEN
4093	mass_solu_aero(:,:) = ccm(:,:,1)
4094	mass_solu_aero_pi(:,:) = ccm(:,:,2)
4095	ENDIF
4096
4097	IF (ok_newmicro) then
4098	! AI IF (iflag_rrtm.NE.0) THEN
4099	IF (iflag_rrtm.EQ.1) THEN
4100	#ifdef CPP_RRTM
4101	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
4102	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
4103	// 'pour ok_cdnc'
4104	CALL abort_physic(modname,abort_message,1)
4105	ENDIF
4106	#else
4107
4108	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
4109	CALL abort_physic(modname,abort_message,1)
4110	#endif
4111	ENDIF
4112	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
4113	paprs, pplay, t_seri, cldliq, picefra, cldfra, &
4114	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4115	flwp, fiwp, flwc, fiwc, &
4116	mass_solu_aero, mass_solu_aero_pi, &
4117	cldtaupi, re, fl, ref_liq, ref_ice, &
4118	ref_liq_pi, ref_ice_pi)
4119	ELSE
4120	CALL nuage (paprs, pplay, &
4121	t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
4122	cldh, cldl, cldm, cldt, cldq, &
4123	ok_aie, &
4124	mass_solu_aero, mass_solu_aero_pi, &
4125	bl95_b0, bl95_b1, &
4126	cldtaupi, re, fl)
4127	ENDIF
4128	!
4129	!IM betaCRF
4130	!
4131	cldtaurad = cldtau
4132	cldtaupirad = cldtaupi
4133	cldemirad = cldemi
4134	cldfrarad = cldfra
4135
4136	!
4137	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4138	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4139	!
4140	! global
4141	!
4142	!IM 251017 begin
4143	! print*,'physiq betaCRF global zdtime=',zdtime
4144	!IM 251017 end
4145	DO k=1, klev
4146	DO i=1, klon
4147	IF (pplay(i,k).GE.pfree) THEN
4148	beta(i,k) = beta_pbl
4149	ELSE
4150	beta(i,k) = beta_free
4151	ENDIF
4152	IF (mskocean_beta) THEN
4153	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4154	ENDIF
4155	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4156	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4157	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4158	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4159	ENDDO
4160	ENDDO
4161	!
4162	ELSE
4163	!
4164	! regional
4165	!
4166	DO k=1, klev
4167	DO i=1,klon
4168	!
4169	IF (longitude_deg(i).ge.lon1_beta.AND. &
4170	longitude_deg(i).le.lon2_beta.AND. &
4171	latitude_deg(i).le.lat1_beta.AND. &
4172	latitude_deg(i).ge.lat2_beta) THEN
4173	IF (pplay(i,k).GE.pfree) THEN
4174	beta(i,k) = beta_pbl
4175	ELSE
4176	beta(i,k) = beta_free
4177	ENDIF
4178	IF (mskocean_beta) THEN
4179	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4180	ENDIF
4181	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4182	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4183	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4184	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4185	ENDIF
4186	!
4187	ENDDO
4188	ENDDO
4189	!
4190	ENDIF
4191
4192	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4193	IF (ok_chlorophyll) THEN
4194	print*,"-- reading chlorophyll"
4195	CALL readchlorophyll(debut)
4196	ENDIF
4197
4198	!--if ok_suntime_rrtm we use ancillay data for RSUN
4199	!--previous values are therefore overwritten
4200	!--this is needed for CMIP6 runs
4201	!--and only possible for new radiation scheme
4202	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4203	#ifdef CPP_RRTM
4204	CALL read_rsun_rrtm(debut)
4205	#endif
4206	ENDIF
4207
4208	IF (mydebug) THEN
4209	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4210	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4211	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4212	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4213	ENDIF
4214
4215	!
4216	!sonia : If Iflag_radia >=2, pertubation of some variables
4217	!input to radiation (DICE)
4218	!
4219	IF (iflag_radia .ge. 2) THEN
4220	zsav_tsol (:) = zxtsol(:)
4221	CALL perturb_radlwsw(zxtsol,iflag_radia)
4222	ENDIF
4223
4224	IF (aerosol_couple.AND.config_inca=='aero') THEN
4225	#ifdef INCA
4226	CALL radlwsw_inca &
4227	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4228	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4229	size(wo,3), wo, &
4230	cldfrarad, cldemirad, cldtaurad, &
4231	heat,heat0,cool,cool0,albpla, &
4232	topsw,toplw,solsw,sollw, &
4233	sollwdown, &
4234	topsw0,toplw0,solsw0,sollw0, &
4235	lwdn0, lwdn, lwup0, lwup, &
4236	swdn0, swdn, swup0, swup, &
4237	ok_ade, ok_aie, &
4238	tau_aero, piz_aero, cg_aero, &
4239	topswad_aero, solswad_aero, &
4240	topswad0_aero, solswad0_aero, &
4241	topsw_aero, topsw0_aero, &
4242	solsw_aero, solsw0_aero, &
4243	cldtaupirad, &
4244	topswai_aero, solswai_aero)
4245	#endif
4246	ELSE
4247	!
4248	!IM calcul radiatif pour le cas actuel
4249	!
4250	RCO2 = RCO2_act
4251	RCH4 = RCH4_act
4252	RN2O = RN2O_act
4253	RCFC11 = RCFC11_act
4254	RCFC12 = RCFC12_act
4255	!
4256	!--interactive CO2 in ppm from carbon cycle
4257	IF (carbon_cycle_rad.AND..NOT.debut) THEN
4258	RCO2=RCO2_glo
4259	ENDIF
4260	!
4261	IF (prt_level .GE.10) THEN
4262	print *,' ->radlwsw, number 1 '
4263	ENDIF
4264	!
4265	CALL radlwsw &
4266	(dist, rmu0, fract, &
4267	!albedo SB >>>
4268	! paprs, pplay,zxtsol,albsol1, albsol2, &
4269	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4270	!albedo SB <<<
4271	t_seri,q_seri,wo, &
4272	cldfrarad, cldemirad, cldtaurad, &
4273	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4274	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4275	tau_aero, piz_aero, cg_aero, &
4276	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4277	! Rajoute par OB pour RRTM
4278	tau_aero_lw_rrtm, &
4279	cldtaupirad, &
4280	! zqsat, flwcrad, fiwcrad, &
4281	zqsat, flwc, fiwc, &
4282	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4283	heat,heat0,cool,cool0,albpla, &
4284	heat_volc,cool_volc, &
4285	topsw,toplw,solsw,solswfdiff,sollw, &
4286	sollwdown, &
4287	topsw0,toplw0,solsw0,sollw0, &
4288	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4289	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4290	topswad_aero, solswad_aero, &
4291	topswai_aero, solswai_aero, &
4292	topswad0_aero, solswad0_aero, &
4293	topsw_aero, topsw0_aero, &
4294	solsw_aero, solsw0_aero, &
4295	topswcf_aero, solswcf_aero, &
4296	!-C. Kleinschmitt for LW diagnostics
4297	toplwad_aero, sollwad_aero,&
4298	toplwai_aero, sollwai_aero, &
4299	toplwad0_aero, sollwad0_aero,&
4300	!-end
4301	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4302	ZSWFT0_i, ZFSDN0, ZFSUP0)
4303
4304	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4305	!schemes
4306	toplw = toplw + betalwoff * (toplw0 - toplw)
4307	sollw = sollw + betalwoff * (sollw0 - sollw)
4308	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4309	lwup = lwup + betalwoff * (lwup0 - lwup)
4310	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4311	sollwdown(:))
4312	cool = cool + betalwoff * (cool0 - cool)
4313
4314	#ifndef CPP_XIOS
4315	!--OB 30/05/2016 modified 21/10/2016
4316	!--here we return swaero_diag and dryaod_diag to FALSE
4317	!--and histdef will switch it back to TRUE if necessary
4318	!--this is necessary to get the right swaero at first step
4319	!--but only in the case of no XIOS as XIOS is covered elsewhere
4320	IF (debut) swaerofree_diag = .FALSE.
4321	IF (debut) swaero_diag = .FALSE.
4322	IF (debut) dryaod_diag = .FALSE.
4323	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4324	!--as for swaero_diag, see above
4325	IF (debut) ok_4xCO2atm = .FALSE.
4326
4327	!
4328	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4329	!IM des taux doit etre different du taux actuel
4330	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4331	!
4332	IF (RCO2_per.NE.RCO2_act.OR. &
4333	RCH4_per.NE.RCH4_act.OR. &
4334	RN2O_per.NE.RN2O_act.OR. &
4335	RCFC11_per.NE.RCFC11_act.OR. &
4336	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4337	#endif
4338	!
4339	IF (ok_4xCO2atm) THEN
4340	!
4341	RCO2 = RCO2_per
4342	RCH4 = RCH4_per
4343	RN2O = RN2O_per
4344	RCFC11 = RCFC11_per
4345	RCFC12 = RCFC12_per
4346	!
4347	IF (prt_level .GE.10) THEN
4348	print *,' ->radlwsw, number 2 '
4349	ENDIF
4350	!
4351	CALL radlwsw &
4352	(dist, rmu0, fract, &
4353	!albedo SB >>>
4354	! paprs, pplay,zxtsol,albsol1, albsol2, &
4355	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4356	!albedo SB <<<
4357	t_seri,q_seri,wo, &
4358	cldfrarad, cldemirad, cldtaurad, &
4359	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4360	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4361	tau_aero, piz_aero, cg_aero, &
4362	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4363	! Rajoute par OB pour RRTM
4364	tau_aero_lw_rrtm, &
4365	cldtaupi, &
4366	! zqsat, flwcrad, fiwcrad, &
4367	zqsat, flwc, fiwc, &
4368	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4369	heatp,heat0p,coolp,cool0p,albplap, &
4370	heat_volc,cool_volc, &
4371	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4372	sollwdownp, &
4373	topsw0p,toplw0p,solsw0p,sollw0p, &
4374	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4375	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4376	topswad_aerop, solswad_aerop, &
4377	topswai_aerop, solswai_aerop, &
4378	topswad0_aerop, solswad0_aerop, &
4379	topsw_aerop, topsw0_aerop, &
4380	solsw_aerop, solsw0_aerop, &
4381	topswcf_aerop, solswcf_aerop, &
4382	!-C. Kleinschmitt for LW diagnostics
4383	toplwad_aerop, sollwad_aerop,&
4384	toplwai_aerop, sollwai_aerop, &
4385	toplwad0_aerop, sollwad0_aerop,&
4386	!-end
4387	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4388	ZSWFT0_i, ZFSDN0, ZFSUP0)
4389	ENDIF !ok_4xCO2atm
4390	ENDIF ! aerosol_couple
4391	itaprad = 0
4392	!
4393	! If Iflag_radia >=2, reset pertubed variables
4394	!
4395	IF (iflag_radia .ge. 2) THEN
4396	zxtsol(:) = zsav_tsol (:)
4397	ENDIF
4398	ENDIF ! MOD(itaprad,radpas)
4399	itaprad = itaprad + 1
4400
4401	IF (iflag_radia.eq.0) THEN
4402	IF (prt_level.ge.9) THEN
4403	PRINT *,'--------------------------------------------------'
4404	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4405	PRINT *,'>>>> heat et cool mis a zero '
4406	PRINT *,'--------------------------------------------------'
4407	ENDIF
4408	heat=0.
4409	cool=0.
4410	sollw=0. ! MPL 01032011
4411	solsw=0.
4412	radsol=0.
4413	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4414	swup0=0.
4415	lwup=0.
4416	lwup0=0.
4417	lwdn=0.
4418	lwdn0=0.
4419	ENDIF
4420
4421	!
4422	! Calculer radsol a l'exterieur de radlwsw
4423	! pour prendre en compte le cycle diurne
4424	! recode par Olivier Boucher en sept 2015
4425	!
4426	radsol=solsw*swradcorr+sollw
4427
4428	IF (ok_4xCO2atm) THEN
4429	radsolp=solswp*swradcorr+sollwp
4430	ENDIF
4431
4432	!
4433	! Ajouter la tendance des rayonnements (tous les pas)
4434	! avec une correction pour le cycle diurne dans le SW
4435	!
4436
4437	DO k=1, klev
4438	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4439	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4440	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4441	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4442	ENDDO
4443
4444	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4445	CALL prt_enerbil('SW',itap)
4446	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4447	CALL prt_enerbil('LW',itap)
4448
4449	!
4450	IF (mydebug) THEN
4451	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4452	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4453	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4454	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4455	ENDIF
4456
4457	! Calculer l'hydrologie de la surface
4458	!
4459	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4460	! . agesno, ftsol,fqsurf,fsnow, ruis)
4461	!
4462
4463	!
4464	! Calculer le bilan du sol et la derive de temperature (couplage)
4465	!
4466	DO i = 1, klon
4467	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4468	! a la demande de JLD
4469	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4470	ENDDO
4471	!
4472	!moddeblott(jan95)
4473	! Appeler le programme de parametrisation de l'orographie
4474	! a l'echelle sous-maille:
4475	!
4476	IF (prt_level .GE.10) THEN
4477	print *,' call orography ? ', ok_orodr
4478	ENDIF
4479	!
4480	IF (ok_orodr) THEN
4481	!
4482	! selection des points pour lesquels le shema est actif:
4483	igwd=0
4484	DO i=1,klon
4485	itest(i)=0
4486	! IF ((zstd(i).gt.10.0)) THEN
4487	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4488	itest(i)=1
4489	igwd=igwd+1
4490	idx(igwd)=i
4491	ENDIF
4492	ENDDO
4493	! igwdim=MAX(1,igwd)
4494	!
4495	IF (ok_strato) THEN
4496
4497	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4498	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4499	igwd,idx,itest, &
4500	t_seri, u_seri, v_seri, &
4501	zulow, zvlow, zustrdr, zvstrdr, &
4502	d_t_oro, d_u_oro, d_v_oro)
4503
4504	ELSE
4505	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4506	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4507	igwd,idx,itest, &
4508	t_seri, u_seri, v_seri, &
4509	zulow, zvlow, zustrdr, zvstrdr, &
4510	d_t_oro, d_u_oro, d_v_oro)
4511	ENDIF
4512	!
4513	! ajout des tendances
4514	!-----------------------------------------------------------------------
4515	! ajout des tendances de la trainee de l'orographie
4516	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4517	abortphy,flag_inhib_tend,itap,0)
4518	CALL prt_enerbil('oro',itap)
4519	!----------------------------------------------------------------------
4520	!
4521	ENDIF ! fin de test sur ok_orodr
4522	!
4523	IF (mydebug) THEN
4524	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4525	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4526	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4527	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4528	ENDIF
4529
4530	IF (ok_orolf) THEN
4531	!
4532	! selection des points pour lesquels le shema est actif:
4533	igwd=0
4534	DO i=1,klon
4535	itest(i)=0
4536	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4537	itest(i)=1
4538	igwd=igwd+1
4539	idx(igwd)=i
4540	ENDIF
4541	ENDDO
4542	! igwdim=MAX(1,igwd)
4543	!
4544	IF (ok_strato) THEN
4545
4546	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4547	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4548	igwd,idx,itest, &
4549	t_seri, u_seri, v_seri, &
4550	zulow, zvlow, zustrli, zvstrli, &
4551	d_t_lif, d_u_lif, d_v_lif )
4552
4553	ELSE
4554	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4555	latitude_deg,zmea,zstd,zpic, &
4556	itest, &
4557	t_seri, u_seri, v_seri, &
4558	zulow, zvlow, zustrli, zvstrli, &
4559	d_t_lif, d_u_lif, d_v_lif)
4560	ENDIF
4561
4562	! ajout des tendances de la portance de l'orographie
4563	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4564	'lif', abortphy,flag_inhib_tend,itap,0)
4565	CALL prt_enerbil('lif',itap)
4566	ENDIF ! fin de test sur ok_orolf
4567
4568	IF (ok_hines) then
4569	! HINES GWD PARAMETRIZATION
4570	east_gwstress=0.
4571	west_gwstress=0.
4572	du_gwd_hines=0.
4573	dv_gwd_hines=0.
4574	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4575	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4576	du_gwd_hines, dv_gwd_hines)
4577	zustr_gwd_hines=0.
4578	zvstr_gwd_hines=0.
4579	DO k = 1, klev
4580	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4581	* (paprs(:, k)-paprs(:, k+1))/rg
4582	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4583	* (paprs(:, k)-paprs(:, k+1))/rg
4584	ENDDO
4585
4586	d_t_hin(:, :)=0.
4587	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4588	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4589	CALL prt_enerbil('hin',itap)
4590	ENDIF
4591
4592	IF (.not. ok_hines .and. ok_gwd_rando) then
4593	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4594	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4595	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4596	dv_gwd_front, east_gwstress, west_gwstress)
4597	zustr_gwd_front=0.
4598	zvstr_gwd_front=0.
4599	DO k = 1, klev
4600	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4601	* (paprs(:, k)-paprs(:, k+1))/rg
4602	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4603	* (paprs(:, k)-paprs(:, k+1))/rg
4604	ENDDO
4605
4606	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4607	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4608	CALL prt_enerbil('front_gwd_rando',itap)
4609	ENDIF
4610
4611	IF (ok_gwd_rando) THEN
4612	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4613	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4614	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4615	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4616	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4617	CALL prt_enerbil('flott_gwd_rando',itap)
4618	zustr_gwd_rando=0.
4619	zvstr_gwd_rando=0.
4620	DO k = 1, klev
4621	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4622	* (paprs(:, k)-paprs(:, k+1))/rg
4623	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4624	* (paprs(:, k)-paprs(:, k+1))/rg
4625	ENDDO
4626	ENDIF
4627
4628	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4629
4630	IF (mydebug) THEN
4631	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4632	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4633	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4634	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4635	ENDIF
4636
4637	DO i = 1, klon
4638	zustrph(i)=0.
4639	zvstrph(i)=0.
4640	ENDDO
4641	DO k = 1, klev
4642	DO i = 1, klon
4643	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4644	(paprs(i,k)-paprs(i,k+1))/rg
4645	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4646	(paprs(i,k)-paprs(i,k+1))/rg
4647	ENDDO
4648	ENDDO
4649	!
4650	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4651	!
4652	IF (is_sequential .and. ok_orodr) THEN
4653	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4654	ra,rg,romega, &
4655	latitude_deg,longitude_deg,pphis, &
4656	zustrdr,zustrli,zustrph, &
4657	zvstrdr,zvstrli,zvstrph, &
4658	paprs,u,v, &
4659	aam, torsfc)
4660	ENDIF
4661	!IM cf. FLott END
4662	!DC Calcul de la tendance due au methane
4663	IF (ok_qch4) THEN
4664	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4665	! ajout de la tendance d'humidite due au methane
4666	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4667	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4668	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4669	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4670	ENDIF
4671	!
4672	!
4673
4674	!===============================================================
4675	! Additional tendency of TKE due to orography
4676	!===============================================================
4677	!
4678	! Inititialization
4679	!------------------
4680
4681	addtkeoro=0
4682	CALL getin_p('addtkeoro',addtkeoro)
4683
4684	IF (prt_level.ge.5) &
4685	print*,'addtkeoro', addtkeoro
4686
4687	alphatkeoro=1.
4688	CALL getin_p('alphatkeoro',alphatkeoro)
4689	alphatkeoro=min(max(0.,alphatkeoro),1.)
4690
4691	smallscales_tkeoro=.FALSE.
4692	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4693
4694
4695	dtadd(:,:)=0.
4696	duadd(:,:)=0.
4697	dvadd(:,:)=0.
4698
4699	! Choices for addtkeoro:
4700	! ** 0 no TKE tendency from orography
4701	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4702	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4703	!
4704
4705	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4706	! -------------------------------------------
4707
4708
4709	! selection des points pour lesquels le schema est actif:
4710
4711
4712	IF (addtkeoro .EQ. 1 ) THEN
4713
4714	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4715	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4716
4717	ELSE IF (addtkeoro .EQ. 2) THEN
4718
4719	IF (smallscales_tkeoro) THEN
4720	igwd=0
4721	DO i=1,klon
4722	itest(i)=0
4723	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4724	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4725	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4726	IF (zstd(i).GT.1.0) THEN
4727	itest(i)=1
4728	igwd=igwd+1
4729	idx(igwd)=i
4730	ENDIF
4731	ENDDO
4732
4733	ELSE
4734
4735	igwd=0
4736	DO i=1,klon
4737	itest(i)=0
4738	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4739	itest(i)=1
4740	igwd=igwd+1
4741	idx(igwd)=i
4742	ENDIF
4743	ENDDO
4744
4745	ENDIF
4746
4747	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4748	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4749	igwd,idx,itest, &
4750	t_seri, u_seri, v_seri, &
4751	zulow, zvlow, zustrdr, zvstrdr, &
4752	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4753
4754	zustrdr(:)=0.
4755	zvstrdr(:)=0.
4756	zulow(:)=0.
4757	zvlow(:)=0.
4758
4759	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4760	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4761	ENDIF
4762
4763
4764	! TKE update from subgrid temperature and wind tendencies
4765	!----------------------------------------------------------
4766	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4767
4768
4769	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4770	!
4771	! Prevent pbl_tke_w from becoming negative
4772	wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
4773	!
4774
4775	ENDIF
4776	! -----
4777	!===============================================================
4778
4779
4780	!====================================================================
4781	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4782	!====================================================================
4783	! Abderrahmane 24.08.09
4784
4785	IF (ok_cosp) THEN
4786	! adeclarer
4787	#ifdef CPP_COSP
4788	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4789
4790	IF (prt_level .GE.10) THEN
4791	print*,'freq_cosp',freq_cosp
4792	ENDIF
4793	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4794	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4795	! s ref_liq,ref_ice
4796	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4797	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4798	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4799	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4800	JrNt,ref_liq,ref_ice, &
4801	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4802	zu10m,zv10m,pphis, &
4803	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4804	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4805	prfl(:,1:klev),psfl(:,1:klev), &
4806	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4807	mr_ozone,cldtau, cldemi)
4808
4809	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4810	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4811	! M clMISR,
4812	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4813	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4814
4815	ENDIF
4816	#endif
4817
4818	#ifdef CPP_COSP2
4819	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4820
4821	IF (prt_level .GE.10) THEN
4822	print*,'freq_cosp',freq_cosp
4823	ENDIF
4824	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4825	print*,'Dans physiq.F avant appel '
4826	! s ref_liq,ref_ice
4827	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4828	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4829	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4830	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4831	JrNt,ref_liq,ref_ice, &
4832	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4833	zu10m,zv10m,pphis, &
4834	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4835	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4836	prfl(:,1:klev),psfl(:,1:klev), &
4837	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4838	mr_ozone,cldtau, cldemi)
4839	ENDIF
4840	#endif
4841
4842	#ifdef CPP_COSPV2
4843	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4844	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4845
4846	IF (prt_level .GE.10) THEN
4847	print*,'freq_cosp',freq_cosp
4848	ENDIF
4849	DO k = 1, klev
4850	DO i = 1, klon
4851	phicosp(i,k) = pphi(i,k) + pphis(i)
4852	ENDDO
4853	ENDDO
4854	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4855	print*,'Dans physiq.F avant appel '
4856	! s ref_liq,ref_ice
4857	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4858	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4859	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4860	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4861	JrNt,ref_liq,ref_ice, &
4862	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4863	zu10m,zv10m,pphis, &
4864	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4865	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4866	prfl(:,1:klev),psfl(:,1:klev), &
4867	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4868	mr_ozone,cldtau, cldemi)
4869	ENDIF
4870	#endif
4871
4872	ENDIF !ok_cosp
4873
4874
4875	! Marine
4876
4877	IF (ok_airs) then
4878
4879	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4880	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4881	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4882	& map_prop_hc,map_prop_hist,&
4883	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4884	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4885	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4886	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4887	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4888	& map_ntot,map_hc,map_hist,&
4889	& map_Cb,map_ThCi,map_Anv,&
4890	& alt_tropo )
4891	ENDIF
4892
4893	ENDIF ! ok_airs
4894
4895
4896	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4897	!AA
4898	!AA Installation de l'interface online-offline pour traceurs
4899	!AA
4900	!====================================================================
4901	! Calcul des tendances traceurs
4902	!====================================================================
4903	!
4904
4905	IF (type_trac=='repr') THEN
4906	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
4907	!MM dans Reprobus
4908	sh_in(:,:) = q_seri(:,:)
4909	#ifdef REPROBUS
4910	d_q_rep(:,:) = 0.
4911	d_ql_rep(:,:) = 0.
4912	d_qi_rep(:,:) = 0.
4913	#endif
4914	ELSE
4915	sh_in(:,:) = qx(:,:,ivap)
4916	IF (nqo .EQ. 3) THEN
4917	ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
4918	ELSE
4919	ch_in(:,:) = qx(:,:,iliq)
4920	ENDIF
4921	ENDIF
4922
4923	#ifdef CPP_Dust
4924	! Avec SPLA, iflag_phytrac est forcé =1
4925	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4926	pdtphys,ftsol, & ! I
4927	t,q_seri,paprs,pplay,RHcl, & ! I
4928	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4929	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4930	u_seri, v_seri, latitude_deg, longitude_deg, &
4931	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4932	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4933	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4934	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4935	fm_therm, entr_therm, rneb, & ! I
4936	beta_prec_fisrt,beta_prec, & !I
4937	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4938	d_tr_dyn,tr_seri)
4939
4940	#else
4941	IF (iflag_phytrac == 1 ) THEN
4942	CALL phytrac ( &
4943	itap, days_elapsed+1, jH_cur, debut, &
4944	lafin, phys_tstep, u, v, t, &
4945	paprs, pplay, pmfu, pmfd, &
4946	pen_u, pde_u, pen_d, pde_d, &
4947	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4948	u1, v1, ftsol, pctsrf, &
4949	zustar, zu10m, zv10m, &
4950	wstar(:,is_ave), ale_bl, ale_wake, &
4951	latitude_deg, longitude_deg, &
4952	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4953	presnivs, pphis, pphi, albsol1, &
4954	sh_in, ch_in, rhcl, cldfra, rneb, &
4955	diafra, cldliq, itop_con, ibas_con, &
4956	pmflxr, pmflxs, prfl, psfl, &
4957	da, phi, mp, upwd, &
4958	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4959	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4960	ev, ep, epmlmMm, eplaMm, & !<<RomP
4961	dnwd, aerosol_couple, flxmass_w, &
4962	tau_aero, piz_aero, cg_aero, ccm, &
4963	rfname, &
4964	d_tr_dyn, & !<<RomP
4965	tr_seri, init_source)
4966	#ifdef REPROBUS
4967
4968
4969	print*,'avt add phys rep',abortphy
4970
4971	CALL add_phys_tend &
4972	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
4973	'rep',abortphy,flag_inhib_tend,itap,0)
4974	IF (abortphy==1) Print*,'ERROR ABORT REP'
4975
4976	print*,'apr add phys rep',abortphy
4977
4978	#endif
4979	ENDIF ! (iflag_phytrac=1)
4980
4981	#endif
4982	!ENDIF ! (iflag_phytrac=1)
4983
4984	IF (offline) THEN
4985
4986	IF (prt_level.ge.9) &
4987	print*,'Attention on met a 0 les thermiques pour phystoke'
4988	CALL phystokenc ( &
4989	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4990	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4991	fm_therm,entr_therm, &
4992	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4993	frac_impa, frac_nucl, &
4994	pphis,cell_area,phys_tstep,itap, &
4995	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4996
4997
4998	ENDIF
4999
5000	!
5001	! Calculer le transport de l'eau et de l'energie (diagnostique)
5002	!
5003	CALL transp (paprs,zxtsol, &
5004	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5005	ve, vq, ue, uq, vwat, uwat)
5006	!
5007	!IM global posePB BEG
5008	IF(1.EQ.0) THEN
5009	!
5010	CALL transp_lay (paprs,zxtsol, &
5011	t_seri, q_seri, u_seri, v_seri, zphi, &
5012	ve_lay, vq_lay, ue_lay, uq_lay)
5013	!
5014	ENDIF !(1.EQ.0) THEN
5015	!IM global posePB END
5016	! Accumuler les variables a stocker dans les fichiers histoire:
5017	!
5018
5019	!================================================================
5020	! Conversion of kinetic and potential energy into heat, for
5021	! parameterisation of subgrid-scale motions
5022	!================================================================
5023
5024	d_t_ec(:,:)=0.
5025	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5026	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
5027	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
5028	zmasse,exner,d_t_ec)
5029	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
5030
5031	!=======================================================================
5032	! SORTIES
5033	!=======================================================================
5034	!
5035	!IM initialisation + calculs divers diag AMIP2
5036	!
5037	include "calcul_divers.h"
5038	!
5039	!IM Interpolation sur les niveaux de pression du NMC
5040	! -------------------------------------------------
5041	!
5042	include "calcul_STDlev.h"
5043	!
5044	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5045	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5046	!
5047	!cc prw = eau precipitable
5048	! prlw = colonne eau liquide
5049	! prlw = colonne eau solide
5050	prw(:) = 0.
5051	prlw(:) = 0.
5052	prsw(:) = 0.
5053	DO k = 1, klev
5054	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
5055	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5056	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
5057	ENDDO
5058	!
5059	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
5060	#ifdef INCA
5061	CALL VTe(VTphysiq)
5062	CALL VTb(VTinca)
5063
5064	CALL chemhook_end ( &
5065	phys_tstep, &
5066	pplay, &
5067	t_seri, &
5068	tr_seri(:,:,1+nqCO2:nbtr), &
5069	nbtr, &
5070	paprs, &
5071	q_seri, &
5072	cell_area, &
5073	pphi, &
5074	pphis, &
5075	zx_rh, &
5076	aps, bps, ap, bp)
5077
5078	CALL VTe(VTinca)
5079	CALL VTb(VTphysiq)
5080	#endif
5081	ENDIF
5082
5083
5084	!
5085	! Convertir les incrementations en tendances
5086	!
5087	IF (prt_level .GE.10) THEN
5088	print *,'Convertir les incrementations en tendances '
5089	ENDIF
5090	!
5091	IF (mydebug) THEN
5092	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5093	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5094	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5095	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5096	ENDIF
5097
5098	DO k = 1, klev
5099	DO i = 1, klon
5100	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5101	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5102	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5103	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5104	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
5105	!CR: on ajoute le contenu en glace
5106	IF (nqo.gt.3) THEN
5107	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5108	ENDIF
5109	!--ice_sursat: nqo=4, on ajoute rneb
5110	IF (nqo.eq.4) THEN
5111	d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
5112	ENDIF
5113	ENDDO
5114	ENDDO
5115	!
5116	IF (nqtot > nqo) THEN
5117	itr = 0
5118	DO iq = 1, nqtot
5119	IF(tracers(iq)%isH2Ofamily) CYCLE
5120	itr = itr+1
5121	DO k = 1, klev
5122	DO i = 1, klon
5123	d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
5124	ENDDO
5125	ENDDO
5126	ENDDO
5127	ENDIF
5128	!
5129	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5130	!IM global posePB include "write_bilKP_ins.h"
5131	!IM global posePB include "write_bilKP_ave.h"
5132	!
5133
5134	!--OB mass fixer
5135	!--profile is corrected to force mass conservation of water
5136	IF (mass_fixer) THEN
5137	qql2(:)=0.0
5138	DO k = 1, klev
5139	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
5140	ENDDO
5141	DO i = 1, klon
5142	!--compute ratio of what q+ql should be with conservation to what it is
5143	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5144	DO k = 1, klev
5145	q_seri(i,k) =q_seri(i,k)*corrqql
5146	ql_seri(i,k)=ql_seri(i,k)*corrqql
5147	ENDDO
5148	ENDDO
5149	ENDIF
5150	!--fin mass fixer
5151
5152	! Sauvegarder les valeurs de t et q a la fin de la physique:
5153	!
5154	u_ancien(:,:) = u_seri(:,:)
5155	v_ancien(:,:) = v_seri(:,:)
5156	t_ancien(:,:) = t_seri(:,:)
5157	q_ancien(:,:) = q_seri(:,:)
5158	ql_ancien(:,:) = ql_seri(:,:)
5159	qs_ancien(:,:) = qs_seri(:,:)
5160	rneb_ancien(:,:) = rneb_seri(:,:)
5161	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5162	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5163	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5164	! !! RomP >>>
5165	IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
5166	! !! RomP <<<
5167	!==========================================================================
5168	! Sorties des tendances pour un point particulier
5169	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5170	! pour le debug
5171	! La valeur de igout est attribuee plus haut dans le programme
5172	!==========================================================================
5173
5174	IF (prt_level.ge.1) THEN
5175	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5176	write(lunout,*) &
5177	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5178	write(lunout,*) &
5179	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5180	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5181	pctsrf(igout,is_sic)
5182	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5183	DO k=1,klev
5184	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5185	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5186	d_t_eva(igout,k)
5187	ENDDO
5188	write(lunout,*) 'cool,heat'
5189	DO k=1,klev
5190	write(lunout,*) cool(igout,k),heat(igout,k)
5191	ENDDO
5192
5193	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5194	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5195	!jyg! do k=1,klev
5196	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5197	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5198	!jyg! enddo
5199	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5200	DO k=1,klev
5201	write(lunout,*) d_t_vdf(igout,k), &
5202	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5203	ENDDO
5204	!>jyg
5205
5206	write(lunout,*) 'd_ps ',d_ps(igout)
5207	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5208	DO k=1,klev
5209	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5210	d_qx(igout,k,1),d_qx(igout,k,2)
5211	ENDDO
5212	ENDIF
5213
5214	!============================================================
5215	! Calcul de la temperature potentielle
5216	!============================================================
5217	DO k = 1, klev
5218	DO i = 1, klon
5219	!JYG/IM theta en debut du pas de temps
5220	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5221	!JYG/IM theta en fin de pas de temps de physique
5222	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5223	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5224	! MPL 20130625
5225	! fth_fonctions.F90 et parkind1.F90
5226	! sinon thetal=theta
5227	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5228	! : ql_seri(i,k))
5229	thetal(i,k)=theta(i,k)
5230	ENDDO
5231	ENDDO
5232	!
5233
5234	! 22.03.04 BEG
5235	!=============================================================
5236	! Ecriture des sorties
5237	!=============================================================
5238	#ifdef CPP_IOIPSL
5239
5240	! Recupere des varibles calcule dans differents modules
5241	! pour ecriture dans histxxx.nc
5242
5243	! Get some variables from module fonte_neige_mod
5244	CALL fonte_neige_get_vars(pctsrf, &
5245	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5246
5247
5248	!=============================================================
5249	! Separation entre thermiques et non thermiques dans les sorties
5250	! de fisrtilp
5251	!=============================================================
5252
5253	IF (iflag_thermals>=1) THEN
5254	d_t_lscth=0.
5255	d_t_lscst=0.
5256	d_q_lscth=0.
5257	d_q_lscst=0.
5258	DO k=1,klev
5259	DO i=1,klon
5260	IF (ptconvth(i,k)) THEN
5261	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5262	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5263	ELSE
5264	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5265	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5266	ENDIF
5267	ENDDO
5268	ENDDO
5269
5270	DO i=1,klon
5271	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5272	plul_th(i)=prfl(i,1)+psfl(i,1)
5273	ENDDO
5274	ENDIF
5275
5276	!On effectue les sorties:
5277
5278	#ifdef CPP_Dust
5279	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5280	pplay, lmax_th, aerosol_couple, &
5281	ok_ade, ok_aie, ivap, ok_sync, &
5282	ptconv, read_climoz, clevSTD, &
5283	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5284	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5285	#else
5286	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5287	pplay, lmax_th, aerosol_couple, &
5288	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, &
5289	ok_sync, ptconv, read_climoz, clevSTD, &
5290	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5291	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5292	#endif
5293
5294	#ifndef CPP_XIOS
5295	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5296	#endif
5297
5298	#endif
5299
5300	!====================================================================
5301	! Arret du modele apres hgardfou en cas de detection d'un
5302	! plantage par hgardfou
5303	!====================================================================
5304
5305	IF (abortphy==1) THEN
5306	abort_message ='Plantage hgardfou'
5307	CALL abort_physic (modname,abort_message,1)
5308	ENDIF
5309
5310	! 22.03.04 END
5311	!
5312	!====================================================================
5313	! Si c'est la fin, il faut conserver l'etat de redemarrage
5314	!====================================================================
5315	!
5316
5317	! Disabling calls to the prt_alerte function
5318	alert_first_call = .FALSE.
5319
5320	IF (lafin) THEN
5321	itau_phy = itau_phy + itap
5322	CALL phyredem ("restartphy.nc")
5323	! open(97,form="unformatted",file="finbin")
5324	! write(97) u_seri,v_seri,t_seri,q_seri
5325	! close(97)
5326
5327	IF (is_omp_master) THEN
5328
5329	IF (read_climoz >= 1) THEN
5330	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5331	DEALLOCATE(press_edg_climoz) ! pointer
5332	DEALLOCATE(press_cen_climoz) ! pointer
5333	ENDIF
5334
5335	ENDIF
5336	#ifdef CPP_XIOS
5337	IF (is_omp_master) CALL xios_context_finalize
5338	#endif
5339	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5340	ENDIF
5341
5342	! first=.false.
5343
5344	END SUBROUTINE physiq
5345
5346	END MODULE physiq_mod

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

Download in other formats:

Original Format