Context Navigation

physiq.F @ 5521

Last change on this file since 5521 was 1928, checked in by idelkadi, 11 years ago

Correction du bug dans l'interface avec la physique : On ne veut pas que la distinction entre les nuages convectifs et stratiformes soit prise en compte dans les calculs Cosp. Dans l'interface avec la physique, sont passees en entree pour Cosp, les quantites totales (stratiforme + convective) du contenus en eau et d'autres variables. La fraction nuageuse convective calculee dans la physique est passee en entree pour Cosp dans la version buggee. Elle est remise a 0 dans cette version corrigee.

Mise a jour pour ISCCP : la fraction d'ensoleillement calculee par LMDZ est passee en argument pour Isccp

Rajouts d'autre variables diagnostiques : sunlit (=0 si nuit et =1 si jour) parasol_crefl (reflectance integree)

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

File size: 127.2 KB

Rev	Line
[1279]	1	! $Id: physiq.F 1928 2014-01-16 12:43:34Z fhourdin $
[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
[1717]	1188	integer iostat
[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',
[1717]	1476	$ form='formatted',iostat=iostat)
	1477	if (iostat == 0) then
[1534]	1478	READ(98,*,end=998) nCFMIP
	1479	998 CONTINUE
	1480	CLOSE(98)
[1536]	1481	CONTINUE
[1534]	1482	IF(nCFMIP.GT.npCFMIP) THEN
	1483	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
	1484	CALL abort
	1485	else
	1486	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
	1487	ENDIF
	1488	c
	1489	ALLOCATE(tabCFMIP(nCFMIP))
	1490	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
	1491	ALLOCATE(tabijGCM(nCFMIP))
	1492	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
	1493	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
	1494	c
	1495	c lecture des nCFMIP stations CFMIP, de leur numero
	1496	c et des coordonnees geographiques lonCFMIP, latCFMIP
	1497	c
	1498	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,
	1499	$lonCFMIP, latCFMIP)
	1500	c
	1501	c identification des
	1502	c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
	1503	c 2) indices points tabijGCM de la grille physique 1d sur klon points
	1504	c 3) indices iGCM, jGCM de la grille physique 2d
	1505	c
	1506	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP,
	1507	$tabijGCM, lonGCM, latGCM, iGCM, jGCM)
	1508	c
[1717]	1509	else
	1510	ALLOCATE(tabijGCM(0))
	1511	ALLOCATE(lonGCM(0), latGCM(0))
	1512	ALLOCATE(iGCM(0), jGCM(0))
	1513	end if
[1398]	1514	ENDIF !debut
[1534]	1515
[1279]	1516	DO i=1,klon
	1517	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1518	ENDDO
	1519
[524]	1520	c34EK
	1521	IF (ok_orodr) THEN
[878]	1522
	1523	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1524	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
	1525	! justement quand ok_orodr = false.
	1526	! ce rugoro est utilise par la couche limite et fait double emploi
	1527	! avec les paramétrisations spécifiques de Francois Lott.
	1528	! DO i=1,klon
	1529	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1530	! ENDDO
	1531	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1001]	1532	IF (ok_strato) THEN
	1533	CALL SUGWD_strato(klon,klev,paprs,pplay)
	1534	ELSE
	1535	CALL SUGWD(klon,klev,paprs,pplay)
	1536	ENDIF
	1537
[782]	1538	DO i=1,klon
	1539	zuthe(i)=0.
	1540	zvthe(i)=0.
	1541	if(zstd(i).gt.10.)then
	1542	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1543	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1544	endif
	1545	ENDDO
[524]	1546	ENDIF
	1547	c
	1548	c
	1549	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1550	WRITE(lunout,*)'La frequence de lecture surface est de ',
	1551	. lmt_pas
	1552	c
	1553	capemaxcels = 't_max(X)'
	1554	t2mincels = 't_min(X)'
	1555	t2maxcels = 't_max(X)'
[644]	1556	tinst = 'inst(X)'
	1557	tave = 'ave(X)'
	1558	cIM cf. AM 081204 BEG
	1559	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1560	cIM cf. AM 081204 END
[524]	1561	c
	1562	c=============================================================
	1563	c Initialisation des sorties
	1564	c=============================================================
	1565
	1566	#ifdef CPP_IOIPSL
	1567
[987]	1568	c$OMP MASTER
[1534]	1569	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM,
	1570	& iGCM,jGCM,lonGCM,latGCM,
	1571	& jjmp1,nlevSTD,clevSTD,
	1572	& nbteta, ctetaSTD, dtime,ok_veget,
	1573	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
	1574	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
	1575	& read_climoz, phys_out_filestations,
	1576	& new_aod, aerosol_couple
	1577	& )
[987]	1578	c$OMP END MASTER
	1579	c$OMP BARRIER
[909]	1580
[524]	1581	#ifdef histISCCP
	1582	#include "ini_histISCCP.h"
	1583	#endif
	1584
[1352]	1585	#ifdef histNMC
	1586	#include "ini_histhfNMC.h"
	1587	#include "ini_histdayNMC.h"
[524]	1588	#include "ini_histmthNMC.h"
	1589	#endif
	1590
[687]	1591	#include "ini_histday_seri.h"
[524]	1592
[687]	1593	#include "ini_paramLMDZ_phy.h"
[524]	1594
[644]	1595	#endif
	1596
[1279]	1597	ecrit_hf2mth = ecrit_mth/ecrit_hf
	1598
	1599	ecrit_hf = ecrit_hf * un_jour
[1352]	1600	cIM
[1279]	1601	IF(ecrit_day.LE.1.) THEN
	1602	ecrit_day = ecrit_day * un_jour !en secondes
	1603	ENDIF
[1352]	1604	cIM
[1279]	1605	ecrit_mth = ecrit_mth * un_jour
	1606	ecrit_ins = ecrit_ins * un_jour
	1607	ecrit_reg = ecrit_reg * un_jour
	1608	ecrit_tra = ecrit_tra * un_jour
	1609	ecrit_LES = ecrit_LES * un_jour
	1610	c
	1611	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
	1612	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
	1613	. ecrit_hf2mth
	1614
[524]	1615	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
[1279]	1616	date0 = jD_ref
[524]	1617	WRITE(,) 'physiq date0 : ',date0
	1618	c
	1619	c
	1620	c
	1621	c Prescrire l'ozone dans l'atmosphere
	1622	c
	1623	c
	1624	cc DO i = 1, klon
	1625	cc DO k = 1, klev
	1626	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1627	cc ENDDO
	1628	cc ENDDO
	1629	c
[959]	1630	IF (config_inca /= 'none') THEN
[524]	1631	#ifdef INCA
[959]	1632	CALL VTe(VTphysiq)
	1633	CALL VTb(VTinca)
[1279]	1634	! iii = MOD(NINT(xjour),360)
	1635	! calday = FLOAT(iii) + jH_cur
	1636	calday = FLOAT(days_elapsed) + jH_cur
	1637	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
[959]	1638
	1639	CALL chemini(
[524]	1640	$ rg,
	1641	$ ra,
	1642	$ airephy,
	1643	$ rlat,
	1644	$ rlon,
	1645	$ presnivs,
	1646	$ calday,
	1647	$ klon,
[1146]	1648	$ nqtot,
[524]	1649	$ pdtphys,
[567]	1650	$ annee_ref,
[1279]	1651	$ day_ref,
	1652	$ itau_phy)
[959]	1653
	1654	CALL VTe(VTinca)
	1655	CALL VTb(VTphysiq)
[524]	1656	#endif
[959]	1657	END IF
[524]	1658	c
[998]	1659	c
	1660	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1661	! Nouvelle initialisation pour le rayonnement RRTM
	1662	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1663
	1664	call iniradia(klon,klev,paprs(1,1:klev+1))
	1665
[1279]	1666	C$omp single
	1667	if (read_climoz >= 1) then
	1668	call open_climoz(ncid_climoz, press_climoz)
	1669	END IF
	1670	C$omp end single
[1534]	1671	c
	1672	cIM betaCRF
	1673	pfree=70000. !Pa
	1674	beta_pbl=1.
	1675	beta_free=1.
	1676	lon1_beta=-180.
	1677	lon2_beta=+180.
	1678	lat1_beta=90.
	1679	lat2_beta=-90.
	1680	mskocean_beta=.FALSE.
	1681
	1682	OPEN(99,file='beta_crf.data',status='old',
	1683	$ form='formatted',err=9999)
	1684	READ(99,*,end=9998) pfree
	1685	READ(99,*,end=9998) beta_pbl
	1686	READ(99,*,end=9998) beta_free
	1687	READ(99,*,end=9998) lon1_beta
	1688	READ(99,*,end=9998) lon2_beta
	1689	READ(99,*,end=9998) lat1_beta
	1690	READ(99,*,end=9998) lat2_beta
	1691	READ(99,*,end=9998) mskocean_beta
	1692	9998 Continue
	1693	CLOSE(99)
	1694	9999 Continue
	1695	WRITE(,)'pfree=',pfree
	1696	WRITE(,)'beta_pbl=',beta_pbl
	1697	WRITE(,)'beta_free=',beta_free
	1698	WRITE(,)'lon1_beta=',lon1_beta
	1699	WRITE(,)'lon2_beta=',lon2_beta
	1700	WRITE(,)'lat1_beta=',lat1_beta
	1701	WRITE(,)'lat2_beta=',lat2_beta
	1702	WRITE(,)'mskocean_beta=',mskocean_beta
[524]	1703	ENDIF
[996]	1704	!
	1705	! ************** Fin de IF ( debut ) *************
	1706	!
	1707	!
	1708	! Incrementer le compteur de la physique
	1709	!
	1710	itap = itap + 1
	1711	!
	1712	! Update fraction of the sub-surfaces (pctsrf) and
	1713	! initialize, where a new fraction has appeared, all variables depending
	1714	! on the surface fraction.
	1715	!
[1279]	1716	CALL change_srf_frac(itap, dtime, days_elapsed+1,
[996]	1717	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
	1718
[904]	1719	! Tendances bidons pour les processus qui n'affectent pas certaines
	1720	! variables.
	1721	du0(:,:)=0.
	1722	dv0(:,:)=0.
	1723	dq0(:,:)=0.
	1724	dql0(:,:)=0.
[524]	1725	c
	1726	c Mettre a zero des variables de sortie (pour securite)
	1727	c
	1728	DO i = 1, klon
	1729	d_ps(i) = 0.0
	1730	ENDDO
	1731	DO k = 1, klev
	1732	DO i = 1, klon
	1733	d_t(i,k) = 0.0
	1734	d_u(i,k) = 0.0
	1735	d_v(i,k) = 0.0
	1736	ENDDO
	1737	ENDDO
[1146]	1738	DO iq = 1, nqtot
[524]	1739	DO k = 1, klev
	1740	DO i = 1, klon
	1741	d_qx(i,k,iq) = 0.0
	1742	ENDDO
	1743	ENDDO
	1744	ENDDO
[660]	1745	da(:,:)=0.
	1746	mp(:,:)=0.
	1747	phi(:,:,:)=0.
[524]	1748	c
	1749	c Ne pas affecter les valeurs entrees de u, v, h, et q
	1750	c
	1751	DO k = 1, klev
	1752	DO i = 1, klon
	1753	t_seri(i,k) = t(i,k)
	1754	u_seri(i,k) = u(i,k)
	1755	v_seri(i,k) = v(i,k)
	1756	q_seri(i,k) = qx(i,k,ivap)
	1757	ql_seri(i,k) = qx(i,k,iliq)
	1758	qs_seri(i,k) = 0.
	1759	ENDDO
	1760	ENDDO
[1146]	1761	IF (nqtot.GE.3) THEN
	1762	DO iq = 3, nqtot
[524]	1763	DO k = 1, klev
	1764	DO i = 1, klon
	1765	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1766	ENDDO
	1767	ENDDO
	1768	ENDDO
	1769	ELSE
	1770	DO k = 1, klev
	1771	DO i = 1, klon
	1772	tr_seri(i,k,1) = 0.0
	1773	ENDDO
	1774	ENDDO
	1775	ENDIF
	1776	C
	1777	DO i = 1, klon
	1778	ztsol(i) = 0.
	1779	ENDDO
	1780	DO nsrf = 1, nbsrf
	1781	DO i = 1, klon
	1782	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1783	ENDDO
	1784	ENDDO
[687]	1785	cIM
	1786	IF (ip_ebil_phy.ge.1) THEN
[524]	1787	ztit='after dynamic'
[687]	1788	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	1789	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1790	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1791	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1792	C on devrait avoir que la variation d'entalpie par la dynamique
	1793	C est egale a la variation de la physique au pas de temps precedent.
	1794	C Donc la somme de ces 2 variations devrait etre nulle.
[687]	1795	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1796	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1797	e , zero_v, zero_v, zero_v, ztsol
	1798	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	1799	s , fs_bound, fq_bound )
	1800	END IF
	1801
	1802	c Diagnostiquer la tendance dynamique
	1803	c
	1804	IF (ancien_ok) THEN
	1805	DO k = 1, klev
	1806	DO i = 1, klon
[1054]	1807	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	1808	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
[524]	1809	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1810	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1811	ENDDO
	1812	ENDDO
	1813	ELSE
	1814	DO k = 1, klev
	1815	DO i = 1, klon
[1054]	1816	d_u_dyn(i,k) = 0.0
	1817	d_v_dyn(i,k) = 0.0
[524]	1818	d_t_dyn(i,k) = 0.0
	1819	d_q_dyn(i,k) = 0.0
	1820	ENDDO
	1821	ENDDO
	1822	ancien_ok = .TRUE.
	1823	ENDIF
	1824	c
	1825	c Ajouter le geopotentiel du sol:
	1826	c
	1827	DO k = 1, klev
	1828	DO i = 1, klon
	1829	zphi(i,k) = pphi(i,k) + pphis(i)
	1830	ENDDO
	1831	ENDDO
	1832	c
	1833	c Verifier les temperatures
	1834	c
[687]	1835	cIM BEG
	1836	IF (check) THEN
	1837	amn=MIN(ftsol(1,is_ter),1000.)
	1838	amx=MAX(ftsol(1,is_ter),-1000.)
	1839	DO i=2, klon
	1840	amn=MIN(ftsol(i,is_ter),amn)
	1841	amx=MAX(ftsol(i,is_ter),amx)
	1842	ENDDO
	1843	c
	1844	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	1845	ENDIF !(check) THEN
	1846	cIM END
	1847	c
[524]	1848	CALL hgardfou(t_seri,ftsol,'debutphy')
	1849	c
[687]	1850	cIM BEG
	1851	IF (check) THEN
	1852	amn=MIN(ftsol(1,is_ter),1000.)
	1853	amx=MAX(ftsol(1,is_ter),-1000.)
	1854	DO i=2, klon
	1855	amn=MIN(ftsol(i,is_ter),amn)
	1856	amx=MAX(ftsol(i,is_ter),amx)
	1857	ENDDO
	1858	c
	1859	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	1860	ENDIF !(check) THEN
	1861	cIM END
	1862	c
[524]	1863	c Mettre en action les conditions aux limites (albedo, sst, etc.).
	1864	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1865	c
[1279]	1866	if (read_climoz >= 1) then
	1867	C Ozone from a file
	1868	! Update required ozone index:
	1869	ro3i = int((days_elapsed + jh_cur - jh_1jan)
	1870	$ / ioget_year_len(year_cur) * 360.) + 1
	1871	if (ro3i == 361) ro3i = 360
	1872	C (This should never occur, except perhaps because of roundup
	1873	C error. See documentation.)
	1874	if (ro3i /= co3i) then
	1875	C Update ozone field:
	1876	if (read_climoz == 1) then
	1877	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
	1878	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
	1879	else
	1880	C read_climoz == 2
	1881	call regr_pr_av(ncid_climoz,
	1882	$ (/"tro3 ", "tro3_daylight"/),
	1883	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
	1884	$ v3=wo)
	1885	end if
	1886	! Convert from mole fraction of ozone to column density of ozone in a
	1887	! cell, in kDU:
	1888	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
	1889	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
	1890	C (By regridding ozone values for LMDZ only once every 360th of
	1891	C year, we have already neglected the variation of pressure in one
	1892	C 360th of year. So do not recompute "wo" at each time step even if
	1893	C "zmasse" changes a little.)
	1894	co3i = ro3i
	1895	end if
	1896	elseif (MOD(itap-1,lmt_pas) == 0) THEN
	1897	C Once per day, update ozone from Royer:
	1898	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
[524]	1899	ENDIF
	1900	c
	1901	c Re-evaporer l'eau liquide nuageuse
	1902	c
	1903	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1904	DO i = 1, klon
	1905	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1906	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1907	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1908	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1909	zb = MAX(0.0,ql_seri(i,k))
	1910	za = - MAX(0.0,ql_seri(i,k))
	1911	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1912	t_seri(i,k) = t_seri(i,k) + za
	1913	q_seri(i,k) = q_seri(i,k) + zb
	1914	ql_seri(i,k) = 0.0
	1915	d_t_eva(i,k) = za
	1916	d_q_eva(i,k) = zb
	1917	ENDDO
	1918	ENDDO
[687]	1919	cIM
	1920	IF (ip_ebil_phy.ge.2) THEN
[524]	1921	ztit='after reevap'
[687]	1922	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
[524]	1923	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1924	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	1925	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1926	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1927	e , zero_v, zero_v, zero_v, ztsol
	1928	e , d_h_vcol, d_qt, d_ec
	1929	s , fs_bound, fq_bound )
	1930	C
	1931	END IF
[782]	1932
[524]	1933	c
[883]	1934	c=========================================================================
	1935	! Calculs de l'orbite.
	1936	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	1937	! doit donc etre placé avant radlwsw et pbl_surface
	1938
[1279]	1939	! calcul selon la routine utilisee pour les planetes
	1940	if (new_orbit) then
	1941	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
	1942	day_since_equinox = (jD_cur + jH_cur) - jD_eq
	1943	! day_since_equinox = (jD_cur) - jD_eq
	1944	call solarlong(day_since_equinox, zlongi, dist)
	1945	else
	1946	! calcul selon la routine utilisee pour l'AR4
[883]	1947	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	1948	! solarlong0
[1279]	1949	if (solarlong0<-999.) then
	1950	CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
	1951	else
	1952	zlongi=solarlong0 ! longitude solaire vraie
	1953	dist=1. ! distance au soleil / moyenne
	1954	endif
[883]	1955	endif
[1279]	1956	if(prt_level.ge.1) &
	1957	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
[883]	1958
	1959	! Avec ou sans cycle diurne
[524]	1960	IF (cycle_diurne) THEN
	1961	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
[1279]	1962	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
[524]	1963	ELSE
[1068]	1964	CALL angle(zlongi, rlat, fract, rmu0)
[524]	1965	ENDIF
	1966
[766]	1967	if (mydebug) then
	1968	call writefield_phy('u_seri',u_seri,llm)
	1969	call writefield_phy('v_seri',v_seri,llm)
	1970	call writefield_phy('t_seri',t_seri,llm)
	1971	call writefield_phy('q_seri',q_seri,llm)
	1972	endif
[782]	1973
	1974	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1975	c Appel au pbl_surface : Planetary Boudary Layer et Surface
	1976	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
	1977	c turbulent du couche limit.
	1978	c
	1979	c Certains varibales de sorties de pbl_surface sont utiliser que pour
	1980	c ecriture des fihiers hist_XXXX.nc, ces sont :
	1981	c qsol, zq2m, s_pblh, s_lcl,
	1982	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	1983	c s_therm, s_trmb1, s_trmb2, s_trmb3,
	1984	c zxrugs, zu10m, zv10m, fder,
	1985	c zxqsurf, rh2m, zxfluxu, zxfluxv,
	1986	c frugs, agesno, fsollw, fsolsw,
	1987	c d_ts, fevap, fluxlat, t2m,
	1988	c wfbils, wfbilo, fluxt, fluxu, fluxv,
[687]	1989	c
[782]	1990	c Certains ne sont pas utiliser du tout :
	1991	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
[687]	1992	c
[883]	1993
[782]	1994	CALL pbl_surface(
[1279]	1995	e dtime, date0, itap, days_elapsed+1,
[782]	1996	e debut, lafin,
	1997	e rlon, rlat, rugoro, rmu0,
	1998	e rain_fall, snow_fall, solsw, sollw,
	1999	e t_seri, q_seri, u_seri, v_seri,
	2000	e pplay, paprs, pctsrf,
[888]	2001	+ ftsol, falb1, falb2, u10m, v10m,
[1067]	2002	s sollwdown, cdragh, cdragm, u1, v1,
[888]	2003	s albsol1, albsol2, sens, evap,
[782]	2004	s zxtsol, zxfluxlat, zt2m, qsat2m,
	2005	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
[1534]	2006	s coefh, coefm, slab_wfbils,
[782]	2007	d qsol, zq2m, s_pblh, s_lcl,
	2008	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	2009	d s_therm, s_trmb1, s_trmb2, s_trmb3,
	2010	d zxrugs, zu10m, zv10m, fder,
	2011	d zxqsurf, rh2m, zxfluxu, zxfluxv,
	2012	d frugs, agesno, fsollw, fsolsw,
	2013	d d_ts, fevap, fluxlat, t2m,
	2014	d wfbils, wfbilo, fluxt, fluxu, fluxv,
	2015	- dsens, devap, zxsnow,
[878]	2016	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
[996]	2017
[1067]	2018
[904]	2019	!-----------------------------------------------------------------------------------------
	2020	! ajout des tendances de la diffusion turbulente
	2021	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
	2022	!-----------------------------------------------------------------------------------------
[766]	2023
	2024	if (mydebug) then
	2025	call writefield_phy('u_seri',u_seri,llm)
	2026	call writefield_phy('v_seri',v_seri,llm)
	2027	call writefield_phy('t_seri',t_seri,llm)
	2028	call writefield_phy('q_seri',q_seri,llm)
	2029	endif
	2030
	2031
[687]	2032	IF (ip_ebil_phy.ge.2) THEN
[782]	2033	ztit='after surface_main'
[687]	2034	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2035	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2036	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2037	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2038	e , zero_v, zero_v, zero_v, zero_v, sens
	2039	e , evap , zero_v, zero_v, ztsol
	2040	e , d_h_vcol, d_qt, d_ec
	2041	s , fs_bound, fq_bound )
	2042	END IF
	2043
[881]	2044	c =================================================================== c
	2045	c Calcul de Qsat
	2046
	2047	DO k = 1, klev
	2048	DO i = 1, klon
	2049	zx_t = t_seri(i,k)
	2050	IF (thermcep) THEN
	2051	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2052	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2053	zx_qs = MIN(0.5,zx_qs)
	2054	zcor = 1./(1.-retv*zx_qs)
	2055	zx_qs = zx_qs*zcor
	2056	ELSE
	2057	IF (zx_t.LT.t_coup) THEN
	2058	zx_qs = qsats(zx_t)/pplay(i,k)
	2059	ELSE
	2060	zx_qs = qsatl(zx_t)/pplay(i,k)
	2061	ENDIF
	2062	ENDIF
	2063	zqsat(i,k)=zx_qs
	2064	ENDDO
	2065	ENDDO
	2066
[942]	2067	if (prt_level.ge.1) then
[881]	2068	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	2069	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	2070	endif
[524]	2071	c
	2072	c Appeler la convection (au choix)
	2073	c
	2074	DO k = 1, klev
	2075	DO i = 1, klon
	2076	conv_q(i,k) = d_q_dyn(i,k)
	2077	. + d_q_vdf(i,k)/dtime
	2078	conv_t(i,k) = d_t_dyn(i,k)
	2079	. + d_t_vdf(i,k)/dtime
	2080	ENDDO
	2081	ENDDO
	2082	IF (check) THEN
	2083	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2084	WRITE(lunout,*) "avantcon=", za
	2085	ENDIF
	2086	zx_ajustq = .FALSE.
	2087	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	2088	IF (zx_ajustq) THEN
	2089	DO i = 1, klon
	2090	z_avant(i) = 0.0
	2091	ENDDO
	2092	DO k = 1, klev
	2093	DO i = 1, klon
	2094	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
	2095	. *(paprs(i,k)-paprs(i,k+1))/RG
	2096	ENDDO
	2097	ENDDO
	2098	ENDIF
[959]	2099
	2100	c Calcule de vitesse verticale a partir de flux de masse verticale
	2101	DO k = 1, klev
	2102	DO i = 1, klon
	2103	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
	2104	END DO
	2105	END DO
[1279]	2106	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
	2107	$ omega(igout, :)
[959]	2108
[524]	2109	IF (iflag_con.EQ.1) THEN
	2110	stop'reactiver le call conlmd dans physiq.F'
	2111	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	2112	c . d_t_con, d_q_con,
	2113	c . rain_con, snow_con, ibas_con, itop_con)
	2114	ELSE IF (iflag_con.EQ.2) THEN
	2115	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
[782]	2116	e conv_t, conv_q, -evap, omega,
[524]	2117	s d_t_con, d_q_con, rain_con, snow_con,
	2118	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	2119	s kcbot, kctop, kdtop, pmflxr, pmflxs)
[1015]	2120	d_u_con = 0.
	2121	d_v_con = 0.
	2122
[524]	2123	WHERE (rain_con < 0.) rain_con = 0.
	2124	WHERE (snow_con < 0.) snow_con = 0.
	2125	DO i = 1, klon
	2126	ibas_con(i) = klev+1 - kcbot(i)
	2127	itop_con(i) = klev+1 - kctop(i)
	2128	ENDDO
	2129	ELSE IF (iflag_con.GE.3) THEN
	2130	c nb of tracers for the KE convection:
[619]	2131	c MAF la partie traceurs est faite dans phytrac
	2132	c on met ntra=1 pour limiter les appels mais on peut
	2133	c supprimer les calculs / ftra.
	2134	ntra = 1
[879]	2135
	2136	c=====================================================================================
	2137	cajout pour la parametrisation des poches froides:
	2138	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	2139	do k=1,klev
	2140	do i=1,klon
	2141	if (iflag_wake.eq.1) then
	2142	t_wake(i,k) = t_seri(i,k)
	2143	. +(1-wake_s(i))*wake_deltat(i,k)
	2144	q_wake(i,k) = q_seri(i,k)
	2145	. +(1-wake_s(i))*wake_deltaq(i,k)
	2146	t_undi(i,k) = t_seri(i,k)
	2147	. -wake_s(i)*wake_deltat(i,k)
	2148	q_undi(i,k) = q_seri(i,k)
	2149	. -wake_s(i)*wake_deltaq(i,k)
	2150	else
	2151	t_wake(i,k) = t_seri(i,k)
	2152	q_wake(i,k) = q_seri(i,k)
	2153	t_undi(i,k) = t_seri(i,k)
	2154	q_undi(i,k) = q_seri(i,k)
	2155	endif
	2156	enddo
	2157	enddo
	2158
	2159	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
	2160	cc-- pour le soulevement des particules dans le modele convectif
	2161	c
	2162	do i = 1,klon
	2163	ALE(i) = 0.
	2164	ALP(i) = 0.
	2165	enddo
	2166	c
	2167	ccalcul de ale_wake et alp_wake
	2168	do i = 1,klon
	2169	if (iflag_wake.eq.1) then
	2170	ale_wake(i) = 0.5wake_cstar(i)*2
	2171	alp_wake(i) = wake_fip(i)
	2172	else
	2173	ale_wake(i) = 0.
	2174	alp_wake(i) = 0.
	2175	endif
	2176	enddo
	2177	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
	2178	cdans le thermique sinon
	2179	if (iflag_coupl.eq.0) then
[942]	2180	if (debut) print*,'ALE et ALP imposes'
[879]	2181	do i = 1,klon
	2182	con ne couple que ale
	2183	c ALE(i) = max(ale_wake(i),Ale_bl(i))
	2184	ALE(i) = max(ale_wake(i),ale_bl_prescr)
	2185	con ne couple que alp
	2186	c ALP(i) = alp_wake(i) + Alp_bl(i)
	2187	ALP(i) = alp_wake(i) + alp_bl_prescr
	2188	enddo
	2189	else
[965]	2190	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
[879]	2191	do i = 1,klon
	2192	ALE(i) = max(ale_wake(i),Ale_bl(i))
	2193	ALP(i) = alp_wake(i) + Alp_bl(i)
	2194	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
	2195	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
	2196	enddo
	2197	endif
[979]	2198	do i=1,klon
	2199	if (alp(i)>alp_max) then
[1146]	2200	IF(prt_level>9)WRITE(lunout,*) &
	2201	& 'WARNING SUPER ALP (seuil=',alp_max,
[979]	2202	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
	2203	alp(i)=alp_max
	2204	endif
	2205	if (ale(i)>ale_max) then
[1146]	2206	IF(prt_level>9)WRITE(lunout,*) &
	2207	& 'WARNING SUPER ALE (seuil=',ale_max,
[979]	2208	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
	2209	ale(i)=ale_max
	2210	endif
	2211	enddo
[879]	2212
	2213	cfin calcul ale et alp
	2214	c=================================================================================================
	2215
	2216
[524]	2217	c sb, oct02:
	2218	c Schema de convection modularise et vectorise:
	2219	c (driver commun aux versions 3 et 4)
	2220	c
	2221	IF (ok_cvl) THEN ! new driver for convectL
	2222
[879]	2223	CALL concvl (iflag_con,iflag_clos,
	2224	. dtime,paprs,pplay,t_undi,q_undi,
[1146]	2225	. t_wake,q_wake,wake_s,
[879]	2226	. u_seri,v_seri,tr_seri,nbtr,
	2227	. ALE,ALP,
[524]	2228	. ema_work1,ema_work2,
	2229	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
[879]	2230	. rain_con, snow_con, ibas_con, itop_con, sigd,
[1398]	2231	. ema_cbmf,upwd,dnwd,dnwd0,
[879]	2232	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
[619]	2233	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
[879]	2234	. pmflxr,pmflxs,da,phi,mp,
	2235	. ftd,fqd,lalim_conv,wght_th)
[619]	2236
[973]	2237	cIM begin
[1045]	2238	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2239	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
[973]	2240	cIM end
[524]	2241	cIM cf. FH
	2242	clwcon0=qcondc
[619]	2243	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[524]	2244
[1334]	2245	do i = 1, klon
	2246	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2247	enddo
	2248
[524]	2249	ELSE ! ok_cvl
[1368]	2250
[619]	2251	c MAF conema3 ne contient pas les traceurs
[524]	2252	CALL conema3 (dtime,
	2253	. paprs,pplay,t_seri,q_seri,
[619]	2254	. u_seri,v_seri,tr_seri,ntra,
[524]	2255	. ema_work1,ema_work2,
	2256	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	2257	. rain_con, snow_con, ibas_con, itop_con,
	2258	. upwd,dnwd,dnwd0,bas,top,
	2259	. Ma,cape,tvp,rflag,
	2260	. pbase
	2261	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
	2262	. ,clwcon0)
	2263
	2264	ENDIF ! ok_cvl
	2265
[766]	2266	c
	2267	c Correction precip
	2268	rain_con = rain_con * cvl_corr
	2269	snow_con = snow_con * cvl_corr
	2270	c
	2271
[524]	2272	IF (.NOT. ok_gust) THEN
	2273	do i = 1, klon
	2274	wd(i)=0.0
	2275	enddo
	2276	ENDIF
	2277
	2278	c =================================================================== c
	2279	c Calcul des proprietes des nuages convectifs
	2280	c
	2281
	2282	c calcul des proprietes des nuages convectifs
	2283	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2284	call clouds_gno
	2285	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2286
	2287	c =================================================================== c
	2288
	2289	DO i = 1, klon
[1334]	2290	ema_pcb(i) = paprs(i,ibas_con(i))
[524]	2291	ENDDO
	2292	DO i = 1, klon
[879]	2293	! L'idicage de itop_con peut cacher un pb potentiel
	2294	! FH sous la dictee de JYG, CR
	2295	ema_pct(i) = paprs(i,itop_con(i)+1)
	2296
[878]	2297	if (itop_con(i).gt.klev-3) then
[1279]	2298	if(prt_level >= 9) then
	2299	write(lunout,*)'La convection monte trop haut '
	2300	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
	2301	endif
[878]	2302	endif
[1398]	2303	ENDDO
[881]	2304	ELSE IF (iflag_con.eq.0) THEN
	2305	write(lunout,*) 'On n appelle pas la convection'
	2306	clwcon0=0.
	2307	rnebcon0=0.
	2308	d_t_con=0.
	2309	d_q_con=0.
	2310	d_u_con=0.
	2311	d_v_con=0.
	2312	rain_con=0.
	2313	snow_con=0.
	2314	bas=1
	2315	top=1
[524]	2316	ELSE
	2317	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2318	CALL abort
	2319	ENDIF
	2320
	2321	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2322	c . d_u_con, d_v_con)
	2323
[904]	2324	!-----------------------------------------------------------------------------------------
	2325	! ajout des tendances de la diffusion turbulente
	2326	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
	2327	!-----------------------------------------------------------------------------------------
[766]	2328
	2329	if (mydebug) then
	2330	call writefield_phy('u_seri',u_seri,llm)
	2331	call writefield_phy('v_seri',v_seri,llm)
	2332	call writefield_phy('t_seri',t_seri,llm)
	2333	call writefield_phy('q_seri',q_seri,llm)
	2334	endif
	2335
[687]	2336	cIM
	2337	IF (ip_ebil_phy.ge.2) THEN
[524]	2338	ztit='after convect'
[687]	2339	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2340	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2341	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2342	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2343	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2344	e , zero_v, rain_con, snow_con, ztsol
	2345	e , d_h_vcol, d_qt, d_ec
	2346	s , fs_bound, fq_bound )
	2347	END IF
	2348	C
	2349	IF (check) THEN
	2350	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2351	WRITE(lunout,*)"aprescon=", za
	2352	zx_t = 0.0
	2353	za = 0.0
	2354	DO i = 1, klon
	2355	za = za + airephy(i)/FLOAT(klon)
	2356	zx_t = zx_t + (rain_con(i)+
	2357	. snow_con(i))*airephy(i)/FLOAT(klon)
	2358	ENDDO
	2359	zx_t = zx_t/za*dtime
	2360	WRITE(lunout,*)"Precip=", zx_t
	2361	ENDIF
	2362	IF (zx_ajustq) THEN
	2363	DO i = 1, klon
	2364	z_apres(i) = 0.0
	2365	ENDDO
	2366	DO k = 1, klev
	2367	DO i = 1, klon
	2368	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
	2369	. *(paprs(i,k)-paprs(i,k+1))/RG
	2370	ENDDO
	2371	ENDDO
	2372	DO i = 1, klon
	2373	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
	2374	. /z_apres(i)
	2375	ENDDO
	2376	DO k = 1, klev
	2377	DO i = 1, klon
	2378	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
	2379	. z_factor(i).LT.(1.0-1.0E-08)) THEN
	2380	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2381	ENDIF
	2382	ENDDO
	2383	ENDDO
	2384	ENDIF
	2385	zx_ajustq=.FALSE.
[879]	2386
[524]	2387	c
[879]	2388	c=============================================================================
	2389	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2390	cpour la couche limite diffuse pour l instant
	2391	c
	2392	if (iflag_wake.eq.1) then
	2393	DO k=1,klev
	2394	DO i=1,klon
	2395	dt_dwn(i,k) = ftd(i,k)
[973]	2396	wdt_PBL(i,k) = 0.
[879]	2397	dq_dwn(i,k) = fqd(i,k)
[973]	2398	wdq_PBL(i,k) = 0.
[879]	2399	M_dwn(i,k) = dnwd0(i,k)
	2400	M_up(i,k) = upwd(i,k)
	2401	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
[973]	2402	udt_PBL(i,k) = 0.
[879]	2403	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
[973]	2404	udq_PBL(i,k) = 0.
[879]	2405	ENDDO
	2406	ENDDO
	2407	c
	2408	ccalcul caracteristiques de la poche froide
	2409	call calWAKE (paprs,pplay,dtime
[953]	2410	: ,t_seri,q_seri,omega
[879]	2411	: ,dt_dwn,dq_dwn,M_dwn,M_up
	2412	: ,dt_a,dq_a,sigd
	2413	: ,wdt_PBL,wdq_PBL
	2414	: ,udt_PBL,udq_PBL
	2415	o ,wake_deltat,wake_deltaq,wake_dth
	2416	o ,wake_h,wake_s,wake_dens
	2417	o ,wake_pe,wake_fip,wake_gfl
	2418	o ,dt_wake,dq_wake
	2419	o ,wake_k, t_undi,q_undi
	2420	o ,wake_omgbdth,wake_dp_omgb
	2421	o ,wake_dtKE,wake_dqKE
	2422	o ,wake_dtPBL,wake_dqPBL
	2423	o ,wake_omg,wake_dp_deltomg
	2424	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
	2425	o ,wake_ddeltat,wake_ddeltaq)
	2426	c
[904]	2427	!-----------------------------------------------------------------------------------------
	2428	! ajout des tendances des poches froides
	2429	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2430	! coherent avec les autres d_t_...
	2431	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2432	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2433	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
	2434	!-----------------------------------------------------------------------------------------
[879]	2435
	2436	endif
	2437	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2438	c
[541]	2439	c===================================================================
	2440	c Convection seche (thermiques ou ajustement)
	2441	c===================================================================
[524]	2442	c
[878]	2443	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
	2444	s ,seuil_inversion,weak_inversion,dthmin)
	2445
	2446
	2447
	2448	d_t_ajsb(:,:)=0.
	2449	d_q_ajsb(:,:)=0.
[541]	2450	d_t_ajs(:,:)=0.
	2451	d_u_ajs(:,:)=0.
	2452	d_v_ajs(:,:)=0.
	2453	d_q_ajs(:,:)=0.
[878]	2454	clwcon0th(:,:)=0.
[541]	2455	c
[973]	2456	fm_therm(:,:)=0.
	2457	entr_therm(:,:)=0.
	2458	detr_therm(:,:)=0.
	2459	c
[557]	2460	IF(prt_level>9)WRITE(lunout,*)
	2461	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
[541]	2462	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2463	if(iflag_thermals.lt.0) then
	2464	c Rien
	2465	c ====
[557]	2466	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
[541]	2467
[878]	2468
[541]	2469	else
[878]	2470
[541]	2471	c Thermiques
	2472	c ==========
[557]	2473	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
[541]	2474	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[878]	2475
	2476
	2477	if (iflag_thermals.gt.1) then
[541]	2478	call calltherm(pdtphys
[878]	2479	s ,pplay,paprs,pphi,weak_inversion
	2480	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
[541]	2481	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
[973]	2482	s ,fm_therm,entr_therm,detr_therm
	2483	s ,zqasc,clwcon0th,lmax_th,ratqscth
[879]	2484	s ,ratqsdiff,zqsatth
	2485	con rajoute ale et alp, et les caracteristiques de la couche alim
[1032]	2486	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
[878]	2487	endif
	2488
	2489
	2490	c Ajustement sec
	2491	c ==============
	2492
[1391]	2493	! Dans le cas oÃ¹ on active les thermiques, on fait partir l'ajustement
[878]	2494	! a partir du sommet des thermiques.
	2495	! Dans le cas contraire, on demarre au niveau 1.
	2496
	2497	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
	2498
	2499	if(iflag_thermals.eq.0) then
	2500	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	2501	limbas(:)=1
	2502	else
	2503	limbas(:)=lmax_th(:)
	2504	endif
	2505
	2506	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	2507	! pour des test de convergence numerique.
	2508	! Le nouveau ajsec est a priori mieux, meme pour le cas
	2509	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	2510	! non nulles numeriquement pour des mailles non concernees.
	2511
	2512	if (iflag_thermals.eq.0) then
	2513	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
	2514	s , d_t_ajsb, d_q_ajsb)
	2515	else
	2516	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
	2517	s , d_t_ajsb, d_q_ajsb)
	2518	endif
	2519
[904]	2520	!-----------------------------------------------------------------------------------------
	2521	! ajout des tendances de l'ajustement sec ou des thermiques
	2522	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
[878]	2523	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2524	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2525
[904]	2526	!-----------------------------------------------------------------------------------------
	2527
[878]	2528	endif
	2529
[541]	2530	endif
	2531	c
	2532	c===================================================================
[687]	2533	cIM
	2534	IF (ip_ebil_phy.ge.2) THEN
[524]	2535	ztit='after dry_adjust'
[687]	2536	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2537	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2538	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2539	END IF
	2540
	2541
	2542	c-------------------------------------------------------------------------
	2543	c Caclul des ratqs
	2544	c-------------------------------------------------------------------------
	2545
	2546	c print*,'calcul des ratqs'
	2547	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
	2548	c ----------------
	2549	c on ecrase le tableau ratqsc calcule par clouds_gno
	2550	if (iflag_cldcon.eq.1) then
	2551	do k=1,klev
	2552	do i=1,klon
	2553	if(ptconv(i,k)) then
	2554	ratqsc(i,k)=ratqsbas
	2555	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
	2556	else
	2557	ratqsc(i,k)=0.
	2558	endif
	2559	enddo
	2560	enddo
[1032]	2561
	2562	c-----------------------------------------------------------------------
	2563	c par nversion de la fonction log normale
	2564	c-----------------------------------------------------------------------
[878]	2565	else if (iflag_cldcon.eq.4) then
	2566	ptconvth(:,:)=.false.
	2567	ratqsc(:,:)=0.
[942]	2568	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
[878]	2569	call clouds_gno
	2570	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
[942]	2571	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
[1032]	2572
	2573	c-----------------------------------------------------------------------
	2574	c par calcul direct de l'ecart-type
	2575	c-----------------------------------------------------------------------
	2576
	2577	else if (iflag_cldcon>=5) then
	2578	wmax_th(:)=0.
	2579	zmax_th(:)=0.
	2580	do k=1,klev
	2581	do i=1,klon
	2582	wmax_th(i)=max(wmax_th(i),zw2(i,k))
	2583	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
	2584	enddo
	2585	enddo
	2586	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
	2587	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
	2588
	2589	c On impose que l'air autour de la fraction couverte par le thermique
	2590	c plus son air detraine durant tau_overturning_th soit superieur
	2591	c a 0.1 q_seri
	2592	zz=0.1
	2593	do k=1,klev
	2594	do i=1,klon
	2595	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
	2596	s tau_overturning_th(i)*detr_therm(i,k)
	2597	s *rg/(paprs(i,k)-paprs(i,k+1))
	2598	znum=(1.-zz)*q_seri(i,k)
	2599	zden=zqasc(i,k)-zz*q_seri(i,k)
	2600	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
	2601	lambda_th(i,k)=min(lambda_th(i,k),0.9)
	2602	enddo
	2603	enddo
	2604
	2605	if(iflag_cldcon==5) then
	2606	do k=1,klev
	2607	do i=1,klon
	2608	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
	2609	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
	2610	enddo
	2611	enddo
	2612	else if(iflag_cldcon==6) then
	2613	do k=1,klev
	2614	do i=1,klon
	2615	ratqsc(i,k)=sqrt(lambda_th(i,k))*
	2616	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
	2617	enddo
	2618	enddo
	2619	endif
	2620
[524]	2621	endif
	2622
	2623	c ratqs stables
	2624	c -------------
	2625
[878]	2626	if (iflag_ratqs.eq.0) then
[524]	2627
[878]	2628	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
	2629	do k=1,klev
	2630	do i=1, klon
	2631	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2632	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
	2633	enddo
	2634	enddo
	2635
	2636	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
	2637	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
	2638	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
	2639	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
	2640	! Il s'agit de differents tests dans la phase de reglage du modele
	2641	! avec thermiques.
	2642
	2643	else if (iflag_ratqs.eq.1) then
	2644
	2645	do k=1,klev
	2646	do i=1, klon
	2647	if (pplay(i,k).ge.60000.) then
	2648	ratqss(i,k)=ratqsbas
	2649	else if ((pplay(i,k).ge.30000.).and.
	2650	s (pplay(i,k).lt.60000.)) then
	2651	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2652	s (60000.-pplay(i,k))/(60000.-30000.)
	2653	else
	2654	ratqss(i,k)=ratqshaut
	2655	endif
	2656	enddo
	2657	enddo
	2658
	2659	else
	2660
	2661	do k=1,klev
	2662	do i=1, klon
	2663	if (pplay(i,k).ge.60000.) then
	2664	ratqss(i,k)=ratqsbas
	2665	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
	2666	else if ((pplay(i,k).ge.30000.).and.
	2667	s (pplay(i,k).lt.60000.)) then
	2668	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2669	s (60000.-pplay(i,k))/(60000.-30000.)
	2670	else
	2671	ratqss(i,k)=ratqshaut
	2672	endif
	2673	enddo
	2674	enddo
	2675	endif
	2676
	2677
	2678
	2679
[524]	2680	c ratqs final
	2681	c -----------
[878]	2682
	2683	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
[1032]	2684	s .or.iflag_cldcon.ge.4) then
[878]	2685
	2686	! On ajoute une constante au ratqsc*2 pour tenir compte de
	2687	! fluctuations turbulentes de petite echelle
	2688
	2689	do k=1,klev
	2690	do i=1,klon
	2691	if ((fm_therm(i,k).gt.1.e-10)) then
	2692	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
	2693	endif
	2694	enddo
	2695	enddo
	2696
[1279]	2697	! les ratqs sont une combinaison de ratqss et ratqsc
[1319]	2698	! print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
[878]	2699
[1279]	2700	if (tau_ratqs>1.e-10) then
	2701	facteur=exp(-pdtphys/tau_ratqs)
	2702	else
	2703	facteur=0.
	2704	endif
[878]	2705	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
[1279]	2706	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2707	! FH 22/09/2009
	2708	! La ligne ci-dessous faisait osciller le modele et donnait une solution
	2709	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
	2710	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
	2711	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2712	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
[524]	2713	else
[878]	2714	! on ne prend que le ratqs stable pour fisrtilp
[524]	2715	ratqs(:,:)=ratqss(:,:)
	2716	endif
	2717
	2718
	2719	c
	2720	c Appeler le processus de condensation a grande echelle
	2721	c et le processus de precipitation
	2722	c-------------------------------------------------------------------------
	2723	CALL fisrtilp(dtime,paprs,pplay,
	2724	. t_seri, q_seri,ptconv,ratqs,
	2725	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
	2726	. rain_lsc, snow_lsc,
	2727	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
	2728	. frac_impa, frac_nucl,
	2729	. prfl, psfl, rhcl)
	2730
	2731	WHERE (rain_lsc < 0) rain_lsc = 0.
	2732	WHERE (snow_lsc < 0) snow_lsc = 0.
[904]	2733	!-----------------------------------------------------------------------------------------
	2734	! ajout des tendances de la diffusion turbulente
	2735	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
	2736	!-----------------------------------------------------------------------------------------
[524]	2737	DO k = 1, klev
	2738	DO i = 1, klon
	2739	cldfra(i,k) = rneb(i,k)
	2740	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	2741	ENDDO
	2742	ENDDO
	2743	IF (check) THEN
	2744	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2745	WRITE(lunout,*)"apresilp=", za
	2746	zx_t = 0.0
	2747	za = 0.0
	2748	DO i = 1, klon
	2749	za = za + airephy(i)/FLOAT(klon)
	2750	zx_t = zx_t + (rain_lsc(i)
	2751	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
	2752	ENDDO
	2753	zx_t = zx_t/za*dtime
	2754	WRITE(lunout,*)"Precip=", zx_t
	2755	ENDIF
[687]	2756	cIM
	2757	IF (ip_ebil_phy.ge.2) THEN
[524]	2758	ztit='after fisrt'
[687]	2759	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2760	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2761	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2762	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2763	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2764	e , zero_v, rain_lsc, snow_lsc, ztsol
	2765	e , d_h_vcol, d_qt, d_ec
	2766	s , fs_bound, fq_bound )
	2767	END IF
[766]	2768
	2769	if (mydebug) then
	2770	call writefield_phy('u_seri',u_seri,llm)
	2771	call writefield_phy('v_seri',v_seri,llm)
	2772	call writefield_phy('t_seri',t_seri,llm)
	2773	call writefield_phy('q_seri',q_seri,llm)
	2774	endif
	2775
[524]	2776	c
	2777	c-------------------------------------------------------------------
	2778	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	2779	c-------------------------------------------------------------------
	2780
	2781	c 1. NUAGES CONVECTIFS
	2782	c
[644]	2783	cIM cf FH
	2784	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
[878]	2785	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
[644]	2786	snow_tiedtke=0.
	2787	c print*,'avant calcul de la pseudo precip '
	2788	c print*,'iflag_cldcon',iflag_cldcon
	2789	if (iflag_cldcon.eq.-1) then
	2790	rain_tiedtke=rain_con
	2791	else
	2792	c print*,'calcul de la pseudo precip '
	2793	rain_tiedtke=0.
	2794	c print*,'calcul de la pseudo precip 0'
	2795	do k=1,klev
	2796	do i=1,klon
	2797	if (d_q_con(i,k).lt.0.) then
	2798	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
	2799	s *(paprs(i,k)-paprs(i,k+1))/rg
	2800	endif
	2801	enddo
	2802	enddo
	2803	endif
	2804	c
	2805	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	2806	c
[524]	2807
	2808	c Nuages diagnostiques pour Tiedtke
	2809	CALL diagcld1(paprs,pplay,
[644]	2810	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
	2811	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
[524]	2812	. diafra,dialiq)
	2813	DO k = 1, klev
	2814	DO i = 1, klon
	2815	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2816	cldliq(i,k) = dialiq(i,k)
	2817	cldfra(i,k) = diafra(i,k)
	2818	ENDIF
	2819	ENDDO
	2820	ENDDO
	2821
[878]	2822	ELSE IF (iflag_cldcon.ge.3) THEN
[524]	2823	c On prend pour les nuages convectifs le max du calcul de la
[766]	2824	c convection et du calcul du pas de temps precedent diminue d'un facteur
[524]	2825	c facttemps
	2826	facteur = pdtphys *facttemps
	2827	do k=1,klev
	2828	do i=1,klon
	2829	rnebcon(i,k)=rnebcon(i,k)*facteur
	2830	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
	2831	s then
	2832	rnebcon(i,k)=rnebcon0(i,k)
	2833	clwcon(i,k)=clwcon0(i,k)
	2834	endif
	2835	enddo
	2836	enddo
	2837
[644]	2838	c
[766]	2839	cjq - introduce the aerosol direct and first indirect radiative forcings
	2840	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	2841	IF (ok_ade.OR.ok_aie) THEN
[1279]	2842	IF (.NOT. aerosol_couple)
	2843	& CALL readaerosol_optic(
	2844	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
	2845	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
	2846	& mass_solu_aero, mass_solu_aero_pi,
	2847	& tau_aero, piz_aero, cg_aero,
	2848	& tausum_aero, tau3d_aero)
[766]	2849	ELSE
[1279]	2850	tau_aero(:,:,:,:) = 0.
	2851	piz_aero(:,:,:,:) = 0.
	2852	cg_aero(:,:,:,:) = 0.
[766]	2853	ENDIF
	2854
[524]	2855	cIM calcul nuages par le simulateur ISCCP
[644]	2856	c
[839]	2857	#ifdef histISCCP
[524]	2858	IF (ok_isccp) THEN
[1035]	2859	c
[1045]	2860	cIM lecture invtau, tautab des fichiers formattes
[1035]	2861	c
[1045]	2862	IF (debut) THEN
	2863	c$OMP MASTER
	2864	c
	2865	open(99,file='tautab.formatted', FORM='FORMATTED')
	2866	read(99,'(f30.20)') tautab_omp
	2867	close(99)
	2868	c
	2869	open(99,file='invtau.formatted',form='FORMATTED')
	2870	read(99,'(i10)') invtau_omp
	2871
	2872	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
	2873	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
	2874
	2875	close(99)
	2876	c$OMP END MASTER
	2877	c$OMP BARRIER
	2878	tautab=tautab_omp
	2879	invtau=invtau_omp
	2880	c
	2881	ENDIF !debut
	2882	c
[828]	2883	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	2884	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
[644]	2885	#include "calcul_simulISCCP.h"
[828]	2886	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
[524]	2887	ENDIF !ok_isccp
[839]	2888	#endif
[524]	2889
	2890	c On prend la somme des fractions nuageuses et des contenus en eau
	2891	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2892	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	2893
	2894	ENDIF
	2895	c
	2896	c 2. NUAGES STARTIFORMES
	2897	c
	2898	IF (ok_stratus) THEN
	2899	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	2900	DO k = 1, klev
	2901	DO i = 1, klon
	2902	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2903	cldliq(i,k) = dialiq(i,k)
	2904	cldfra(i,k) = diafra(i,k)
	2905	ENDIF
	2906	ENDDO
	2907	ENDDO
	2908	ENDIF
	2909	c
	2910	c Precipitation totale
	2911	c
	2912	DO i = 1, klon
	2913	rain_fall(i) = rain_con(i) + rain_lsc(i)
	2914	snow_fall(i) = snow_con(i) + snow_lsc(i)
	2915	ENDDO
[687]	2916	cIM
	2917	IF (ip_ebil_phy.ge.2) THEN
[524]	2918	ztit="after diagcld"
[687]	2919	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2920	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2921	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2922	END IF
	2923	c
	2924	c Calculer l'humidite relative pour diagnostique
	2925	c
	2926	DO k = 1, klev
	2927	DO i = 1, klon
	2928	zx_t = t_seri(i,k)
	2929	IF (thermcep) THEN
	2930	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2931	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2932	zx_qs = MIN(0.5,zx_qs)
	2933	zcor = 1./(1.-retv*zx_qs)
	2934	zx_qs = zx_qs*zcor
	2935	ELSE
	2936	IF (zx_t.LT.t_coup) THEN
	2937	zx_qs = qsats(zx_t)/pplay(i,k)
	2938	ELSE
	2939	zx_qs = qsatl(zx_t)/pplay(i,k)
	2940	ENDIF
	2941	ENDIF
	2942	zx_rh(i,k) = q_seri(i,k)/zx_qs
	2943	zqsat(i,k)=zx_qs
	2944	ENDDO
	2945	ENDDO
[782]	2946
[687]	2947	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	2948	c equivalente a 2m (tpote) pour diagnostique
	2949	c
	2950	DO i = 1, klon
	2951	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	2952	IF (thermcep) THEN
	2953	IF(zt2m(i).LT.RTT) then
	2954	Lheat=RLSTT
	2955	ELSE
	2956	Lheat=RLVTT
	2957	ENDIF
	2958	ELSE
	2959	IF (zt2m(i).LT.RTT) THEN
	2960	Lheat=RLSTT
	2961	ELSE
	2962	Lheat=RLVTT
	2963	ENDIF
	2964	ENDIF
	2965	tpote(i) = tpot(i)*
	2966	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	2967	ENDDO
[524]	2968
[959]	2969	IF (config_inca /= 'none') THEN
[524]	2970	#ifdef INCA
[959]	2971	CALL VTe(VTphysiq)
	2972	CALL VTb(VTinca)
[1279]	2973	calday = FLOAT(days_elapsed + 1) + jH_cur
[524]	2974
[1287]	2975	call chemtime(itap+itau_phy-1, date0, dtime)
[959]	2976	IF (config_inca == 'aero') THEN
[1015]	2977	CALL AEROSOL_METEO_CALC(
	2978	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
[1279]	2979	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
[959]	2980	END IF
[524]	2981
[959]	2982	zxsnow_dummy(:) = 0.0
[625]	2983
[959]	2984	CALL chemhook_begin (calday,
[1279]	2985	$ days_elapsed+1,
	2986	$ jH_cur,
[593]	2987	$ pctsrf(1,1),
[524]	2988	$ rlat,
	2989	$ rlon,
	2990	$ airephy,
	2991	$ paprs,
	2992	$ pplay,
[1067]	2993	$ coefh,
[524]	2994	$ pphi,
	2995	$ t_seri,
	2996	$ u,
	2997	$ v,
[1279]	2998	$ wo(:, :, 1),
[524]	2999	$ q_seri,
	3000	$ zxtsol,
[782]	3001	$ zxsnow_dummy,
[524]	3002	$ solsw,
[888]	3003	$ albsol1,
[524]	3004	$ rain_fall,
	3005	$ snow_fall,
	3006	$ itop_con,
	3007	$ ibas_con,
	3008	$ cldfra,
	3009	$ iim,
	3010	$ jjm,
[616]	3011	$ tr_seri,
	3012	$ ftsol,
	3013	$ paprs,
	3014	$ cdragh,
	3015	$ cdragm,
	3016	$ pctsrf,
	3017	$ pdtphys,
	3018	$ itap)
	3019
[959]	3020	CALL VTe(VTinca)
	3021	CALL VTb(VTphysiq)
	3022	#endif
	3023	END IF !config_inca /= 'none'
[524]	3024	c
	3025	c Calculer les parametres optiques des nuages et quelques
	3026	c parametres pour diagnostiques:
	3027	c
[959]	3028
	3029	IF (aerosol_couple) THEN
[1279]	3030	mass_solu_aero(:,:) = ccm(:,:,1)
	3031	mass_solu_aero_pi(:,:) = ccm(:,:,2)
	3032	END IF
[1534]	3033	c
[524]	3034	if (ok_newmicro) then
	3035	CALL newmicro (paprs, pplay,ok_newmicro,
	3036	. t_seri, cldliq, cldfra, cldtau, cldemi,
	3037	. cldh, cldl, cldm, cldt, cldq,
	3038	. flwp, fiwp, flwc, fiwc,
	3039	e ok_aie,
[1279]	3040	e mass_solu_aero, mass_solu_aero_pi,
[524]	3041	e bl95_b0, bl95_b1,
[1279]	3042	s cldtaupi, re, fl, ref_liq, ref_ice)
[524]	3043	else
	3044	CALL nuage (paprs, pplay,
	3045	. t_seri, cldliq, cldfra, cldtau, cldemi,
	3046	. cldh, cldl, cldm, cldt, cldq,
	3047	e ok_aie,
[1279]	3048	e mass_solu_aero, mass_solu_aero_pi,
[524]	3049	e bl95_b0, bl95_b1,
	3050	s cldtaupi, re, fl)
	3051
	3052	endif
	3053	c
[1534]	3054	cIM betaCRF
	3055	c
	3056	cldtaurad = cldtau
	3057	cldemirad = cldemi
	3058	c
	3059	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
	3060	$lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
	3061	c
	3062	c global
	3063	c
	3064	DO k=1, klev
	3065	DO i=1, klon
	3066	if (pplay(i,k).GE.pfree) THEN
	3067	beta(i,k) = beta_pbl
	3068	else
	3069	beta(i,k) = beta_free
	3070	endif
	3071	if (mskocean_beta) THEN
	3072	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3073	endif
	3074	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3075	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3076	ENDDO
	3077	ENDDO
	3078	c
	3079	else
	3080	c
	3081	c regional
	3082	c
	3083	DO k=1, klev
	3084	DO i=1,klon
	3085	c
	3086	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.
	3087	$ rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
	3088	if (pplay(i,k).GE.pfree) THEN
	3089	beta(i,k) = beta_pbl
	3090	else
	3091	beta(i,k) = beta_free
	3092	endif
	3093	if (mskocean_beta) THEN
	3094	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3095	endif
	3096	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3097	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3098	endif
	3099	c
	3100	ENDDO
	3101	ENDDO
	3102	c
	3103	endif
	3104	c
[524]	3105	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	3106	c
	3107	IF (MOD(itaprad,radpas).EQ.0) THEN
[782]	3108
[524]	3109	DO i = 1, klon
[888]	3110	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
	3111	. + falb1(i,is_lic) * pctsrf(i,is_lic)
	3112	. + falb1(i,is_ter) * pctsrf(i,is_ter)
	3113	. + falb1(i,is_sic) * pctsrf(i,is_sic)
	3114	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
	3115	. + falb2(i,is_lic) * pctsrf(i,is_lic)
	3116	. + falb2(i,is_ter) * pctsrf(i,is_ter)
	3117	. + falb2(i,is_sic) * pctsrf(i,is_sic)
[524]	3118	ENDDO
[766]	3119
	3120	if (mydebug) then
	3121	call writefield_phy('u_seri',u_seri,llm)
	3122	call writefield_phy('v_seri',v_seri,llm)
	3123	call writefield_phy('t_seri',t_seri,llm)
	3124	call writefield_phy('q_seri',q_seri,llm)
	3125	endif
	3126
[955]	3127	IF (aerosol_couple) THEN
[959]	3128	#ifdef INCA
[1279]	3129	CALL radlwsw_inca
	3130	e (kdlon,kflev,dist, rmu0, fract, solaire,
	3131	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
	3132	e wo(:, :, 1),
[1534]	3133	e cldfra, cldemirad, cldtaurad,
[1279]	3134	s heat,heat0,cool,cool0,radsol,albpla,
	3135	s topsw,toplw,solsw,sollw,
	3136	s sollwdown,
	3137	s topsw0,toplw0,solsw0,sollw0,
	3138	s lwdn0, lwdn, lwup0, lwup,
	3139	s swdn0, swdn, swup0, swup,
	3140	e ok_ade, ok_aie,
	3141	e tau_aero, piz_aero, cg_aero,
	3142	s topswad_aero, solswad_aero,
	3143	s topswad0_aero, solswad0_aero,
	3144	s topsw_aero, topsw0_aero,
	3145	s solsw_aero, solsw0_aero,
	3146	e cldtaupi,
	3147	s topswai_aero, solswai_aero)
	3148
[955]	3149	#endif
	3150	ELSE
[1534]	3151	c
	3152	cIM calcul radiatif pour le cas actuel
	3153	c
	3154	RCO2 = RCO2_act
	3155	RCH4 = RCH4_act
	3156	RN2O = RN2O_act
	3157	RCFC11 = RCFC11_act
	3158	RCFC12 = RCFC12_act
	3159	c
[1279]	3160	CALL radlwsw
	3161	e (dist, rmu0, fract,
	3162	e paprs, pplay,zxtsol,albsol1, albsol2,
	3163	e t_seri,q_seri,wo,
[1534]	3164	e cldfra, cldemirad, cldtaurad,
[1279]	3165	e ok_ade, ok_aie,
	3166	e tau_aero, piz_aero, cg_aero,
	3167	e cldtaupi,new_aod,
	3168	e zqsat, flwc, fiwc,
	3169	s heat,heat0,cool,cool0,radsol,albpla,
	3170	s topsw,toplw,solsw,sollw,
	3171	s sollwdown,
	3172	s topsw0,toplw0,solsw0,sollw0,
	3173	s lwdn0, lwdn, lwup0, lwup,
	3174	s swdn0, swdn, swup0, swup,
	3175	s topswad_aero, solswad_aero,
	3176	s topswai_aero, solswai_aero,
	3177	o topswad0_aero, solswad0_aero,
	3178	o topsw_aero, topsw0_aero,
	3179	o solsw_aero, solsw0_aero,
	3180	o topswcf_aero, solswcf_aero)
	3181
[1534]	3182	c
	3183	cIM 2eme calcul radiatif pour le cas perturbe ou au moins un
	3184	cIM des taux doit etre different du taux actuel
	3185	cIM Par defaut on a les taux perturbes egaux aux taux actuels
	3186	c
	3187	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.
	3188	$RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.
	3189	$RCFC12_per.NE.RCFC12_act) THEN
	3190	c
	3191	RCO2 = RCO2_per
	3192	RCH4 = RCH4_per
	3193	RN2O = RN2O_per
	3194	RCFC11 = RCFC11_per
	3195	RCFC12 = RCFC12_per
	3196	c
	3197	CALL radlwsw
	3198	e (dist, rmu0, fract,
	3199	e paprs, pplay,zxtsol,albsol1, albsol2,
	3200	e t_seri,q_seri,wo,
	3201	e cldfra, cldemi, cldtau,
	3202	e ok_ade, ok_aie,
	3203	e tau_aero, piz_aero, cg_aero,
	3204	e cldtaupi,new_aod,
	3205	e zqsat, flwc, fiwc,
	3206	s heatp,heat0p,coolp,cool0p,radsolp,albplap,
	3207	s topswp,toplwp,solswp,sollwp,
	3208	s sollwdownp,
	3209	s topsw0p,toplw0p,solsw0p,sollw0p,
	3210	s lwdn0p, lwdnp, lwup0p, lwupp,
	3211	s swdn0p, swdnp, swup0p, swupp,
	3212	s topswad_aerop, solswad_aerop,
	3213	s topswai_aerop, solswai_aerop,
	3214	o topswad0_aerop, solswad0_aerop,
	3215	o topsw_aerop, topsw0_aerop,
	3216	o solsw_aerop, solsw0_aerop,
	3217	o topswcf_aerop, solswcf_aerop)
	3218	endif
	3219	c
[1279]	3220	ENDIF ! aerosol_couple
[524]	3221	itaprad = 0
[1279]	3222	ENDIF ! MOD(itaprad,radpas)
[524]	3223	itaprad = itaprad + 1
[879]	3224
	3225	if (iflag_radia.eq.0) then
	3226	print *,'--------------------------------------------------'
	3227	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	3228	print *,'>>>> heat et cool mis a zero '
	3229	print *,'--------------------------------------------------'
	3230	heat=0.
	3231	cool=0.
	3232	endif
	3233
[524]	3234	c
	3235	c Ajouter la tendance des rayonnements (tous les pas)
	3236	c
	3237	DO k = 1, klev
	3238	DO i = 1, klon
	3239	t_seri(i,k) = t_seri(i,k)
[1035]	3240	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
[524]	3241	ENDDO
	3242	ENDDO
[766]	3243	c
	3244	if (mydebug) then
	3245	call writefield_phy('u_seri',u_seri,llm)
	3246	call writefield_phy('v_seri',v_seri,llm)
	3247	call writefield_phy('t_seri',t_seri,llm)
	3248	call writefield_phy('q_seri',q_seri,llm)
	3249	endif
	3250
[687]	3251	cIM
	3252	IF (ip_ebil_phy.ge.2) THEN
[524]	3253	ztit='after rad'
[687]	3254	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	3255	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3256	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	3257	call diagphy(airephy,ztit,ip_ebil_phy
[524]	3258	e , topsw, toplw, solsw, sollw, zero_v
	3259	e , zero_v, zero_v, zero_v, ztsol
	3260	e , d_h_vcol, d_qt, d_ec
	3261	s , fs_bound, fq_bound )
	3262	END IF
	3263	c
	3264	c
	3265	c Calculer l'hydrologie de la surface
	3266	c
	3267	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	3268	c . agesno, ftsol,fqsurf,fsnow, ruis)
	3269	c
[782]	3270
[524]	3271	c
	3272	c Calculer le bilan du sol et la derive de temperature (couplage)
	3273	c
	3274	DO i = 1, klon
	3275	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	3276	c a la demande de JLD
	3277	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	3278	ENDDO
	3279	c
	3280	cmoddeblott(jan95)
	3281	c Appeler le programme de parametrisation de l'orographie
	3282	c a l'echelle sous-maille:
	3283	c
	3284	IF (ok_orodr) THEN
	3285	c
	3286	c selection des points pour lesquels le shema est actif:
	3287	igwd=0
	3288	DO i=1,klon
	3289	itest(i)=0
	3290	c IF ((zstd(i).gt.10.0)) THEN
	3291	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	3292	itest(i)=1
	3293	igwd=igwd+1
	3294	idx(igwd)=i
	3295	ENDIF
	3296	ENDDO
	3297	c igwdim=MAX(1,igwd)
	3298	c
[1001]	3299	IF (ok_strato) THEN
	3300
	3301	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
	3302	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	3303	e igwd,idx,itest,
	3304	e t_seri, u_seri, v_seri,
	3305	s zulow, zvlow, zustrdr, zvstrdr,
	3306	s d_t_oro, d_u_oro, d_v_oro)
	3307
	3308	ELSE
[524]	3309	CALL drag_noro(klon,klev,dtime,paprs,pplay,
	3310	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	3311	e igwd,idx,itest,
	3312	e t_seri, u_seri, v_seri,
[644]	3313	s zulow, zvlow, zustrdr, zvstrdr,
[524]	3314	s d_t_oro, d_u_oro, d_v_oro)
[1001]	3315	ENDIF
[524]	3316	c
	3317	c ajout des tendances
[904]	3318	!-----------------------------------------------------------------------------------------
	3319	! ajout des tendances de la trainee de l'orographie
	3320	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
	3321	!-----------------------------------------------------------------------------------------
[524]	3322	c
	3323	ENDIF ! fin de test sur ok_orodr
	3324	c
[766]	3325	if (mydebug) then
	3326	call writefield_phy('u_seri',u_seri,llm)
	3327	call writefield_phy('v_seri',v_seri,llm)
	3328	call writefield_phy('t_seri',t_seri,llm)
	3329	call writefield_phy('q_seri',q_seri,llm)
	3330	endif
	3331
[524]	3332	IF (ok_orolf) THEN
	3333	c
	3334	c selection des points pour lesquels le shema est actif:
	3335	igwd=0
	3336	DO i=1,klon
	3337	itest(i)=0
	3338	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	3339	itest(i)=1
	3340	igwd=igwd+1
	3341	idx(igwd)=i
	3342	ENDIF
	3343	ENDDO
	3344	c igwdim=MAX(1,igwd)
	3345	c
[1001]	3346	IF (ok_strato) THEN
	3347
	3348	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
	3349	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
	3350	e igwd,idx,itest,
	3351	e t_seri, u_seri, v_seri,
	3352	s zulow, zvlow, zustrli, zvstrli,
	3353	s d_t_lif, d_u_lif, d_v_lif )
	3354
	3355	ELSE
	3356	CALL lift_noro(klon,klev,dtime,paprs,pplay,
[524]	3357	e rlat,zmea,zstd,zpic,
	3358	e itest,
	3359	e t_seri, u_seri, v_seri,
[644]	3360	s zulow, zvlow, zustrli, zvstrli,
[524]	3361	s d_t_lif, d_u_lif, d_v_lif)
[1001]	3362	ENDIF
	3363	c
[904]	3364	!-----------------------------------------------------------------------------------------
	3365	! ajout des tendances de la portance de l'orographie
	3366	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
	3367	!-----------------------------------------------------------------------------------------
[524]	3368	c
	3369	ENDIF ! fin de test sur ok_orolf
[1001]	3370	C HINES GWD PARAMETRIZATION
	3371
	3372	IF (ok_hines) then
	3373
	3374	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
	3375	i rlat,t_seri,u_seri,v_seri,
	3376	o zustrhi,zvstrhi,
	3377	o d_t_hin, d_u_hin, d_v_hin)
[524]	3378	c
[1001]	3379	c ajout des tendances
[1534]	3380	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
[1001]	3381
	3382	ENDIF
	3383	c
	3384
	3385	c
[644]	3386	cIM cf. FLott BEG
	3387	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	3388
[766]	3389	if (mydebug) then
	3390	call writefield_phy('u_seri',u_seri,llm)
	3391	call writefield_phy('v_seri',v_seri,llm)
	3392	call writefield_phy('t_seri',t_seri,llm)
	3393	call writefield_phy('q_seri',q_seri,llm)
	3394	endif
	3395
[644]	3396	DO i = 1, klon
	3397	zustrph(i)=0.
	3398	zvstrph(i)=0.
	3399	ENDDO
	3400	DO k = 1, klev
	3401	DO i = 1, klon
	3402	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	3403	c (paprs(i,k)-paprs(i,k+1))/rg
	3404	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	3405	c (paprs(i,k)-paprs(i,k+1))/rg
	3406	ENDDO
	3407	ENDDO
	3408	c
	3409	cIM calcul composantes axiales du moment angulaire et couple des montagnes
	3410	c
[1279]	3411	IF (is_sequential) THEN
[766]	3412
[1279]	3413	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
[766]	3414	C ra,rg,romega,
	3415	C rlat,rlon,pphis,
	3416	C zustrdr,zustrli,zustrph,
	3417	C zvstrdr,zvstrli,zvstrph,
	3418	C paprs,u,v,
	3419	C aam, torsfc)
	3420	ENDIF
[644]	3421	cIM cf. FLott END
[687]	3422	cIM
	3423	IF (ip_ebil_phy.ge.2) THEN
[524]	3424	ztit='after orography'
[687]	3425	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	3426	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3427	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3428	END IF
	3429	c
	3430	c
[1279]	3431	!====================================================================
	3432	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
	3433	!====================================================================
	3434	! Abderrahmane 24.08.09
	3435
	3436	IF (ok_cosp) THEN
	3437	! adeclarer
	3438	#ifdef CPP_COSP
	3439	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
	3440
	3441	print*,'freq_cosp',freq_cosp
	3442	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
	3443	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
	3444	! s ref_liq,ref_ice
	3445	call phys_cosp(itap,dtime,freq_cosp,
[1368]	3446	$ ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
	3447	$ ecrit_mth,ecrit_day,ecrit_hf,
	3448	$ klon,klev,rlon,rlat,presnivs,overlap,
[1928]	3449	$ fract,ref_liq,ref_ice,
[1279]	3450	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
[1528]	3451	$ zu10m,zv10m,pphis,
[1279]	3452	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
	3453	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
	3454	$ prfl(:,1:klev),psfl(:,1:klev),
	3455	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
[1534]	3456	$ mr_ozone,cldtaurad, cldemirad)
[1334]	3457
[1279]	3458	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
	3459	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
	3460	! M clMISR,
	3461	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
	3462	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
	3463
	3464	ENDIF
	3465
	3466	#endif
	3467	ENDIF !ok_cosp
	3468	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[524]	3469	cAA
	3470	cAA Installation de l'interface online-offline pour traceurs
	3471	cAA
	3472	c====================================================================
	3473	c Calcul des tendances traceurs
	3474	c====================================================================
	3475	C
[959]	3476
[1279]	3477	call phytrac (
	3478	I itap, days_elapsed+1, jH_cur, debut,
	3479	I lafin, dtime, u, v, t,
	3480	I paprs, pplay, pmfu, pmfd,
	3481	I pen_u, pde_u, pen_d, pde_d,
	3482	I cdragh, coefh, fm_therm, entr_therm,
	3483	I u1, v1, ftsol, pctsrf,
	3484	I rlat, frac_impa, frac_nucl,rlon,
	3485	I presnivs, pphis, pphi, albsol1,
	3486	I qx(:,:,ivap),rhcl, cldfra, rneb,
	3487	I diafra, cldliq, itop_con, ibas_con,
	3488	I pmflxr, pmflxs, prfl, psfl,
	3489	I da, phi, mp, upwd,
	3490	I dnwd, aerosol_couple, flxmass_w,
	3491	I tau_aero, piz_aero, cg_aero, ccm,
	3492	I rfname,
	3493	O tr_seri)
[524]	3494
	3495	IF (offline) THEN
	3496
[541]	3497	print*,'Attention on met a 0 les thermiques pour phystoke'
[524]	3498	call phystokenc (
[1279]	3499	I nlon,klev,pdtphys,rlon,rlat,
[524]	3500	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[541]	3501	I fm_therm,entr_therm,
[1067]	3502	I cdragh,coefh,u1,v1,ftsol,pctsrf,
[524]	3503	I frac_impa, frac_nucl,
	3504	I pphis,airephy,dtime,itap)
	3505
	3506
	3507	ENDIF
	3508
	3509	c
	3510	c Calculer le transport de l'eau et de l'energie (diagnostique)
	3511	c
	3512	CALL transp (paprs,zxtsol,
	3513	e t_seri, q_seri, u_seri, v_seri, zphi,
	3514	s ve, vq, ue, uq)
	3515	c
[687]	3516	cIM global posePB BEG
	3517	IF(1.EQ.0) THEN
[524]	3518	c
[644]	3519	CALL transp_lay (paprs,zxtsol,
	3520	e t_seri, q_seri, u_seri, v_seri, zphi,
	3521	s ve_lay, vq_lay, ue_lay, uq_lay)
[524]	3522	c
[687]	3523	ENDIF !(1.EQ.0) THEN
	3524	cIM global posePB END
[644]	3525	c Accumuler les variables a stocker dans les fichiers histoire:
[524]	3526	c
	3527	c+jld ec_conser
	3528	DO k = 1, klev
	3529	DO i = 1, klon
	3530	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
	3531	d_t_ec(i,k)=0.5/ZRCPD
	3532	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
[1279]	3533	ENDDO
	3534	ENDDO
	3535
	3536	DO k = 1, klev
	3537	DO i = 1, klon
[524]	3538	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	3539	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	3540	END DO
	3541	END DO
	3542	c-jld ec_conser
[687]	3543	cIM
	3544	IF (ip_ebil_phy.ge.1) THEN
[524]	3545	ztit='after physic'
[687]	3546	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	3547	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3548	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3549	C Comme les tendances de la physique sont ajoute dans la dynamique,
	3550	C on devrait avoir que la variation d'entalpie par la dynamique
	3551	C est egale a la variation de la physique au pas de temps precedent.
	3552	C Donc la somme de ces 2 variations devrait etre nulle.
[1279]	3553
[687]	3554	call diagphy(airephy,ztit,ip_ebil_phy
[524]	3555	e , topsw, toplw, solsw, sollw, sens
	3556	e , evap, rain_fall, snow_fall, ztsol
	3557	e , d_h_vcol, d_qt, d_ec
	3558	s , fs_bound, fq_bound )
	3559	C
	3560	d_h_vcol_phy=d_h_vcol
	3561	C
	3562	END IF
	3563	C
	3564	c=======================================================================
	3565	c SORTIES
	3566	c=======================================================================
	3567
[644]	3568	cIM Interpolation sur les niveaux de pression du NMC
	3569	c -------------------------------------------------
[524]	3570	c
[644]	3571	#include "calcul_STDlev.h"
[1352]	3572	twriteSTD(:,:,1)=tsumSTD(:,:,1)
	3573	qwriteSTD(:,:,1)=qsumSTD(:,:,1)
	3574	rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
	3575	phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
	3576	uwriteSTD(:,:,1)=usumSTD(:,:,1)
	3577	vwriteSTD(:,:,1)=vsumSTD(:,:,1)
	3578	wwriteSTD(:,:,1)=wsumSTD(:,:,1)
[1055]	3579
[1352]	3580	twriteSTD(:,:,2)=tsumSTD(:,:,2)
	3581	qwriteSTD(:,:,2)=qsumSTD(:,:,2)
	3582	rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
	3583	phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
	3584	uwriteSTD(:,:,2)=usumSTD(:,:,2)
	3585	vwriteSTD(:,:,2)=vsumSTD(:,:,2)
	3586	wwriteSTD(:,:,2)=wsumSTD(:,:,2)
[1055]	3587
	3588	twriteSTD(:,:,3)=tlevSTD(:,:)
	3589	qwriteSTD(:,:,3)=qlevSTD(:,:)
	3590	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
	3591	phiwriteSTD(:,:,3)=philevSTD(:,:)
	3592	uwriteSTD(:,:,3)=ulevSTD(:,:)
	3593	vwriteSTD(:,:,3)=vlevSTD(:,:)
	3594	wwriteSTD(:,:,3)=wlevSTD(:,:)
	3595
	3596	twriteSTD(:,:,4)=tlevSTD(:,:)
	3597	qwriteSTD(:,:,4)=qlevSTD(:,:)
	3598	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
	3599	phiwriteSTD(:,:,4)=philevSTD(:,:)
	3600	uwriteSTD(:,:,4)=ulevSTD(:,:)
	3601	vwriteSTD(:,:,4)=vlevSTD(:,:)
	3602	wwriteSTD(:,:,4)=wlevSTD(:,:)
[524]	3603	c
[1352]	3604	cIM initialisation 5eme fichier de sortie
	3605	twriteSTD(:,:,5)=tlevSTD(:,:)
	3606	qwriteSTD(:,:,5)=qlevSTD(:,:)
	3607	rhwriteSTD(:,:,5)=rhlevSTD(:,:)
	3608	phiwriteSTD(:,:,5)=philevSTD(:,:)
	3609	uwriteSTD(:,:,5)=ulevSTD(:,:)
	3610	vwriteSTD(:,:,5)=vlevSTD(:,:)
	3611	wwriteSTD(:,:,5)=wlevSTD(:,:)
[1534]	3612	c
	3613	cIM initialisation 6eme fichier de sortie
	3614	twriteSTD(:,:,6)=tlevSTD(:,:)
	3615	qwriteSTD(:,:,6)=qlevSTD(:,:)
	3616	rhwriteSTD(:,:,6)=rhlevSTD(:,:)
	3617	phiwriteSTD(:,:,6)=philevSTD(:,:)
	3618	uwriteSTD(:,:,6)=ulevSTD(:,:)
	3619	vwriteSTD(:,:,6)=vlevSTD(:,:)
	3620	wwriteSTD(:,:,6)=wlevSTD(:,:)
[1352]	3621	cIM for NMC files
	3622	DO n=1, nlevSTD3
	3623	DO k=1, nlevSTD
	3624	if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
	3625	twriteSTD3(:,n)=tlevSTD(:,k)
	3626	qwriteSTD3(:,n)=qlevSTD(:,k)
	3627	rhwriteSTD3(:,n)=rhlevSTD(:,k)
	3628	phiwriteSTD3(:,n)=philevSTD(:,k)
	3629	uwriteSTD3(:,n)=ulevSTD(:,k)
	3630	vwriteSTD3(:,n)=vlevSTD(:,k)
	3631	wwriteSTD3(:,n)=wlevSTD(:,k)
	3632	endif !rlevSTD3(n).EQ.rlevSTD(k)
	3633	ENDDO
	3634	ENDDO
	3635	c
	3636	DO n=1, nlevSTD8
	3637	DO k=1, nlevSTD
	3638	if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
	3639	tnondefSTD8(:,n)=tnondef(:,k,2)
	3640	twriteSTD8(:,n)=tsumSTD(:,k,2)
	3641	qwriteSTD8(:,n)=qsumSTD(:,k,2)
	3642	rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
	3643	phiwriteSTD8(:,n)=phisumSTD(:,k,2)
	3644	uwriteSTD8(:,n)=usumSTD(:,k,2)
	3645	vwriteSTD8(:,n)=vsumSTD(:,k,2)
	3646	wwriteSTD8(:,n)=wsumSTD(:,k,2)
	3647	endif !rlevSTD8(n).EQ.rlevSTD(k)
	3648	ENDDO
	3649	ENDDO
	3650	c
[524]	3651	c slp sea level pressure
	3652	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	3653	c
	3654	ccc prw = eau precipitable
	3655	DO i = 1, klon
	3656	prw(i) = 0.
	3657	DO k = 1, klev
	3658	prw(i) = prw(i) +
	3659	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	3660	ENDDO
	3661	ENDDO
	3662	c
[644]	3663	cIM initialisation + calculs divers diag AMIP2
[524]	3664	c
[644]	3665	#include "calcul_divers.h"
	3666	c
[959]	3667	IF (config_inca /= 'none') THEN
[655]	3668	#ifdef INCA
[959]	3669	CALL VTe(VTphysiq)
	3670	CALL VTb(VTinca)
	3671
[1287]	3672	CALL chemhook_end (
[655]	3673	$ dtime,
	3674	$ pplay,
	3675	$ t_seri,
	3676	$ tr_seri,
	3677	$ nbtr,
	3678	$ paprs,
	3679	$ q_seri,
[791]	3680	$ airephy,
[655]	3681	$ pphi,
	3682	$ pphis,
[766]	3683	$ zx_rh)
[959]	3684
	3685	CALL VTe(VTinca)
	3686	CALL VTb(VTphysiq)
[655]	3687	#endif
[959]	3688	END IF
[655]	3689
[524]	3690	c=============================================================
	3691	c
	3692	c Convertir les incrementations en tendances
	3693	c
[766]	3694	if (mydebug) then
	3695	call writefield_phy('u_seri',u_seri,llm)
	3696	call writefield_phy('v_seri',v_seri,llm)
	3697	call writefield_phy('t_seri',t_seri,llm)
	3698	call writefield_phy('q_seri',q_seri,llm)
	3699	endif
	3700
[524]	3701	DO k = 1, klev
	3702	DO i = 1, klon
	3703	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	3704	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	3705	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	3706	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	3707	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	3708	ENDDO
	3709	ENDDO
	3710	c
[1146]	3711	IF (nqtot.GE.3) THEN
	3712	DO iq = 3, nqtot
[524]	3713	DO k = 1, klev
	3714	DO i = 1, klon
	3715	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	3716	ENDDO
	3717	ENDDO
	3718	ENDDO
	3719	ENDIF
	3720	c
[644]	3721	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
[687]	3722	cIM global posePB#include "write_bilKP_ins.h"
	3723	cIM global posePB#include "write_bilKP_ave.h"
[644]	3724	c
[1334]	3725
[524]	3726	c Sauvegarder les valeurs de t et q a la fin de la physique:
	3727	c
	3728	DO k = 1, klev
	3729	DO i = 1, klon
[1054]	3730	u_ancien(i,k) = u_seri(i,k)
	3731	v_ancien(i,k) = v_seri(i,k)
[524]	3732	t_ancien(i,k) = t_seri(i,k)
	3733	q_ancien(i,k) = q_seri(i,k)
	3734	ENDDO
	3735	ENDDO
	3736	c
[879]	3737	!==========================================================================
	3738	! Sorties des tendances pour un point particulier
	3739	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	3740	! pour le debug
	3741	! La valeur de igout est attribuee plus haut dans le programme
	3742	!==========================================================================
	3743
[942]	3744	if (prt_level.ge.1) then
[879]	3745	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	3746	write(lunout,*)
[1279]	3747	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
[879]	3748	write(lunout,*)
[1279]	3749	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
[930]	3750	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
[929]	3751	s pctsrf(igout,is_sic)
[879]	3752	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
[1279]	3753	do k=1,klev
[879]	3754	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
	3755	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
	3756	s d_t_eva(igout,k)
	3757	enddo
	3758	write(lunout,*) 'cool,heat'
[1279]	3759	do k=1,klev
[879]	3760	write(lunout,*) cool(igout,k),heat(igout,k)
	3761	enddo
	3762
	3763	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
[1279]	3764	do k=1,klev
[879]	3765	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
	3766	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	3767	enddo
	3768
	3769	write(lunout,*) 'd_ps ',d_ps(igout)
	3770	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
[1279]	3771	do k=1,klev
[879]	3772	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
	3773	s d_qx(igout,k,1),d_qx(igout,k,2)
	3774	enddo
	3775	endif
	3776
	3777	!==========================================================================
	3778
	3779	c============================================================
	3780	c Calcul de la temperature potentielle
	3781	c============================================================
	3782	DO k = 1, klev
	3783	DO i = 1, klon
[1398]	3784	cJYG/IM theta en debut du pas de temps
	3785	cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	3786	cJYG/IM theta en fin de pas de temps de physique
	3787	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
[879]	3788	ENDDO
	3789	ENDDO
	3790	c
	3791
[524]	3792	c 22.03.04 BEG
	3793	c=============================================================
	3794	c Ecriture des sorties
	3795	c=============================================================
	3796	#ifdef CPP_IOIPSL
[782]	3797
	3798	c Recupere des varibles calcule dans differents modules
	3799	c pour ecriture dans histxxx.nc
[524]	3800
[888]	3801	! Get some variables from module fonte_neige_mod
[782]	3802	CALL fonte_neige_get_vars(pctsrf,
	3803	. zxfqcalving, zxfqfonte, zxffonte)
	3804
	3805
[909]	3806	#include "phys_output_write.h"
	3807
[524]	3808	#ifdef histISCCP
	3809	#include "write_histISCCP.h"
	3810	#endif
	3811
[1352]	3812	#ifdef histNMC
	3813	#include "write_histhfNMC.h"
	3814	#include "write_histdayNMC.h"
[524]	3815	#include "write_histmthNMC.h"
	3816	#endif
	3817
[687]	3818	#include "write_histday_seri.h"
	3819
	3820	#include "write_paramLMDZ_phy.h"
	3821
[524]	3822	#endif
	3823
	3824	c 22.03.04 END
	3825	c
	3826	c====================================================================
	3827	c Si c'est la fin, il faut conserver l'etat de redemarrage
	3828	c====================================================================
	3829	c
[782]	3830
	3831
[524]	3832	IF (lafin) THEN
	3833	itau_phy = itau_phy + itap
[967]	3834	CALL phyredem ("restartphy.nc")
[1001]	3835	! open(97,form="unformatted",file="finbin")
	3836	! write(97) u_seri,v_seri,t_seri,q_seri
	3837	! close(97)
[1279]	3838	C$OMP MASTER
	3839	if (read_climoz >= 1) then
	3840	if (is_mpi_root) then
	3841	call nf95_close(ncid_climoz)
	3842	end if
	3843	deallocate(press_climoz) ! pointer
	3844	end if
	3845	C$OMP END MASTER
[524]	3846	ENDIF
	3847
[1279]	3848	! first=.false.
[524]	3849
	3850	RETURN
	3851	END
	3852	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	3853	IMPLICIT none
	3854	c
	3855	c Calculer et imprimer l'eau totale. A utiliser pour verifier
	3856	c la conservation de l'eau
	3857	c
	3858	#include "YOMCST.h"
	3859	INTEGER klon,klev
	3860	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	3861	REAL aire(klon)
	3862	REAL qtotal, zx, qcheck
	3863	INTEGER i, k
	3864	c
	3865	zx = 0.0
	3866	DO i = 1, klon
	3867	zx = zx + aire(i)
	3868	ENDDO
	3869	qtotal = 0.0
	3870	DO k = 1, klev
	3871	DO i = 1, klon
	3872	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
	3873	. *(paprs(i,k)-paprs(i,k+1))/RG
	3874	ENDDO
	3875	ENDDO
	3876	c
	3877	qcheck = qtotal/zx
	3878	c
	3879	RETURN
	3880	END
	3881	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	3882	IMPLICIT none
	3883	c
	3884	c Tranformer une variable de la grille physique a
	3885	c la grille d'ecriture
	3886	c
	3887	INTEGER nfield,nlon,iim,jjmp1, jjm
	3888	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	3889	c
	3890	INTEGER i, n, ig
	3891	c
	3892	jjm = jjmp1 - 1
	3893	DO n = 1, nfield
	3894	DO i=1,iim
	3895	ecrit(i,n) = fi(1,n)
	3896	ecrit(i+jjm*iim,n) = fi(nlon,n)
	3897	ENDDO
	3898	DO ig = 1, nlon - 2
	3899	ecrit(iim+ig,n) = fi(1+ig,n)
	3900	ENDDO
	3901	ENDDO
	3902	RETURN
	3903	END

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

Download in other formats:

Original Format