Context Navigation

physiq.F @ 1251

Last change on this file since 1251 was 556, checked in by lmdzadmin, 20 years ago
Initialisations diverses necessaires au portage Prism AC LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 87.1 KB

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

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

Download in other formats:

Original Format