Context Navigation

physiq_mod.F90 @ 3893

Last change on this file since 3893 was 3888, checked in by jyg, 4 years ago

New provisional version of the splitting of the
diffusive boundary layer into inwake and offwake
PBLs. The splitting of the diffuse BL should NOT
be activated yet for general purpose simulations.

The splitting is activated by:
mod(iflag_pbl_split,10)=1 for the option with
fixed surface temperature and
mod(iflag_pbl_split,10)=2 for the option with
coupled surface temperature.

iflag_pbl_split=0 ==> no splittingat all.
iflag_pbl_split=10 ==> splitting of thermals.
iflag_pbl_split=11 ==> splitting of thermals and
of vertical diffusion (fixed surf. temp.).
iflag_pbl_split=12 ==> splitting of thermals and
of vertical diffusion (coupled surf. temp.).

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: 186.0 KB

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

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

Download in other formats:

Original Format