Context Navigation

physiq.F @ 772

Last change on this file since 772 was 766, checked in by Laurent Fairhead, 18 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 116.9 KB

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

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

Download in other formats:

Original Format