Context Navigation

physiq.F @ 1536

Last change on this file since 1536 was 1536, checked in by musat, 13 years ago
Petites corrections pour le cas ou l'on veut pas sortir de fichier histstn (et qu'on a pas les fichiers initiaux correspondants). IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 127.0 KB

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

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

Download in other formats:

Original Format