Context Navigation

physiq.F90 @ 2201

Last change on this file since 2201 was 2200, checked in by acozic, 9 years ago
Aco : add some variable in inca routine call to correct a bug on time_axis with xios output
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: 133.3 KB

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

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

Download in other formats:

Original Format