Context Navigation

physiq.F @ 1230

Last change on this file since 1230 was 1229, checked in by idelkadi, 16 years ago

Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 114.7 KB

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

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

Download in other formats:

Original Format