Context Navigation

physiq_mod.F90 @ 3352

Last change on this file since 3352 was 3352, checked in by acozic, 6 years ago
merge with rev 3351 on the trunk (bug fixed on a call to Inca)
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: 170.2 KB

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

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

Download in other formats:

Original Format