Context Navigation

physiq.F @ 1533

Last change on this file since 1533 was 1533, checked in by musat, 13 years ago

Ajout fichier 1D mensuel paramLMDZ_phy.nc contenant

les parametres orbitaux
les taux des GES
les moyennes globales ((bils, evap, evap_land, flat, nettop0, nettop, precip, tsol, t2m, prw)

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

File size: 120.6 KB

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

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

Download in other formats:

Original Format