Context Navigation

physiq_mod.F90 @ 2812

Last change on this file since 2812 was 2812, checked in by jyg, 7 years ago
Addition of a diagnostic mode for add_phys_tend, with a new argument diag_mode (=0 for normal use; =1 for diagnostic use).
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: 155.9 KB

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

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

Download in other formats:

Original Format