Context Navigation

physiq.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 123.1 KB

Rev	Line
[1]	1	! $Id: physiq.F 1428 2010-09-13 08:43:37Z fairhead $
	2	c#define IO_DEBUG
	3
	4	SUBROUTINE physiq (nlon,nlev,
	5	. debut,lafin,jD_cur, jH_cur,pdtphys,
	6	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
	7	. u,v,t,qx,
	8	. flxmass_w,
	9	. d_u, d_v, d_t, d_qx, d_ps
	10	. , dudyn
	11	. , PVteta)
	12
	13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
	14	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
	15	USE comgeomphy
	16	USE phys_cal_mod
	17	USE write_field_phy
	18	USE dimphy
	19	USE infotrac
	20	USE mod_grid_phy_lmdz
	21	USE mod_phys_lmdz_para
	22	USE iophy
	23	USE misc_mod, mydebug=>debug
	24	USE vampir
	25	USE pbl_surface_mod, ONLY : pbl_surface
	26	USE change_srf_frac_mod
	27	USE surface_data, ONLY : type_ocean, ok_veget
	28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
	29	USE phys_state_var_mod ! Variables sauvegardees de la physique
	30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
	31	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
	32	USE phys_output_mod
	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
	36	cIM for NMC files
	37	use netcdf, only: nf90_fill_real
	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
	43	USE control_mod
	44
	45
	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 =====================
	59	#define histNMC
	60	c#define histISCCP
	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
	67	c nlev----input-I-nombre de couches verticales, doit etre egale a klev
	68	c debut---input-L-variable logique indiquant le premier passage
	69	c lafin---input-L-variable logique indiquant le dernier passage
	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)
	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)
	84	c flxmass_w -input-R- flux de masse verticale
	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
	90	cIM
	91	c PVteta--output-R-vorticite potentielle a des thetas constantes
	92	c======================================================================
	93	#include "dimensions.h"
	94	integer jjmp1
	95	parameter (jjmp1=jjm+1-1/jjm)
	96	integer iip1
	97	parameter (iip1=iim+1)
	98
	99	#include "regdim.h"
	100	#include "indicesol.h"
	101	#include "dimsoil.h"
	102	#include "clesphys.h"
	103	#include "temps.h"
	104	#include "iniprint.h"
	105	#include "thermcell.h"
	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.)
	111	integer iflag_radia ! active ou non le rayonnement (MPL)
	112	save iflag_radia
	113	c$OMP THREADPRIVATE(iflag_radia)
	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======================================================================
	120	REAL amn, amx
	121	INTEGER igout
	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
	147	c$OMP THREADPRIVATE(ok_journe)
	148	c
	149	LOGICAL ok_mensuel ! sortir le fichier mensuel
	150	save ok_mensuel
	151	c$OMP THREADPRIVATE(ok_mensuel)
	152	c
	153	LOGICAL ok_instan ! sortir le fichier instantane
	154	save ok_instan
	155	c$OMP THREADPRIVATE(ok_instan)
	156	c
	157	LOGICAL ok_LES ! sortir le fichier LES
	158	save ok_LES
	159	c$OMP THREADPRIVATE(ok_LES)
	160	c
	161	LOGICAL ok_region ! sortir le fichier regional
	162	PARAMETER (ok_region=.FALSE.)
	163	c======================================================================
	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)
	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)
	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
	182	c======================================================================
	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
	195	REAL, intent(in):: jD_cur, jH_cur
	196
	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)
	205	real pir
	206
	207	REAL u(klon,klev)
	208	REAL v(klon,klev)
	209	REAL t(klon,klev),theta(klon,klev)
	210	REAL qx(klon,klev,nqtot)
	211	REAL flxmass_w(klon,klev)
	212	REAL omega(klon,klev) ! vitesse verticale en Pa/s
	213	REAL d_u(klon,klev)
	214	REAL d_v(klon,klev)
	215	REAL d_t(klon,klev)
	216	REAL d_qx(klon,klev,nqtot)
	217	REAL d_ps(klon)
	218	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
	219	!IM definition dynamique o_trac dans phys_output_open
	220	! type(ctrl_out) :: o_trac(nqtot)
	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'/
	228	SAVE ctetaSTD
	229	c$OMP THREADPRIVATE(ctetaSTD)
	230	REAL rtetaSTD(nbteta)
	231	DATA rtetaSTD/350., 380., 405./
	232	SAVE rtetaSTD
	233	c$OMP THREADPRIVATE(rtetaSTD)
	234	c
	235	REAL PVteta(klon,nbteta)
	236	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
	237	c
	238	cMI Amip2 PV a theta constante
	239
	240	cym INTEGER klevp1, klevm1
	241	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
	242	cym#include "raddim.h"
	243	c
	244	c
	245	cIM Amip2
	246	c variables a une pression donnee
	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./
	252	SAVE rlevstd
	253	c$OMP THREADPRIVATE(rlevstd)
	254	CHARACTER*4 clevSTD(nlevSTD)
	255	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
	256	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
	257	.'70 ','50 ','30 ','20 ','10 '/
	258	SAVE clevSTD
	259	c$OMP THREADPRIVATE(clevSTD)
	260	c
	261	CHARACTER*4 bb2
	262	CHARACTER*2 bb3
	263
	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)
	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)
	289	c
	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
	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
	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
	322	INTEGER linv, kp1
	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
	332	cIM ISCCP simulator v3.4
	333	c dans clesphys.h top_height, overlap
	334	cv3.4
	335	INTEGER debug, debugcol
	336	cym INTEGER npoints
	337	cym PARAMETER(npoints=klon)
	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
	349	cIM parametres ISCCP BEG
	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
	354	INTEGER n
	355	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
	356	DATA ifreq_isccp/3/
	357	SAVE ifreq_isccp
	358	c$OMP THREADPRIVATE(ifreq_isccp)
	359	CHARACTER*5 typinout(napisccp)
	360	DATA typinout/'i3od'/
	361	SAVE typinout
	362	c$OMP THREADPRIVATE(typinout)
	363	cIM verif boxptop BEG
	364	CHARACTER*1 verticaxe(napisccp)
	365	DATA verticaxe/'1'/
	366	SAVE verticaxe
	367	c$OMP THREADPRIVATE(verticaxe)
	368	cIM verif boxptop END
	369	INTEGER nvlev(napisccp)
	370	c INTEGER nvlev
	371	REAL t1, aa
	372	REAL seed_re(klon,napisccp)
	373	cym !!!! A voir plus tard
	374	cym INTEGER iphy(iim,jjmp1)
	375	cIM parametres ISCCP END
	376	c
	377	c ncol = nb. de sous-colonnes pour chaque maille du GCM
	378	c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
	379	c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
	380	INTEGER,SAVE :: ncol(napisccp)
	381	c$OMP THREADPRIVATE(ncol)
	382	INTEGER ncolmx, seed(klon,napisccp)
	383	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
	384	c PARAMETER(ncolmx=1500)
	385	PARAMETER(ncolmx=300)
	386	c
	387	cIM verif boxptop BEG
	388	REAL vertlev(ncolmx,napisccp)
	389	cIM verif boxptop END
	390	c
	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)
	394	REAL emsfc_lw
	395	PARAMETER(emsfc_lw=0.99)
	396	c REAL ran0 ! type for random number fuction
	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
	417	INTEGER kmax, lmax, lmax3
	418	PARAMETER(kmax=8, lmax=8, lmax3=3)
	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
	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
	439	INTEGER iw, iwmax
	440	REAL wmin, pas_w
	441	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
	442	cIM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
	443	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
	444	REAL o500(klon)
	445	c
	446
	447	c sorties ISCCP
	448
	449	integer nid_isccp
	450	save nid_isccp
	451	c$OMP THREADPRIVATE(nid_isccp)
	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./
	455	SAVE zx_tau
	456	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
	457	SAVE zx_pc
	458	c$OMP THREADPRIVATE(zx_tau,zx_pc)
	459	c cldtopres pression au sommet des nuages
	460	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
	461	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
	462	DATA cldtopres3/440., 680., 1000./
	463	SAVE cldtopres,cldtopres3
	464	c$OMP THREADPRIVATE(cldtopres,cldtopres3)
	465	cIM 051005 BEG
	466	INTEGER komega, nhoriRD
	467
	468	c taulev: numero du niveau de tau dans les sorties ISCCP
	469	CHARACTER *4 taulev(kmaxm1)
	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'/
	474	SAVE taulev,pclev
	475	c$OMP THREADPRIVATE(taulev,pclev)
	476	c
	477	c cnameisccp
	478	CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
	479	cIM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
	480	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
	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',
	486	. 'pc= 800-1000hPa, tau< 0.3',
	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',
	493	. 'pc= 800-1000hPa, tau= 0.3-1.3',
	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',
	500	. 'pc= 800-1000hPa, tau= 1.3-3.6',
	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',
	507	. 'pc= 800-1000hPa, tau= 3.6-9.4',
	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',
	514	. 'pc= 800-1000hPa, tau= 9.4-23',
	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',
	521	. 'pc= 800-1000hPa, tau= 23-60',
	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.',
	527	. 'pc= 680-800hPa, tau> 60.',
	528	. 'pc= 800-1000hPa, tau> 60.'/
	529	SAVE cnameisccp
	530	c$OMP THREADPRIVATE(cnameisccp)
	531	c
	532	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
	533	c INTEGER nhorix7
	534	cIM: region='3d' <==> sorties en global
	535	CHARACTER*3 region
	536	PARAMETER(region='3d')
	537	c
	538	cIM ISCCP simulator v3.4
	539	c
	540	logical ok_hf
	541	c
	542	integer nid_hf, nid_hf3d
	543	save ok_hf, nid_hf, nid_hf3d
	544	c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
	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
	554	c$OMP THREADPRIVATE(itap)
	555	c
	556	real slp(klon) ! sea level pressure
	557	c
	558	REAL fevap(klon,nbsrf)
	559	REAL fluxlat(klon,nbsrf)
	560	c
	561	REAL qsol(klon)
	562	REAL,save :: solarlong0
	563	c$OMP THREADPRIVATE(solarlong0)
	564
	565	c
	566	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	567	c
	568	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	569	REAL zulow(klon),zvlow(klon)
	570	c
	571	INTEGER igwd,idx(klon),itest(klon)
	572	c
	573	REAL agesno(klon,nbsrf)
	574	c
	575	c REAL,allocatable,save :: run_off_lic_0(:)
	576	cc$OMP THREADPRIVATE(run_off_lic_0)
	577	cym SAVE run_off_lic_0
	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
	584	c$OMP THREADPRIVATE(bas, top)
	585
	586	REAL wdn(klon), tdn(klon), qdn(klon)
	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
	593	REAL Vprecip(klon,klev+1) ! precipitation vertical profile
	594	c
	595	REAL wape_prescr, fip_prescr
	596	INTEGER it_wape_prescr
	597	SAVE wape_prescr, fip_prescr, it_wape_prescr
	598	c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
	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)
	604
	605	real, save :: alp_bl_prescr=0.
	606	real, save :: ale_bl_prescr=0.
	607
	608	real, save :: ale_max=100.
	609	real, save :: alp_max=2.
	610
	611	c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
	612	c$OMP THREADPRIVATE(ale_max,alp_max)
	613
	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
	635	c
	636	cpourquoi y'a pas de save??
	637	REAL wake_h(klon) ! Wake : hauteur de la poche froide
	638	c
	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
	644	cjyg
	645	ccc REAL wake_pe(klon) ! Wake potential energy - WAPE
	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)
	661	REAL, SAVE :: alp_offset
	662	c$OMP THREADPRIVATE(alp_offset)
	663
	664	c
	665	cRR:fin declarations poches froides
	666	c=======================================================================================================
	667
	668	REAL zw2(klon,klev+1)
	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)
	674
	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
	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
	691
	692	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
	693
	694	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
	695	c@$$ PARAMETER (offline=.false.)
	696	c@$$ INTEGER physid
	697	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	698	REAL frac_nucl(klon,klev) ! idem (nucleation)
	699	INTEGER :: iii
	700	REAL :: calday
	701
	702	cIM cf FH pour Tiedtke 080604
	703	REAL rain_tiedtke(klon),snow_tiedtke(klon)
	704	c
	705	cIM 050204 END
	706	REAL evap(klon), devap(klon) ! evaporation et sa derivee
	707	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	708
	709	REAL bils(klon) ! bilan de chaleur au sol
	710	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
	711	C ! type de sous-surface et pondere par la fraction
	712	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
	713	C ! type de sous-surface et pondere par la fraction
	714	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
	715
	716	REAL fder(klon)
	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
	722	REAL frugs(klon,nbsrf)
	723	REAL zxrugs(klon) ! longueur de rugosite
	724	c
	725	c Conditions aux limites
	726	c
	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
	735	c
	736	INTEGER lmt_pas
	737	SAVE lmt_pas ! frequence de mise a jour
	738	c$OMP THREADPRIVATE(lmt_pas)
	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
	742
	743	cIM sorties
	744	REAL un_jour
	745	PARAMETER(un_jour=86400.)
	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
	763	CC EXTERNAL o3cm ! initialiser l'ozone
	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
	767	EXTERNAL suphel ! initialiser certaines constantes
	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
	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)
	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
	805	c
	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
	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
	813	c
	814	INTEGER itaprad
	815	SAVE itaprad
	816	c$OMP THREADPRIVATE(itaprad)
	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)
	825	REAL zxsnow_dummy(klon)
	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)
	840	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
	841	REAL t_coup
	842	PARAMETER (t_coup=234.0)
	843	c
	844	REAL zphi(klon,klev)
	845	cym A voir plus tard !!
	846	cym REAL zx_relief(iim,jjmp1)
	847	cym REAL zx_aire(iim,jjmp1)
	848	c
	849	c Grandeurs de sorties
	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)
	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
	865
	866	c -- convect43:
	867	INTEGER ntra ! nb traceurs pour convect4.3
	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)
	886
	887	! tendance nulles
	888	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
	889
	890	c
	891	*********************************************************
	892	* declarations
	893
	894	*********************************************************
	895	cIM 081204 END
	896	c
	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
	904	REAL rain_lsc(klon)
	905	REAL snow_lsc(klon)
	906	c
	907	REAL ratqss(klon,klev),ratqsc(klon,klev)
	908	real ratqsbas,ratqshaut,tau_ratqs
	909	save ratqsbas,ratqshaut,tau_ratqs
	910	c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
	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
	918	real ref_liq(klon,klev), ref_ice(klon,klev)
	919	c$OMP THREADPRIVATE(ok_newmicro)
	920	save fact_cldcon,facttemps
	921	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
	922	real facteur
	923
	924	integer iflag_cldcon
	925	save iflag_cldcon
	926	c$OMP THREADPRIVATE(iflag_cldcon)
	927	logical ptconv(klon,klev)
	928	cIM cf. AM 081204 BEG
	929	logical ptconvth(klon,klev)
	930	cIM cf. AM 081204 END
	931	c
	932	c Variables liees a l'ecriture de la bande histoire physique
	933	c
	934	c======================================================================
	935	c
	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
	940	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
	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
	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)
	958	cIM RH a 2m (la surface)
	959	REAL rh2m(klon), qsat2m(klon)
	960	REAL tpot(klon), tpote(klon)
	961	REAL Lheat
	962
	963	INTEGER length
	964	PARAMETER ( length = 100 )
	965	REAL tabcntr0( length )
	966	c
	967	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
	968	cIM
	969	INTEGER ndex2d1(iwmax)
	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)
	977	REAL zustrhi(klon), zvstrhi(klon)
	978	REAL aam, torsfc
	979	cIM 141004 END
	980	cIM 190504 BEG
	981	INTEGER ij, imp1jmp1
	982	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
	983	cym A voir plus tard
	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
	993	cym A voir plus tard
	994	cym REAL zm_wo(jjmp1, klev)
	995	cIM AMIP2 END
	996	c
	997	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	998	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
	999	c#ifdef histNMC
	1000	cym A voir plus tard !!!!
	1001	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
	1002	REAL zx_tmp_fiNC(klon,nlevSTD)
	1003	c#endif
	1004	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
	1005	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
	1006	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
	1007	cIM for NMC files
	1008	REAL missing_val
	1009	REAL, SAVE :: freq_moyNMC(nout)
	1010	c$OMP THREADPRIVATE(freq_moyNMC)
	1011	c
	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)
	1018	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
	1019	c
	1020	cIM 280405 BEG
	1021	INTEGER nid_bilKPins, nid_bilKPave
	1022	SAVE nid_bilKPins, nid_bilKPave
	1023	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
	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
	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
	1038	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
	1039
	1040	real zjulian
	1041	save zjulian
	1042	c$OMP THREADPRIVATE(zjulian)
	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
	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)
	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
	1072	c$OMP THREADPRIVATE(d_h_vcol_phy)
	1073	REAL zero_v(klon)
	1074	CHARACTER*15 ztit
	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)
	1082	c+jld ec_conser
	1083	REAL ZRCPD
	1084	c-jld ec_conser
	1085	REAL t2m(klon,nbsrf) ! temperature a 2m
	1086	REAL q2m(klon,nbsrf) ! humidite a 2m
	1087
	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
	1092	CHARACTER*40 tinst, tave, typeval
	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
	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
	1106
	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)
	1110	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
	1111	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
	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)
	1115	INTEGER, SAVE :: flag_aerosol
	1116	c$OMP THREADPRIVATE(flag_aerosol)
	1117	LOGICAL, SAVE :: new_aod
	1118	c$OMP THREADPRIVATE(new_aod)
	1119
	1120	c
	1121	c Declaration des constantes et des fonctions thermodynamiques
	1122	c
	1123	LOGICAL,SAVE :: first=.true.
	1124	c$OMP THREADPRIVATE(first)
	1125
	1126	integer iunit
	1127
	1128	integer, save:: read_climoz ! read ozone climatology
	1129	C (let it keep the default OpenMP shared attribute)
	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
	1137	C (let it keep the default OpenMP shared attribute)
	1138
	1139	real, pointer, save:: press_climoz(:)
	1140	C (let it keep the default OpenMP shared attribute)
	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
	1152	#include "YOMCST.h"
	1153	#include "YOETHF.h"
	1154	#include "FCTTRE.h"
	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
	1162	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1163	c Declarations pour Simulateur COSP
	1164	c============================================================
	1165	real :: mr_ozone(klon,klev)
	1166	cIM for NMC files
	1167	missing_val=nf90_fill_real
	1168	c======================================================================
	1169	! Gestion calendrier : mise a jour du module phys_cal_mod
	1170	!
	1171	CALL phys_cal_update(jD_cur,jH_cur)
	1172
	1173	c======================================================================
	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======================================================================
	1178
	1179	if (prt_level.ge.1) then
	1180	igout=klon/2+1/klon
	1181	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1182	write(lunout,*)
	1183	s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
	1184	write(lunout,*)
	1185	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
	1186
	1187	write(lunout,*) 'paprs, play, phi, u, v, t'
	1188	do k=1,klev
	1189	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
	1190	s u(igout,k),v(igout,k),t(igout,k)
	1191	enddo
	1192	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
	1193	do k=1,klev
	1194	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
	1195	enddo
	1196	endif
	1197
	1198	c======================================================================
	1199
	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.
	1207
	1208	torsfc=0.
	1209	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
	1210
	1211	if (first) then
	1212
	1213	cCR:nvelles variables convection/poches froides
	1214
	1215	print*, '================================================='
	1216	print*, 'Allocation des variables locales et sauvegardees'
	1217	call phys_local_var_init
	1218	c
	1219	pasphys=pdtphys
	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
	1233	. iflag_coupl,iflag_clos,iflag_wake, read_climoz,
	1234	& alp_offset)
	1235	call phys_state_var_init(read_climoz)
	1236	call phys_output_var_init
	1237	print*, '================================================='
	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
	1245	cIM beg
	1246	dnwd0=0.0
	1247	ftd=0.0
	1248	fqd=0.0
	1249	cin=0.
	1250	cym Attention pbase pas initialise dans concvl !!!!
	1251	pbase=0
	1252	paire_ter(:)=0.
	1253	cIM 180608
	1254	c pmflxr=0.
	1255	c pmflxs=0.
	1256
	1257	itau_con=0
	1258	first=.false.
	1259
	1260	endif ! first
	1261
	1262	modname = 'physiq'
	1263	cIM
	1264	IF (ip_ebil_phy.ge.1) THEN
	1265	DO i=1,klon
	1266	zero_v(i)=0.
	1267	END DO
	1268	END IF
	1269	ok_sync=.TRUE.
	1270
	1271	IF (debut) THEN
	1272	CALL suphel ! initialiser constantes et parametres phys.
	1273	ENDIF
	1274
	1275	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
	1276
	1277
	1278	c======================================================================
	1279	! Gestion calendrier : mise a jour du module phys_cal_mod
	1280	!
	1281	c CALL phys_cal_update(jD_cur,jH_cur)
	1282
	1283	c
	1284	c Si c'est le debut, il faut initialiser plusieurs choses
	1285	c ********
	1286	c
	1287	IF (debut) THEN
	1288	!rv
	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
	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.
	1302
	1303	lambda_th(:,:)=0.
	1304	wmax_th(:)=0.
	1305	tau_overturning_th(:)=0.
	1306
	1307	IF (config_inca /= 'none') THEN
	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.
	1313	END IF
	1314
	1315	rnebcon0(:,:) = 0.0
	1316	clwcon0(:,:) = 0.0
	1317	rnebcon(:,:) = 0.0
	1318	clwcon(:,:) = 0.0
	1319
	1320	cIM
	1321	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
	1322	c
	1323	print*,'iflag_coupl,iflag_clos,iflag_wake',
	1324	. iflag_coupl,iflag_clos,iflag_wake
	1325	print*,'CYCLE_DIURNE', cycle_diurne
	1326	c
	1327	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
	1328	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
	1329	CALL abort_gcm (modname,abort_message,1)
	1330	ENDIF
	1331	c
	1332	IF(ok_isccp.AND.iflag_con.LE.2) THEN
	1333	abort_message = 'ISCCP-like outputs may be available for KE
	1334	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
	1335	CALL abort_gcm (modname,abort_message,1)
	1336	ENDIF
	1337	c
	1338	c Initialiser les compteurs:
	1339	c
	1340	itap = 0
	1341	itaprad = 0
	1342
	1343	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1344	!! Un petit travail à faire ici.
	1345	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1346
	1347	if (iflag_pbl>1) then
	1348	PRINT*, "Using method MELLOR&YAMADA"
	1349	endif
	1350
	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
	1361	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
	1362	cIM begin
	1363	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
	1364	$,ratqs(1,1)
	1365	cIM end
	1366
	1367
	1368
	1369	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	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
	1377	cIM cf. AM 081204 BEG
	1378	PRINT*,'cycle_diurne3 =',cycle_diurne
	1379	cIM cf. AM 081204 END
	1380	c
	1381	CALL printflag( tabcntr0,radpas,ok_journe,
	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 '
	1388	! call abort_gcm(modname,abort_message,1)
	1389	dtime=pdtphys
	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
	1404	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
	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 "
	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"
	1423
	1424	DO i = 1, klon
	1425	ema_cbmf(i) = 0.
	1426	ema_pcb(i) = 0.
	1427	ema_pct(i) = 0.
	1428	c ema_workcbmf(i) = 0.
	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
	1433	itop_con(i) = 1
	1434	ENDDO
	1435	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	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
	1445
	1446	do i = 1,klon
	1447	Ale_bl(i)=0.
	1448	Alp_bl(i)=0.
	1449	enddo
	1450
	1451	c================================================================================
	1452
	1453	ENDIF !debut
	1454
	1455	DO i=1,klon
	1456	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1457	ENDDO
	1458
	1459	c34EK
	1460	IF (ok_orodr) THEN
	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	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
	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
	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
	1492	cIM 030306 END
	1493
	1494	capemaxcels = 't_max(X)'
	1495	t2mincels = 't_min(X)'
	1496	t2maxcels = 't_max(X)'
	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
	1502	c
	1503	c=============================================================
	1504	c Initialisation des sorties
	1505	c=============================================================
	1506
	1507	#ifdef CPP_IOIPSL
	1508
	1509	c$OMP MASTER
	1510	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
	1511	& ctetaSTD,dtime,ok_veget,
	1512	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
	1513	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
	1514	& read_climoz, new_aod, aerosol_couple
	1515	& )
	1516	c$OMP END MASTER
	1517	c$OMP BARRIER
	1518
	1519	#ifdef histISCCP
	1520	#include "ini_histISCCP.h"
	1521	#endif
	1522
	1523	#ifdef histNMC
	1524	#include "ini_histhfNMC.h"
	1525	#include "ini_histdayNMC.h"
	1526	#include "ini_histmthNMC.h"
	1527	#endif
	1528
	1529	#include "ini_histday_seri.h"
	1530
	1531	#include "ini_paramLMDZ_phy.h"
	1532
	1533	#endif
	1534
	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
	1539	cIM
	1540	IF(ecrit_day.LE.1.) THEN
	1541	ecrit_day = ecrit_day * un_jour !en secondes
	1542	ENDIF
	1543	cIM
	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
	1557	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1558	date0 = jD_ref
	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
	1572	IF (config_inca /= 'none') THEN
	1573	#ifdef INCA
	1574	CALL VTe(VTphysiq)
	1575	CALL VTb(VTinca)
	1576	! iii = MOD(NINT(xjour),360)
	1577	! calday = REAL(iii) + jH_cur
	1578	calday = REAL(days_elapsed) + jH_cur
	1579	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
	1580
	1581	CALL chemini(
	1582	$ rg,
	1583	$ ra,
	1584	$ airephy,
	1585	$ rlat,
	1586	$ rlon,
	1587	$ presnivs,
	1588	$ calday,
	1589	$ klon,
	1590	$ nqtot,
	1591	$ pdtphys,
	1592	$ annee_ref,
	1593	$ day_ref,
	1594	$ itau_phy)
	1595
	1596	CALL VTe(VTinca)
	1597	CALL VTb(VTphysiq)
	1598	#endif
	1599	END IF
	1600	c
	1601	c
	1602	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1603	! Nouvelle initialisation pour le rayonnement RRTM
	1604	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1605
	1606	call iniradia(klon,klev,paprs(1,1:klev+1))
	1607
	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
	1613	ENDIF
	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	!
	1626	CALL change_srf_frac(itap, dtime, days_elapsed+1,
	1627	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
	1628
	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.
	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
	1648	DO iq = 1, nqtot
	1649	DO k = 1, klev
	1650	DO i = 1, klon
	1651	d_qx(i,k,iq) = 0.0
	1652	ENDDO
	1653	ENDDO
	1654	ENDDO
	1655	da(:,:)=0.
	1656	mp(:,:)=0.
	1657	phi(:,:,:)=0.
	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
	1671	IF (nqtot.GE.3) THEN
	1672	DO iq = 3, nqtot
	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
	1695	cIM
	1696	IF (ip_ebil_phy.ge.1) THEN
	1697	ztit='after dynamic'
	1698	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
	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.
	1705	call diagphy(airephy,ztit,ip_ebil_phy
	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
	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
	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
	1726	d_u_dyn(i,k) = 0.0
	1727	d_v_dyn(i,k) = 0.0
	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
	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
	1758	CALL hgardfou(t_seri,ftsol,'debutphy')
	1759	c
	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
	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
	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))
	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
	1829	cIM
	1830	IF (ip_ebil_phy.ge.2) THEN
	1831	ztit='after reevap'
	1832	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
	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)
	1835	call diagphy(airephy,ztit,ip_ebil_phy
	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
	1842
	1843	c
	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
	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
	1857	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	1858	! solarlong0
	1859	if (solarlong0<-999.) then
	1860	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
	1861	else
	1862	zlongi=solarlong0 ! longitude solaire vraie
	1863	dist=1. ! distance au soleil / moyenne
	1864	endif
	1865	endif
	1866	if(prt_level.ge.1) &
	1867	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
	1868
	1869	! Avec ou sans cycle diurne
	1870	IF (cycle_diurne) THEN
	1871	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
	1872	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
	1873	ELSE
	1874	CALL angle(zlongi, rlat, fract, rmu0)
	1875	ENDIF
	1876
	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)
	1881	call writefield_phy('q_seri',q_seri,llm)
	1882	endif
	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,
	1899	c
	1900	c Certains ne sont pas utiliser du tout :
	1901	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
	1902	c
	1903
	1904	CALL pbl_surface(
	1905	e dtime, date0, itap, days_elapsed+1,
	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,
	1911	+ ftsol, falb1, falb2, u10m, v10m,
	1912	s sollwdown, cdragh, cdragm, u1, v1,
	1913	s albsol1, albsol2, sens, evap,
	1914	s zxtsol, zxfluxlat, zt2m, qsat2m,
	1915	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
	1916	s coefh, slab_wfbils,
	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,
	1926	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
	1927
	1928
	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	!-----------------------------------------------------------------------------------------
	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)
	1938	call writefield_phy('q_seri',q_seri,llm)
	1939	endif
	1940
	1941
	1942	IF (ip_ebil_phy.ge.2) THEN
	1943	ztit='after surface_main'
	1944	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	1947	call diagphy(airephy,ztit,ip_ebil_phy
	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
	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
	1977	if (prt_level.ge.1) then
	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
	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
	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
	2016	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
	2017	$ omega(igout, :)
	2018
	2019	IF (iflag_con.EQ.1) THEN
	2020	abort_message ='reactiver le call conlmd dans physiq.F'
	2021	CALL abort_gcm (modname,abort_message,1)
	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,
	2027	e conv_t, conv_q, -evap, omega,
	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)
	2031	d_u_con = 0.
	2032	d_v_con = 0.
	2033
	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:
	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
	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
	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
	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
	2102	if (debut.and.prt_level.gt.9)
	2103	$ WRITE(lunout,*)'ALE et ALP imposes'
	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
	2113	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
	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
	2141	endif
	2142	do i=1,klon
	2143	if (alp(i)>alp_max) then
	2144	IF(prt_level>9)WRITE(lunout,*) &
	2145	& 'WARNING SUPER ALP (seuil=',alp_max,
	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
	2150	IF(prt_level>9)WRITE(lunout,*) &
	2151	& 'WARNING SUPER ALE (seuil=',ale_max,
	2152	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
	2153	ale(i)=ale_max
	2154	endif
	2155	enddo
	2156
	2157	cfin calcul ale et alp
	2158	c=================================================================================================
	2159
	2160
	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
	2167	CALL concvl (iflag_con,iflag_clos,
	2168	. dtime,paprs,pplay,t_undi,q_undi,
	2169	. t_wake,q_wake,wake_s,
	2170	. u_seri,v_seri,tr_seri,nbtr,
	2171	. ALE,ALP,
	2172	. ema_work1,ema_work2,
	2173	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	2174	. rain_con, snow_con, ibas_con, itop_con, sigd,
	2175	. ema_cbmf,upwd,dnwd,dnwd0,
	2176	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
	2177	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
	2178	. pmflxr,pmflxs,da,phi,mp,
	2179	. ftd,fqd,lalim_conv,wght_th)
	2180
	2181	cIM begin
	2182	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2183	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
	2184	cIM end
	2185	cIM cf. FH
	2186	clwcon0=qcondc
	2187	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
	2188
	2189	do i = 1, klon
	2190	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2191	enddo
	2192
	2193	ELSE ! ok_cvl
	2194
	2195	c MAF conema3 ne contient pas les traceurs
	2196	CALL conema3 (dtime,
	2197	. paprs,pplay,t_seri,q_seri,
	2198	. u_seri,v_seri,tr_seri,ntra,
	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
	2210	c
	2211	c Correction precip
	2212	rain_con = rain_con * cvl_corr
	2213	snow_con = snow_con * cvl_corr
	2214	c
	2215
	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
	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
	2239	ema_pcb(i) = paprs(i,ibas_con(i))
	2240	ENDDO
	2241	DO i = 1, klon
	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
	2246	if (itop_con(i).gt.klev-3) then
	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
	2251	endif
	2252	ENDDO
	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
	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
	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	!-----------------------------------------------------------------------------------------
	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)
	2282	call writefield_phy('q_seri',q_seri,llm)
	2283	endif
	2284
	2285	cIM
	2286	IF (ip_ebil_phy.ge.2) THEN
	2287	ztit='after convect'
	2288	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	2291	call diagphy(airephy,ztit,ip_ebil_phy
	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
	2304	za = za + airephy(i)/REAL(klon)
	2305	zx_t = zx_t + (rain_con(i)+
	2306	. snow_con(i))*airephy(i)/REAL(klon)
	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.
	2335
	2336	c
	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)
	2345	wdt_PBL(i,k) = 0.
	2346	dq_dwn(i,k) = fqd(i,k)
	2347	wdq_PBL(i,k) = 0.
	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)
	2351	udt_PBL(i,k) = 0.
	2352	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2353	udq_PBL(i,k) = 0.
	2354	ENDDO
	2355	ENDDO
	2356	c
	2357	ccalcul caracteristiques de la poche froide
	2358	call calWAKE (paprs,pplay,dtime
	2359	: ,t_seri,q_seri,omega
	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
	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	!-----------------------------------------------------------------------------------------
	2384
	2385	endif
	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
	2401	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2402	c
	2403	c===================================================================
	2404	c Convection seche (thermiques ou ajustement)
	2405	c===================================================================
	2406	c
	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.
	2414	d_t_ajs(:,:)=0.
	2415	d_u_ajs(:,:)=0.
	2416	d_v_ajs(:,:)=0.
	2417	d_q_ajs(:,:)=0.
	2418	clwcon0th(:,:)=0.
	2419	c
	2420	c fm_therm(:,:)=0.
	2421	c entr_therm(:,:)=0.
	2422	c detr_therm(:,:)=0.
	2423	c
	2424	IF(prt_level>9)WRITE(lunout,*)
	2425	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
	2426	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2427	if(iflag_thermals.lt.0) then
	2428	c Rien
	2429	c ====
	2430	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	2431
	2432
	2433	else
	2434
	2435	c Thermiques
	2436	c ==========
	2437	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
	2438	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2439
	2440
	2441	if (iflag_thermals.gt.1) then
	2442	call calltherm(pdtphys
	2443	s ,pplay,paprs,pphi,weak_inversion
	2444	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
	2445	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
	2446	s ,fm_therm,entr_therm,detr_therm
	2447	s ,zqasc,clwcon0th,lmax_th,ratqscth
	2448	s ,ratqsdiff,zqsatth
	2449	con rajoute ale et alp, et les caracteristiques de la couche alim
	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
	2462	endif
	2463
	2464
	2465
	2466	c Ajustement sec
	2467	c ==============
	2468
	2469	! Dans le cas oÃ¹ on active les thermiques, on fait partir l'ajustement
	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
	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')
	2499	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2500	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2501
	2502	!-----------------------------------------------------------------------------------------
	2503
	2504	endif
	2505
	2506	endif
	2507	c
	2508	c===================================================================
	2509	cIM
	2510	IF (ip_ebil_phy.ge.2) THEN
	2511	ztit='after dry_adjust'
	2512	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	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 )
	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
	2542
	2543	c-----------------------------------------------------------------------
	2544	c par nversion de la fonction log normale
	2545	c-----------------------------------------------------------------------
	2546	else if (iflag_cldcon.eq.4) then
	2547	ptconvth(:,:)=.false.
	2548	ratqsc(:,:)=0.
	2549	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
	2550	call clouds_gno
	2551	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
	2552	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
	2553
	2554	endif
	2555
	2556	c ratqs stables
	2557	c -------------
	2558
	2559	if (iflag_ratqs.eq.0) then
	2560
	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
	2613	c ratqs final
	2614	c -----------
	2615
	2616	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
	2617	s .or.iflag_cldcon.eq.4) then
	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
	2630	! les ratqs sont une combinaison de ratqss et ratqsc
	2631	if(prt_level.ge.9)
	2632	$ write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
	2633
	2634	if (tau_ratqs>1.e-10) then
	2635	facteur=exp(-pdtphys/tau_ratqs)
	2636	else
	2637	facteur=0.
	2638	endif
	2639	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
	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(:,:))
	2647	else
	2648	! on ne prend que le ratqs stable pour fisrtilp
	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,
	2663	. prfl, psfl, rhcl,
	2664	. zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
	2665
	2666	WHERE (rain_lsc < 0) rain_lsc = 0.
	2667	WHERE (snow_lsc < 0) snow_lsc = 0.
	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	!-----------------------------------------------------------------------------------------
	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
	2684	za = za + airephy(i)/REAL(klon)
	2685	zx_t = zx_t + (rain_lsc(i)
	2686	. + snow_lsc(i))*airephy(i)/REAL(klon)
	2687	ENDDO
	2688	zx_t = zx_t/za*dtime
	2689	WRITE(lunout,*)"Precip=", zx_t
	2690	ENDIF
	2691	cIM
	2692	IF (ip_ebil_phy.ge.2) THEN
	2693	ztit='after fisrt'
	2694	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	2697	call diagphy(airephy,ztit,ip_ebil_phy
	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
	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)
	2708	call writefield_phy('q_seri',q_seri,llm)
	2709	endif
	2710
	2711	c
	2712	c-------------------------------------------------------------------
	2713	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	2714	c-------------------------------------------------------------------
	2715
	2716	c 1. NUAGES CONVECTIFS
	2717	c
	2718	cIM cf FH
	2719	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
	2720	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
	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
	2742
	2743	c Nuages diagnostiques pour Tiedtke
	2744	CALL diagcld1(paprs,pplay,
	2745	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
	2746	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
	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
	2757	ELSE IF (iflag_cldcon.ge.3) THEN
	2758	c On prend pour les nuages convectifs le max du calcul de la
	2759	c convection et du calcul du pas de temps precedent diminue d'un facteur
	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
	2773	c
	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
	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)
	2784	ELSE
	2785	cIM 170310 BEG
	2786	tausum_aero(:,:,:) = 0.
	2787	cIM 170310 END
	2788	tau_aero(:,:,:,:) = 0.
	2789	piz_aero(:,:,:,:) = 0.
	2790	cg_aero(:,:,:,:) = 0.
	2791	ENDIF
	2792
	2793	cIM calcul nuages par le simulateur ISCCP
	2794	c
	2795	#ifdef histISCCP
	2796	IF (ok_isccp) THEN
	2797	c
	2798	cIM lecture invtau, tautab des fichiers formattes
	2799	c
	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
	2821	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	2822	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
	2823	#include "calcul_simulISCCP.h"
	2824	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
	2825	ENDIF !ok_isccp
	2826	#endif
	2827
	2828	c On prend la somme des fractions nuageuses et des contenus en eau
	2829
	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
	2847	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2848	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	2849	endif
	2850
	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
	2873	cIM
	2874	IF (ip_ebil_phy.ge.2) THEN
	2875	ztit="after diagcld"
	2876	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	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 )
	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
	2908
	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
	2916	Lheat=RLSTT
	2917	ELSE
	2918	Lheat=RLVTT
	2919	ENDIF
	2920	ELSE
	2921	IF (zt2m(i).LT.RTT) THEN
	2922	Lheat=RLSTT
	2923	ELSE
	2924	Lheat=RLVTT
	2925	ENDIF
	2926	ENDIF
	2927	tpote(i) = tpot(i)*
	2928	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	2929	ENDDO
	2930
	2931	IF (config_inca /= 'none') THEN
	2932	#ifdef INCA
	2933	CALL VTe(VTphysiq)
	2934	CALL VTb(VTinca)
	2935	calday = REAL(days_elapsed + 1) + jH_cur
	2936
	2937	call chemtime(itap+itau_phy-1, date0, dtime)
	2938	IF (config_inca == 'aero') THEN
	2939	CALL AEROSOL_METEO_CALC(
	2940	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
	2941	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
	2942	END IF
	2943
	2944	zxsnow_dummy(:) = 0.0
	2945
	2946	CALL chemhook_begin (calday,
	2947	$ days_elapsed+1,
	2948	$ jH_cur,
	2949	$ pctsrf(1,1),
	2950	$ rlat,
	2951	$ rlon,
	2952	$ airephy,
	2953	$ paprs,
	2954	$ pplay,
	2955	$ coefh,
	2956	$ pphi,
	2957	$ t_seri,
	2958	$ u,
	2959	$ v,
	2960	$ wo(:, :, 1),
	2961	$ q_seri,
	2962	$ zxtsol,
	2963	$ zxsnow_dummy,
	2964	$ solsw,
	2965	$ albsol1,
	2966	$ rain_fall,
	2967	$ snow_fall,
	2968	$ itop_con,
	2969	$ ibas_con,
	2970	$ cldfra,
	2971	$ iim,
	2972	$ jjm,
	2973	$ tr_seri,
	2974	$ ftsol,
	2975	$ paprs,
	2976	$ cdragh,
	2977	$ cdragm,
	2978	$ pctsrf,
	2979	$ pdtphys,
	2980	$ itap)
	2981
	2982	CALL VTe(VTinca)
	2983	CALL VTb(VTphysiq)
	2984	#endif
	2985	END IF !config_inca /= 'none'
	2986	c
	2987	c Calculer les parametres optiques des nuages et quelques
	2988	c parametres pour diagnostiques:
	2989	c
	2990
	2991	IF (aerosol_couple) THEN
	2992	mass_solu_aero(:,:) = ccm(:,:,1)
	2993	mass_solu_aero_pi(:,:) = ccm(:,:,2)
	2994	END IF
	2995
	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,
	3002	e mass_solu_aero, mass_solu_aero_pi,
	3003	e bl95_b0, bl95_b1,
	3004	s cldtaupi, re, fl, ref_liq, ref_ice)
	3005	else
	3006	CALL nuage (paprs, pplay,
	3007	. t_seri, cldliq, cldfra, cldtau, cldemi,
	3008	. cldh, cldl, cldm, cldt, cldq,
	3009	e ok_aie,
	3010	e mass_solu_aero, mass_solu_aero_pi,
	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
	3019
	3020	DO i = 1, klon
	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)
	3029	ENDDO
	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)
	3035	call writefield_phy('q_seri',q_seri,llm)
	3036	endif
	3037
	3038	IF (aerosol_couple) THEN
	3039	#ifdef INCA
	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
	3060	#endif
	3061	ELSE
	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
	3087	itaprad = 0
	3088	ENDIF ! MOD(itaprad,radpas)
	3089	itaprad = itaprad + 1
	3090
	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
	3101
	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)
	3108	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
	3109	ENDDO
	3110	ENDDO
	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)
	3116	call writefield_phy('q_seri',q_seri,llm)
	3117	endif
	3118
	3119	cIM
	3120	IF (ip_ebil_phy.ge.2) THEN
	3121	ztit='after rad'
	3122	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	3125	call diagphy(airephy,ztit,ip_ebil_phy
	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
	3138
	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
	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
	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,
	3181	s zulow, zvlow, zustrdr, zvstrdr,
	3182	s d_t_oro, d_u_oro, d_v_oro)
	3183	ENDIF
	3184	c
	3185	c ajout des tendances
	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	!-----------------------------------------------------------------------------------------
	3190	c
	3191	ENDIF ! fin de test sur ok_orodr
	3192	c
	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)
	3197	call writefield_phy('q_seri',q_seri,llm)
	3198	endif
	3199
	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
	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,
	3225	e rlat,zmea,zstd,zpic,
	3226	e itest,
	3227	e t_seri, u_seri, v_seri,
	3228	s zulow, zvlow, zustrli, zvstrli,
	3229	s d_t_lif, d_u_lif, d_v_lif)
	3230	ENDIF
	3231	c
	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	!-----------------------------------------------------------------------------------------
	3236	c
	3237	ENDIF ! fin de test sur ok_orolf
	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)
	3246	c
	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
	3254	cIM cf. FLott BEG
	3255	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	3256
	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)
	3261	call writefield_phy('q_seri',q_seri,llm)
	3262	endif
	3263
	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
	3279	IF (is_sequential) THEN
	3280
	3281	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
	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
	3289	cIM cf. FLott END
	3290	cIM
	3291	IF (ip_ebil_phy.ge.2) THEN
	3292	ztit='after orography'
	3293	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
	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)
	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 )
	3301	END IF
	3302	c
	3303	c
	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,
	3319	$ ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
	3320	$ ecrit_mth,ecrit_day,ecrit_hf,
	3321	$ klon,klev,rlon,rlat,presnivs,overlap,
	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)
	3330
	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	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
	3349
	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)
	3367
	3368	IF (offline) THEN
	3369
	3370	IF (prt_level.ge.9)
	3371	$ print*,'Attention on met a 0 les thermiques pour phystoke'
	3372	call phystokenc (
	3373	I nlon,klev,pdtphys,rlon,rlat,
	3374	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	3375	I fm_therm,entr_therm,
	3376	I cdragh,coefh,u1,v1,ftsol,pctsrf,
	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
	3390	cIM global posePB BEG
	3391	IF(1.EQ.0) THEN
	3392	c
	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)
	3396	c
	3397	ENDIF !(1.EQ.0) THEN
	3398	cIM global posePB END
	3399	c Accumuler les variables a stocker dans les fichiers histoire:
	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)
	3407	ENDDO
	3408	ENDDO
	3409
	3410	DO k = 1, klev
	3411	DO i = 1, klon
	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
	3417	cIM
	3418	IF (ip_ebil_phy.ge.1) THEN
	3419	ztit='after physic'
	3420	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
	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.
	3427
	3428	call diagphy(airephy,ztit,ip_ebil_phy
	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
	3442	cIM Interpolation sur les niveaux de pression du NMC
	3443	c -------------------------------------------------
	3444	c
	3445	#include "calcul_STDlev.h"
	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)
	3453
	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)
	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(:,:)
	3477	c
	3478	cIM initialisation 5eme fichier de sortie
	3479	cIM ajoute 5eme niveau 170310 BEG
	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
	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
	3529	cIM initialisation + calculs divers diag AMIP2
	3530	c
	3531	#include "calcul_divers.h"
	3532	c
	3533	IF (config_inca /= 'none') THEN
	3534	#ifdef INCA
	3535	CALL VTe(VTphysiq)
	3536	CALL VTb(VTinca)
	3537
	3538	CALL chemhook_end (
	3539	$ dtime,
	3540	$ pplay,
	3541	$ t_seri,
	3542	$ tr_seri,
	3543	$ nbtr,
	3544	$ paprs,
	3545	$ q_seri,
	3546	$ airephy,
	3547	$ pphi,
	3548	$ pphis,
	3549	$ zx_rh)
	3550
	3551	CALL VTe(VTinca)
	3552	CALL VTb(VTphysiq)
	3553	#endif
	3554	END IF
	3555
	3556	c=============================================================
	3557	c
	3558	c Convertir les incrementations en tendances
	3559	c
	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)
	3564	call writefield_phy('q_seri',q_seri,llm)
	3565	endif
	3566
	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
	3577	IF (nqtot.GE.3) THEN
	3578	DO iq = 3, nqtot
	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
	3587	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	3588	cIM global posePB#include "write_bilKP_ins.h"
	3589	cIM global posePB#include "write_bilKP_ave.h"
	3590	c
	3591
	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
	3596	u_ancien(i,k) = u_seri(i,k)
	3597	v_ancien(i,k) = v_seri(i,k)
	3598	t_ancien(i,k) = t_seri(i,k)
	3599	q_ancien(i,k) = q_seri(i,k)
	3600	ENDDO
	3601	ENDDO
	3602	c
	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
	3610	if (prt_level.ge.1) then
	3611	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	3612	write(lunout,*)
	3613	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
	3614	write(lunout,*)
	3615	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
	3616	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
	3617	s pctsrf(igout,is_sic)
	3618	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	3619	do k=1,klev
	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'
	3625	do k=1,klev
	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'
	3630	do k=1,klev
	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 '
	3637	do k=1,klev
	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
	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)
	3654	ENDDO
	3655	ENDDO
	3656	c
	3657
	3658	c 22.03.04 BEG
	3659	c=============================================================
	3660	c Ecriture des sorties
	3661	c=============================================================
	3662	#ifdef CPP_IOIPSL
	3663
	3664	c Recupere des varibles calcule dans differents modules
	3665	c pour ecriture dans histxxx.nc
	3666
	3667	! Get some variables from module fonte_neige_mod
	3668	CALL fonte_neige_get_vars(pctsrf,
	3669	. zxfqcalving, zxfqfonte, zxffonte)
	3670
	3671
	3672	#include "phys_output_write.h"
	3673
	3674	#ifdef histISCCP
	3675	#include "write_histISCCP.h"
	3676	#endif
	3677
	3678	#ifdef histNMC
	3679	#include "write_histhfNMC.h"
	3680	#include "write_histdayNMC.h"
	3681	#include "write_histmthNMC.h"
	3682	#endif
	3683
	3684	#include "write_histday_seri.h"
	3685
	3686	#include "write_paramLMDZ_phy.h"
	3687
	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
	3696
	3697
	3698	IF (lafin) THEN
	3699	itau_phy = itau_phy + itap
	3700	CALL phyredem ("restartphy.nc")
	3701	! open(97,form="unformatted",file="finbin")
	3702	! write(97) u_seri,v_seri,t_seri,q_seri
	3703	! close(97)
	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
	3712	ENDIF
	3713
	3714	! first=.false.
	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