Context Navigation

physiq.F @ 5095

Last change on this file since 5095 was 710, checked in by Laurent Fairhead, 18 years ago
Integration des modifications d'OM et JLD pour une meilleure prise en compte de la conservation du flux d'eau LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 89.0 KB

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

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

Download in other formats:

Original Format