Context Navigation

physiq.F @ 444

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

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

Download in other formats:

Original Format