Context Navigation

physiq_mod.F90 @ 2795

Last change on this file since 2795 was 2795, checked in by Laurent Fairhead, 7 years ago
Upgrading time counter in the physics to a count of seconds from the start of the run. The old counter made the model diverge when runs made in 4 months chunks were compared to runs made in 1 month chunks
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 154.0 KB

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

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

Download in other formats:

Original Format