Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

physiq.F @ 1274

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

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

Download in other formats:

Original Format