Context Navigation

physiq.F @ 786

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

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

Download in other formats:

Original Format