Context Navigation

physiq_mod.F90 @ 2952

Last change on this file since 2952 was 2914, checked in by acozic, 7 years ago
Add an argument to chemini call ( = rev 2906 on trunk LMDZ5)
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: 160.8 KB

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

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

Download in other formats:

Original Format