Context Navigation

physiq_mod.F90 @ 3607

Last change on this file since 3607 was 3607, checked in by jghattas, 4 years ago
Revert commit [3604]. This is only done temporary because we need to make a correction for Jean-Zay and to have a revision without these changes to go into IPSLCM6.1.10.
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: 172.4 KB

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

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

Download in other formats:

Original Format