Context Navigation

physiq_mod.F90 @ 3179

Last change on this file since 3179 was 3179, checked in by jyg, 6 years ago

New structure for the representation of vdf
splitting in pbl_surface. Computations are
gathered in two subroutines included in the module
wx_pbl_mod: wx_pbl0_fuse, called before the
sub-surfaces, determines the single column
equivalent to the (w) and (x) columns.
wx_pbl0_split, called after the subsurfaces,
determines the distinct (w) and (x) surface
fluxes.

This is a first version with no surface
temperature difference between (w) and (x) (hence
the index 0).

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

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

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

Download in other formats:

Original Format