Context Navigation

physiq.F @ 1189

Last change on this file since 1189 was 705, checked in by Laurent Fairhead, 18 years ago
Pour la conservation du flux d'eau OM, JLD LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 89.1 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
	4	C
	5	c $Header$
	6	c
	7	SUBROUTINE physiq (nlon,nlev,nqmax,
	8	. debut,lafin,rjourvrai,gmtime,pdtphys,
	9	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
	10	. u,v,t,qx,
	11	. omega,
	12	#ifdef INCA_CH4
	13	. flxmass_w,
	14	#endif
	15	. d_u, d_v, d_t, d_qx, d_ps)
	16
	17	USE ioipsl
	18	USE histcom
	19	#ifdef INCA
	20	USE chemshut
	21	#ifdef INCA_CH4
	22	! USE obs_pos
	23	#endif
	24	#endif
	25	IMPLICIT none
	26	c======================================================================
	27	c
	28	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	29	c
	30	c Objet: Moniteur general de la physique du modele
	31	cAA Modifications quant aux traceurs :
	32	cAA - uniformisation des parametrisations ds phytrac
	33	cAA - stockage des moyennes des champs necessaires
	34	cAA en mode traceur off-line
	35	c======================================================================
	36	c CLEFS CPP POUR LES IO
	37	c =====================
	38	c #define histhf
	39	#define histday
	40	#define histmth
	41	#define histins
	42	c #define histISCCP
	43	c #define histREGDYN
	44	c #define histmthNMC
	45	c======================================================================
	46	c modif ( P. Le Van , 12/10/98 )
	47	c
	48	c Arguments:
	49	c
	50	c nlon----input-I-nombre de points horizontaux
	51	c nlev----input-I-nombre de couches verticales
	52	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
	53	c debut---input-L-variable logique indiquant le premier passage
	54	c lafin---input-L-variable logique indiquant le dernier passage
	55	c rjour---input-R-numero du jour de l'experience
	56	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
	57	c pdtphys-input-R-pas d'integration pour la physique (seconde)
	58	c paprs---input-R-pression pour chaque inter-couche (en Pa)
	59	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	60	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	61	c pphis---input-R-geopotentiel du sol
	62	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	63	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
	64	c v-------input-R-vitesse Y (de S a N) en m/s
	65	c t-------input-R-temperature (K)
	66	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	67	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	68	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
	69	c omega---input-R-vitesse verticale en Pa/s
	70	c
	71	c d_u-----output-R-tendance physique de "u" (m/s/s)
	72	c d_v-----output-R-tendance physique de "v" (m/s/s)
	73	c d_t-----output-R-tendance physique de "t" (K/s)
	74	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	75	c d_ps----output-R-tendance physique de la pression au sol
	76	c======================================================================
	77	#include "dimensions.h"
	78	integer jjmp1
	79	parameter (jjmp1=jjm+1-1/jjm)
	80	#include "dimphy.h"
	81	#include "regdim.h"
	82	#include "indicesol.h"
	83	#include "dimsoil.h"
	84	#include "clesphys.h"
	85	#include "control.h"
	86	#include "temps.h"
	87	#include "comgeomphy.h"
	88	#include "advtrac.h"
	89	#include "iniprint.h"
	90	c======================================================================
	91	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	92	PARAMETER (ok_cvl=.TRUE.)
	93	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	94	PARAMETER (ok_gust=.FALSE.)
	95	c======================================================================
	96	LOGICAL check ! Verifier la conservation du modele en eau
	97	PARAMETER (check=.FALSE.)
	98	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	99	PARAMETER (ok_stratus=.FALSE.)
	100	c======================================================================
	101	c Parametres lies au coupleur OASIS:
	102	#include "oasis.h"
	103	INTEGER,SAVE :: npas, nexca
	104	logical rnpb
	105	#ifdef INCA
	106	parameter(rnpb=.false.)
	107	#else
	108	parameter(rnpb=.true.)
	109	#endif
	110	c ocean = type de modele ocean a utiliser: force, slab, couple
	111	character*6 ocean
	112	SAVE ocean
	113
	114	c parameter (ocean = 'force ')
	115	c parameter (ocean = 'couple')
	116	logical ok_ocean
	117	c======================================================================
	118	c Clef controlant l'activation du cycle diurne:
	119	ccc LOGICAL cycle_diurne
	120	ccc PARAMETER (cycle_diurne=.FALSE.)
	121	c======================================================================
	122	c Modele thermique du sol, a activer pour le cycle diurne:
	123	ccc LOGICAL soil_model
	124	ccc PARAMETER (soil_model=.FALSE.)
	125	logical ok_veget
	126	save ok_veget
	127	c parameter (ok_veget = .true.)
	128	c parameter (ok_veget = .false.)
	129	c======================================================================
	130	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	131	c le calcul du rayonnement est celle apres la precipitation des nuages.
	132	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	133	c la condensation et la precipitation. Cette cle augmente les impacts
	134	c radiatifs des nuages.
	135	ccc LOGICAL new_oliq
	136	ccc PARAMETER (new_oliq=.FALSE.)
	137	c======================================================================
	138	c Clefs controlant deux parametrisations de l'orographie:
	139	cc LOGICAL ok_orodr
	140	ccc PARAMETER (ok_orodr=.FALSE.)
	141	ccc LOGICAL ok_orolf
	142	ccc PARAMETER (ok_orolf=.FALSE.)
	143	c======================================================================
	144	LOGICAL ok_journe ! sortir le fichier journalier
	145	save ok_journe
	146	c PARAMETER (ok_journe=.true.)
	147	c
	148	LOGICAL ok_mensuel ! sortir le fichier mensuel
	149	save ok_mensuel
	150	c PARAMETER (ok_mensuel=.true.)
	151	c
	152	LOGICAL ok_instan ! sortir le fichier instantane
	153	save ok_instan
	154	c PARAMETER (ok_instan=.true.)
	155	c
	156	LOGICAL ok_region ! sortir le fichier regional
	157	PARAMETER (ok_region=.FALSE.)
	158	c======================================================================
	159	c
	160	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	161	PARAMETER (ivap=1)
	162	INTEGER iliq ! indice de traceurs pour eau liquide
	163	PARAMETER (iliq=2)
	164
	165	c
	166	c
	167	c Variables argument:
	168	c
	169	INTEGER nlon
	170	INTEGER nlev
	171	INTEGER nqmax
	172	REAL rjourvrai
	173	REAL gmtime
	174	REAL pdtphys
	175	LOGICAL debut, lafin
	176	REAL paprs(klon,klev+1)
	177	REAL pplay(klon,klev)
	178	REAL pphi(klon,klev)
	179	REAL pphis(klon)
	180	REAL presnivs(klev)
	181	REAL znivsig(klev)
	182	REAL zsurf(nbsrf)
	183
	184	REAL u(klon,klev)
	185	REAL v(klon,klev)
	186	REAL t(klon,klev)
	187	REAL qx(klon,klev,nqmax)
	188
	189	REAL t_ancien(klon,klev), q_ancien(klon,klev)
	190	SAVE t_ancien, q_ancien
	191	LOGICAL ancien_ok
	192	SAVE ancien_ok
	193
	194	REAL d_t_dyn(klon,klev)
	195	REAL d_q_dyn(klon,klev)
	196
	197	REAL omega(klon,klev)
	198
	199	#ifdef INCA_CH4
	200	REAL flxmass_w(klon,klev)
	201	#endif
	202	REAL d_u(klon,klev)
	203	REAL d_v(klon,klev)
	204	REAL d_t(klon,klev)
	205	REAL d_qx(klon,klev,nqmax)
	206	REAL d_ps(klon)
	207
	208	INTEGER klevp1, klevm1
	209	PARAMETER(klevp1=klev+1,klevm1=klev-1)
	210	#include "raddim.h"
	211	c
	212	cIM 080304 REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
	213	REAL swdn0(klon,klevp1), swdn(klon,klevp1)
	214	REAL swup0(klon,klevp1), swup(klon,klevp1)
	215	SAVE swdn0 , swdn, swup0, swup
	216	c
	217	REAL SWdn200clr(klon), SWdn200(klon)
	218	REAL SWup200clr(klon), SWup200(klon)
	219	SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
	220	c
	221	REAL lwdn0(klon,klevp1), lwdn(klon,klevp1)
	222	REAL lwup0(klon,klevp1), lwup(klon,klevp1)
	223	SAVE lwdn0 , lwdn, lwup0, lwup
	224	c
	225	REAL LWdn200clr(klon), LWdn200(klon)
	226	REAL LWup200clr(klon), LWup200(klon)
	227	SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
	228	c
	229	REAL LWdnTOA(klon), LWdnTOAclr(klon)
	230	SAVE LWdnTOA, LWdnTOAclr
	231	c
	232	c vents meridien et zonal a un niveau de pression
	233	c
	234	integer nlevSTD
	235	PARAMETER(nlevSTD=17)
	236	real rlevSTD(nlevSTD)
	237	DATA rlevSTD/100000., 92500., 85000., 70000.,
	238	.60000., 50000., 40000., 30000., 25000., 20000.,
	239	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
	240	CHARACTER*5 clevSTD(nlevSTD), aa, bb
	241	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
	242	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
	243	.'70 ','50 ','30 ','20 ','10 '/
	244	c
	245	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
	246	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
	247	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
	248	c
	249	cIM ENSEMBLES BEG
	250	c
	251	integer nlevENS
	252	PARAMETER(nlevENS=4)
	253	integer indENS(nlevENS)
	254	save indENS
	255	real rlevENS(nlevENS)
	256	DATA rlevENS/85000., 70000., 50000., 20000./
	257	CHARACTER*3 clev(nlevENS)
	258	DATA clev/'850','700','500','200'/
	259
	260	real tlev(klon,nlevENS), qlev(klon,nlevENS), rhlev(klon,nlevENS)
	261	real ulev(klon,nlevENS), vlev(klon,nlevENS), philev(klon,nlevENS)
	262	real wlev(klon,nlevENS)
	263	cIM ENSEMBLES END
	264	c
	265	c prw: precipitable water
	266	real prw(klon)
	267
	268	REAL convliq(klon,klev) ! eau liquide nuageuse convective
	269	REAL convfra(klon,klev) ! fraction nuageuse convective
	270
	271	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
	272	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
	273	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
	274	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
	275
	276	INTEGER linv, kp1
	277	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
	278	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
	279	REAL flwp(klon), fiwp(klon)
	280	REAL flwc(klon,klev), fiwc(klon,klev)
	281	REAL flwp_c(klon), fiwp_c(klon)
	282	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
	283	REAL flwp_s(klon), fiwp_s(klon)
	284	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
	285
	286	c ISCCP simulator v3.4
	287	c dans clesphys.h top_height, overlap
	288	cv3.4
	289	INTEGER debug, debugcol
	290	INTEGER npoints
	291	PARAMETER(npoints=klon)
	292	c
	293	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
	294	INTEGER nregISCtot
	295	PARAMETER(nregISCtot=1)
	296	c
	297	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
	298	c y compris pour 1 point
	299	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
	300	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
	301	INTEGER imin_debut, nbpti
	302	INTEGER jmin_debut, nbptj
	303	c
	304	REAL nbsunlit(nregISCtot,klon) !nbsunlit : moyenne de sunlit
	305	INTEGER ncol, seed(klon)
	306
	307	c ncol = nb. de sous-colonnes pour chaque maille du GCM
	308	c PARAMETER(ncol=100)
	309	c PARAMETER(ncol=625)
	310	c PARAMETER(ncol=10)
	311	PARAMETER(ncol=25)
	312	REAL tautab(0:255)
	313	INTEGER invtau(-20:45000)
	314	REAL emsfc_lw
	315	PARAMETER(emsfc_lw=0.99)
	316	REAL ran0 ! type for random number fuction
	317	c
	318	REAL cldtot(klon,klev)
	319	c variables de haut en bas pour le simulateur ISCCP
	320	REAL dtau_s(klon,klev) !tau nuages startiformes
	321	REAL dtau_c(klon,klev) !tau nuages convectifs
	322	REAL dem_s(klon,klev) !emissivite nuages startiformes
	323	REAL dem_c(klon,klev) !emissivite nuages convectifs
	324	c
	325	c variables de haut en bas pour le simulateur ISCCP
	326	REAL pfull(klon,klev)
	327	REAL phalf(klon,klev+1)
	328	REAL qv(klon,klev)
	329	REAL cc(klon,klev)
	330	REAL conv(klon,klev)
	331	REAL dtau_sH2B(klon,klev)
	332	REAL dtau_cH2B(klon,klev)
	333	REAL at(klon,klev)
	334	REAL dem_sH2B(klon,klev)
	335	REAL dem_cH2B(klon,klev)
	336
	337	c output from ISCCP simulator
	338	REAL fq_isccp(klon,7,7)
	339	REAL totalcldarea(klon)
	340	REAL meanptop(klon)
	341	REAL meantaucld(klon)
	342	REAL boxtau(klon,ncol)
	343	REAL boxptop(klon,ncol)
	344	c
	345	INTEGER l, ni, nj, kmax, lmax
	346	PARAMETER(kmax=8, lmax=8)
	347	INTEGER kmaxm1, lmaxm1
	348	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
	349	INTEGER iimx7, jjmx7, jjmp1x7
	350	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
	351	.jjmp1x7=jjmp1*lmaxm1)
	352	REAL fq4d(iim,jjmp1,kmaxm1,lmaxm1)
	353	REAL fq3d(iimx7, jjmp1x7)
	354	c
	355	INTEGER iw, iwmax
	356	REAL wmin, pas_w
	357	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
	358	PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
	359	REAL o500(klon)
	360	c
	361	cIM: nbregdyn = nbre regions pour calculs statistiques sur output du ISCCP
	362	cIM: dynamiques
	363	INTEGER nreg, nbregdyn
	364	PARAMETER(nbregdyn=5)
	365	REAL histoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
	366	REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
	367	REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbregdyn)
	368	SAVE nhistoWt
	369
	370	INTEGER linv
	371	INTEGER pct_ocean(klon,nbregdyn)
	372	REAL rlonPOS(klon)
	373
	374	c sorties ISCCP
	375
	376	logical ok_isccp
	377	real ecrit_isccp
	378	integer nid_isccp
	379	save ok_isccp, ecrit_isccp, nid_isccp
	380
	381	#ifdef histISCCP
	382	data ok_isccp/.true./
	383	#else
	384	data ok_isccp/.false./
	385	#endif
	386
	387	c sorties statistiques regime dynamique
	388	logical ok_regdyn
	389	real ecrit_regdyn
	390	integer nid_regdyn
	391	save ok_regdyn, ecrit_regdyn, nid_regdyn
	392
	393	#ifdef histREGDYN
	394	c data ok_regdyn,ecrit_regdyn/.true.,0.125/
	395	c data ok_regdyn,ecrit_regdyn/.true.,1./
	396	data ok_regdyn/.true./
	397	#else
	398	data ok_regdyn/.false./
	399	#endif
	400
	401	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
	402	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
	403	DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
	404
	405	c cldtopres pression au sommet des nuages
	406	REAL cldtopres(lmaxm1)
	407	DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
	408
	409	INTEGER komega, nhoriRD
	410
	411	c taulev: numero du niveau de tau dans les sorties ISCCP
	412	CHARACTER *4 taulev(kmaxm1)
	413	DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
	414
	415	REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
	416	INTEGER nhorix7
	417	cIM: region='3d' <==> sorties en global
	418	CHARACTER*3 region
	419	PARAMETER(region='3d')
	420	c
	421	logical ok_hf
	422	real ecrit_hf
	423	integer nid_hf, nid_hf3d
	424	save ok_hf, ecrit_hf, nid_hf, nid_hf3d
	425
	426	c QUESTION : noms de variables ?
	427
	428	#ifdef histhf
	429	data ok_hf,ecrit_hf/.true.,0.25/
	430	#else
	431	data ok_hf/.false./
	432	#endif
	433
	434	INTEGER longcles
	435	PARAMETER ( longcles = 20 )
	436	REAL clesphy0( longcles )
	437	c
	438	c Variables quasi-arguments
	439	c
	440	REAL xjour
	441	SAVE xjour
	442	c
	443	c
	444	c Variables propres a la physique
	445	c
	446	REAL dtime
	447	SAVE dtime ! pas temporel de la physique
	448	c
	449	INTEGER radpas
	450	SAVE radpas ! frequence d'appel rayonnement
	451	c
	452	REAL radsol(klon)
	453	SAVE radsol ! bilan radiatif au sol calcule par code radiatif
	454	c
	455	REAL rlat(klon)
	456	SAVE rlat ! latitude pour chaque point
	457	c
	458	REAL rlon(klon)
	459	SAVE rlon ! longitude pour chaque point
	460	c
	461	cc INTEGER iflag_con
	462	cc SAVE iflag_con ! indicateur de la convection
	463	c
	464	INTEGER itap
	465	SAVE itap ! compteur pour la physique
	466	c
	467	REAL co2_ppm_etat0
	468	c
	469	REAL solaire_etat0
	470	c
	471	real slp(klon) ! sea level pressure
	472
	473	REAL ftsol(klon,nbsrf)
	474	SAVE ftsol ! temperature du sol
	475	c
	476	REAL ftsoil(klon,nsoilmx,nbsrf)
	477	SAVE ftsoil ! temperature dans le sol
	478	c
	479	REAL fevap(klon,nbsrf)
	480	SAVE fevap ! evaporation
	481	REAL fluxlat(klon,nbsrf)
	482	SAVE fluxlat
	483	c
	484	REAL deltat(klon)
	485	SAVE deltat ! ecart avec la SST de reference
	486	c
	487	REAL fqsurf(klon,nbsrf)
	488	SAVE fqsurf ! humidite de l'air au contact de la surface
	489	c
	490	REAL qsol(klon)
	491	SAVE qsol ! hauteur d'eau dans le sol
	492	c
	493	REAL fsnow(klon,nbsrf)
	494	SAVE fsnow ! epaisseur neigeuse
	495	c
	496	REAL falbe(klon,nbsrf)
	497	SAVE falbe ! albedo par type de surface
	498	REAL falblw(klon,nbsrf)
	499	SAVE falblw ! albedo par type de surface
	500
	501	c
	502	c
	503	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	504	c
	505	REAL zmea(klon)
	506	SAVE zmea ! orographie moyenne
	507	c
	508	REAL zstd(klon)
	509	SAVE zstd ! deviation standard de l'OESM
	510	c
	511	REAL zsig(klon)
	512	SAVE zsig ! pente de l'OESM
	513	c
	514	REAL zgam(klon)
	515	save zgam ! anisotropie de l'OESM
	516	c
	517	REAL zthe(klon)
	518	SAVE zthe ! orientation de l'OESM
	519	c
	520	REAL zpic(klon)
	521	SAVE zpic ! Maximum de l'OESM
	522	c
	523	REAL zval(klon)
	524	SAVE zval ! Minimum de l'OESM
	525	c
	526	REAL rugoro(klon)
	527	SAVE rugoro ! longueur de rugosite de l'OESM
	528	c
	529	REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	530	c
	531	REAL zuthe(klon),zvthe(klon)
	532	SAVE zuthe
	533	SAVE zvthe
	534	INTEGER igwd,idx(klon),itest(klon)
	535	c
	536	REAL agesno(klon,nbsrf)
	537	SAVE agesno ! age de la neige
	538	c
	539	REAL alb_neig(klon)
	540	SAVE alb_neig ! albedo de la neige
	541	c
	542	REAL run_off_lic_0(klon)
	543	SAVE run_off_lic_0
	544	cKE43
	545	c Variables liees a la convection de K. Emanuel (sb):
	546	c
	547	REAL ema_workcbmf(klon) ! cloud base mass flux
	548	SAVE ema_workcbmf
	549
	550	REAL ema_cbmf(klon) ! cloud base mass flux
	551	SAVE ema_cbmf
	552
	553	REAL ema_pcb(klon) ! cloud base pressure
	554	SAVE ema_pcb
	555
	556	REAL ema_pct(klon) ! cloud top pressure
	557	SAVE ema_pct
	558
	559	REAL bas, top ! cloud base and top levels
	560	SAVE bas
	561	SAVE top
	562
	563	REAL Ma(klon,klev) ! undilute upward mass flux
	564	SAVE Ma
	565	REAL qcondc(klon,klev) ! in-cld water content from convect
	566	SAVE qcondc
	567	REAL ema_work1(klon, klev), ema_work2(klon, klev)
	568	SAVE ema_work1, ema_work2
	569	REAL wdn(klon), tdn(klon), qdn(klon)
	570
	571	REAL wd(klon) ! sb
	572	SAVE wd ! sb
	573
	574	c Variables locales pour la couche limite (al1):
	575	c
	576	cAl1 REAL pblh(klon) ! Hauteur de couche limite
	577	cAl1 SAVE pblh
	578	c34EK
	579	c
	580	c Variables locales:
	581	c
	582	REAL cdragh(klon) ! drag coefficient pour T and Q
	583	REAL cdragm(klon) ! drag coefficient pour vent
	584	cAA
	585	cAA Pour phytrac
	586	cAA
	587	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	588	REAL yu1(klon) ! vents dans la premiere couche U
	589	REAL yv1(klon) ! vents dans la premiere couche V
	590	REAL ffonte(klon,nbsrf) !Flux thermique utilise pour fondre la neige
	591	REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface
	592	c !et necessaire pour limiter la
	593	c !hauteur de neige, en kg/m2/s
[705]	594	REAL fqfonte(klon,nbsrf) !Quantite d'eau de fonte des glaciers
	595	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
[524]	596
	597	LOGICAL offline ! Controle du stockage ds "physique"
	598	PARAMETER (offline=.false.)
	599	INTEGER physid
	600	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
	601	save pfrac_impa
	602	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
	603	save pfrac_nucl
	604	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
	605	save pfrac_1nucl
	606	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	607	REAL frac_nucl(klon,klev) ! idem (nucleation)
	608	#ifdef INCA
	609	REAL :: calday
	610	#endif
	611
	612	cAA
	613	REAL rain_fall(klon) ! pluie
	614	REAL snow_fall(klon) ! neige
	615	save snow_fall, rain_fall
	616	cIM 050204 BEG
	617	REAL total_rain(klon), nday_rain(klon)
	618	save total_rain, nday_rain
	619	cIM 050204 END
	620	REAL evap(klon), devap(klon) ! evaporation et sa derivee
	621	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	622	REAL dlw(klon) ! derivee infra rouge
	623	REAL bils(klon) ! bilan de chaleur au sol
	624	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
	625	C ! type de sous-surface et pondere par la fraction
	626	REAL fder(klon) ! Derive de flux (sensible et latente)
	627	save fder
	628	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	629	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	630	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	631	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	632	c
	633	REAL frugs(klon,nbsrf) ! longueur de rugosite
	634	save frugs
	635	REAL zxrugs(klon) ! longueur de rugosite
	636	c
	637	c Conditions aux limites
	638	c
	639	INTEGER julien
	640	c
	641	INTEGER lmt_pas
	642	SAVE lmt_pas ! frequence de mise a jour
	643	REAL pctsrf(klon,nbsrf)
	644	cIM
	645	REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
	646	REAL paire_ter(klon) !surfaces terre
	647	cIM
	648	SAVE pctsrf ! sous-fraction du sol
	649	REAL albsol(klon)
	650	SAVE albsol ! albedo du sol total
	651	REAL albsollw(klon)
	652	SAVE albsollw ! albedo du sol total
	653
	654	REAL wo(klon,klev)
	655	SAVE wo ! ozone
	656	c======================================================================
	657	c
	658	c Declaration des procedures appelees
	659	c
	660	EXTERNAL angle ! calculer angle zenithal du soleil
	661	EXTERNAL alboc ! calculer l'albedo sur ocean
	662	EXTERNAL ajsec ! ajustement sec
	663	EXTERNAL clmain ! couche limite
	664	EXTERNAL conlmd ! convection (schema LMD)
	665	cKE43
	666	EXTERNAL conema3 ! convect4.3
	667	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	668	cAA
	669	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	670	c ! stockage des coefficients necessaires au
	671	c ! lessivage OFF-LINE et ON-LINE
	672	EXTERNAL hgardfou ! verifier les temperatures
	673	EXTERNAL nuage ! calculer les proprietes radiatives
	674	EXTERNAL o3cm ! initialiser l'ozone
	675	EXTERNAL orbite ! calculer l'orbite terrestre
	676	EXTERNAL ozonecm ! prescrire l'ozone
	677	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	678	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	679	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	680	EXTERNAL suphec ! initialiser certaines constantes
	681	EXTERNAL transp ! transport total de l'eau et de l'energie
	682	EXTERNAL ecribina ! ecrire le fichier binaire global
	683	EXTERNAL ecribins ! ecrire le fichier binaire global
	684	EXTERNAL ecrirega ! ecrire le fichier binaire regional
	685	EXTERNAL ecriregs ! ecrire le fichier binaire regional
	686	cIM
	687	EXTERNAL haut2bas !variables de haut en bas
	688	INTEGER lnblnk1
	689	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
	690	!caracter
	691	c
	692	c Variables locales
	693	c
	694	real clwcon(klon,klev),rnebcon(klon,klev)
	695	real clwcon0(klon,klev),rnebcon0(klon,klev)
	696	save rnebcon, clwcon
	697
	698	REAL rhcl(klon,klev) ! humiditi relative ciel clair
	699	REAL dialiq(klon,klev) ! eau liquide nuageuse
	700	REAL diafra(klon,klev) ! fraction nuageuse
	701	REAL cldliq(klon,klev) ! eau liquide nuageuse
	702	REAL cldfra(klon,klev) ! fraction nuageuse
	703	REAL cldtau(klon,klev) ! epaisseur optique
	704	REAL cldemi(klon,klev) ! emissivite infrarouge
	705	c
	706	CXXX PB
	707	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	708	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
	709	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
	710	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
	711	c
	712	REAL zxfluxt(klon, klev)
	713	REAL zxfluxq(klon, klev)
	714	REAL zxfluxu(klon, klev)
	715	REAL zxfluxv(klon, klev)
	716	CXXX
	717	REAL heat(klon,klev) ! chauffage solaire
	718	REAL heat0(klon,klev) ! chauffage solaire ciel clair
	719	REAL cool(klon,klev) ! refroidissement infrarouge
	720	REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair
	721	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
	722	real sollwdown(klon) ! downward LW flux at surface
	723	cIM BEG
	724	real sollwdownclr(klon) ! downward CS LW flux at surface
	725	real toplwdown(klon) ! downward CS LW flux at TOA
	726	real toplwdownclr(klon) ! downward CS LW flux at TOA
	727	cIM END
	728	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
	729	REAL albpla(klon)
	730	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
	731	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
	732	c Le rayonnement n'est pas calcule tous les pas, il faut donc
	733	c sauvegarder les sorties du rayonnement
	734	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
	735	SAVE sollwdownclr, toplwdown, toplwdownclr
	736	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
	737	c
	738	INTEGER itaprad
	739	SAVE itaprad
	740	c
	741	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	742	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
	743	c
	744	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
	745	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
	746	c
	747	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
	748	c
	749	REAL dist, rmu0(klon), fract(klon)
	750	REAL zdtime, zlongi
	751	c
	752	CHARACTER*2 str2
	753	CHARACTER*2 iqn
	754	c
	755	REAL qcheck
	756	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	757	LOGICAL zx_ajustq
	758	c
	759	REAL za, zb
	760	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
	761	real zqsat(klon,klev)
	762	INTEGER i, k, iq, ig, j, nsrf, ll
	763	REAL t_coup
	764	PARAMETER (t_coup=234.0)
	765	c
	766	REAL zphi(klon,klev)
	767	REAL zx_tmp_x(iim), zx_tmp_yjjmp1
	768	REAL zx_relief(iim,jjmp1)
	769	REAL zx_aire(iim,jjmp1)
	770	cKE43
	771	c Variables locales pour la convection de K. Emanuel (sb):
	772	c
	773	REAL upwd(klon,klev) ! saturated updraft mass flux
	774	REAL dnwd(klon,klev) ! saturated downdraft mass flux
	775	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
	776	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	777	REAL cape(klon) ! CAPE
	778	SAVE cape
	779	CHARACTER*40 capemaxcels !max(CAPE)
	780
	781	REAL pbase(klon) ! cloud base pressure
	782	SAVE pbase
	783	REAL bbase(klon) ! cloud base buoyancy
	784	SAVE bbase
	785	REAL rflag(klon) ! flag fonctionnement de convect
	786	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
	787	c -- convect43:
	788	INTEGER ntra ! nb traceurs pour convect4.3
	789	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
	790	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
	791	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	792	REAL dplcldt(klon), dplcldr(klon)
	793	c? . condm_con(klon,klev),conda_con(klon,klev),
	794	c? . mr_con(klon,klev),ep_con(klon,klev)
	795	c? . ,sadiab(klon,klev),wadiab(klon,klev)
	796	c --
	797	c34EK
	798	c
	799	c Variables du changement
	800	c
	801	c con: convection
	802	c lsc: condensation a grande echelle (Large-Scale-Condensation)
	803	c ajs: ajustement sec
	804	c eva: evaporation de l'eau liquide nuageuse
	805	c vdf: couche limite (Vertical DiFfusion)
	806	REAL d_t_con(klon,klev),d_q_con(klon,klev)
	807	REAL d_u_con(klon,klev),d_v_con(klon,klev)
	808	REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
	809	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
	810	REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
	811	REAL rneb(klon,klev)
	812	c
	813	REAL pmfu(klon,klev), pmfd(klon,klev)
	814	REAL pen_u(klon,klev), pen_d(klon,klev)
	815	REAL pde_u(klon,klev), pde_d(klon,klev)
	816	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	817	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
	818	REAL prfl(klon,klev+1), psfl(klon,klev+1)
	819	c
	820	INTEGER ibas_con(klon), itop_con(klon)
	821	REAL rain_con(klon), rain_lsc(klon)
	822	REAL snow_con(klon), snow_lsc(klon)
	823	REAL d_ts(klon,nbsrf)
	824	c
	825	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	826	REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
	827	c
	828	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	829	REAL d_t_oro(klon,klev)
	830	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	831	REAL d_t_lif(klon,klev)
	832	REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif
	833
	834	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
	835	real ratqsbas,ratqshaut
	836	save ratqsbas,ratqshaut, ratqs
	837	real zpt_conv(klon,klev)
	838
	839	c Parametres lies au nouveau schema de nuages (SB, PDF)
	840	real fact_cldcon
	841	real facttemps
	842	logical ok_newmicro
	843	save ok_newmicro
	844	save fact_cldcon,facttemps
	845	real facteur
	846
	847	integer iflag_cldcon
	848	save iflag_cldcon
	849
	850	logical ptconv(klon,klev)
	851
	852	c
	853	c Variables liees a l'ecriture de la bande histoire physique
	854	c
	855	INTEGER ecrit_mth
	856	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	857	c
	858	INTEGER ecrit_day
	859	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	860	c
	861	INTEGER ecrit_ins
	862	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	863	c
	864	INTEGER ecrit_reg
	865	SAVE ecrit_reg ! frequence d'ecriture
	866	c
	867	integer itau_w ! pas de temps ecriture = itap + itau_phy
	868	c
	869	c
	870	c Variables locales pour effectuer les appels en serie
	871	c
	872	REAL t_seri(klon,klev), q_seri(klon,klev)
	873	REAL ql_seri(klon,klev),qs_seri(klon,klev)
	874	REAL u_seri(klon,klev), v_seri(klon,klev)
	875	c
	876	REAL tr_seri(klon,klev,nbtr)
	877	REAL d_tr(klon,klev,nbtr)
	878
	879	REAL zx_rh(klon,klev)
	880
	881	INTEGER length
	882	PARAMETER ( length = 100 )
	883	REAL tabcntr0( length )
	884	c
	885	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
	886	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	887	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
	888	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
	889	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
	890	c
	891	INTEGER nid_day, nid_mth, nid_ins, nid_nmc
	892	SAVE nid_day, nid_mth, nid_ins, nid_nmc
	893	c
	894	INTEGER nhori, nvert
	895	REAL zsto, zout, zsto1, zsto2
	896	real zjulian
	897	save zjulian
	898
	899	character*20 modname
	900	character*80 abort_message
	901	logical ok_sync
	902	real date0
	903	integer idayref
	904
	905	C essai writephys
	906	integer fid_day, fid_mth, fid_ins
	907	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	908	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	909	parameter (prof2d_on = 1, prof3d_on = 2,
	910	. prof2d_av = 3, prof3d_av = 4)
	911	character*30 nom_fichier
	912	character*10 varname
	913	character*40 vartitle
	914	character*20 varunits
	915	C Variables liees au bilan d'energie et d'enthalpi
	916	REAL ztsol(klon)
	917	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	918	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	919	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	920	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	921	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	922	REAL d_h_vcol_phy
	923	REAL fs_bound, fq_bound
	924	SAVE d_h_vcol_phy
	925	REAL zero_v(klon)
	926	CHARACTER*15 ztit
	927	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	928	SAVE ip_ebil
	929	DATA ip_ebil/0/
	930	INTEGER if_ebil ! level for energy conserv. dignostics
	931	SAVE if_ebil
	932	c+jld ec_conser
	933	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
	934	REAL ZRCPD
	935	c-jld ec_conser
	936	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
	937	REAL t2m(klon,nbsrf), q2m(klon,nbsrf) !temperature, humidite a 2m
	938	REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
	939	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
	940	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
	941	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
	942	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
	943	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
	944	REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
	945	SAVE sulfate_pi
	946
	947	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
	948
	949	REAL re(klon, klev) ! Cloud droplet effective radius
	950	REAL fl(klon, klev) ! denominator of re
	951
	952	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
	953
	954	! Aerosol optical properties
	955	REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
	956	REAL cg_ae(klon,klev,2)
	957
	958	REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
	959	! ok_ade=T -ADE=topswad-topsw
	960
	961	REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
	962	! ok_aie=T ->
	963	! ok_ade=T -AIE=topswai-topswad
	964	! ok_ade=F -AIE=topswai-topsw
	965
	966	REAL aerindex(klon) ! POLDER aerosol index
	967
	968	! Parameters
	969	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
	970	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
	971	cjq-end
	972	c
	973	c Declaration des constantes et des fonctions thermodynamiques
	974	c
[705]	975	SAVE fqfonte
[524]	976	#include "YOMCST.h"
	977	#include "YOETHF.h"
	978	#include "FCTTRE.h"
	979	c======================================================================
	980	modname = 'physiq'
	981	IF (if_ebil.ge.1) THEN
	982	DO i=1,klon
	983	zero_v(i)=0.
	984	END DO
	985	END IF
	986	ok_sync=.TRUE.
	987	IF (nqmax .LT. 2) THEN
	988	abort_message = 'eaux vapeur et liquide sont indispensables'
	989	CALL abort_gcm (modname,abort_message,1)
	990	ENDIF
	991	IF (debut) THEN
	992	CALL suphec ! initialiser constantes et parametres phys.
	993	c
	994	cIM 050204 BEG
	995	DO i=1, klon
	996	nday_rain(i)=0.
	997	ENDDO
	998	cIM 050204 END
	999	c
	1000	c======================================================================
	1001	cIM BEG
	1002	DO k=1, nlev
	1003	DO l=1, nlevSTD
	1004	c
	1005	bb=clevSTD(l)
	1006	c
	1007	IF(l.GE.2) THEN
	1008	aa=clevSTD(l)
	1009	bb=aa(1:lnblnk1(aa))
	1010	ENDIF
	1011	c
	1012	IF(bb.EQ.clev(k)) THEN
	1013	c print*,'k=',k,'l=',l,'clev=',clev(k)
	1014	indENS(k)=l
	1015	c print*,'k=',k,'l=',l,'clev=',clev(k),'indENS=',indENS(k)
	1016	ENDIF
	1017	c
	1018	ENDDO
	1019	ENDDO
	1020	c
	1021	ENDIF !debut
	1022	cIM END
	1023	xjour = rjourvrai
	1024	c
	1025	c Si c'est le debut, il faut initialiser plusieurs choses
	1026	c ********
	1027	c
	1028	IF (debut) THEN
	1029	C
[705]	1030	fqfonte(:,:)=0.
[524]	1031	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	1032	c
	1033	c appel a la lecture du run.def physique
	1034	c
	1035	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
	1036	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
	1037	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
	1038	. ok_ade, ok_aie,
	1039	. bl95_b0, bl95_b1)
	1040	cIM . , RI0)
	1041
	1042	c
	1043	c
	1044	c Initialiser les compteurs:
	1045	c
	1046
	1047	frugs = 0.
	1048	itap = 0
	1049	itaprad = 0
	1050	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
	1051	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
	1052	. falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
	1053	. dlw,radsol,frugs,agesno,clesphy0,
	1054	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
	1055	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
	1056	. run_off_lic_0)
	1057
	1058	c
	1059	radpas = NINT( 86400./dtime/nbapp_rad)
	1060	c
	1061	C on remet le calendrier a zero
	1062	c
	1063	IF (raz_date .eq. 1) THEN
	1064	itau_phy = 0
	1065	ENDIF
	1066
	1067	c
	1068	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
	1069	, ok_instan, ok_region )
	1070	c
	1071	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1072	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
	1073	. pdtphys
	1074	abort_message='Pas physique n est pas correct '
	1075	call abort_gcm(modname,abort_message,1)
	1076	ENDIF
	1077	IF (nlon .NE. klon) THEN
	1078	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	1079	. klon
	1080	abort_message='nlon et klon ne sont pas coherents'
	1081	call abort_gcm(modname,abort_message,1)
	1082	ENDIF
	1083	IF (nlev .NE. klev) THEN
	1084	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	1085	. klev
	1086	abort_message='nlev et klev ne sont pas coherents'
	1087	call abort_gcm(modname,abort_message,1)
	1088	ENDIF
	1089	c
	1090	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
	1091	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1092	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1093	abort_message='Nbre d appels au rayonnement insuffisant'
	1094	call abort_gcm(modname,abort_message,1)
	1095	ENDIF
	1096	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1097	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
	1098	. ok_cvl
	1099	c
	1100	cKE43
	1101	c Initialisation pour la convection de K.E. (sb):
	1102	IF (iflag_con.GE.3) THEN
	1103
	1104	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
	1105	WRITE(lunout,*)
	1106	. "On va utiliser le melange convectif des traceurs qui"
	1107	WRITE(lunout,*)"est calcule dans convect4.3"
	1108	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
	1109
	1110	DO i = 1, klon
	1111	ema_cbmf(i) = 0.
	1112	ema_pcb(i) = 0.
	1113	ema_pct(i) = 0.
	1114	ema_workcbmf(i) = 0.
	1115	ENDDO
	1116	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1117	DO i = 1, klon
	1118	ibas_con(i) = 1
	1119	itop_con(i) = klev+1
	1120	ENDDO
	1121	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	1122
	1123	ENDIF
	1124
	1125	c34EK
	1126	IF (ok_orodr) THEN
	1127	DO i=1,klon
	1128	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1129	ENDDO
	1130	CALL SUGWD(klon,klev,paprs,pplay)
	1131	DO i=1,klon
	1132	zuthe(i)=0.
	1133	zvthe(i)=0.
	1134	if(zstd(i).gt.10.)then
	1135	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1136	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1137	endif
	1138	ENDDO
	1139	ENDIF
	1140	c
	1141	c
	1142	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1143	WRITE(lunout,*)'La frequence de lecture surface est de ',
	1144	. lmt_pas
	1145	c
	1146	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
	1147	IF (ok_mensuel) THEN
	1148	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	1149	. ecrit_mth
	1150	ENDIF
	1151	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
	1152	IF (ok_journe) THEN
	1153	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	1154	. ecrit_day
	1155	ENDIF
	1156	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
	1157	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
	1158	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
	1159	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
	1160	IF (ok_instan) THEN
	1161	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	1162	. ecrit_ins
	1163	ENDIF
	1164	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
	1165	IF (ok_region) THEN
	1166	WRITE(lunout,*)'La frequence de sortie region est de ',
	1167	. ecrit_reg
	1168	ENDIF
	1169
	1170	c
	1171	c Initialiser le couplage si necessaire
	1172	c
	1173	npas = 0
	1174	nexca = 0
	1175	if (ocean == 'couple') then
	1176	npas = itaufin/ iphysiq
	1177	nexca = 86400 / dtime
	1178	write(lunout,*)' ##### Ocean couple #####'
	1179	write(lunout,*)' Valeurs des pas de temps'
	1180	write(lunout,*)' npas = ', npas
	1181	write(lunout,*)' nexca = ', nexca
	1182	endif
	1183	c
	1184	c
	1185	cIM
	1186	capemaxcels = 't_max(X)'
	1187	t2mincels = 't_min(X)'
	1188	t2maxcels = 't_max(X)'
	1189
	1190	c
	1191	c=============================================================
	1192	c Initialisation des sorties
	1193	c=============================================================
	1194
	1195	#ifdef CPP_IOIPSL
	1196
	1197	#ifdef histhf
	1198	#include "ini_histhf.h"
	1199	#endif
	1200
	1201	#ifdef histday
	1202	#include "ini_histday.h"
	1203	#endif
	1204
	1205	#ifdef histmth
	1206	#include "ini_histmth.h"
	1207	#endif
	1208
	1209	#ifdef histins
	1210	#include "ini_histins.h"
	1211	#endif
	1212
	1213	#ifdef histISCCP
	1214	#include "ini_histISCCP.h"
	1215	#endif
	1216
	1217	#ifdef histmthNMC
	1218	#include "ini_histmthNMC.h"
	1219	#endif
	1220
	1221	#ifdef histREGDYN
	1222	#include "ini_histREGDYN.h"
	1223	#endif
	1224
	1225	#ifdef histISCCP
	1226	#include "ini_histISCCP.h"
	1227	#endif
	1228	#endif
	1229
	1230	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1231	date0 = zjulian
	1232	C date0 = day_ini
	1233	WRITE(,) 'physiq date0 : ',date0
	1234	c
	1235	c
	1236	c
	1237	c Prescrire l'ozone dans l'atmosphere
	1238	c
	1239	c
	1240	cc DO i = 1, klon
	1241	cc DO k = 1, klev
	1242	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1243	cc ENDDO
	1244	cc ENDDO
	1245	c
	1246	#ifdef INCA
	1247	iii = MOD(NINT(xjour),360)
	1248	calday = FLOAT(iii) + gmtime
	1249	WRITE(lunout,*) 'initial time ', xjour, calday
	1250	#ifdef INCAINFO
	1251	WRITE(lunout,*) 'Appel CHEMINI ...'
	1252	#endif
	1253	CALL chemini( rpi,
	1254	$ rg,
	1255	$ ra,
	1256	$ airephy,
	1257	$ rlat,
	1258	$ rlon,
	1259	$ presnivs,
	1260	$ calday,
	1261	$ tracnam,
	1262	$ natsnam,
	1263	c $ mxoutflds,
	1264	c $ outinst,
	1265	c $ outtimav,
	1266	$ klon,
	1267	$ nqmax,
	1268	$ pdtphys,
	1269	$ anne_ref,
	1270	$ day_ini)
	1271	#ifdef INCAINFO
	1272	WRITE(lunout,*) 'OK.'
	1273	#endif
	1274	#endif
	1275	c
	1276	ENDIF
	1277	c
	1278	c ************** Fin de IF ( debut ) *************
	1279	c
	1280	c
	1281	c Mettre a zero des variables de sortie (pour securite)
	1282	c
	1283	DO i = 1, klon
	1284	d_ps(i) = 0.0
	1285	ENDDO
	1286	DO k = 1, klev
	1287	DO i = 1, klon
	1288	d_t(i,k) = 0.0
	1289	d_u(i,k) = 0.0
	1290	d_v(i,k) = 0.0
	1291	ENDDO
	1292	ENDDO
	1293	DO iq = 1, nqmax
	1294	DO k = 1, klev
	1295	DO i = 1, klon
	1296	d_qx(i,k,iq) = 0.0
	1297	ENDDO
	1298	ENDDO
	1299	ENDDO
	1300	c
	1301	c Ne pas affecter les valeurs entrees de u, v, h, et q
	1302	c
	1303	DO k = 1, klev
	1304	DO i = 1, klon
	1305	t_seri(i,k) = t(i,k)
	1306	u_seri(i,k) = u(i,k)
	1307	v_seri(i,k) = v(i,k)
	1308	q_seri(i,k) = qx(i,k,ivap)
	1309	ql_seri(i,k) = qx(i,k,iliq)
	1310	qs_seri(i,k) = 0.
	1311	ENDDO
	1312	ENDDO
	1313	IF (nqmax.GE.3) THEN
	1314	DO iq = 3, nqmax
	1315	DO k = 1, klev
	1316	DO i = 1, klon
	1317	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1318	ENDDO
	1319	ENDDO
	1320	ENDDO
	1321	ELSE
	1322	DO k = 1, klev
	1323	DO i = 1, klon
	1324	tr_seri(i,k,1) = 0.0
	1325	ENDDO
	1326	ENDDO
	1327	ENDIF
	1328	C
	1329	DO i = 1, klon
	1330	ztsol(i) = 0.
	1331	ENDDO
	1332	DO nsrf = 1, nbsrf
	1333	DO i = 1, klon
	1334	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1335	ENDDO
	1336	ENDDO
	1337	C
	1338	IF (if_ebil.ge.1) THEN
	1339	ztit='after dynamic'
	1340	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	1341	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1342	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1343	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1344	C on devrait avoir que la variation d'entalpie par la dynamique
	1345	C est egale a la variation de la physique au pas de temps precedent.
	1346	C Donc la somme de ces 2 variations devrait etre nulle.
	1347	call diagphy(airephy,ztit,ip_ebil
	1348	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1349	e , zero_v, zero_v, zero_v, ztsol
	1350	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	1351	s , fs_bound, fq_bound )
	1352	END IF
	1353
	1354	c Diagnostiquer la tendance dynamique
	1355	c
	1356	IF (ancien_ok) THEN
	1357	DO k = 1, klev
	1358	DO i = 1, klon
	1359	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1360	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1361	ENDDO
	1362	ENDDO
	1363	ELSE
	1364	DO k = 1, klev
	1365	DO i = 1, klon
	1366	d_t_dyn(i,k) = 0.0
	1367	d_q_dyn(i,k) = 0.0
	1368	ENDDO
	1369	ENDDO
	1370	ancien_ok = .TRUE.
	1371	ENDIF
	1372	c
	1373	c Ajouter le geopotentiel du sol:
	1374	c
	1375	DO k = 1, klev
	1376	DO i = 1, klon
	1377	zphi(i,k) = pphi(i,k) + pphis(i)
	1378	ENDDO
	1379	ENDDO
	1380	c
	1381	c Verifier les temperatures
	1382	c
	1383	CALL hgardfou(t_seri,ftsol,'debutphy')
	1384	c
	1385	c Incrementer le compteur de la physique
	1386	c
	1387	itap = itap + 1
	1388	julien = MOD(NINT(xjour),360)
	1389	if (julien .eq. 0) julien = 360
	1390	c
	1391	c Mettre en action les conditions aux limites (albedo, sst, etc.).
	1392	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1393	c
	1394	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
	1395	WRITE(lunout,*)' PHYS cond julien ',julien
	1396	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
	1397	ENDIF
	1398	c
	1399	c Re-evaporer l'eau liquide nuageuse
	1400	c
	1401	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1402	DO i = 1, klon
	1403	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1404	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1405	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1406	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1407	zb = MAX(0.0,ql_seri(i,k))
	1408	za = - MAX(0.0,ql_seri(i,k))
	1409	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1410	t_seri(i,k) = t_seri(i,k) + za
	1411	q_seri(i,k) = q_seri(i,k) + zb
	1412	ql_seri(i,k) = 0.0
	1413	d_t_eva(i,k) = za
	1414	d_q_eva(i,k) = zb
	1415	ENDDO
	1416	ENDDO
	1417	c
	1418	IF (if_ebil.ge.2) THEN
	1419	ztit='after reevap'
	1420	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
	1421	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1422	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1423	call diagphy(airephy,ztit,ip_ebil
	1424	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1425	e , zero_v, zero_v, zero_v, ztsol
	1426	e , d_h_vcol, d_qt, d_ec
	1427	s , fs_bound, fq_bound )
	1428	C
	1429	END IF
	1430	C
	1431	c
	1432	c Appeler la diffusion verticale (programme de couche limite)
	1433	c
	1434	DO i = 1, klon
	1435	c if (.not. ok_veget) then
	1436	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
	1437	c endif
	1438	c frugs(i,is_lic) = rugoro(i)
	1439	c frugs(i,is_oce) = rugmer(i)
	1440	c frugs(i,is_sic) = 0.001
	1441	zxrugs(i) = 0.0
	1442	ENDDO
	1443	DO nsrf = 1, nbsrf
	1444	DO i = 1, klon
	1445	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
	1446	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
	1447	ENDDO
	1448	ENDDO
	1449	DO nsrf = 1, nbsrf
	1450	DO i = 1, klon
	1451	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
	1452	ENDDO
	1453	ENDDO
	1454	c
	1455	C calculs necessaires au calcul de l'albedo dans l'interface
	1456	c
	1457	CALL orbite(FLOAT(julien),zlongi,dist)
	1458	IF (cycle_diurne) THEN
	1459	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
	1460	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
	1461	ELSE
	1462	rmu0 = -999.999
	1463	ENDIF
	1464	cIM BEG
	1465	DO i=1, klon
	1466	sunlit(i)=1
	1467	IF(rmu0(i).EQ.0.) sunlit(i)=0
	1468	nbsunlit(1,i)=FLOAT(sunlit(i))
	1469	ENDDO
	1470	cIM END
	1471	C Calcul de l'abedo moyen par maille
	1472	albsol(:)=0.
	1473	albsollw(:)=0.
	1474	DO nsrf = 1, nbsrf
	1475	DO i = 1, klon
	1476	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
	1477	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
	1478	ENDDO
	1479	ENDDO
	1480	C
	1481	C Repartition sous maille des flux LW et SW
	1482	C Modif OM+PASB+JLD
	1483	C Repartition du longwave par sous-surface linearisee
	1484	Cn
	1485
	1486	DO nsrf = 1, nbsrf
	1487	DO i = 1, klon
	1488	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
	1489	c$$$ fsollw(i,nsrf) = sollw(i)
	1490	fsollw(i,nsrf) = sollw(i)
	1491	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
	1492	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
	1493	ENDDO
	1494	ENDDO
	1495
	1496	fder = dlw
	1497
	1498
	1499	CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
	1500	e t_seri,q_seri,u_seri,v_seri,
	1501	e julien, rmu0, co2_ppm,
	1502	e ok_veget, ocean, npas, nexca, ftsol,
	1503	$ soil_model,cdmmax, cdhmax,
	1504	$ ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
	1505	$ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
	1506	$ fluxlat,
	1507	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
	1508	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
	1509	e rlon, rlat, cuphy, cvphy, frugs,
	1510	e debut, lafin, agesno,rugoro ,
	1511	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
	1512	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
	1513	s dsens, devap,
	1514	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
[705]	1515	s fqcalving, fqfonte,ffonte, run_off_lic_0)
[524]	1516	c
	1517	CXXX PB
	1518	CXXX Incrementation des flux
	1519	CXXX
	1520
	1521	zxfluxt=0.
	1522	zxfluxq=0.
	1523	zxfluxu=0.
	1524	zxfluxv=0.
	1525	DO nsrf = 1, nbsrf
	1526	DO k = 1, klev
	1527	DO i = 1, klon
	1528	zxfluxt(i,k) = zxfluxt(i,k) +
	1529	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
	1530	zxfluxq(i,k) = zxfluxq(i,k) +
	1531	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
	1532	zxfluxu(i,k) = zxfluxu(i,k) +
	1533	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
	1534	zxfluxv(i,k) = zxfluxv(i,k) +
	1535	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
	1536	END DO
	1537	END DO
	1538	END DO
	1539	DO i = 1, klon
	1540	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
	1541	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
	1542	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
	1543	fder(i) = dlw(i) + dsens(i) + devap(i)
	1544	ENDDO
	1545
	1546
	1547	DO k = 1, klev
	1548	DO i = 1, klon
	1549	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	1550	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
	1551	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	1552	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	1553	ENDDO
	1554	ENDDO
	1555	c
	1556	IF (if_ebil.ge.2) THEN
	1557	ztit='after clmain'
	1558	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1559	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1560	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1561	call diagphy(airephy,ztit,ip_ebil
	1562	e , zero_v, zero_v, zero_v, zero_v, sens
	1563	e , evap , zero_v, zero_v, ztsol
	1564	e , d_h_vcol, d_qt, d_ec
	1565	s , fs_bound, fq_bound )
	1566	END IF
	1567	C
	1568	c
	1569	c Incrementer la temperature du sol
	1570	c
	1571	DO i = 1, klon
	1572	zxtsol(i) = 0.0
	1573	zxfluxlat(i) = 0.0
	1574	c
	1575	zt2m(i) = 0.0
	1576	zq2m(i) = 0.0
	1577	zu10m(i) = 0.0
	1578	zv10m(i) = 0.0
	1579	cIM cf JLD ??
	1580	zxffonte(i) = 0.0
	1581	zxfqcalving(i) = 0.0
[705]	1582	zxfqfonte(i) = 0.0
[524]	1583	c
	1584	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
	1585	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
	1586	$ THEN
	1587	WRITE(,) 'physiq : pb sous surface au point ', i,
	1588	$ pctsrf(i, 1 : nbsrf)
	1589	ENDIF
	1590	ENDDO
	1591	DO nsrf = 1, nbsrf
	1592	DO i = 1, klon
	1593	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
	1594	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
	1595	cIM cf. JLD
	1596	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
	1597	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
	1598	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1599	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
	1600	cccIM
	1601	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
	1602	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
	1603	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
	1604	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
	1605	cIM cf JLD ??
	1606	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
	1607	zxfqcalving(i) = zxfqcalving(i) +
	1608	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
[705]	1609	zxfqfonte(i) = zxfqfonte(i) +
	1610	. fqfonte(i,nsrf)*pctsrf(i,nsrf)
[524]	1611	c ENDIF
	1612	ENDDO
	1613	ENDDO
	1614
	1615	c
	1616	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
	1617	c
	1618	DO nsrf = 1, nbsrf
	1619	DO i = 1, klon
	1620	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
	1621	cccIM
	1622	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
	1623	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
	1624	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
	1625	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
	1626	cIM cf JLD ??
	1627	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
	1628	IF (pctsrf(i,nsrf) .LT. epsfra)
	1629	. fqcalving(i,nsrf) = zxfqcalving(i)
[705]	1630	IF (pctsrf(i,nsrf) .LT. epsfra)
	1631	. fqfonte(i,nsrf) = zxfqfonte(i)
[524]	1632	ENDDO
	1633	ENDDO
	1634	c
	1635	c
	1636	c Calculer la derive du flux infrarouge
	1637	c
	1638	cXXX DO nsrf = 1, nbsrf
	1639	DO i = 1, klon
	1640	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
	1641	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
	1642	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
	1643	cXXX . *pctsrf(i,nsrf)
	1644	cXXX ENDIF
	1645	cXXX ENDDO
	1646	ENDDO
	1647	c
	1648	c Appeler la convection (au choix)
	1649	c
	1650	DO k = 1, klev
	1651	DO i = 1, klon
	1652	conv_q(i,k) = d_q_dyn(i,k)
	1653	. + d_q_vdf(i,k)/dtime
	1654	conv_t(i,k) = d_t_dyn(i,k)
	1655	. + d_t_vdf(i,k)/dtime
	1656	ENDDO
	1657	ENDDO
	1658	IF (check) THEN
	1659	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	1660	WRITE(lunout,*) "avantcon=", za
	1661	ENDIF
	1662	zx_ajustq = .FALSE.
	1663	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	1664	IF (zx_ajustq) THEN
	1665	DO i = 1, klon
	1666	z_avant(i) = 0.0
	1667	ENDDO
	1668	DO k = 1, klev
	1669	DO i = 1, klon
	1670	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
	1671	. *(paprs(i,k)-paprs(i,k+1))/RG
	1672	ENDDO
	1673	ENDDO
	1674	ENDIF
	1675	IF (iflag_con.EQ.1) THEN
	1676	stop'reactiver le call conlmd dans physiq.F'
	1677	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	1678	c . d_t_con, d_q_con,
	1679	c . rain_con, snow_con, ibas_con, itop_con)
	1680	ELSE IF (iflag_con.EQ.2) THEN
	1681	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
	1682	e conv_t, conv_q, zxfluxq(1,1), omega,
	1683	s d_t_con, d_q_con, rain_con, snow_con,
	1684	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	1685	s kcbot, kctop, kdtop, pmflxr, pmflxs)
	1686	WHERE (rain_con < 0.) rain_con = 0.
	1687	WHERE (snow_con < 0.) snow_con = 0.
	1688	DO i = 1, klon
	1689	ibas_con(i) = klev+1 - kcbot(i)
	1690	itop_con(i) = klev+1 - kctop(i)
	1691	ENDDO
	1692	ELSE IF (iflag_con.GE.3) THEN
	1693	c nb of tracers for the KE convection:
	1694	if (nqmax .GE. 4) then
	1695	ntra = nbtr
	1696	else
	1697	ntra = 1
	1698	endif
	1699	c
	1700	c sb, oct02:
	1701	c Schema de convection modularise et vectorise:
	1702	c (driver commun aux versions 3 et 4)
	1703	c
	1704	IF (ok_cvl) THEN ! new driver for convectL
	1705
	1706	CALL concvl (iflag_con,
	1707	. dtime,paprs,pplay,t_seri,q_seri,
	1708	. u_seri,v_seri,tr_seri,nbtr,
	1709	. ema_work1,ema_work2,
	1710	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	1711	. rain_con, snow_con, ibas_con, itop_con,
	1712	. upwd,dnwd,dnwd0,
	1713	. Ma,cape,tvp,iflagctrl,
	1714	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
	1715	cIM cf. FH
	1716	clwcon0=qcondc
	1717
	1718	ELSE ! ok_cvl
	1719
	1720	CALL conema3 (dtime,
	1721	. paprs,pplay,t_seri,q_seri,
	1722	. u_seri,v_seri,tr_seri,nbtr,
	1723	. ema_work1,ema_work2,
	1724	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	1725	. rain_con, snow_con, ibas_con, itop_con,
	1726	. upwd,dnwd,dnwd0,bas,top,
	1727	. Ma,cape,tvp,rflag,
	1728	. pbase
	1729	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
	1730	. ,clwcon0)
	1731
	1732	ENDIF ! ok_cvl
	1733
	1734	IF (.NOT. ok_gust) THEN
	1735	do i = 1, klon
	1736	wd(i)=0.0
	1737	enddo
	1738	ENDIF
	1739
	1740	c =================================================================== c
	1741	c Calcul des proprietes des nuages convectifs
	1742	c
	1743	DO k = 1, klev
	1744	DO i = 1, klon
	1745	zx_t = t_seri(i,k)
	1746	IF (thermcep) THEN
	1747	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	1748	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	1749	zx_qs = MIN(0.5,zx_qs)
	1750	zcor = 1./(1.-retv*zx_qs)
	1751	zx_qs = zx_qs*zcor
	1752	ELSE
	1753	IF (zx_t.LT.t_coup) THEN
	1754	zx_qs = qsats(zx_t)/pplay(i,k)
	1755	ELSE
	1756	zx_qs = qsatl(zx_t)/pplay(i,k)
	1757	ENDIF
	1758	ENDIF
	1759	zqsat(i,k)=zx_qs
	1760	ENDDO
	1761	ENDDO
	1762
	1763	c calcul des proprietes des nuages convectifs
	1764	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	1765	call clouds_gno
	1766	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	1767
	1768	c =================================================================== c
	1769
	1770	DO i = 1, klon
	1771	ema_pcb(i) = pbase(i)
	1772	ENDDO
	1773	DO i = 1, klon
	1774	ema_pct(i) = paprs(i,itop_con(i))
	1775	ENDDO
	1776	DO i = 1, klon
	1777	ema_cbmf(i) = ema_workcbmf(i)
	1778	ENDDO
	1779	ELSE
	1780	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	1781	CALL abort
	1782	ENDIF
	1783
	1784	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	1785	c . d_u_con, d_v_con)
	1786
	1787	DO k = 1, klev
	1788	DO i = 1, klon
	1789	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
	1790	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
	1791	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
	1792	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
	1793	ENDDO
	1794	ENDDO
	1795	c
	1796	IF (if_ebil.ge.2) THEN
	1797	ztit='after convect'
	1798	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1799	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1800	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1801	call diagphy(airephy,ztit,ip_ebil
	1802	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1803	e , zero_v, rain_con, snow_con, ztsol
	1804	e , d_h_vcol, d_qt, d_ec
	1805	s , fs_bound, fq_bound )
	1806	END IF
	1807	C
	1808	IF (check) THEN
	1809	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	1810	WRITE(lunout,*)"aprescon=", za
	1811	zx_t = 0.0
	1812	za = 0.0
	1813	DO i = 1, klon
	1814	za = za + airephy(i)/FLOAT(klon)
	1815	zx_t = zx_t + (rain_con(i)+
	1816	. snow_con(i))*airephy(i)/FLOAT(klon)
	1817	ENDDO
	1818	zx_t = zx_t/za*dtime
	1819	WRITE(lunout,*)"Precip=", zx_t
	1820	ENDIF
	1821	IF (zx_ajustq) THEN
	1822	DO i = 1, klon
	1823	z_apres(i) = 0.0
	1824	ENDDO
	1825	DO k = 1, klev
	1826	DO i = 1, klon
	1827	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
	1828	. *(paprs(i,k)-paprs(i,k+1))/RG
	1829	ENDDO
	1830	ENDDO
	1831	DO i = 1, klon
	1832	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
	1833	. /z_apres(i)
	1834	ENDDO
	1835	DO k = 1, klev
	1836	DO i = 1, klon
	1837	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
	1838	. z_factor(i).LT.(1.0-1.0E-08)) THEN
	1839	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	1840	ENDIF
	1841	ENDDO
	1842	ENDDO
	1843	ENDIF
	1844	zx_ajustq=.FALSE.
	1845	c
	1846	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
	1847	c
	1848	IF (iflag_con .NE. 2 .AND. debut) THEN
	1849	WRITE(lunout,*)'Pour l instant, seul conflx fonctionne ',
	1850	$ 'avec traceurs', iflag_con
	1851	WRITE(lunout,*)' Mettre iflag_con',
	1852	$ ' = 2 dans run.def et repasser'
	1853	c CALL abort
	1854	ENDIF
	1855	c
	1856	ENDIF !--nqmax.GT.2
	1857	c
	1858	c Appeler l'ajustement sec
	1859	c
	1860	CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
	1861	DO k = 1, klev
	1862	DO i = 1, klon
	1863	t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
	1864	q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
	1865	ENDDO
	1866	ENDDO
	1867	c
	1868	IF (if_ebil.ge.2) THEN
	1869	ztit='after dry_adjust'
	1870	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1871	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1872	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1873	END IF
	1874
	1875
	1876	c-------------------------------------------------------------------------
	1877	c Caclul des ratqs
	1878	c-------------------------------------------------------------------------
	1879
	1880	c print*,'calcul des ratqs'
	1881	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
	1882	c ----------------
	1883	c on ecrase le tableau ratqsc calcule par clouds_gno
	1884	if (iflag_cldcon.eq.1) then
	1885	do k=1,klev
	1886	do i=1,klon
	1887	if(ptconv(i,k)) then
	1888	ratqsc(i,k)=ratqsbas
	1889	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
	1890	else
	1891	ratqsc(i,k)=0.
	1892	endif
	1893	enddo
	1894	enddo
	1895	endif
	1896
	1897	c ratqs stables
	1898	c -------------
	1899	do k=1,klev
	1900	cIM RAJOUT boucle do=i
	1901	do i=1, klon
	1902	cIM ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
	1903	cIM s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
	1904	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	1905	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
	1906	cIM print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
	1907	enddo
	1908	enddo
	1909
	1910
	1911	c ratqs final
	1912	c -----------
	1913	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
	1914	c les ratqs sont une conbinaison de ratqss et ratqsc
	1915	c ratqs final
	1916	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
	1917	c relaxation des ratqs
	1918	c facttemps=exp(-pdtphys/1.e4)
	1919	facteur=exp(-pdtphys*facttemps)
	1920	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
	1921	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
	1922	c print*,'calcul des ratqs fini'
	1923	else
	1924	c on ne prend que le ratqs stable pour fisrtilp
	1925	ratqs(:,:)=ratqss(:,:)
	1926	endif
	1927
	1928
	1929	c
	1930	c Appeler le processus de condensation a grande echelle
	1931	c et le processus de precipitation
	1932	c-------------------------------------------------------------------------
	1933	CALL fisrtilp(dtime,paprs,pplay,
	1934	. t_seri, q_seri,ptconv,ratqs,
	1935	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
	1936	. rain_lsc, snow_lsc,
	1937	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
	1938	. frac_impa, frac_nucl,
	1939	. prfl, psfl, rhcl)
	1940
	1941	WHERE (rain_lsc < 0) rain_lsc = 0.
	1942	WHERE (snow_lsc < 0) snow_lsc = 0.
	1943	DO k = 1, klev
	1944	DO i = 1, klon
	1945	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
	1946	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
	1947	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
	1948	cldfra(i,k) = rneb(i,k)
	1949	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	1950	ENDDO
	1951	ENDDO
	1952	IF (check) THEN
	1953	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	1954	WRITE(lunout,*)"apresilp=", za
	1955	zx_t = 0.0
	1956	za = 0.0
	1957	DO i = 1, klon
	1958	za = za + airephy(i)/FLOAT(klon)
	1959	zx_t = zx_t + (rain_lsc(i)
	1960	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
	1961	ENDDO
	1962	zx_t = zx_t/za*dtime
	1963	WRITE(lunout,*)"Precip=", zx_t
	1964	ENDIF
	1965	c
	1966	IF (if_ebil.ge.2) THEN
	1967	ztit='after fisrt'
	1968	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1969	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1970	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1971	call diagphy(airephy,ztit,ip_ebil
	1972	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1973	e , zero_v, rain_lsc, snow_lsc, ztsol
	1974	e , d_h_vcol, d_qt, d_ec
	1975	s , fs_bound, fq_bound )
	1976	END IF
	1977	c
	1978	c-------------------------------------------------------------------
	1979	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	1980	c-------------------------------------------------------------------
	1981
	1982	c 1. NUAGES CONVECTIFS
	1983	c
	1984	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
	1985
	1986	c Nuages diagnostiques pour Tiedtke
	1987	CALL diagcld1(paprs,pplay,
	1988	. rain_con,snow_con,ibas_con,itop_con,
	1989	. diafra,dialiq)
	1990	DO k = 1, klev
	1991	DO i = 1, klon
	1992	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	1993	cldliq(i,k) = dialiq(i,k)
	1994	cldfra(i,k) = diafra(i,k)
	1995	ENDIF
	1996	ENDDO
	1997	ENDDO
	1998
	1999	ELSE IF (iflag_cldcon.eq.3) THEN
	2000	c On prend pour les nuages convectifs le max du calcul de la
	2001	c convection et du calcul du pas de temps précédent diminué d'un facteur
	2002	c facttemps
	2003	c facttemps=pdtphys/1.e4
	2004	facteur = pdtphys *facttemps
	2005	do k=1,klev
	2006	do i=1,klon
	2007	rnebcon(i,k)=rnebcon(i,k)*facteur
	2008	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
	2009	s then
	2010	rnebcon(i,k)=rnebcon0(i,k)
	2011	clwcon(i,k)=clwcon0(i,k)
	2012	endif
	2013	enddo
	2014	enddo
	2015
	2016	cIM calcul nuages par le simulateur ISCCP
	2017	IF (ok_isccp) THEN
	2018	cIM calcul tau. emi nuages convectifs
	2019	convfra(:,:)=rnebcon(:,:)
	2020	convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
	2021	CALL newmicro (paprs, pplay,ok_newmicro,
	2022	. t_seri, convliq, convfra, dtau_c, dem_c,
	2023	. cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
	2024	. flwp_c, fiwp_c, flwc_c, fiwc_c,
	2025	e ok_aie,
	2026	e sulfate, sulfate_pi,
	2027	e bl95_b0, bl95_b1,
	2028	s cldtaupi, re, fl)
	2029	c
	2030	cIM calcul tau. emi nuages startiformes
	2031	CALL newmicro (paprs, pplay,ok_newmicro,
	2032	. t_seri, cldliq, cldfra, dtau_s, dem_s,
	2033	. cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
	2034	. flwp_s, fiwp_s, flwc_s, fiwc_s,
	2035	e ok_aie,
	2036	e sulfate, sulfate_pi,
	2037	e bl95_b0, bl95_b1,
	2038	s cldtaupi, re, fl)
	2039	c
	2040	cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2041
	2042	cIM inversion des niveaux de pression ==> de haut en bas
	2043	CALL haut2bas(klon, klev, pplay, pfull)
	2044	CALL haut2bas(klon, klev, q_seri, qv)
	2045	CALL haut2bas(klon, klev, cldtot, cc)
	2046	CALL haut2bas(klon, klev, rnebcon, conv)
	2047	CALL haut2bas(klon, klev, dtau_s, dtau_sH2B)
	2048	CALL haut2bas(klon, klev, dtau_c, dtau_cH2B)
	2049	CALL haut2bas(klon, klev, t_seri, at)
	2050	CALL haut2bas(klon, klev, dem_s, dem_sH2B)
	2051	CALL haut2bas(klon, klev, dem_c, dem_cH2B)
	2052	CALL haut2bas(klon, klevp1, paprs, phalf)
	2053
	2054	c open(99,file='tautab.bin',access='sequential',
	2055	c $ form='unformatted',status='old')
	2056	c read(99) tautab
	2057
	2058	cIM210503
	2059	IF (debut) THEN
	2060	open(99,file='tautab.formatted', FORM='FORMATTED')
	2061	read(99,'(f30.20)') tautab
	2062	close(99)
	2063	c
	2064	open(99,file='invtau.formatted',form='FORMATTED')
	2065	read(99,'(i10)') invtau
	2066	close(99)
	2067	c
	2068	cIM: calcul coordonnees regions pour statistiques distribution
	2069	cIM: nuages en ftion du regime dynamique pour regions oceaniques
	2070	IF (ok_regdyn) THEN !histREGDYN
	2071	nsrf=3
	2072	DO nreg=1, nbregdyn
	2073	DO i=1, klon
	2074
	2075	c IF (debut) THEN
	2076	IF(rlon(i).LT.0.) THEN
	2077	rlonPOS(i)=rlon(i)+360.
	2078	ELSE
	2079	rlonPOS(i)=rlon(i)
	2080	ENDIF
	2081	c ENDIF
	2082
	2083	pct_ocean(i,nreg)=0
	2084
	2085	c test si c'est 1 point d'ocean
	2086	IF(pctsrf(i,nsrf).EQ.1.) THEN
	2087
	2088	IF(nreg.EQ.1) THEN
	2089
	2090	c TROP
	2091	IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
	2092	pct_ocean(i,nreg)=1
	2093	ENDIF
	2094
	2095	c PACIFIQUE NORD
	2096	ELSEIF(nreg.EQ.2) THEN
	2097	IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
	2098	IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
	2099	pct_ocean(i,nreg)=1
	2100	ENDIF
	2101	ENDIF
	2102	c CALIFORNIE ST-CU
	2103	ELSEIF(nreg.EQ.3) THEN
	2104	IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
	2105	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
	2106	pct_ocean(i,nreg)=1
	2107	ENDIF
	2108	ENDIF
	2109	c HAWAI
	2110	ELSEIF(nreg.EQ.4) THEN
	2111	IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
	2112	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
	2113	pct_ocean(i,nreg)=1
	2114	ENDIF
	2115	ENDIF
	2116	c WARM POOL
	2117	ELSEIF(nreg.EQ.5) THEN
	2118	IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
	2119	IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
	2120	pct_ocean(i,nreg)=1
	2121	ENDIF
	2122	ENDIF
	2123	ENDIF !nbregdyn
	2124	c TROP
	2125	c IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
	2126	c pct_ocean(i)=.TRUE.
	2127	c WRITE(,) 'pct_ocean =',i, rlon(i), rlat(i)
	2128	c ENDIF !lon
	2129	c ENDIF !lat
	2130
	2131	ENDIF !pctsrf
	2132	ENDDO !klon
	2133	ENDDO !nbregdyn
	2134	ENDIF !ok_regdyn
	2135
	2136	cIM somme de toutes les nhistoW BEG
	2137	DO nreg = 1, nbregdyn
	2138	DO k = 1, kmaxm1
	2139	DO l = 1, lmaxm1
	2140	DO iw = 1, iwmax
	2141	nhistoWt(k,l,iw,nreg)=0.
	2142	ENDDO !iw
	2143	ENDDO !l
	2144	ENDDO !k
	2145	ENDDO !nreg
	2146	cIM somme de toutes les nhistoW END
	2147	ENDIF
	2148	cIM: initialisation de seed
	2149	DO i=1, klon
	2150	seed(i)=i+100
	2151	ENDDO
	2152
	2153	cIM: pas de debug, debugcol
	2154	debug=0
	2155	debugcol=0
	2156	cIM260503
	2157	c o500 ==> distribution nuage ftion du regime dynamique a 500 hPa
	2158	DO k=1, klevm1
	2159	kp1=k+1
	2160	c PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1)
	2161	if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN
	2162	DO i=1, klon
	2163	o500(i)=omega(i,k)*RDAY/100.
	2164	c if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1)
	2165	ENDDO
	2166	GOTO 1000
	2167	endif
	2168	1000 continue
	2169	ENDDO
	2170
	2171	CALL ISCCP_CLOUD_TYPES(
	2172	& debug,
	2173	& debugcol,
	2174	& klon,
	2175	& sunlit,
	2176	& klev,
	2177	& ncol,
	2178	& seed,
	2179	& pfull,
	2180	& phalf,
	2181	& qv, cc, conv, dtau_sH2B, dtau_cH2B,
	2182	& top_height,
	2183	& overlap,
	2184	& tautab,
	2185	& invtau,
	2186	& ztsol,
	2187	& emsfc_lw,
	2188	& at, dem_sH2B, dem_cH2B,
	2189	& fq_isccp,
	2190	& totalcldarea,
	2191	& meanptop,
	2192	& meantaucld,
	2193	& boxtau,
	2194	& boxptop)
	2195
	2196
	2197	c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
	2198	DO l=1, lmaxm1
	2199	DO k=1, kmaxm1
	2200	DO i=1, iim
	2201	fq4d(i,1,k,l)=fq_isccp(1,k,l)
	2202	ENDDO
	2203	DO j=2, jjm
	2204	DO i=1, iim
	2205	ig=i+1+(j-2)*iim
	2206	fq4d(i,j,k,l)=fq_isccp(ig,k,l)
	2207	ENDDO
	2208	ENDDO
	2209	DO i=1, iim
	2210	fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
	2211	ENDDO
	2212	ENDDO
	2213	ENDDO
	2214	c
	2215	DO l=1, lmaxm1
	2216	DO k=1, kmaxm1
	2217	DO j=1, jjmp1
	2218	DO i=1, iim
	2219	ni=(i-1)*lmaxm1+l
	2220	nj=(j-1)*kmaxm1+k
	2221	fq3d(ni,nj)=fq4d(i,j,k,l)
	2222	ENDDO
	2223	ENDDO
	2224	ENDDO
	2225	ENDDO
	2226
	2227	c
	2228	c calculs statistiques distribution nuage ftion du regime dynamique
	2229	c
	2230	c Ce calcul doit etre fait a partir de valeurs mensuelles ??
	2231	CALL histo_o500_pctau(nbregdyn,pct_ocean,o500,fq_isccp,
	2232	&histoW,nhistoW)
	2233	c
	2234	c nhistoWt = somme de toutes les nhistoW
	2235	DO nreg=1, nbregdyn
	2236	DO k = 1, kmaxm1
	2237	DO l = 1, lmaxm1
	2238	DO iw = 1, iwmax
	2239	nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
	2240	& nhistoW(k,l,iw,nreg)
	2241	ENDDO
	2242	ENDDO
	2243	ENDDO
	2244	ENDDO
	2245	c
	2246	ENDIF !ok_isccp
	2247
	2248	c On prend la somme des fractions nuageuses et des contenus en eau
	2249	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2250	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	2251
	2252	ENDIF
	2253
	2254	c
	2255	c 2. NUAGES STARTIFORMES
	2256	c
	2257	IF (ok_stratus) THEN
	2258	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	2259	DO k = 1, klev
	2260	DO i = 1, klon
	2261	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2262	cldliq(i,k) = dialiq(i,k)
	2263	cldfra(i,k) = diafra(i,k)
	2264	ENDIF
	2265	ENDDO
	2266	ENDDO
	2267	ENDIF
	2268	c
	2269	c Precipitation totale
	2270	c
	2271	DO i = 1, klon
	2272	rain_fall(i) = rain_con(i) + rain_lsc(i)
	2273	snow_fall(i) = snow_con(i) + snow_lsc(i)
	2274	ENDDO
	2275	c
	2276	IF (if_ebil.ge.2) THEN
	2277	ztit="after diagcld"
	2278	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	2279	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2280	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2281	END IF
	2282	c
	2283	c Calculer l'humidite relative pour diagnostique
	2284	c
	2285	DO k = 1, klev
	2286	DO i = 1, klon
	2287	zx_t = t_seri(i,k)
	2288	IF (thermcep) THEN
	2289	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2290	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2291	zx_qs = MIN(0.5,zx_qs)
	2292	zcor = 1./(1.-retv*zx_qs)
	2293	zx_qs = zx_qs*zcor
	2294	ELSE
	2295	IF (zx_t.LT.t_coup) THEN
	2296	zx_qs = qsats(zx_t)/pplay(i,k)
	2297	ELSE
	2298	zx_qs = qsatl(zx_t)/pplay(i,k)
	2299	ENDIF
	2300	ENDIF
	2301	zx_rh(i,k) = q_seri(i,k)/zx_qs
	2302	zqsat(i,k)=zx_qs
	2303	ENDDO
	2304	ENDDO
	2305	cjq - introduce the aerosol direct and first indirect radiative forcings
	2306	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	2307	IF (ok_ade.OR.ok_aie) THEN
	2308	! Get sulfate aerosol distribution
	2309	CALL readsulfate(rjourvrai, debut, sulfate)
	2310	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
	2311
	2312	! Calculate aerosol optical properties (Olivier Boucher)
	2313	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
	2314	. tau_ae, piz_ae, cg_ae, aerindex)
	2315	ENDIF
	2316
	2317	#ifdef INCA
	2318	calday = FLOAT(julien) + gmtime
	2319
	2320	#ifdef INCA_AER
	2321	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
	2322	& prfl,psfl,pctsrf,airephy,xjour,rlat,rlon)
	2323	#endif
	2324
	2325	#ifdef INCAINFO
	2326	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
	2327	#endif
	2328	CALL chemhook_begin (calday,
	2329	$ pctsrf(1,3),
	2330	$ rlat,
	2331	$ rlon,
	2332	$ airephy,
	2333	$ paprs,
	2334	$ pplay,
	2335	$ ycoefh,
	2336	$ pphi,
	2337	$ t_seri,
	2338	$ u,
	2339	$ v,
	2340	$ wo,
	2341	$ q_seri,
	2342	$ zxtsol,
	2343	$ zxsnow,
	2344	$ solsw,
	2345	$ albsol,
	2346	$ rain_fall,
	2347	$ snow_fall,
	2348	$ itop_con,
	2349	$ ibas_con,
	2350	$ cldfra,
	2351	$ iim,
	2352	$ jjm,
	2353	$ tr_seri)
	2354	#ifdef INCAINFO
	2355	WRITE(lunout,*)'OK.'
	2356	#endif
	2357	#endif
	2358	c
	2359	c Calculer les parametres optiques des nuages et quelques
	2360	c parametres pour diagnostiques:
	2361	c
	2362	if (ok_newmicro) then
	2363	CALL newmicro (paprs, pplay,ok_newmicro,
	2364	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2365	. cldh, cldl, cldm, cldt, cldq,
	2366	. flwp, fiwp, flwc, fiwc,
	2367	e ok_aie,
	2368	e sulfate, sulfate_pi,
	2369	e bl95_b0, bl95_b1,
	2370	s cldtaupi, re, fl)
	2371	else
	2372	CALL nuage (paprs, pplay,
	2373	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2374	. cldh, cldl, cldm, cldt, cldq,
	2375	e ok_aie,
	2376	e sulfate, sulfate_pi,
	2377	e bl95_b0, bl95_b1,
	2378	s cldtaupi, re, fl)
	2379
	2380	endif
	2381	c
	2382	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	2383	c
	2384	IF (MOD(itaprad,radpas).EQ.0) THEN
	2385	DO i = 1, klon
	2386	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
	2387	. + falbe(i,is_lic) * pctsrf(i,is_lic)
	2388	. + falbe(i,is_ter) * pctsrf(i,is_ter)
	2389	. + falbe(i,is_sic) * pctsrf(i,is_sic)
	2390	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
	2391	. + falblw(i,is_lic) * pctsrf(i,is_lic)
	2392	. + falblw(i,is_ter) * pctsrf(i,is_ter)
	2393	. + falblw(i,is_sic) * pctsrf(i,is_sic)
	2394	ENDDO
	2395	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
	2396	e (dist, rmu0, fract,
	2397	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
	2398	e wo,
	2399	e cldfra, cldemi, cldtau,
	2400	s heat,heat0,cool,cool0,radsol,albpla,
	2401	s topsw,toplw,solsw,sollw,
	2402	s sollwdown,
	2403	s topsw0,toplw0,solsw0,sollw0,
	2404	s lwdn0, lwdn, lwup0, lwup,
	2405	s swdn0, swdn, swup0, swup,
	2406	e ok_ade, ok_aie, ! new for aerosol radiative effects
	2407	e tau_ae, piz_ae, cg_ae, ! ="=
	2408	s topswad, solswad, ! ="=
	2409	e cldtaupi, ! ="=
	2410	s topswai, solswai) ! ="=
	2411	itaprad = 0
	2412	ENDIF
	2413	itaprad = itaprad + 1
	2414	c
	2415	c Ajouter la tendance des rayonnements (tous les pas)
	2416	c
	2417	DO k = 1, klev
	2418	DO i = 1, klon
	2419	t_seri(i,k) = t_seri(i,k)
	2420	. + (heat(i,k)-cool(i,k)) * dtime/86400.
	2421	ENDDO
	2422	ENDDO
	2423	c
	2424	IF (if_ebil.ge.2) THEN
	2425	ztit='after rad'
	2426	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	2427	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2428	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2429	call diagphy(airephy,ztit,ip_ebil
	2430	e , topsw, toplw, solsw, sollw, zero_v
	2431	e , zero_v, zero_v, zero_v, ztsol
	2432	e , d_h_vcol, d_qt, d_ec
	2433	s , fs_bound, fq_bound )
	2434	END IF
	2435	c
	2436	c
	2437	c Calculer l'hydrologie de la surface
	2438	c
	2439	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	2440	c . agesno, ftsol,fqsurf,fsnow, ruis)
	2441	c
	2442	DO i = 1, klon
	2443	zxqsurf(i) = 0.0
	2444	zxsnow(i) = 0.0
	2445	ENDDO
	2446	DO nsrf = 1, nbsrf
	2447	DO i = 1, klon
	2448	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
	2449	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
	2450	ENDDO
	2451	ENDDO
	2452	c
	2453	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
	2454	c
	2455	cXXX DO nsrf = 1, nbsrf
	2456	cXXX DO i = 1, klon
	2457	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
	2458	cXXX fqsurf(i,nsrf) = zxqsurf(i)
	2459	cXXX fsnow(i,nsrf) = zxsnow(i)
	2460	cXXX ENDIF
	2461	cXXX ENDDO
	2462	cXXX ENDDO
	2463	c
	2464	c Calculer le bilan du sol et la derive de temperature (couplage)
	2465	c
	2466	DO i = 1, klon
	2467	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	2468	c a la demande de JLD
	2469	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	2470	ENDDO
	2471	c
	2472	cmoddeblott(jan95)
	2473	c Appeler le programme de parametrisation de l'orographie
	2474	c a l'echelle sous-maille:
	2475	c
	2476	IF (ok_orodr) THEN
	2477	c
	2478	c selection des points pour lesquels le shema est actif:
	2479	igwd=0
	2480	DO i=1,klon
	2481	itest(i)=0
	2482	c IF ((zstd(i).gt.10.0)) THEN
	2483	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	2484	itest(i)=1
	2485	igwd=igwd+1
	2486	idx(igwd)=i
	2487	ENDIF
	2488	ENDDO
	2489	c igwdim=MAX(1,igwd)
	2490	c
	2491	CALL drag_noro(klon,klev,dtime,paprs,pplay,
	2492	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	2493	e igwd,idx,itest,
	2494	e t_seri, u_seri, v_seri,
	2495	s zulow, zvlow, zustr, zvstr,
	2496	s d_t_oro, d_u_oro, d_v_oro)
	2497	c
	2498	c ajout des tendances
	2499	DO k = 1, klev
	2500	DO i = 1, klon
	2501	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
	2502	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
	2503	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
	2504	ENDDO
	2505	ENDDO
	2506	c
	2507	ENDIF ! fin de test sur ok_orodr
	2508	c
	2509	IF (ok_orolf) THEN
	2510	c
	2511	c selection des points pour lesquels le shema est actif:
	2512	igwd=0
	2513	DO i=1,klon
	2514	itest(i)=0
	2515	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	2516	itest(i)=1
	2517	igwd=igwd+1
	2518	idx(igwd)=i
	2519	ENDIF
	2520	ENDDO
	2521	c igwdim=MAX(1,igwd)
	2522	c
	2523	CALL lift_noro(klon,klev,dtime,paprs,pplay,
	2524	e rlat,zmea,zstd,zpic,
	2525	e itest,
	2526	e t_seri, u_seri, v_seri,
	2527	s zulow, zvlow, zustr, zvstr,
	2528	s d_t_lif, d_u_lif, d_v_lif)
	2529	c
	2530	c ajout des tendances
	2531	DO k = 1, klev
	2532	DO i = 1, klon
	2533	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
	2534	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
	2535	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
	2536	ENDDO
	2537	ENDDO
	2538	c
	2539	ENDIF ! fin de test sur ok_orolf
	2540	c
	2541	IF (if_ebil.ge.2) THEN
	2542	ztit='after orography'
	2543	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	2544	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2545	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2546	END IF
	2547	c
	2548	c
	2549	cAA
	2550	cAA Installation de l'interface online-offline pour traceurs
	2551	cAA
	2552	c====================================================================
	2553	c Calcul des tendances traceurs
	2554	c====================================================================
	2555	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
	2556	cKE43 des traceurs => il faut donc mettre des flags a .false.
	2557	IF (iflag_con.GE.3) THEN
	2558	c on ajoute les tendances calculees par KE43
	2559	cXXX OM on onhibe la convection sur les traceurs
	2560	DO iq=1, nqmax-2 ! Sandrine a -3 ???
	2561	cXXX OM on inhibe la convection sur les traceur
	2562	cXXX DO k = 1, nlev
	2563	cXXX DO i = 1, klon
	2564	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
	2565	cXXX ENDDO
	2566	cXXX ENDDO
	2567	WRITE(iqn,'(i2.2)') iq
	2568	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
	2569	ENDDO
	2570	CMAF modif pour garder info du nombre de traceurs auxquels
	2571	C la physique s'applique
	2572	ELSE
	2573	CMAF modif pour garder info du nombre de traceurs auxquels
	2574	C la physique s'applique
	2575	C
	2576	call phytrac (rnpb,
	2577	I itap, julien, gmtime,
	2578	I debut,lafin,
	2579	I nqmax-2,
	2580	I nlon,nlev,dtime,
	2581	I u,v,t,paprs,pplay,
	2582	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	2583	I ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
	2584	I frac_impa, frac_nucl,
	2585	I rlon,presnivs,pphis,pphi,
	2586	I albsol,
	2587	I qx(1,1,1), rhcl,
	2588	I cldfra, rneb, diafra, cldliq, itop_con,
	2589	I ibas_con,
	2590	I pmflxr,pmflxs,prfl,psfl,
	2591	#ifdef INCA_CH4
	2592	I flxmass_w,
	2593	#endif
	2594	O tr_seri)
	2595	ENDIF
	2596
	2597	IF (offline) THEN
	2598
	2599	call phystokenc (
	2600	I nlon,nlev,pdtphys,rlon,rlat,
	2601	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	2602	I ycoefh,yu1,yv1,ftsol,pctsrf,
	2603	I frac_impa, frac_nucl,
	2604	I pphis,airephy,dtime,itap)
	2605
	2606
	2607	ENDIF
	2608
	2609	c
	2610	c Calculer le transport de l'eau et de l'energie (diagnostique)
	2611	c
	2612	CALL transp (paprs,zxtsol,
	2613	e t_seri, q_seri, u_seri, v_seri, zphi,
	2614	s ve, vq, ue, uq)
	2615	c
	2616	c Accumuler les variables a stocker dans les fichiers histoire:
	2617	c
	2618	c
	2619	c
	2620	c+jld ec_conser
	2621	DO k = 1, klev
	2622	DO i = 1, klon
	2623	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
	2624	d_t_ec(i,k)=0.5/ZRCPD
	2625	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	2626	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	2627	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	2628	END DO
	2629	END DO
	2630	c-jld ec_conser
	2631	IF (if_ebil.ge.1) THEN
	2632	ztit='after physic'
	2633	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	2634	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2635	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2636	C Comme les tendances de la physique sont ajoute dans la dynamique,
	2637	C on devrait avoir que la variation d'entalpie par la dynamique
	2638	C est egale a la variation de la physique au pas de temps precedent.
	2639	C Donc la somme de ces 2 variations devrait etre nulle.
	2640	call diagphy(airephy,ztit,ip_ebil
	2641	e , topsw, toplw, solsw, sollw, sens
	2642	e , evap, rain_fall, snow_fall, ztsol
	2643	e , d_h_vcol, d_qt, d_ec
	2644	s , fs_bound, fq_bound )
	2645	C
	2646	d_h_vcol_phy=d_h_vcol
	2647	C
	2648	END IF
	2649	C
	2650	c=======================================================================
	2651	c SORTIES
	2652	c=======================================================================
	2653
	2654	c Interpollation sur quelques niveaux de pression
	2655	c -----------------------------------------------
	2656	c
	2657	c on moyenne mensuellement les champs 3D et on les interpole sur les niveaux STD
	2658	c if(itap.EQ.1.OR.itap.EQ.13.OR.itap.EQ.25.OR.itap.EQ.37) THEN
	2659	c if(MOD(itap,12).EQ.1) THEN
	2660	cIM 120304 END
	2661	DO k=1, nlevSTD
	2662	call plevel(klon,klev,.true.,pplay,rlevSTD(k),
	2663	. t_seri,tlevSTD(:,k))
	2664	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
	2665	. u_seri,ulevSTD(:,k))
	2666	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
	2667	. v_seri,vlevSTD(:,k))
	2668	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
	2669	. zphi,philevSTD(:,k))
	2670	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
	2671	. qx(:,:,ivap),qlevSTD(:,k))
	2672	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
	2673	. zx_rh,rhlevSTD(:,k))
	2674	ENDDO !nlevSTD
	2675	c ENSEMBLES BEG
	2676	DO k=1, nlevENS
	2677	cIM 170304
	2678	tlev(:,k)=tlevSTD(:,indENS(k))
	2679	ulev(:,k)=ulevSTD(:,indENS(k))
	2680	vlev(:,k)=vlevSTD(:,indENS(k))
	2681	philev(:,k)=philevSTD(:,indENS(k))
	2682	qlev(:,k)=qlevSTD(:,indENS(k))
	2683	rhlev(:,k)=rhlevSTD(:,indENS(k))
	2684	c
	2685	call plevel(klon,klevp1,.true.,paprs,rlevENS(k),
	2686	. omega,wlev(:,k))
	2687	c
	2688	ENDDO !k=1, nlevENS
	2689	cIM 100304 BEG
	2690	cIM interpolation a chaque pas de temps du SWup(clr) et SWdn(clr) a 200 hPa
	2691	call plevel(klon,klevp1,.true.,paprs,20000.,
	2692	$ swdn0,SWdn200clr)
	2693	call plevel(klon,klevp1,.false.,paprs,20000.,
	2694	$ swdn,SWdn200)
	2695	call plevel(klon,klevp1,.false.,paprs,20000.,
	2696	$ swup0,SWup200clr)
	2697	call plevel(klon,klevp1,.false.,paprs,20000.,
	2698	$ swup,SWup200)
	2699	c
	2700	call plevel(klon,klevp1,.false.,paprs,20000.,
	2701	$ lwdn0,LWdn200clr)
	2702	call plevel(klon,klevp1,.false.,paprs,20000.,
	2703	$ lwdn,LWdn200)
	2704	call plevel(klon,klevp1,.false.,paprs,20000.,
	2705	$ lwup0,LWup200clr)
	2706	call plevel(klon,klevp1,.false.,paprs,20000.,
	2707	$ lwup,LWup200)
	2708	c
	2709	cIM 100304 END
	2710	c
	2711	c ENSEMBLES END
	2712	c
	2713	c slp sea level pressure
	2714	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	2715	c
	2716	ccc prw = eau precipitable
	2717	DO i = 1, klon
	2718	prw(i) = 0.
	2719	DO k = 1, klev
	2720	prw(i) = prw(i) +
	2721	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	2722	ENDDO
	2723	ENDDO
	2724	c
	2725	cIM sorties bilans energie cinetique et potentielle MJO
	2726	DO k = 1, klev
	2727	DO i = 1, klon
	2728	d_u_oli(i,k) = d_u_oro(i,k) + d_u_lif(i,k)
	2729	d_v_oli(i,k) = d_v_oro(i,k) + d_v_lif(i,k)
	2730	ENDDO
	2731	ENDDO
	2732	c
	2733	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
	2734	cIM PRINT *,' PHYS cond julien ',julien
	2735	c CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
	2736	DO i = 1, klon
	2737	total_rain(i)=rain_fall(i)+snow_fall(i)
	2738	IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
	2739	ENDDO
	2740	c
	2741	ENDIF
	2742	c surface terre
	2743	IF (debut) THEN
	2744	DO i=1, klon
	2745	IF(pctsrf_new(i,is_ter).GT.0.) THEN
	2746	paire_ter(i)=airephy(i)*pctsrf_new(i,is_ter)
	2747	ENDIF
	2748	ENDDO
	2749	ENDIF
	2750	cIM 050204 END
	2751
	2752	c=============================================================
	2753	c
	2754	c Convertir les incrementations en tendances
	2755	c
	2756	DO k = 1, klev
	2757	DO i = 1, klon
	2758	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	2759	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	2760	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	2761	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	2762	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	2763	ENDDO
	2764	ENDDO
	2765	c
	2766	IF (nqmax.GE.3) THEN
	2767	DO iq = 3, nqmax
	2768	DO k = 1, klev
	2769	DO i = 1, klon
	2770	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	2771	ENDDO
	2772	ENDDO
	2773	ENDDO
	2774	ENDIF
	2775	c
	2776	c Sauvegarder les valeurs de t et q a la fin de la physique:
	2777	c
	2778	DO k = 1, klev
	2779	DO i = 1, klon
	2780	t_ancien(i,k) = t_seri(i,k)
	2781	q_ancien(i,k) = q_seri(i,k)
	2782	ENDDO
	2783	ENDDO
	2784	c
	2785	c 22.03.04 BEG
	2786	c=============================================================
	2787	c Ecriture des sorties
	2788	c=============================================================
	2789	#ifdef CPP_IOIPSL
	2790
	2791	#ifdef histhf
	2792	#include "write_histhf.h"
	2793	#endif
	2794
	2795	#ifdef histday
	2796	#include "write_histday.h"
	2797	#endif
	2798
	2799	#ifdef histmth
	2800	#include "write_histmth.h"
	2801	#endif
	2802
	2803	#ifdef histins
	2804	#include "write_histins.h"
	2805	#endif
	2806
	2807	#ifdef histREGDYN
	2808	#include "write_histREGDYN.h"
	2809	#endif
	2810
	2811	#ifdef histISCCP
	2812	#include "write_histISCCP.h"
	2813	#endif
	2814
	2815	#ifdef histmthNMC
	2816	#include "write_histmthNMC.h"
	2817	#endif
	2818
	2819	#endif
	2820
	2821	#ifdef INCA
	2822	#ifdef INCAINFO
	2823	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
	2824	#endif
	2825	CALL chemhook_end (calday,
	2826	$ dtime,
	2827	$ pplay,
	2828	$ t_seri,
	2829	$ tr_seri,
	2830	$ nbtr,
	2831	$ paprs,
	2832	$ q_seri,
	2833	$ anne_ini,
	2834	$ day_ini,
	2835	$ xjour)
	2836	#ifdef INCAINFO
	2837	WRITE(lunout,*)'OK.'
	2838	#endif
	2839	#endif
	2840
	2841	c 22.03.04 END
	2842	c
	2843	c====================================================================
	2844	c Si c'est la fin, il faut conserver l'etat de redemarrage
	2845	c====================================================================
	2846	c
	2847	IF (lafin) THEN
	2848	itau_phy = itau_phy + itap
	2849	ccc IF (ok_oasis) CALL quitcpl
	2850	CALL phyredem ("restartphy.nc",dtime,radpas,
	2851	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
	2852	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
	2853	. solsw, sollwdown,dlw,
	2854	. radsol,frugs,agesno,
	2855	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
	2856	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
	2857	ENDIF
	2858
	2859
	2860	RETURN
	2861	END
	2862	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	2863	IMPLICIT none
	2864	c
	2865	c Calculer et imprimer l'eau totale. A utiliser pour verifier
	2866	c la conservation de l'eau
	2867	c
	2868	#include "YOMCST.h"
	2869	INTEGER klon,klev
	2870	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	2871	REAL aire(klon)
	2872	REAL qtotal, zx, qcheck
	2873	INTEGER i, k
	2874	c
	2875	zx = 0.0
	2876	DO i = 1, klon
	2877	zx = zx + aire(i)
	2878	ENDDO
	2879	qtotal = 0.0
	2880	DO k = 1, klev
	2881	DO i = 1, klon
	2882	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
	2883	. *(paprs(i,k)-paprs(i,k+1))/RG
	2884	ENDDO
	2885	ENDDO
	2886	c
	2887	qcheck = qtotal/zx
	2888	c
	2889	RETURN
	2890	END
	2891	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	2892	IMPLICIT none
	2893	c
	2894	c Tranformer une variable de la grille physique a
	2895	c la grille d'ecriture
	2896	c
	2897	INTEGER nfield,nlon,iim,jjmp1, jjm
	2898	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	2899	c
	2900	INTEGER i, n, ig
	2901	c
	2902	jjm = jjmp1 - 1
	2903	DO n = 1, nfield
	2904	DO i=1,iim
	2905	ecrit(i,n) = fi(1,n)
	2906	ecrit(i+jjm*iim,n) = fi(nlon,n)
	2907	ENDDO
	2908	DO ig = 1, nlon - 2
	2909	ecrit(iim+ig,n) = fi(1+ig,n)
	2910	ENDDO
	2911	ENDDO
	2912	RETURN
	2913	END

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

Download in other formats:

Original Format