Context Navigation

physiq_mod.F90 @ 2793

Last change on this file since 2793 was 2788, checked in by dcugnet, 8 years ago

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

records 2-13 are obtained with records 1-12 of "climoz.nc".
records 1 and 14 are obtained respectively with:
- record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
- record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.

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: 153.9 KB

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

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

Download in other formats:

Original Format