Context Navigation

physiq.F90 @ 2101

Last change on this file since 2101 was 2100, checked in by lguez, 10 years ago

Removed "on rentre dans guide_main" from guide_main in dyn3dpar, was
already commented out in the dyn3dmem version.

Keeping length of lines under 80 characters in physiq (for
readability). Removed wrong comments "ajout des tendances de la
diffusion turbulente". Replaced "con" by "convection" as an argument
of add_phys_tend.

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: 125.6 KB

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

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

Download in other formats:

Original Format