Context Navigation

physiq.F @ 2238

Last change on this file since 2238 was 1428, checked in by Laurent Fairhead, 14 years ago

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG

Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

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

File size: 123.0 KB

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

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

Download in other formats:

Original Format