Context Navigation

physiq_mod.F90 @ 3992

Last change on this file since 3992 was 3989, checked in by oboucher, 3 years ago
addition of flag_volc_surfstrat required for VOLMIP This flag can select either surface or the atmospheric effects of volcanic aerosols to separate the effects.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 186.0 KB

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

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

Download in other formats:

Original Format