Context Navigation

physiq.F90 @ 2031

Last change on this file since 2031 was 2009, checked in by Laurent Fairhead, 11 years ago

Rajout de la clé CPP_RRTM pour la compilation avec RRTM

Added CPP_RRTM key for RRTM compilation

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

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

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

Download in other formats:

Original Format