Context Navigation

physiq_mod.F90 @ 3356

Last change on this file since 3356 was 3356, checked in by Laurent Fairhead, 6 years ago
First attempt at merging with trunk
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 173.3 KB

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

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

Download in other formats:

Original Format