Context Navigation

physiq.F @ 517

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

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

Download in other formats:

Original Format