Context Navigation

physiq.F @ 523

Last change on this file since 523 was 523, checked in by lmdzadmin, 20 years ago

Dernieres modifs (ha, ha) sur le tag IPSL-CM4_IPCC, MED & JLD:

avant 1930 on force la lecture des sulfates naturels
bug de debordement de tableau dans physiq.F
taucalv passe à 10 ans

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format