Context Navigation

physiq.F @ 434

Last change on this file since 434 was 434, checked in by lmdzadmin, 22 years ago
Suppression du tableau veget() qui ne sert plus à rien ici LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 63.8 KB

Rev	Line
[179]	1	c
	2	c $Header$
	3	c
[411]	4	SUBROUTINE physiq (nlon,nlev,nqmax,
[2]	5	. debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
	6	. paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
	7	. u,v,t,qx,
[177]	8	. omega, cufi, cvfi,
[2]	9	. d_u, d_v, d_t, d_qx, d_ps)
	10	USE ioipsl
[177]	11	USE histcom
[271]	12	USE writephys
[177]	13
[2]	14	IMPLICIT none
	15	c======================================================================
	16	c
	17	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	18	c
	19	c Objet: Moniteur general de la physique du modele
	20	cAA Modifications quant aux traceurs :
	21	cAA - uniformisation des parametrisations ds phytrac
	22	cAA - stockage des moyennes des champs necessaires
	23	cAA en mode traceur off-line
	24	c======================================================================
	25	c modif ( P. Le Van , 12/10/98 )
	26	c
	27	c Arguments:
	28	c
	29	c nlon----input-I-nombre de points horizontaux
	30	c nlev----input-I-nombre de couches verticales
	31	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
	32	c debut---input-L-variable logique indiquant le premier passage
	33	c lafin---input-L-variable logique indiquant le dernier passage
	34	c rjour---input-R-numero du jour de l'experience
	35	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
	36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
	37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
	38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	39	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	40	c pphis---input-R-geopotentiel du sol
	41	c paire---input-R-aire de chaque maille
	42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
	44	c v-------input-R-vitesse Y (de S a N) en m/s
	45	c t-------input-R-temperature (K)
	46	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
[46]	48	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
[2]	49	c omega---input-R-vitesse verticale en Pa/s
[171]	50	c cufi----input-R-resolution des mailles en x (m)
	51	c cvfi----input-R-resolution des mailles en y (m)
[2]	52	c
	53	c d_u-----output-R-tendance physique de "u" (m/s/s)
	54	c d_v-----output-R-tendance physique de "v" (m/s/s)
	55	c d_t-----output-R-tendance physique de "t" (K/s)
	56	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	57	c d_ps----output-R-tendance physique de la pression au sol
	58	c======================================================================
	59	#include "dimensions.h"
[80]	60	integer jjmp1
	61	parameter (jjmp1=jjm+1-1/jjm)
[2]	62	#include "dimphy.h"
	63	#include "regdim.h"
	64	#include "indicesol.h"
	65	#include "dimsoil.h"
	66	#include "clesphys.h"
	67	#include "control.h"
[53]	68	#include "temps.h"
[2]	69	c======================================================================
[411]	70	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	71	PARAMETER (ok_cvl=.TRUE.)
	72	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	73	PARAMETER (ok_gust=.FALSE.)
	74	c======================================================================
[2]	75	LOGICAL check ! Verifier la conservation du modele en eau
	76	PARAMETER (check=.FALSE.)
[46]	77	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	78	PARAMETER (ok_stratus=.FALSE.)
[2]	79	c======================================================================
	80	c Parametres lies au coupleur OASIS:
	81	#include "oasis.h"
[179]	82	INTEGER,SAVE :: npas, nexca
[2]	83	logical rnpb
	84	parameter(rnpb=.true.)
[112]	85	c PARAMETER (npas=1440)
	86	c PARAMETER (nexca=48)
[46]	87	EXTERNAL fromcpl, intocpl, inicma
[132]	88	c ocean = type de modele ocean a utiliser: force, slab, couple
[223]	89	character*6 ocean
	90	SAVE ocean
	91
	92	c parameter (ocean = 'force ')
[178]	93	c parameter (ocean = 'couple')
[158]	94	logical ok_ocean
[2]	95	c======================================================================
	96	c Clef controlant l'activation du cycle diurne:
	97	ccc LOGICAL cycle_diurne
	98	ccc PARAMETER (cycle_diurne=.FALSE.)
	99	c======================================================================
	100	c Modele thermique du sol, a activer pour le cycle diurne:
	101	ccc LOGICAL soil_model
	102	ccc PARAMETER (soil_model=.FALSE.)
[98]	103	logical ok_veget
[223]	104	save ok_veget
	105	c parameter (ok_veget = .true.)
[205]	106	c parameter (ok_veget = .false.)
[2]	107	c======================================================================
	108	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	109	c le calcul du rayonnement est celle apres la precipitation des nuages.
	110	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	111	c la condensation et la precipitation. Cette cle augmente les impacts
	112	c radiatifs des nuages.
	113	ccc LOGICAL new_oliq
	114	ccc PARAMETER (new_oliq=.FALSE.)
	115	c======================================================================
	116	c Clefs controlant deux parametrisations de l'orographie:
	117	cc LOGICAL ok_orodr
	118	ccc PARAMETER (ok_orodr=.FALSE.)
	119	ccc LOGICAL ok_orolf
	120	ccc PARAMETER (ok_orolf=.FALSE.)
	121	c======================================================================
	122	LOGICAL ok_journe ! sortir le fichier journalier
[223]	123	save ok_journe
	124	c PARAMETER (ok_journe=.true.)
[433]	125
	126	integer lev_histday
	127	save lev_histday
	128	data lev_histday/1/
[2]	129	c
	130	LOGICAL ok_mensuel ! sortir le fichier mensuel
[223]	131	save ok_mensuel
	132	c PARAMETER (ok_mensuel=.true.)
[2]	133	c
	134	LOGICAL ok_instan ! sortir le fichier instantane
[223]	135	save ok_instan
	136	c PARAMETER (ok_instan=.true.)
[2]	137	c
	138	LOGICAL ok_region ! sortir le fichier regional
	139	PARAMETER (ok_region=.FALSE.)
	140	c======================================================================
	141	c
	142	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	143	PARAMETER (ivap=1)
	144	INTEGER iliq ! indice de traceurs pour eau liquide
	145	PARAMETER (iliq=2)
[98]	146
[2]	147	c
[98]	148	c
[2]	149	c Variables argument:
	150	c
	151	INTEGER nlon
	152	INTEGER nlev
	153	INTEGER nqmax
	154	REAL rjourvrai, rjour_ecri
	155	REAL gmtime
	156	REAL pdtphys
	157	LOGICAL debut, lafin
	158	REAL paprs(klon,klev+1)
	159	REAL pplay(klon,klev)
	160	REAL pphi(klon,klev)
	161	REAL pphis(klon)
	162	REAL paire(klon)
	163	REAL presnivs(klev)
	164	REAL znivsig(klev)
[171]	165	REAL zsurf(nbsrf)
	166	real cufi(klon), cvfi(klon)
[2]	167
	168	REAL u(klon,klev)
	169	REAL v(klon,klev)
	170	REAL t(klon,klev)
	171	REAL qx(klon,klev,nqmax)
	172
[46]	173	REAL t_ancien(klon,klev), q_ancien(klon,klev)
	174	SAVE t_ancien, q_ancien
	175	LOGICAL ancien_ok
	176	SAVE ancien_ok
	177
[2]	178	REAL d_t_dyn(klon,klev)
[46]	179	REAL d_q_dyn(klon,klev)
[2]	180
	181	REAL omega(klon,klev)
	182
	183	REAL d_u(klon,klev)
	184	REAL d_v(klon,klev)
	185	REAL d_t(klon,klev)
	186	REAL d_qx(klon,klev,nqmax)
	187	REAL d_ps(klon)
	188
[411]	189	cccIM
	190	INTEGER klevp1
	191	PARAMETER(klevp1=klev+1)
	192	#include "raddim.h"
	193	REAL*8 ZFSUP(KDLON,KFLEV+1)
	194	REAL*8 ZFSDN(KDLON,KFLEV+1)
	195	REAL*8 ZFSUP0(KDLON,KFLEV+1)
	196	REAL*8 ZFSDN0(KDLON,KFLEV+1)
	197
	198	cccIM cf. FH
	199	real u850(klon),v850(klon),u200(klon),v200(klon)
[433]	200	real u500(klon),v500(klon),phi500(klon),w500(klon)
[411]	201
	202	logical ok_hf
	203	real ecrit_hf
	204	integer nid_hf
	205	save ok_hf, ecrit_hf, nid_hf
	206
	207	c QUESTION : noms de variables ?
	208
	209	#define histhf
	210	#ifdef histhf
	211	data ok_hf,ecrit_hf/.true.,0.25/
	212	#else
	213	data ok_hf/.false./
	214	#endif
	215
[2]	216	INTEGER longcles
	217	PARAMETER ( longcles = 20 )
	218	REAL clesphy0( longcles )
	219	c
	220	c Variables quasi-arguments
	221	c
	222	REAL xjour
	223	SAVE xjour
	224	c
	225	c
	226	c Variables propres a la physique
	227	c
	228	REAL dtime
	229	SAVE dtime ! pas temporel de la physique
	230	c
	231	INTEGER radpas
	232	SAVE radpas ! frequence d'appel rayonnement
	233	c
	234	REAL radsol(klon)
	235	SAVE radsol ! bilan radiatif au sol
	236	c
	237	REAL rlat(klon)
	238	SAVE rlat ! latitude pour chaque point
	239	c
	240	REAL rlon(klon)
	241	SAVE rlon ! longitude pour chaque point
	242	c
	243	cc INTEGER iflag_con
	244	cc SAVE iflag_con ! indicateur de la convection
	245	c
	246	INTEGER itap
	247	SAVE itap ! compteur pour la physique
	248	c
[433]	249	REAL co2_ppm_etat0
[2]	250	c
[433]	251	REAL solaire_etat0
[2]	252	c
[433]	253	real slp(klon) ! sea level pressure
	254
[2]	255	REAL ftsol(klon,nbsrf)
	256	SAVE ftsol ! temperature du sol
	257	c
	258	REAL ftsoil(klon,nsoilmx,nbsrf)
	259	SAVE ftsoil ! temperature dans le sol
	260	c
[98]	261	REAL fevap(klon,nbsrf)
	262	SAVE fevap ! evaporation
[177]	263	REAL fluxlat(klon,nbsrf)
	264	SAVE fluxlat
[98]	265	c
[2]	266	REAL deltat(klon)
	267	SAVE deltat ! ecart avec la SST de reference
	268	c
	269	REAL fqsol(klon,nbsrf)
	270	SAVE fqsol ! humidite du sol
	271	c
	272	REAL fsnow(klon,nbsrf)
	273	SAVE fsnow ! epaisseur neigeuse
	274	c
[98]	275	REAL falbe(klon,nbsrf)
	276	SAVE falbe ! albedo par type de surface
[282]	277	REAL falblw(klon,nbsrf)
	278	SAVE falblw ! albedo par type de surface
	279
[98]	280	c
[2]	281	c
	282	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	283	c
	284	REAL zmea(klon)
	285	SAVE zmea ! orographie moyenne
	286	c
	287	REAL zstd(klon)
	288	SAVE zstd ! deviation standard de l'OESM
	289	c
	290	REAL zsig(klon)
	291	SAVE zsig ! pente de l'OESM
	292	c
	293	REAL zgam(klon)
	294	save zgam ! anisotropie de l'OESM
	295	c
	296	REAL zthe(klon)
	297	SAVE zthe ! orientation de l'OESM
	298	c
	299	REAL zpic(klon)
	300	SAVE zpic ! Maximum de l'OESM
	301	c
	302	REAL zval(klon)
	303	SAVE zval ! Minimum de l'OESM
	304	c
	305	REAL rugoro(klon)
	306	SAVE rugoro ! longueur de rugosite de l'OESM
	307	c
	308	REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	309	c
	310	REAL zuthe(klon),zvthe(klon)
	311	SAVE zuthe
	312	SAVE zvthe
[158]	313	INTEGER igwd,idx(klon),itest(klon)
[2]	314	c
[258]	315	REAL agesno(klon,nbsrf)
[2]	316	SAVE agesno ! age de la neige
	317	c
[230]	318	REAL alb_neig(klon)
	319	SAVE alb_neig ! albedo de la neige
	320	cKE43
	321	c Variables liees a la convection de K. Emanuel (sb):
[2]	322	c
[230]	323	REAL ema_workcbmf(klon) ! cloud base mass flux
	324	SAVE ema_workcbmf
	325
	326	REAL ema_cbmf(klon) ! cloud base mass flux
	327	SAVE ema_cbmf
	328
	329	REAL ema_pcb(klon) ! cloud base pressure
	330	SAVE ema_pcb
	331
	332	REAL ema_pct(klon) ! cloud top pressure
	333	SAVE ema_pct
	334
	335	REAL bas, top ! cloud base and top levels
	336	SAVE bas
	337	SAVE top
	338
	339	REAL Ma(klon,klev) ! undilute upward mass flux
	340	SAVE Ma
[411]	341	REAL qcondc(klon,klev) ! in-cld water content from convect
	342	SAVE qcondc
[230]	343	REAL ema_work1(klon, klev), ema_work2(klon, klev)
	344	SAVE ema_work1, ema_work2
	345	REAL wdn(klon), tdn(klon), qdn(klon)
[411]	346
	347	REAL wd(klon) ! sb
	348	SAVE wd ! sb
	349
[230]	350	c Variables locales pour la couche limite (al1):
	351	c
	352	cAl1 REAL pblh(klon) ! Hauteur de couche limite
	353	cAl1 SAVE pblh
	354	c34EK
	355	c
[2]	356	c Variables locales:
	357	c
	358	REAL cdragh(klon) ! drag coefficient pour T and Q
	359	REAL cdragm(klon) ! drag coefficient pour vent
	360	cAA
	361	cAA Pour phytrac
	362	cAA
	363	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	364	REAL yu1(klon) ! vents dans la premiere couche U
	365	REAL yv1(klon) ! vents dans la premiere couche V
	366	LOGICAL offline ! Controle du stockage ds "physique"
[80]	367	PARAMETER (offline=.false.)
[53]	368	INTEGER physid
[2]	369	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
	370	save pfrac_impa
	371	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
	372	save pfrac_nucl
	373	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
	374	save pfrac_1nucl
	375	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	376	REAL frac_nucl(klon,klev) ! idem (nucleation)
	377	cAA
	378	REAL rain_fall(klon) ! pluie
	379	REAL snow_fall(klon) ! neige
[158]	380	save snow_fall, rain_fall
[2]	381	REAL evap(klon), devap(klon) ! evaporation et sa derivee
	382	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
[258]	383	REAL dlw(klon) ! derivee infra rouge
[2]	384	REAL bils(klon) ! bilan de chaleur au sol
[433]	385	cIM cf. JLD
	386	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
	387	C type de sous-surface et pondere par la fraction
[2]	388	REAL fder(klon) ! Derive de flux (sensible et latente)
[158]	389	save fder
[2]	390	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	391	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	392	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	393	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	394	c
	395	REAL frugs(klon,nbsrf) ! longueur de rugosite
[158]	396	save frugs
[2]	397	REAL zxrugs(klon) ! longueur de rugosite
	398	c
	399	c Conditions aux limites
	400	c
	401	INTEGER julien
	402	c
	403	INTEGER lmt_pas
	404	SAVE lmt_pas ! frequence de mise a jour
	405	REAL pctsrf(klon,nbsrf)
	406	SAVE pctsrf ! sous-fraction du sol
	407	REAL albsol(klon)
	408	SAVE albsol ! albedo du sol total
[282]	409	REAL albsollw(klon)
	410	SAVE albsollw ! albedo du sol total
[295]	411	REAL albsol1(klon)
	412	SAVE albsol1 ! albedo du sol total
	413	REAL albsollw1(klon)
	414	SAVE albsollw1 ! albedo du sol total
[282]	415
[2]	416	REAL wo(klon,klev)
	417	SAVE wo ! ozone
	418	c======================================================================
	419	c
	420	c Declaration des procedures appelees
	421	c
	422	EXTERNAL angle ! calculer angle zenithal du soleil
	423	EXTERNAL alboc ! calculer l'albedo sur ocean
	424	EXTERNAL albsno ! calculer albedo sur neige
	425	EXTERNAL ajsec ! ajustement sec
	426	EXTERNAL clmain ! couche limite
	427	EXTERNAL condsurf ! lire les conditions aux limites
	428	EXTERNAL conlmd ! convection (schema LMD)
[230]	429	cKE43
[373]	430	EXTERNAL conema3 ! convect4.3
[2]	431	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	432	cAA
	433	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	434	c ! stockage des coefficients necessaires au
	435	c ! lessivage OFF-LINE et ON-LINE
	436	cAA
	437	EXTERNAL hgardfou ! verifier les temperatures
	438	EXTERNAL nuage ! calculer les proprietes radiatives
	439	EXTERNAL o3cm ! initialiser l'ozone
	440	EXTERNAL orbite ! calculer l'orbite terrestre
	441	EXTERNAL ozonecm ! prescrire l'ozone
	442	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	443	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	444	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	445	EXTERNAL suphec ! initialiser certaines constantes
	446	EXTERNAL transp ! transport total de l'eau et de l'energie
	447	EXTERNAL ecribina ! ecrire le fichier binaire global
	448	EXTERNAL ecribins ! ecrire le fichier binaire global
	449	EXTERNAL ecrirega ! ecrire le fichier binaire regional
	450	EXTERNAL ecriregs ! ecrire le fichier binaire regional
	451	c
	452	c Variables locales
	453	c
[373]	454	real clwcon(klon,klev),rnebcon(klon,klev)
	455	real clwcon0(klon,klev),rnebcon0(klon,klev)
	456	save rnebcon, clwcon
	457
	458	REAL rhcl(klon,klev) ! humiditi relative ciel clair
[2]	459	REAL dialiq(klon,klev) ! eau liquide nuageuse
	460	REAL diafra(klon,klev) ! fraction nuageuse
	461	REAL cldliq(klon,klev) ! eau liquide nuageuse
	462	REAL cldfra(klon,klev) ! fraction nuageuse
	463	REAL cldtau(klon,klev) ! epaisseur optique
	464	REAL cldemi(klon,klev) ! emissivite infrarouge
	465	c
[411]	466	CXXX PB
[98]	467	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	468	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
	469	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
	470	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
[2]	471	c
[98]	472	REAL zxfluxt(klon, klev)
	473	REAL zxfluxq(klon, klev)
	474	REAL zxfluxu(klon, klev)
	475	REAL zxfluxv(klon, klev)
[411]	476	CXXX
[2]	477	REAL heat(klon,klev) ! chauffage solaire
	478	REAL heat0(klon,klev) ! chauffage solaire ciel clair
	479	REAL cool(klon,klev) ! refroidissement infrarouge
	480	REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair
	481	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
[177]	482	real sollwdown(klon) ! downward LW flux at surface
[2]	483	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
	484	REAL albpla(klon)
[433]	485	cIM cf. JLD
	486	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
	487	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
[2]	488	c Le rayonnement n'est pas calcule tous les pas, il faut donc
	489	c sauvegarder les sorties du rayonnement
[177]	490	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
[2]	491	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
[411]	492	cccIM
	493	SAVE ZFSUP,ZFSDN,ZFSUP0,ZFSDN0
	494
[2]	495	INTEGER itaprad
	496	SAVE itaprad
	497	c
	498	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	499	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
	500	c
	501	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
	502	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
	503	c
[391]	504	REAL zxtsol(klon), zxqsol(klon), zxsnow(klon), zxfluxlat(klon)
[2]	505	c
	506	REAL dist, rmu0(klon), fract(klon)
	507	REAL zdtime, zlongi
	508	c
	509	CHARACTER*2 str2
[235]	510	CHARACTER*2 iqn
[2]	511	c
[46]	512	REAL qcheck
	513	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	514	LOGICAL zx_ajustq
	515	c
[2]	516	REAL za, zb
[373]	517	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
	518	real zqsat(klon,klev)
	519	INTEGER i, k, iq, ig, j, nsrf, ll
[2]	520	REAL t_coup
	521	PARAMETER (t_coup=234.0)
	522	c
	523	REAL zphi(klon,klev)
[230]	524	REAL zx_tmp_x(iim), zx_tmp_yjjmp1
	525	REAL zx_relief(iim,jjmp1)
	526	REAL zx_aire(iim,jjmp1)
	527	cKE43
	528	c Variables locales pour la convection de K. Emanuel (sb):
[2]	529	c
[230]	530	REAL upwd(klon,klev) ! saturated updraft mass flux
	531	REAL dnwd(klon,klev) ! saturated downdraft mass flux
	532	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
	533	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	534	REAL cape(klon) ! CAPE
	535	SAVE cape
[433]	536	cccIM
	537	CHARACTER*40 capemaxcels
	538
[230]	539	REAL pbase(klon) ! cloud base pressure
	540	SAVE pbase
	541	REAL bbase(klon) ! cloud base buoyancy
	542	SAVE bbase
	543	REAL rflag(klon) ! flag fonctionnement de convect
[263]	544	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
[230]	545	c -- convect43:
	546	INTEGER ntra ! nb traceurs pour convect4.3
	547	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
	548	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
	549	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	550	REAL dplcldt(klon), dplcldr(klon)
	551	c? . condm_con(klon,klev),conda_con(klon,klev),
	552	c? . mr_con(klon,klev),ep_con(klon,klev)
	553	c? . ,sadiab(klon,klev),wadiab(klon,klev)
	554	c --
	555	c34EK
	556	c
[2]	557	c Variables du changement
	558	c
	559	c con: convection
	560	c lsc: condensation a grande echelle (Large-Scale-Condensation)
	561	c ajs: ajustement sec
	562	c eva: evaporation de l'eau liquide nuageuse
	563	c vdf: couche limite (Vertical DiFfusion)
	564	REAL d_t_con(klon,klev),d_q_con(klon,klev)
	565	REAL d_u_con(klon,klev),d_v_con(klon,klev)
	566	REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
[46]	567	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
[2]	568	REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
	569	REAL rneb(klon,klev)
	570	c
	571	REAL pmfu(klon,klev), pmfd(klon,klev)
	572	REAL pen_u(klon,klev), pen_d(klon,klev)
	573	REAL pde_u(klon,klev), pde_d(klon,klev)
	574	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	575	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
[23]	576	REAL prfl(klon,klev+1), psfl(klon,klev+1)
[2]	577	c
	578	INTEGER ibas_con(klon), itop_con(klon)
	579	REAL rain_con(klon), rain_lsc(klon)
	580	REAL snow_con(klon), snow_lsc(klon)
	581	REAL d_ts(klon,nbsrf)
	582	c
	583	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	584	REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
	585	c
	586	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	587	REAL d_t_oro(klon,klev)
	588	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	589	REAL d_t_lif(klon,klev)
[230]	590
[373]	591	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
	592	real ratqsbas,ratqshaut
	593	save ratqsbas,ratqshaut, ratqs
[287]	594	real zpt_conv(klon,klev)
[230]	595
[373]	596	c Parametres lies au nouveau schema de nuages (SB, PDF)
	597	real fact_cldcon
	598	real facttemps
	599	logical ok_newmicro
	600	save ok_newmicro
	601	save fact_cldcon,facttemps
[385]	602	real facteur
[373]	603
	604	integer iflag_cldcon
	605	save iflag_cldcon
	606
	607	logical ptconv(klon,klev)
	608
[2]	609	c
	610	c Variables liees a l'ecriture de la bande histoire physique
	611	c
	612	INTEGER ecrit_mth
	613	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	614	c
	615	INTEGER ecrit_day
	616	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	617	c
	618	INTEGER ecrit_ins
	619	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	620	c
	621	INTEGER ecrit_reg
	622	SAVE ecrit_reg ! frequence d'ecriture
	623	c
[353]	624	integer itau_w ! pas de temps ecriture = itap + itau_phy
[2]	625	c
	626	c
	627	c Variables locales pour effectuer les appels en serie
	628	c
	629	REAL t_seri(klon,klev), q_seri(klon,klev)
[385]	630	REAL ql_seri(klon,klev),qs_seri(klon,klev)
[2]	631	REAL u_seri(klon,klev), v_seri(klon,klev)
	632	c
	633	REAL tr_seri(klon,klev,nbtr)
[235]	634	REAL d_tr(klon,klev,nbtr)
[2]	635
	636	REAL zx_rh(klon,klev)
	637
	638	INTEGER length
	639	PARAMETER ( length = 100 )
	640	REAL tabcntr0( length )
	641	c
[80]	642	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
[2]	643	REAL zx_tmp_fi2d(klon)
[80]	644	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
	645	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
[2]	646	c
	647	INTEGER nid_day, nid_mth, nid_ins
	648	SAVE nid_day, nid_mth, nid_ins
	649	c
[152]	650	INTEGER nhori, nvert
[295]	651	REAL zsto, zout
	652	real zjulian
	653	save zjulian
[2]	654
	655	character*20 modname
	656	character*80 abort_message
	657	logical ok_sync
[205]	658	real date0
[353]	659	integer idayref
[2]	660
[271]	661	C essai writephys
	662	integer fid_day, fid_mth, fid_ins
	663	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	664	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	665	parameter (prof2d_on = 1, prof3d_on = 2,
	666	. prof2d_av = 3, prof3d_av = 4)
	667	character*30 nom_fichier
	668	character*10 varname
	669	character*40 vartitle
	670	character*20 varunits
[385]	671	C Variables liees au bilan d'energie et d'enthalpi
	672	REAL ztsol(klon)
	673	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	674	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	675	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	676	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	677	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	678	REAL d_h_vcol_phy
	679	REAL fs_bound, fq_bound
	680	SAVE d_h_vcol_phy
	681	REAL zero_v(klon)
	682	CHARACTER*15 ztit
	683	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	684	SAVE ip_ebil
	685	DATA ip_ebil/2/
[411]	686	INTEGER if_ebil ! level for energy conserv. dignostics
	687	SAVE if_ebil
	688	c+jld ec_conser
	689	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
	690	REAL ZRCPD
	691	c-jld ec_conser
	692	cIM
	693	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
	694	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
	695	REAL zt2m(klon), zq2m(klon)
	696	REAL zu10m(klon), zv10m(klon)
	697	CHARACTER*40 t2mincels, t2maxcels
[2]	698	c
	699	c Declaration des constantes et des fonctions thermodynamiques
	700	c
	701	#include "YOMCST.h"
	702	#include "YOETHF.h"
	703	#include "FCTTRE.h"
	704	c======================================================================
	705	modname = 'physiq'
[385]	706	IF (if_ebil.ge.1) THEN
	707	DO i=1,klon
	708	zero_v(i)=0.
	709	END DO
	710	END IF
[2]	711	ok_sync=.TRUE.
	712	IF (nqmax .LT. 2) THEN
	713	PRINT*, 'eaux vapeur et liquide sont indispensables'
	714	CALL ABORT
	715	ENDIF
	716	IF (debut) THEN
	717	CALL suphec ! initialiser constantes et parametres phys.
	718	ENDIF
[88]	719
	720
[2]	721	c======================================================================
	722	xjour = rjourvrai
	723	c
	724	c Si c'est le debut, il faut initialiser plusieurs choses
	725	c ********
	726	c
	727	IF (debut) THEN
[385]	728	C
	729	IF (if_ebil.ge.1) d_h_vcol_phy=0.
[223]	730	c
	731	c appel a la lecture du run.def physique
	732	c
	733	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
[373]	734	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
[395]	735	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil)
[433]	736	cIM . , RI0)
[223]	737
[2]	738	c
[98]	739	c
[2]	740	c Initialiser les compteurs:
	741	c
	742
[158]	743	frugs = 0.
[2]	744	itap = 0
	745	itaprad = 0
	746	c
[433]	747	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
[98]	748	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
[177]	749	. falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
[258]	750	. dlw,radsol,frugs,agesno,clesphy0,
[46]	751	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
[373]	752	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
[2]	753
	754	c
	755	radpas = NINT( 86400./dtime/nbapp_rad)
	756
	757	c
	758	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
	759	, ok_instan, ok_region )
	760	c
	761	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	762	PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
	763	abort_message=' See above '
	764	call abort_gcm(modname,abort_message,1)
	765	ENDIF
	766	IF (nlon .NE. klon) THEN
	767	PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
	768	abort_message=' See above '
	769	call abort_gcm(modname,abort_message,1)
	770	ENDIF
	771	IF (nlev .NE. klev) THEN
	772	PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
	773	abort_message=' See above '
	774	call abort_gcm(modname,abort_message,1)
	775	ENDIF
	776	c
	777	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
	778	PRINT*, 'Nbre d appels au rayonnement insuffisant'
	779	PRINT*, "Au minimum 4 appels par jour si cycle diurne"
	780	abort_message=' See above '
	781	call abort_gcm(modname,abort_message,1)
	782	ENDIF
	783	PRINT*, "Clef pour la convection, iflag_con=", iflag_con
[411]	784	PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl
[2]	785	c
[230]	786	cKE43
	787	c Initialisation pour la convection de K.E. (sb):
[301]	788	IF (iflag_con.GE.3) THEN
[230]	789
	790	PRINT, "** Convection de Kerry Emanuel 4.3 "
	791	PRINT*, "On va utiliser le melange convectif des traceurs qui"
	792	PRINT*, "est calcule dans convect4.3"
	793	PRINT*, " !!! penser aux logical flags de phytrac"
	794
	795	DO i = 1, klon
	796	ema_cbmf(i) = 0.
	797	ema_pcb(i) = 0.
	798	ema_pct(i) = 0.
	799	ema_workcbmf(i) = 0.
	800	ENDDO
[433]	801
	802	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	803	DO i = 1, klon
	804	ibas_con(i) = 1
	805	itop_con(i) = klev+1
	806	ENDDO
	807	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	808
[230]	809	ENDIF
[433]	810
[230]	811	c34EK
[2]	812	IF (ok_orodr) THEN
	813	DO i=1,klon
	814	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	815	ENDDO
	816	CALL SUGWD(klon,klev,paprs,pplay)
	817	DO i=1,klon
	818	zuthe(i)=0.
	819	zvthe(i)=0.
	820	if(zstd(i).gt.10.)then
	821	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	822	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	823	endif
	824	ENDDO
	825	ENDIF
	826	c
	827	c
	828	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	829	PRINT*,'La frequence de lecture surface est de ', lmt_pas
	830	c
	831	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
	832	IF (ok_mensuel) THEN
	833	PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
	834	ENDIF
	835	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
	836	IF (ok_journe) THEN
	837	PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
	838	ENDIF
	839	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
[80]	840	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
[411]	841	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
[364]	842	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
[2]	843	IF (ok_instan) THEN
	844	PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
	845	ENDIF
	846	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
	847	IF (ok_region) THEN
	848	PRINT*, 'La frequence de sortie region est de ', ecrit_reg
	849	ENDIF
[112]	850
[2]	851	c
[112]	852	c Initialiser le couplage si necessaire
[2]	853	c
[112]	854	npas = 0
	855	nexca = 0
	856	if (ocean == 'couple') then
	857	npas = itaufin/ iphysiq
	858	nexca = 86400 / dtime
	859	write(,)' ##### Ocean couple #####'
	860	write(,)' Valeurs des pas de temps'
	861	write(,)' npas = ', npas
	862	write(,)' nexca = ', nexca
	863	endif
	864	c
	865	c
[411]	866	cccIM
[433]	867	capemaxcels = 't_max(X)'
[411]	868	t2mincels = 't_min(X)'
	869	t2maxcels = 't_max(X)'
[295]	870
[411]	871	cccIM cf. FH
[295]	872	c
[411]	873	c=============================================================
	874	c Initialisation des sorties
	875	c=============================================================
	876	#ifdef histhf
	877	#include "ini_histhf.h"
	878	#endif
[98]	879
[411]	880	#include "ini_histday.h"
	881	#include "ini_histmth.h"
	882	#include "ini_histins.h"
[177]	883
[411]	884	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
[361]	885	date0 = zjulian
	886	C date0 = day_ini
[290]	887	WRITE(,) 'physiq date0 : ',date0
[2]	888	c
	889	c
	890	c
	891	c Prescrire l'ozone dans l'atmosphere
	892	c
	893	c
	894	cc DO i = 1, klon
	895	cc DO k = 1, klev
	896	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	897	cc ENDDO
	898	cc ENDDO
	899	c
	900	c
	901	ENDIF
	902	c
	903	c ************** Fin de IF ( debut ) *************
	904	c
	905	c
	906	c Mettre a zero des variables de sortie (pour securite)
	907	c
	908	DO i = 1, klon
	909	d_ps(i) = 0.0
	910	ENDDO
	911	DO k = 1, klev
	912	DO i = 1, klon
	913	d_t(i,k) = 0.0
	914	d_u(i,k) = 0.0
	915	d_v(i,k) = 0.0
	916	ENDDO
	917	ENDDO
	918	DO iq = 1, nqmax
	919	DO k = 1, klev
	920	DO i = 1, klon
	921	d_qx(i,k,iq) = 0.0
	922	ENDDO
	923	ENDDO
	924	ENDDO
	925	c
	926	c Ne pas affecter les valeurs entrees de u, v, h, et q
	927	c
	928	DO k = 1, klev
	929	DO i = 1, klon
	930	t_seri(i,k) = t(i,k)
	931	u_seri(i,k) = u(i,k)
	932	v_seri(i,k) = v(i,k)
	933	q_seri(i,k) = qx(i,k,ivap)
	934	ql_seri(i,k) = qx(i,k,iliq)
[385]	935	qs_seri(i,k) = 0.
[2]	936	ENDDO
	937	ENDDO
	938	IF (nqmax.GE.3) THEN
	939	DO iq = 3, nqmax
	940	DO k = 1, klev
	941	DO i = 1, klon
	942	tr_seri(i,k,iq-2) = qx(i,k,iq)
	943	ENDDO
	944	ENDDO
	945	ENDDO
	946	ELSE
	947	DO k = 1, klev
	948	DO i = 1, klon
	949	tr_seri(i,k,1) = 0.0
	950	ENDDO
	951	ENDDO
	952	ENDIF
[385]	953	C
	954	IF (if_ebil.ge.1) THEN
	955	DO i = 1, klon
	956	ztsol(i) = 0.
	957	ENDDO
[391]	958	DO nsrf = 1, nbsrf
	959	DO i = 1, klon
	960	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	961	ENDDO
[385]	962	ENDDO
	963	ztit='after dynamic'
	964	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
	965	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	966	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	967	C Comme les tendances de la physique sont ajoute dans la dynamique,
	968	C on devrait avoir que la variation d'entalpie par la dynamique
	969	C est egale a la variation de la physique au pas de temps precedent.
	970	C Donc la somme de ces 2 variations devrait etre nulle.
	971	call diagphy(paire,ztit,ip_ebil
	972	e , zero_v, zero_v, zero_v, zero_v, zero_v
	973	e , zero_v, zero_v, zero_v, ztsol
	974	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	975	s , fs_bound, fq_bound )
	976	END IF
	977
[46]	978	c Diagnostiquer la tendance dynamique
	979	c
	980	IF (ancien_ok) THEN
	981	DO k = 1, klev
	982	DO i = 1, klon
	983	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	984	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	985	ENDDO
	986	ENDDO
	987	ELSE
	988	DO k = 1, klev
	989	DO i = 1, klon
	990	d_t_dyn(i,k) = 0.0
	991	d_q_dyn(i,k) = 0.0
	992	ENDDO
	993	ENDDO
	994	ancien_ok = .TRUE.
	995	ENDIF
	996	c
[2]	997	c Ajouter le geopotentiel du sol:
	998	c
	999	DO k = 1, klev
	1000	DO i = 1, klon
	1001	zphi(i,k) = pphi(i,k) + pphis(i)
	1002	ENDDO
	1003	ENDDO
	1004	c
	1005	c Verifier les temperatures
	1006	c
	1007	CALL hgardfou(t_seri,ftsol,'debutphy')
	1008	c
	1009	c Incrementer le compteur de la physique
	1010	c
	1011	itap = itap + 1
	1012	julien = MOD(NINT(xjour),360)
	1013	c
	1014	c Mettre en action les conditions aux limites (albedo, sst, etc.).
	1015	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1016	c
	1017	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
[353]	1018	PRINT *,' PHYS cond julien ',julien
[2]	1019	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
	1020	ENDIF
	1021	c
	1022	c Re-evaporer l'eau liquide nuageuse
	1023	c
	1024	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1025	DO i = 1, klon
	1026	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
[373]	1027	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1028	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
[2]	1029	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1030	zb = MAX(0.0,ql_seri(i,k))
	1031	za = - MAX(0.0,ql_seri(i,k))
	1032	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1033	t_seri(i,k) = t_seri(i,k) + za
	1034	q_seri(i,k) = q_seri(i,k) + zb
	1035	ql_seri(i,k) = 0.0
	1036	d_t_eva(i,k) = za
	1037	d_q_eva(i,k) = zb
	1038	ENDDO
	1039	ENDDO
	1040	c
[385]	1041	IF (if_ebil.ge.2) THEN
	1042	ztit='after reevap'
[411]	1043	CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime
[385]	1044	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1045	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1046	call diagphy(paire,ztit,ip_ebil
	1047	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1048	e , zero_v, zero_v, zero_v, ztsol
	1049	e , d_h_vcol, d_qt, d_ec
	1050	s , fs_bound, fq_bound )
	1051	C
	1052	END IF
	1053	C
	1054	c
[2]	1055	c Appeler la diffusion verticale (programme de couche limite)
	1056	c
	1057	DO i = 1, klon
[152]	1058	c if (.not. ok_veget) then
	1059	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
	1060	c endif
	1061	c frugs(i,is_lic) = rugoro(i)
	1062	c frugs(i,is_oce) = rugmer(i)
	1063	c frugs(i,is_sic) = 0.001
[2]	1064	zxrugs(i) = 0.0
	1065	ENDDO
	1066	DO nsrf = 1, nbsrf
	1067	DO i = 1, klon
[395]	1068	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
	1069	cccc frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
[2]	1070	ENDDO
	1071	ENDDO
	1072	DO nsrf = 1, nbsrf
	1073	DO i = 1, klon
[177]	1074	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
[2]	1075	ENDDO
	1076	ENDDO
	1077	c
[109]	1078	C calculs necessaires au calcul de l'albedo dans l'interface
	1079	c
	1080	CALL orbite(FLOAT(julien),zlongi,dist)
	1081	IF (cycle_diurne) THEN
	1082	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
	1083	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
	1084	ELSE
	1085	rmu0 = -999.999
	1086	ENDIF
	1087
[433]	1088	DO nsrf = 1, nbsrf
	1089	DO i = 1, klon
	1090	fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
	1091	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
	1092	ENDDO
	1093	ENDDO
	1094
[258]	1095	fder = dlw
[152]	1096
[364]	1097	CALL clmain(dtime,itap,date0,pctsrf,
[109]	1098	e t_seri,q_seri,u_seri,v_seri,
	1099	e julien, rmu0,
[112]	1100	e ok_veget, ocean, npas, nexca, ftsol,
[177]	1101	$ soil_model,ftsoil,
[282]	1102	$ paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,falblw,
	1103	$ fluxlat,
[433]	1104	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
	1105	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
[171]	1106	e rlon, rlat, cufi, cvfi, frugs,
	1107	e debut, lafin, agesno,rugoro ,
[2]	1108	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
[171]	1109	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
[2]	1110	s dsens, devap,
[411]	1111	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m)
[2]	1112	c
[411]	1113	CXXX PB
	1114	CXXX Incrementation des flux
	1115	CXXX
[98]	1116	zxfluxt=0.
	1117	zxfluxq=0.
	1118	zxfluxu=0.
	1119	zxfluxv=0.
	1120	DO nsrf = 1, nbsrf
	1121	DO k = 1, klev
	1122	DO i = 1, klon
	1123	zxfluxt(i,k) = zxfluxt(i,k) +
	1124	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
	1125	zxfluxq(i,k) = zxfluxq(i,k) +
	1126	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
	1127	zxfluxu(i,k) = zxfluxu(i,k) +
	1128	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
	1129	zxfluxv(i,k) = zxfluxv(i,k) +
	1130	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
	1131	END DO
	1132	END DO
	1133	END DO
[2]	1134	DO i = 1, klon
[98]	1135	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
	1136	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
	1137	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
[258]	1138	fder(i) = dlw(i) + dsens(i) + devap(i)
[2]	1139	ENDDO
[80]	1140
[287]	1141
[2]	1142	DO k = 1, klev
	1143	DO i = 1, klon
	1144	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	1145	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
	1146	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	1147	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	1148	ENDDO
	1149	ENDDO
	1150	c
[385]	1151	IF (if_ebil.ge.2) THEN
	1152	ztit='after clmain'
	1153	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1154	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1155	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1156	call diagphy(paire,ztit,ip_ebil
	1157	e , zero_v, zero_v, zero_v, zero_v, sens
	1158	e , evap , zero_v, zero_v, ztsol
	1159	e , d_h_vcol, d_qt, d_ec
	1160	s , fs_bound, fq_bound )
	1161	END IF
	1162	C
	1163	c
[2]	1164	c Incrementer la temperature du sol
	1165	c
	1166	DO i = 1, klon
	1167	zxtsol(i) = 0.0
[391]	1168	zxfluxlat(i) = 0.0
[411]	1169	cccIM
	1170	zt2m(i) = 0.0
	1171	zq2m(i) = 0.0
	1172	zu10m(i) = 0.0
	1173	zv10m(i) = 0.0
	1174	c
[98]	1175	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
	1176	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
	1177	$ THEN
	1178	WRITE(,) 'physiq : pb sous surface au point ', i,
	1179	$ pctsrf(i, 1 : nbsrf)
	1180	ENDIF
[2]	1181	ENDDO
	1182	DO nsrf = 1, nbsrf
[391]	1183	DO i = 1, klon
[411]	1184	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
[177]	1185	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
[433]	1186	cIM cf. JLD
	1187	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
	1188	$ + fluxt(i,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
[177]	1189	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
[391]	1190	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
[411]	1191	cccIM
	1192	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
	1193	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
	1194	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
	1195	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
	1196	c ENDIF
[391]	1197	ENDDO
[2]	1198	ENDDO
	1199
	1200	c
	1201	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
	1202	c
	1203	DO nsrf = 1, nbsrf
[230]	1204	DO i = 1, klon
	1205	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
[411]	1206	cccIM
	1207	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
	1208	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
	1209	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
	1210	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
[230]	1211	ENDDO
[2]	1212	ENDDO
	1213	c
[411]	1214	c
[2]	1215	c Calculer la derive du flux infrarouge
	1216	c
[411]	1217	cXXX DO nsrf = 1, nbsrf
[2]	1218	DO i = 1, klon
[411]	1219	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
[258]	1220	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
[411]	1221	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
	1222	cXXX . *pctsrf(i,nsrf)
	1223	cXXX ENDIF
	1224	cXXX ENDDO
[2]	1225	ENDDO
	1226	c
	1227	c Appeler la convection (au choix)
	1228	c
	1229	DO k = 1, klev
	1230	DO i = 1, klon
[46]	1231	conv_q(i,k) = d_q_dyn(i,k)
[2]	1232	. + d_q_vdf(i,k)/dtime
	1233	conv_t(i,k) = d_t_dyn(i,k)
	1234	. + d_t_vdf(i,k)/dtime
	1235	ENDDO
	1236	ENDDO
	1237	IF (check) THEN
[46]	1238	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
	1239	PRINT*, "avantcon=", za
[2]	1240	ENDIF
[46]	1241	zx_ajustq = .FALSE.
	1242	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	1243	IF (zx_ajustq) THEN
	1244	DO i = 1, klon
	1245	z_avant(i) = 0.0
	1246	ENDDO
	1247	DO k = 1, klev
	1248	DO i = 1, klon
	1249	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
	1250	. *(paprs(i,k)-paprs(i,k+1))/RG
	1251	ENDDO
	1252	ENDDO
	1253	ENDIF
[2]	1254	IF (iflag_con.EQ.1) THEN
	1255	stop'reactiver le call conlmd dans physiq.F'
	1256	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	1257	c . d_t_con, d_q_con,
	1258	c . rain_con, snow_con, ibas_con, itop_con)
	1259	ELSE IF (iflag_con.EQ.2) THEN
	1260	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
[98]	1261	e conv_t, conv_q, zxfluxq(1,1), omega,
[2]	1262	s d_t_con, d_q_con, rain_con, snow_con,
	1263	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	1264	s kcbot, kctop, kdtop, pmflxr, pmflxs)
[177]	1265	WHERE (rain_con < 0.) rain_con = 0.
	1266	WHERE (snow_con < 0.) snow_con = 0.
[2]	1267	DO i = 1, klon
	1268	ibas_con(i) = klev+1 - kcbot(i)
	1269	itop_con(i) = klev+1 - kctop(i)
	1270	ENDDO
[301]	1271	ELSE IF (iflag_con.GE.3) THEN
[230]	1272	c nb of tracers for the KE convection:
	1273	if (nqmax .GE. 4) then
	1274	ntra = nbtr
	1275	else
	1276	ntra = 1
	1277	endif
[411]	1278	c
	1279	c sb, oct02:
	1280	c Schema de convection modularise et vectorise:
	1281	c (driver commun aux versions 3 et 4)
	1282	c
	1283	IF (ok_cvl) THEN ! new driver for convectL
	1284
	1285	CALL concvl (iflag_con,
	1286	. dtime,paprs,pplay,t_seri,q_seri,
	1287	. u_seri,v_seri,tr_seri,nbtr,
	1288	. ema_work1,ema_work2,
	1289	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	1290	. rain_con, snow_con, ibas_con, itop_con,
	1291	. upwd,dnwd,dnwd0,
	1292	. Ma,cape,tvp,iflagctrl,
	1293	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
[433]	1294	cIM cf. FH
	1295	clwcon0=qcondc
[411]	1296
	1297	ELSE ! ok_cvl
	1298
[373]	1299	c print*,'Avant conema OUI'
	1300	CALL conema3 (dtime,
	1301	. paprs,pplay,t_seri,q_seri,
[301]	1302	. u_seri,v_seri,tr_seri,nbtr,
[230]	1303	. ema_work1,ema_work2,
	1304	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	1305	. rain_con, snow_con, ibas_con, itop_con,
[301]	1306	. upwd,dnwd,dnwd0,bas,top,
	1307	. Ma,cape,tvp,rflag,
[373]	1308	. pbase
	1309	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
	1310	. ,clwcon0)
[395]	1311	print*,'Apres conema3 '
[373]	1312
[433]	1313	ENDIF ! ok_cvl
	1314
[411]	1315	IF (.NOT. ok_gust) THEN
	1316	do i = 1, klon
	1317	wd(i)=0.0
	1318	enddo
	1319	ENDIF
	1320
[433]	1321	c =================================================================== c
	1322	c Calcul des proprietes des nuages convectifs
[373]	1323	c
	1324	DO k = 1, klev
	1325	DO i = 1, klon
	1326	zx_t = t_seri(i,k)
	1327	IF (thermcep) THEN
	1328	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	1329	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	1330	zx_qs = MIN(0.5,zx_qs)
	1331	zcor = 1./(1.-retv*zx_qs)
	1332	zx_qs = zx_qs*zcor
	1333	ELSE
	1334	IF (zx_t.LT.t_coup) THEN
	1335	zx_qs = qsats(zx_t)/pplay(i,k)
	1336	ELSE
	1337	zx_qs = qsatl(zx_t)/pplay(i,k)
	1338	ENDIF
	1339	ENDIF
	1340	zqsat(i,k)=zx_qs
	1341	ENDDO
	1342	ENDDO
	1343
[411]	1344	c calcul des proprietes des nuages convectifs
[373]	1345	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	1346	call clouds_gno
	1347	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	1348
[433]	1349	c =================================================================== c
[411]	1350
[230]	1351	DO i = 1, klon
	1352	ema_pcb(i) = pbase(i)
	1353	ENDDO
	1354	DO i = 1, klon
	1355	ema_pct(i) = paprs(i,itop_con(i))
	1356	ENDDO
	1357	DO i = 1, klon
	1358	ema_cbmf(i) = ema_workcbmf(i)
	1359	ENDDO
[2]	1360	ELSE
[230]	1361	PRINT*, "iflag_con non-prevu", iflag_con
	1362	CALL abort
[2]	1363	ENDIF
[205]	1364
[364]	1365	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	1366	c . d_u_con, d_v_con)
[2]	1367	DO k = 1, klev
[205]	1368	DO i = 1, klon
[2]	1369	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
	1370	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
	1371	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
	1372	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
[205]	1373	ENDDO
[2]	1374	ENDDO
[385]	1375	c
	1376	IF (if_ebil.ge.2) THEN
	1377	ztit='after convect'
	1378	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1379	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1380	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1381	call diagphy(paire,ztit,ip_ebil
	1382	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1383	e , zero_v, rain_con, snow_con, ztsol
	1384	e , d_h_vcol, d_qt, d_ec
	1385	s , fs_bound, fq_bound )
	1386	END IF
	1387	C
[2]	1388	IF (check) THEN
[230]	1389	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
	1390	PRINT*, "aprescon=", za
	1391	zx_t = 0.0
	1392	za = 0.0
	1393	DO i = 1, klon
[46]	1394	za = za + paire(i)/FLOAT(klon)
	1395	zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
[230]	1396	ENDDO
	1397	zx_t = zx_t/za*dtime
	1398	PRINT*, "Precip=", zx_t
[2]	1399	ENDIF
[46]	1400	IF (zx_ajustq) THEN
[230]	1401	DO i = 1, klon
[46]	1402	z_apres(i) = 0.0
[230]	1403	ENDDO
	1404	DO k = 1, klev
	1405	DO i = 1, klon
	1406	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
	1407	. *(paprs(i,k)-paprs(i,k+1))/RG
	1408	ENDDO
	1409	ENDDO
	1410	DO i = 1, klon
	1411	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
	1412	. /z_apres(i)
	1413	ENDDO
	1414	DO k = 1, klev
	1415	DO i = 1, klon
	1416	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
	1417	. z_factor(i).LT.(1.0-1.0E-08)) THEN
	1418	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	1419	ENDIF
	1420	ENDDO
	1421	ENDDO
[46]	1422	ENDIF
	1423	zx_ajustq=.FALSE.
[2]	1424	c
	1425	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
	1426	c
[373]	1427	IF (iflag_con .NE. 2 .AND. debut) THEN
[230]	1428	PRINT*, 'Pour l instant, seul conflx fonctionne ',
	1429	$ 'avec traceurs', iflag_con
	1430	PRINT*,' Mettre iflag_con',
[373]	1431	$ ' = 2 dans run.def et repasser'
	1432	c CALL abort
[230]	1433	ENDIF
[2]	1434	c
	1435	ENDIF !--nqmax.GT.2
	1436	c
	1437	c Appeler l'ajustement sec
	1438	c
[34]	1439	CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
	1440	DO k = 1, klev
	1441	DO i = 1, klon
	1442	t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
	1443	q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
	1444	ENDDO
	1445	ENDDO
[385]	1446	c
	1447	IF (if_ebil.ge.2) THEN
	1448	ztit='after dry_adjust'
	1449	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1450	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1451	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1452	END IF
[230]	1453
[373]	1454
	1455	c-------------------------------------------------------------------------
	1456	c Caclul des ratqs
	1457	c-------------------------------------------------------------------------
	1458
	1459	c print*,'calcul des ratqs'
	1460	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
	1461	c ----------------
	1462	c on ecrase le tableau ratqsc calcule par clouds_gno
	1463	if (iflag_cldcon.eq.1) then
	1464	do k=1,klev
	1465	do i=1,klon
	1466	if(ptconv(i,k)) then
	1467	ratqsc(i,k)=ratqsbas
	1468	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
	1469	else
	1470	ratqsc(i,k)=0.
	1471	endif
	1472	enddo
	1473	enddo
	1474	endif
	1475
	1476	c ratqs stables
	1477	c -------------
	1478	do k=1,klev
	1479	ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
	1480	s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
	1481	enddo
	1482
	1483
	1484	c ratqs final
	1485	c -----------
	1486	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
	1487	c les ratqs sont une conbinaison de ratqss et ratqsc
	1488	c ratqs final
	1489	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
	1490	c relaxation des ratqs
[385]	1491	c facttemps=exp(-pdtphys/1.e4)
	1492	facteur=exp(-pdtphys*facttemps)
	1493	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
[373]	1494	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
	1495	c print*,'calcul des ratqs fini'
[304]	1496	else
[373]	1497	c on ne prend que le ratqs stable pour fisrtilp
	1498	ratqs(:,:)=ratqss(:,:)
[304]	1499	endif
[373]	1500
	1501
[2]	1502	c
	1503	c Appeler le processus de condensation a grande echelle
	1504	c et le processus de precipitation
[373]	1505	c-------------------------------------------------------------------------
	1506	CALL fisrtilp(dtime,paprs,pplay,
	1507	. t_seri, q_seri,ptconv,ratqs,
[2]	1508	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
	1509	. rain_lsc, snow_lsc,
	1510	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
[23]	1511	. frac_impa, frac_nucl,
[373]	1512	. prfl, psfl, rhcl)
	1513
[177]	1514	WHERE (rain_lsc < 0) rain_lsc = 0.
	1515	WHERE (snow_lsc < 0) snow_lsc = 0.
[2]	1516	DO k = 1, klev
	1517	DO i = 1, klon
	1518	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
	1519	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
	1520	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
	1521	cldfra(i,k) = rneb(i,k)
	1522	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	1523	ENDDO
	1524	ENDDO
	1525	IF (check) THEN
[46]	1526	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
	1527	PRINT*, "apresilp=", za
[2]	1528	zx_t = 0.0
[46]	1529	za = 0.0
[2]	1530	DO i = 1, klon
[46]	1531	za = za + paire(i)/FLOAT(klon)
	1532	zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
	1533	ENDDO
	1534	zx_t = zx_t/za*dtime
[2]	1535	PRINT*, "Precip=", zx_t
	1536	ENDIF
	1537	c
[385]	1538	IF (if_ebil.ge.2) THEN
	1539	ztit='after fisrt'
	1540	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1541	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1542	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1543	call diagphy(paire,ztit,ip_ebil
	1544	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1545	e , zero_v, rain_lsc, snow_lsc, ztsol
	1546	e , d_h_vcol, d_qt, d_ec
	1547	s , fs_bound, fq_bound )
	1548	END IF
	1549	c
[373]	1550	c-------------------------------------------------------------------
	1551	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	1552	c-------------------------------------------------------------------
	1553
	1554	c 1. NUAGES CONVECTIFS
[2]	1555	c
[373]	1556	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
	1557
	1558	c Nuages diagnostiques pour Tiedtke
[46]	1559	CALL diagcld1(paprs,pplay,
[2]	1560	. rain_con,snow_con,ibas_con,itop_con,
	1561	. diafra,dialiq)
	1562	DO k = 1, klev
	1563	DO i = 1, klon
	1564	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	1565	cldliq(i,k) = dialiq(i,k)
	1566	cldfra(i,k) = diafra(i,k)
	1567	ENDIF
	1568	ENDDO
	1569	ENDDO
[373]	1570
	1571	ELSE IF (iflag_cldcon.eq.3) THEN
	1572	c On prend pour les nuages convectifs le max du calcul de la
	1573	c convection et du calcul du pas de temps précédent diminué d'un facteur
	1574	c facttemps
[385]	1575	c facttemps=pdtphys/1.e4
	1576	facteur = pdtphys *facttemps
[373]	1577	do k=1,klev
	1578	do i=1,klon
[385]	1579	rnebcon(i,k)=rnebcon(i,k)*facteur
[373]	1580	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
	1581	s then
	1582	rnebcon(i,k)=rnebcon0(i,k)
	1583	clwcon(i,k)=clwcon0(i,k)
	1584	endif
	1585	enddo
	1586	enddo
	1587
	1588	c On prend la somme des fractions nuageuses et des contenus en eau
	1589	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	1590	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	1591
	1592
[2]	1593	ENDIF
	1594	c
[373]	1595	c 2. NUAGES STARTIFORMES
[46]	1596	c
	1597	IF (ok_stratus) THEN
	1598	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	1599	DO k = 1, klev
	1600	DO i = 1, klon
	1601	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	1602	cldliq(i,k) = dialiq(i,k)
	1603	cldfra(i,k) = diafra(i,k)
	1604	ENDIF
	1605	ENDDO
	1606	ENDDO
	1607	ENDIF
	1608	c
[2]	1609	c Precipitation totale
	1610	c
	1611	DO i = 1, klon
	1612	rain_fall(i) = rain_con(i) + rain_lsc(i)
	1613	snow_fall(i) = snow_con(i) + snow_lsc(i)
	1614	ENDDO
	1615	c
[385]	1616	IF (if_ebil.ge.2) THEN
	1617	ztit="after diagcld"
	1618	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1619	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1620	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1621	END IF
	1622	c
[2]	1623	c Calculer l'humidite relative pour diagnostique
	1624	c
	1625	DO k = 1, klev
	1626	DO i = 1, klon
	1627	zx_t = t_seri(i,k)
	1628	IF (thermcep) THEN
	1629	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	1630	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	1631	zx_qs = MIN(0.5,zx_qs)
	1632	zcor = 1./(1.-retv*zx_qs)
	1633	zx_qs = zx_qs*zcor
	1634	ELSE
	1635	IF (zx_t.LT.t_coup) THEN
	1636	zx_qs = qsats(zx_t)/pplay(i,k)
	1637	ELSE
	1638	zx_qs = qsatl(zx_t)/pplay(i,k)
	1639	ENDIF
	1640	ENDIF
	1641	zx_rh(i,k) = q_seri(i,k)/zx_qs
[373]	1642	zqsat(i,k)=zx_qs
[2]	1643	ENDDO
	1644	ENDDO
	1645	c
	1646	c Calculer les parametres optiques des nuages et quelques
	1647	c parametres pour diagnostiques:
	1648	c
[373]	1649	if (ok_newmicro) then
	1650	CALL newmicro (paprs, pplay,ok_newmicro,
	1651	. t_seri, cldliq, cldfra, cldtau, cldemi,
	1652	. cldh, cldl, cldm, cldt, cldq)
	1653	else
[2]	1654	CALL nuage (paprs, pplay,
	1655	. t_seri, cldliq, cldfra, cldtau, cldemi,
	1656	. cldh, cldl, cldm, cldt, cldq)
[373]	1657	endif
[2]	1658	c
	1659	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	1660	c
	1661	IF (MOD(itaprad,radpas).EQ.0) THEN
	1662	DO i = 1, klon
[98]	1663	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
	1664	. + falbe(i,is_lic) * pctsrf(i,is_lic)
	1665	. + falbe(i,is_ter) * pctsrf(i,is_ter)
	1666	. + falbe(i,is_sic) * pctsrf(i,is_sic)
[282]	1667	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
	1668	. + falblw(i,is_lic) * pctsrf(i,is_lic)
	1669	. + falblw(i,is_ter) * pctsrf(i,is_ter)
	1670	. + falblw(i,is_sic) * pctsrf(i,is_sic)
[2]	1671	ENDDO
[295]	1672	! if (debut) then
	1673	! albsol1 = albsol
	1674	! albsollw1 = albsollw
	1675	! endif
	1676	! albsol = albsol1
	1677	! albsollw = albsollw1
[2]	1678	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
[433]	1679	cIM e (dist, rmu0, fract, co2_ppm, solaire,
	1680	e (dist, rmu0, fract,
[282]	1681	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
	1682	e wo,
[2]	1683	e cldfra, cldemi, cldtau,
	1684	s heat,heat0,cool,cool0,radsol,albpla,
	1685	s topsw,toplw,solsw,sollw,
[177]	1686	s sollwdown,
[411]	1687	cccIMs topsw0,toplw0,solsw0,sollw0)
	1688	s topsw0,toplw0,solsw0,sollw0,
	1689	s ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
[2]	1690	itaprad = 0
	1691	ENDIF
	1692	itaprad = itaprad + 1
[433]	1693
[2]	1694	c
	1695	c Ajouter la tendance des rayonnements (tous les pas)
	1696	c
	1697	DO k = 1, klev
	1698	DO i = 1, klon
	1699	t_seri(i,k) = t_seri(i,k)
	1700	. + (heat(i,k)-cool(i,k)) * dtime/86400.
	1701	ENDDO
	1702	ENDDO
	1703	c
[385]	1704	IF (if_ebil.ge.2) THEN
	1705	ztit='after rad'
	1706	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1707	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1708	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1709	call diagphy(paire,ztit,ip_ebil
	1710	e , topsw, toplw, solsw, sollw, zero_v
	1711	e , zero_v, zero_v, zero_v, ztsol
	1712	e , d_h_vcol, d_qt, d_ec
	1713	s , fs_bound, fq_bound )
	1714	END IF
	1715	c
	1716	c
[2]	1717	c Calculer l'hydrologie de la surface
	1718	c
[98]	1719	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	1720	c . agesno, ftsol,fqsol,fsnow, ruis)
[2]	1721	c
	1722	DO i = 1, klon
	1723	zxqsol(i) = 0.0
	1724	zxsnow(i) = 0.0
	1725	ENDDO
	1726	DO nsrf = 1, nbsrf
	1727	DO i = 1, klon
	1728	zxqsol(i) = zxqsol(i) + fqsol(i,nsrf)*pctsrf(i,nsrf)
	1729	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
	1730	ENDDO
	1731	ENDDO
	1732	c
	1733	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
	1734	c
[411]	1735	cXXX DO nsrf = 1, nbsrf
	1736	cXXX DO i = 1, klon
	1737	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
	1738	cXXX fqsol(i,nsrf) = zxqsol(i)
	1739	cXXX fsnow(i,nsrf) = zxsnow(i)
	1740	cXXX ENDIF
	1741	cXXX ENDDO
	1742	cXXX ENDDO
[2]	1743	c
	1744	c Calculer le bilan du sol et la derive de temperature (couplage)
	1745	c
	1746	DO i = 1, klon
[391]	1747	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	1748	c a la demande de JLD
	1749	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
[2]	1750	ENDDO
	1751	c
	1752	cmoddeblott(jan95)
	1753	c Appeler le programme de parametrisation de l'orographie
	1754	c a l'echelle sous-maille:
	1755	c
	1756	IF (ok_orodr) THEN
	1757	c
	1758	c selection des points pour lesquels le shema est actif:
	1759	igwd=0
	1760	DO i=1,klon
	1761	itest(i)=0
	1762	c IF ((zstd(i).gt.10.0)) THEN
	1763	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	1764	itest(i)=1
	1765	igwd=igwd+1
	1766	idx(igwd)=i
	1767	ENDIF
	1768	ENDDO
[158]	1769	c igwdim=MAX(1,igwd)
[2]	1770	c
	1771	CALL drag_noro(klon,klev,dtime,paprs,pplay,
	1772	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
[158]	1773	e igwd,idx,itest,
[2]	1774	e t_seri, u_seri, v_seri,
	1775	s zulow, zvlow, zustr, zvstr,
	1776	s d_t_oro, d_u_oro, d_v_oro)
	1777	c
	1778	c ajout des tendances
	1779	DO k = 1, klev
	1780	DO i = 1, klon
	1781	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
	1782	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
	1783	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
	1784	ENDDO
	1785	ENDDO
	1786	c
	1787	ENDIF ! fin de test sur ok_orodr
	1788	c
	1789	IF (ok_orolf) THEN
	1790	c
	1791	c selection des points pour lesquels le shema est actif:
	1792	igwd=0
	1793	DO i=1,klon
	1794	itest(i)=0
	1795	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	1796	itest(i)=1
	1797	igwd=igwd+1
	1798	idx(igwd)=i
	1799	ENDIF
	1800	ENDDO
[158]	1801	c igwdim=MAX(1,igwd)
[2]	1802	c
	1803	CALL lift_noro(klon,klev,dtime,paprs,pplay,
[158]	1804	e rlat,zmea,zstd,zpic,
	1805	e itest,
[2]	1806	e t_seri, u_seri, v_seri,
	1807	s zulow, zvlow, zustr, zvstr,
	1808	s d_t_lif, d_u_lif, d_v_lif)
	1809	c
	1810	c ajout des tendances
	1811	DO k = 1, klev
	1812	DO i = 1, klon
	1813	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
	1814	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
	1815	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
	1816	ENDDO
	1817	ENDDO
	1818	c
	1819	ENDIF ! fin de test sur ok_orolf
	1820	c
[385]	1821	IF (if_ebil.ge.2) THEN
	1822	ztit='after orography'
	1823	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	1824	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1825	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1826	END IF
	1827	c
	1828	c
[2]	1829	cAA
	1830	cAA Installation de l'interface online-offline pour traceurs
	1831	cAA
	1832	c====================================================================
	1833	c Calcul des tendances traceurs
	1834	c====================================================================
[230]	1835	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
	1836	cKE43 des traceurs => il faut donc mettre des flags a .false.
[301]	1837	IF (iflag_con.GE.3) THEN
[230]	1838	c on ajoute les tendances calculees par KE43
[411]	1839	cXXX OM on onhibe la convection sur les traceurs
[230]	1840	DO iq=1, nqmax-2 ! Sandrine a -3 ???
[411]	1841	cXXX OM on inhibe la convection sur les traceur
	1842	cXXX DO k = 1, nlev
	1843	cXXX DO i = 1, klon
	1844	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
	1845	cXXX ENDDO
	1846	cXXX ENDDO
[230]	1847	WRITE(iqn,'(i2.2)') iq
	1848	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
	1849	ENDDO
[2]	1850	CMAF modif pour garder info du nombre de traceurs auxquels
	1851	C la physique s'applique
[230]	1852	ELSE
	1853	CMAF modif pour garder info du nombre de traceurs auxquels
	1854	C la physique s'applique
[2]	1855	C
	1856	call phytrac (rnpb,
[230]	1857	I debut,lafin,
[2]	1858	I nqmax-2,
	1859	I nlon,nlev,dtime,
	1860	I t,paprs,pplay,
	1861	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	1862	I ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
	1863	I frac_impa, frac_nucl,
	1864	I rlon,presnivs,paire,pphis,
	1865	O tr_seri)
[230]	1866	ENDIF
[2]	1867
	1868	IF (offline) THEN
[53]	1869
	1870	call phystokenc (
	1871	I nlon,nlev,pdtphys,rlon,rlat,
[230]	1872	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[2]	1873	I ycoefh,yu1,yv1,ftsol,pctsrf,
[53]	1874	I frac_impa, frac_nucl,
[230]	1875	I pphis,paire,dtime,itap)
[2]	1876
[230]	1877
[2]	1878	ENDIF
	1879
	1880	c
	1881	c Calculer le transport de l'eau et de l'energie (diagnostique)
	1882	c
	1883	CALL transp (paprs,zxtsol,
	1884	e t_seri, q_seri, u_seri, v_seri, zphi,
	1885	s ve, vq, ue, uq)
	1886	c
	1887	c Accumuler les variables a stocker dans les fichiers histoire:
	1888	c
	1889	c
	1890	c
[411]	1891	c+jld ec_conser
	1892	DO k = 1, klev
	1893	DO i = 1, klon
	1894	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
	1895	d_t_ec(i,k)=0.5/ZRCPD
	1896	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	1897	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	1898	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	1899	END DO
	1900	END DO
	1901	c-jld ec_conser
[385]	1902	IF (if_ebil.ge.1) THEN
	1903	ztit='after physic'
	1904	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
	1905	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1906	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1907	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1908	C on devrait avoir que la variation d'entalpie par la dynamique
	1909	C est egale a la variation de la physique au pas de temps precedent.
	1910	C Donc la somme de ces 2 variations devrait etre nulle.
	1911	call diagphy(paire,ztit,ip_ebil
	1912	e , topsw, toplw, solsw, sollw, sens
	1913	e , evap, rain_fall, snow_fall, ztsol
	1914	e , d_h_vcol, d_qt, d_ec
	1915	s , fs_bound, fq_bound )
	1916	C
	1917	d_h_vcol_phy=d_h_vcol
	1918	C
	1919	END IF
	1920	C
[411]	1921	cccIM cf. FH
	1922	c=======================================================================
	1923	c SORTIES
	1924	c=======================================================================
[158]	1925
[411]	1926	c Interpollation sur quelques niveaux de pression
	1927	c -----------------------------------------------
[271]	1928
[411]	1929	call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850)
	1930	call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850)
	1931	call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500)
	1932	call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500)
	1933	call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200)
	1934	call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200)
	1935	call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500)
[433]	1936	call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
[271]	1937
[433]	1938	slp(:) = paprs(:,1)exp(pphis(:)/(289.t_seri(:,1)))
	1939	c
	1940
	1941
[411]	1942	c=============================================================
	1943	c Ecriture des sorties
	1944	c=============================================================
[271]	1945
[411]	1946	#ifdef histhf
	1947	#include "write_histhf.h"
	1948	#endif
[158]	1949
[411]	1950	#include "write_histday.h"
	1951	#include "write_histmth.h"
	1952	#include "write_histins.h"
[29]	1953
[411]	1954	c=============================================================
[2]	1955	c
	1956	c Convertir les incrementations en tendances
	1957	c
	1958	DO k = 1, klev
	1959	DO i = 1, klon
	1960	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	1961	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	1962	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	1963	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	1964	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	1965	ENDDO
	1966	ENDDO
	1967	c
	1968	IF (nqmax.GE.3) THEN
	1969	DO iq = 3, nqmax
	1970	DO k = 1, klev
	1971	DO i = 1, klon
	1972	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	1973	ENDDO
	1974	ENDDO
	1975	ENDDO
	1976	ENDIF
	1977	c
[46]	1978	c Sauvegarder les valeurs de t et q a la fin de la physique:
	1979	c
	1980	DO k = 1, klev
	1981	DO i = 1, klon
	1982	t_ancien(i,k) = t_seri(i,k)
	1983	q_ancien(i,k) = q_seri(i,k)
	1984	ENDDO
	1985	ENDDO
	1986	c
[2]	1987	c====================================================================
	1988	c Si c'est la fin, il faut conserver l'etat de redemarrage
	1989	c====================================================================
	1990	c
	1991	IF (lafin) THEN
[353]	1992	itau_phy = itau_phy + itap
[46]	1993	ccc IF (ok_oasis) CALL quitcpl
[2]	1994	CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
[98]	1995	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
	1996	. falbe, fevap, rain_fall, snow_fall,
[258]	1997	. solsw, sollwdown,dlw,
[112]	1998	. radsol,frugs,agesno,
[98]	1999	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
[373]	2000	. t_ancien, q_ancien, rnebcon, ratqs, clwcon)
[2]	2001	ENDIF
	2002
	2003	RETURN
	2004	END
[46]	2005	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
[2]	2006	IMPLICIT none
	2007	c
	2008	c Calculer et imprimer l'eau totale. A utiliser pour verifier
	2009	c la conservation de l'eau
	2010	c
	2011	#include "YOMCST.h"
	2012	INTEGER klon,klev
	2013	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
[46]	2014	REAL aire(klon)
	2015	REAL qtotal, zx, qcheck
[2]	2016	INTEGER i, k
	2017	c
[46]	2018	zx = 0.0
	2019	DO i = 1, klon
	2020	zx = zx + aire(i)
	2021	ENDDO
[2]	2022	qtotal = 0.0
	2023	DO k = 1, klev
	2024	DO i = 1, klon
[46]	2025	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
[2]	2026	. *(paprs(i,k)-paprs(i,k+1))/RG
	2027	ENDDO
	2028	ENDDO
	2029	c
[46]	2030	qcheck = qtotal/zx
[2]	2031	c
[46]	2032	RETURN
[2]	2033	END
	2034	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	2035	IMPLICIT none
	2036	c
	2037	c Tranformer une variable de la grille physique a
	2038	c la grille d'ecriture
	2039	c
	2040	INTEGER nfield,nlon,iim,jjmp1, jjm
[15]	2041	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
[2]	2042	c
[158]	2043	INTEGER i, n, ig
[2]	2044	c
	2045	jjm = jjmp1 - 1
	2046	DO n = 1, nfield
	2047	DO i=1,iim
[15]	2048	ecrit(i,n) = fi(1,n)
	2049	ecrit(i+jjm*iim,n) = fi(nlon,n)
[2]	2050	ENDDO
[15]	2051	DO ig = 1, nlon - 2
	2052	ecrit(iim+ig,n) = fi(1+ig,n)
[2]	2053	ENDDO
	2054	ENDDO
	2055	RETURN
	2056	END
[15]	2057

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

Download in other formats:

Original Format