Context Navigation

physiq_mod.F90 @ 3817

Last change on this file since 3817 was 3817, checked in by musat, 3 years ago
Nouvelle formulation des calculs a 2m
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 185.4 KB

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

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

Download in other formats:

Original Format