Context Navigation

physiq_mod.F90 @ 2799

Last change on this file since 2799 was 2799, checked in by jyg, 7 years ago

New diagnostics to verify energy conservation.
Some corrections to improve energy conservation.
Some bug correction.

New output variables and new flags are

introduced:

(1)New variables: d_xx_col is the budget over

each column of variable "xx".

(2) fl_ebil: 0 -> no diag; 1 -> diags activated
(3) fl_cor_ebil: 0 -> no correction ensuring

energy conservation; 1 -> first set of
corrections.

+ ok_bad_ecmwf_thermo: if true, a bug setting

RVTMP2 to 0 is fixed; for backward compatibility
the default is "false".

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.0 KB

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

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

Download in other formats:

Original Format