Context Navigation

physiq.F @ 1425

Last change on this file since 1425 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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

File size: 122.9 KB

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

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

Download in other formats:

Original Format