Context Navigation

physiq_mod.F90 @ 3762

Last change on this file since 3762 was 3754, checked in by adurocher, 4 years ago
Fix conditional compilation to compile without XIOS
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 `Author Date Id Revision`
File size: 185.9 KB

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

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

Download in other formats:

Original Format