Context Navigation

physiq_mod.F90 @ 2695

Last change on this file since 2695 was 2692, checked in by oboucher, 8 years ago
Adding a call to stratosphere_mask in the case of StratAer?. A lot of cosmetic changes on if/endif do/enddo and calls
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: 157.2 KB

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

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

Download in other formats:

Original Format