Context Navigation

physiq_mod.F90 @ 2684

Last change on this file since 2684 was 2684, checked in by acozic, 8 years ago
Modification to allow read_climoz=1 or 2 when we run a configuration with INCA (LMDZORINCA)
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: 157.3 KB

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

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

Download in other formats:

Original Format