Context Navigation

physiq.F90 @ 2343

Last change on this file since 2343 was 2343, checked in by Ehouarn Millour, 9 years ago
Another step towards a clean separation between physics and dynamics: adapted read_pstoke.F90, read_pstoke0.F90 initphysto.F90 and phystokenc.F90 (now module phystokenc_mod.F90) to not explicitely include/use file/modules from the dynamics. In the process, added module "time_phylmdz_mod.F90" in the physics, which contains the information otherwise found in "temps.h" (which is in the dynamics) and should be used instead. EM
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: 141.9 KB

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

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

Download in other formats:

Original Format