Context Navigation

physiq.F @ 1224

Last change on this file since 1224 was 1221, checked in by jghattas, 16 years ago

Modification pour les aerosols sels marin.

Nicolas Yan, Yves Balkanski LSCE

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 114.1 KB

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

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

Download in other formats:

Original Format