Context Navigation

physiq_mod.F90 @ 4143

Last change on this file since 4143 was 4143, checked in by dcugnet, 2 years ago
Some variables are renamed or replaced by direct equivalents: iso_indnum -> tracers(:)%iso_iName niso_possibles -> niso iqiso -> iqIsoPha ; index_trac -> itZonIso ok_iso_verif -> isoCheck ntraceurs_zone -> nzone ; ntraciso -> ntiso qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times. Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18). The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
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: 190.7 KB

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

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

Download in other formats:

Original Format