Context Navigation

physiq.F @ 1085

Last change on this file since 1085 was 1068, checked in by Laurent Fairhead, 16 years ago
Correction sans cycle diurne FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.8 KB

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

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

Download in other formats:

Original Format