Context Navigation

physiq_mod.F90 @ 2753

Last change on this file since 2753 was 2753, checked in by oboucher, 7 years ago
Adding aerosol properties to Dust version
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 152.8 KB

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

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

Download in other formats:

Original Format