Context Navigation

physiq.F @ 505

Last change on this file since 505 was 505, checked in by lmdzadmin, 21 years ago
IM: petit oubli sorties 3d niveaux NMC (fis. "histNMC.nc")
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 83.8 KB

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

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

Download in other formats:

Original Format