Context Navigation

physiq_mod.F90 @ 3756

Last change on this file since 3756 was 3756, checked in by oboucher, 4 years ago
The fraction of SW down radiation that is diffuse is extracted from radlwsw.F90 and passed on to pbl_surface for inclusion in fields_out for transfer to ORCHIDEE if it requested in input file coupling_fields.def. Hence there is nothing automatic here. The fraction of SW down radiation is also added as a diagnostic. To satisfy the case when this new quantity is used in ORCHIDEE, then it is added in the restart (and read in case it is there).
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: 185.0 KB

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

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

Download in other formats:

Original Format