Context Navigation

physiq.F @ 5456

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

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

Download in other formats:

Original Format