Context Navigation

physiq.F90 @ 2157

Last change on this file since 2157 was 2073, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r2054:2070 into testing branch
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: 123.9 KB

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

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

Download in other formats:

Original Format