Context Navigation

physiq.F @ 455

Last change on this file since 455 was 449, checked in by lmdzadmin, 21 years ago
Oublie de remplacer fqsol par qsol dans l'appel de phyetat0 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

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE physiq (nlon,nlev,nqmax,
5	. debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
6	. paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
7	. u,v,t,qx,
8	. omega, cufi, cvfi,
9	. d_u, d_v, d_t, d_qx, d_ps)
10	USE ioipsl
11	USE histcom
12	USE writephys
13
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)
48	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
49	c omega---input-R-vitesse verticale en Pa/s
50	c cufi----input-R-resolution des mailles en x (m)
51	c cvfi----input-R-resolution des mailles en y (m)
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"
60	integer jjmp1
61	parameter (jjmp1=jjm+1-1/jjm)
62	#include "dimphy.h"
63	#include "regdim.h"
64	#include "indicesol.h"
65	#include "dimsoil.h"
66	#include "clesphys.h"
67	#include "control.h"
68	#include "temps.h"
69	c======================================================================
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======================================================================
75	LOGICAL check ! Verifier la conservation du modele en eau
76	PARAMETER (check=.FALSE.)
77	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
78	PARAMETER (ok_stratus=.FALSE.)
79	c======================================================================
80	c Parametres lies au coupleur OASIS:
81	#include "oasis.h"
82	INTEGER,SAVE :: npas, nexca
83	logical rnpb
84	parameter(rnpb=.true.)
85	c PARAMETER (npas=1440)
86	c PARAMETER (nexca=48)
87	EXTERNAL fromcpl, intocpl, inicma
88	c ocean = type de modele ocean a utiliser: force, slab, couple
89	character*6 ocean
90	SAVE ocean
91
92	c parameter (ocean = 'force ')
93	c parameter (ocean = 'couple')
94	logical ok_ocean
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.)
103	logical ok_veget
104	save ok_veget
105	c parameter (ok_veget = .true.)
106	c parameter (ok_veget = .false.)
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
123	save ok_journe
124	c PARAMETER (ok_journe=.true.)
125
126	integer lev_histday
127	save lev_histday
128	data lev_histday/1/
129	c
130	LOGICAL ok_mensuel ! sortir le fichier mensuel
131	save ok_mensuel
132	c PARAMETER (ok_mensuel=.true.)
133	c
134	LOGICAL ok_instan ! sortir le fichier instantane
135	save ok_instan
136	c PARAMETER (ok_instan=.true.)
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)
146
147	c
148	c
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)
165	REAL zsurf(nbsrf)
166	real cufi(klon), cvfi(klon)
167
168	REAL u(klon,klev)
169	REAL v(klon,klev)
170	REAL t(klon,klev)
171	REAL qx(klon,klev,nqmax)
172
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
178	REAL d_t_dyn(klon,klev)
179	REAL d_q_dyn(klon,klev)
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
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)
200	real u500(klon),v500(klon),phi500(klon),w500(klon)
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
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
249	REAL co2_ppm_etat0
250	c
251	REAL solaire_etat0
252	c
253	real slp(klon) ! sea level pressure
254
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
261	REAL fevap(klon,nbsrf)
262	SAVE fevap ! evaporation
263	REAL fluxlat(klon,nbsrf)
264	SAVE fluxlat
265	c
266	REAL deltat(klon)
267	SAVE deltat ! ecart avec la SST de reference
268	c
269	REAL fqsurf(klon,nbsrf)
270	SAVE fqsurf ! humidite de l'air au contact de la surface
271	c
272	REAL qsol(klon)
273	SAVE qsol ! hauteur d'eau dans le sol
274	c
275	REAL fsnow(klon,nbsrf)
276	SAVE fsnow ! epaisseur neigeuse
277	c
278	REAL falbe(klon,nbsrf)
279	SAVE falbe ! albedo par type de surface
280	REAL falblw(klon,nbsrf)
281	SAVE falblw ! albedo par type de surface
282
283	c
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
316	INTEGER igwd,idx(klon),itest(klon)
317	c
318	REAL agesno(klon,nbsrf)
319	SAVE agesno ! age de la neige
320	c
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):
325	c
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
344	REAL qcondc(klon,klev) ! in-cld water content from convect
345	SAVE qcondc
346	REAL ema_work1(klon, klev), ema_work2(klon, klev)
347	SAVE ema_work1, ema_work2
348	REAL wdn(klon), tdn(klon), qdn(klon)
349
350	REAL wd(klon) ! sb
351	SAVE wd ! sb
352
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
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"
370	PARAMETER (offline=.false.)
371	INTEGER physid
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
383	save snow_fall, rain_fall
384	REAL evap(klon), devap(klon) ! evaporation et sa derivee
385	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
386	REAL dlw(klon) ! derivee infra rouge
387	REAL bils(klon) ! bilan de chaleur au sol
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
391	REAL fder(klon) ! Derive de flux (sensible et latente)
392	save fder
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
399	save frugs
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
412	REAL albsollw(klon)
413	SAVE albsollw ! albedo du sol total
414	REAL albsol1(klon)
415	SAVE albsol1 ! albedo du sol total
416	REAL albsollw1(klon)
417	SAVE albsollw1 ! albedo du sol total
418
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)
432	cKE43
433	EXTERNAL conema3 ! convect4.3
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
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
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
469	CXXX PB
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
474	c
475	REAL zxfluxt(klon, klev)
476	REAL zxfluxq(klon, klev)
477	REAL zxfluxu(klon, klev)
478	REAL zxfluxv(klon, klev)
479	CXXX
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)
485	real sollwdown(klon) ! downward LW flux at surface
486	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
487	REAL albpla(klon)
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
491	c Le rayonnement n'est pas calcule tous les pas, il faut donc
492	c sauvegarder les sorties du rayonnement
493	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
494	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
495	cccIM
496	SAVE ZFSUP,ZFSDN,ZFSUP0,ZFSDN0
497
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
507	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
508	c
509	REAL dist, rmu0(klon), fract(klon)
510	REAL zdtime, zlongi
511	c
512	CHARACTER*2 str2
513	CHARACTER*2 iqn
514	c
515	REAL qcheck
516	REAL z_avant(klon), z_apres(klon), z_factor(klon)
517	LOGICAL zx_ajustq
518	c
519	REAL za, zb
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
523	REAL t_coup
524	PARAMETER (t_coup=234.0)
525	c
526	REAL zphi(klon,klev)
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):
532	c
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
539	cccIM
540	CHARACTER*40 capemaxcels
541
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
547	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
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
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)
570	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
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)
579	REAL prfl(klon,klev+1), psfl(klon,klev+1)
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)
593
594	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
595	real ratqsbas,ratqshaut
596	save ratqsbas,ratqshaut, ratqs
597	real zpt_conv(klon,klev)
598
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
605	real facteur
606
607	integer iflag_cldcon
608	save iflag_cldcon
609
610	logical ptconv(klon,klev)
611
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
627	integer itau_w ! pas de temps ecriture = itap + itau_phy
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)
633	REAL ql_seri(klon,klev),qs_seri(klon,klev)
634	REAL u_seri(klon,klev), v_seri(klon,klev)
635	c
636	REAL tr_seri(klon,klev,nbtr)
637	REAL d_tr(klon,klev,nbtr)
638
639	REAL zx_rh(klon,klev)
640
641	INTEGER length
642	PARAMETER ( length = 100 )
643	REAL tabcntr0( length )
644	c
645	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
646	REAL zx_tmp_fi2d(klon)
647	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
648	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
649	c
650	INTEGER nid_day, nid_mth, nid_ins
651	SAVE nid_day, nid_mth, nid_ins
652	c
653	INTEGER nhori, nvert
654	REAL zsto, zout
655	real zjulian
656	save zjulian
657
658	character*20 modname
659	character*80 abort_message
660	logical ok_sync
661	real date0
662	integer idayref
663
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
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/
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
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'
709	IF (if_ebil.ge.1) THEN
710	DO i=1,klon
711	zero_v(i)=0.
712	END DO
713	END IF
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
722
723
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
731	C
732	IF (if_ebil.ge.1) d_h_vcol_phy=0.
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,
737	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
738	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil)
739	cIM . , RI0)
740
741	c
742	c
743	c Initialiser les compteurs:
744	c
745
746	frugs = 0.
747	itap = 0
748	itaprad = 0
749	c
750	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
751	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
752	. falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
753	. dlw,radsol,frugs,agesno,clesphy0,
754	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
755	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
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
787	PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl
788	c
789	cKE43
790	c Initialisation pour la convection de K.E. (sb):
791	IF (iflag_con.GE.3) THEN
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
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
812	ENDIF
813
814	c34EK
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
843	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
844	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
845	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
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
853
854	c
855	c Initialiser le couplage si necessaire
856	c
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
869	cccIM
870	capemaxcels = 't_max(X)'
871	t2mincels = 't_min(X)'
872	t2maxcels = 't_max(X)'
873
874	cccIM cf. FH
875	c
876	c=============================================================
877	c Initialisation des sorties
878	c=============================================================
879	#ifdef histhf
880	#include "ini_histhf.h"
881	#endif
882
883	#include "ini_histday.h"
884	#include "ini_histmth.h"
885	#include "ini_histins.h"
886
887	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
888	date0 = zjulian
889	C date0 = day_ini
890	WRITE(,) 'physiq date0 : ',date0
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)
938	qs_seri(i,k) = 0.
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
956	C
957	IF (if_ebil.ge.1) THEN
958	DO i = 1, klon
959	ztsol(i) = 0.
960	ENDDO
961	DO nsrf = 1, nbsrf
962	DO i = 1, klon
963	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
964	ENDDO
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
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
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
1021	PRINT *,' PHYS cond julien ',julien
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))
1030	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1031	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
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
1044	IF (if_ebil.ge.2) THEN
1045	ztit='after reevap'
1046	CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime
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
1058	c Appeler la diffusion verticale (programme de couche limite)
1059	c
1060	DO i = 1, klon
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
1067	zxrugs(i) = 0.0
1068	ENDDO
1069	DO nsrf = 1, nbsrf
1070	DO i = 1, klon
1071	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1072	cccc frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1073	ENDDO
1074	ENDDO
1075	DO nsrf = 1, nbsrf
1076	DO i = 1, klon
1077	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1078	ENDDO
1079	ENDDO
1080	c
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
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
1098	fder = dlw
1099
1100	CALL clmain(dtime,itap,date0,pctsrf,
1101	e t_seri,q_seri,u_seri,v_seri,
1102	e julien, rmu0,
1103	e ok_veget, ocean, npas, nexca, ftsol,
1104	$ soil_model,ftsoil, qsol,
1105	$ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1106	$ fluxlat,
1107	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1108	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1109	e rlon, rlat, cufi, cvfi, frugs,
1110	e debut, lafin, agesno,rugoro ,
1111	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1112	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1113	s dsens, devap,
1114	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m)
1115	c
1116	CXXX PB
1117	CXXX Incrementation des flux
1118	CXXX
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
1137	DO i = 1, klon
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
1141	fder(i) = dlw(i) + dsens(i) + devap(i)
1142	ENDDO
1143
1144
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
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
1167	c Incrementer la temperature du sol
1168	c
1169	DO i = 1, klon
1170	zxtsol(i) = 0.0
1171	zxfluxlat(i) = 0.0
1172	cccIM
1173	zt2m(i) = 0.0
1174	zq2m(i) = 0.0
1175	zu10m(i) = 0.0
1176	zv10m(i) = 0.0
1177	c
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
1184	ENDDO
1185	DO nsrf = 1, nbsrf
1186	DO i = 1, klon
1187	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1188	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1189	cIM cf. JLD
1190	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1191	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1192	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1193	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
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
1200	ENDDO
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
1207	DO i = 1, klon
1208	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
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)
1214	ENDDO
1215	ENDDO
1216	c
1217	c
1218	c Calculer la derive du flux infrarouge
1219	c
1220	cXXX DO nsrf = 1, nbsrf
1221	DO i = 1, klon
1222	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1223	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1224	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1225	cXXX . *pctsrf(i,nsrf)
1226	cXXX ENDIF
1227	cXXX ENDDO
1228	ENDDO
1229	c
1230	c Appeler la convection (au choix)
1231	c
1232	DO k = 1, klev
1233	DO i = 1, klon
1234	conv_q(i,k) = d_q_dyn(i,k)
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
1241	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1242	PRINT*, "avantcon=", za
1243	ENDIF
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
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,
1264	e conv_t, conv_q, zxfluxq(1,1), omega,
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)
1268	WHERE (rain_con < 0.) rain_con = 0.
1269	WHERE (snow_con < 0.) snow_con = 0.
1270	DO i = 1, klon
1271	ibas_con(i) = klev+1 - kcbot(i)
1272	itop_con(i) = klev+1 - kctop(i)
1273	ENDDO
1274	ELSE IF (iflag_con.GE.3) THEN
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
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)
1297	cIM cf. FH
1298	clwcon0=qcondc
1299
1300	ELSE ! ok_cvl
1301
1302	c print*,'Avant conema OUI'
1303	CALL conema3 (dtime,
1304	. paprs,pplay,t_seri,q_seri,
1305	. u_seri,v_seri,tr_seri,nbtr,
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,
1309	. upwd,dnwd,dnwd0,bas,top,
1310	. Ma,cape,tvp,rflag,
1311	. pbase
1312	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1313	. ,clwcon0)
1314	print*,'Apres conema3 '
1315
1316	ENDIF ! ok_cvl
1317
1318	IF (.NOT. ok_gust) THEN
1319	do i = 1, klon
1320	wd(i)=0.0
1321	enddo
1322	ENDIF
1323
1324	c =================================================================== c
1325	c Calcul des proprietes des nuages convectifs
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
1347	c calcul des proprietes des nuages convectifs
1348	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1349	call clouds_gno
1350	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1351
1352	c =================================================================== c
1353
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
1363	ELSE
1364	PRINT*, "iflag_con non-prevu", iflag_con
1365	CALL abort
1366	ENDIF
1367
1368	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1369	c . d_u_con, d_v_con)
1370	DO k = 1, klev
1371	DO i = 1, klon
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)
1376	ENDDO
1377	ENDDO
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
1391	IF (check) THEN
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
1397	za = za + paire(i)/FLOAT(klon)
1398	zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1399	ENDDO
1400	zx_t = zx_t/za*dtime
1401	PRINT*, "Precip=", zx_t
1402	ENDIF
1403	IF (zx_ajustq) THEN
1404	DO i = 1, klon
1405	z_apres(i) = 0.0
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
1425	ENDIF
1426	zx_ajustq=.FALSE.
1427	c
1428	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1429	c
1430	IF (iflag_con .NE. 2 .AND. debut) THEN
1431	PRINT*, 'Pour l instant, seul conflx fonctionne ',
1432	$ 'avec traceurs', iflag_con
1433	PRINT*,' Mettre iflag_con',
1434	$ ' = 2 dans run.def et repasser'
1435	c CALL abort
1436	ENDIF
1437	c
1438	ENDIF !--nqmax.GT.2
1439	c
1440	c Appeler l'ajustement sec
1441	c
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
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
1456
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
1494	c facttemps=exp(-pdtphys/1.e4)
1495	facteur=exp(-pdtphys*facttemps)
1496	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1497	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1498	c print*,'calcul des ratqs fini'
1499	else
1500	c on ne prend que le ratqs stable pour fisrtilp
1501	ratqs(:,:)=ratqss(:,:)
1502	endif
1503
1504
1505	c
1506	c Appeler le processus de condensation a grande echelle
1507	c et le processus de precipitation
1508	c-------------------------------------------------------------------------
1509	CALL fisrtilp(dtime,paprs,pplay,
1510	. t_seri, q_seri,ptconv,ratqs,
1511	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1512	. rain_lsc, snow_lsc,
1513	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
1514	. frac_impa, frac_nucl,
1515	. prfl, psfl, rhcl)
1516
1517	WHERE (rain_lsc < 0) rain_lsc = 0.
1518	WHERE (snow_lsc < 0) snow_lsc = 0.
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
1529	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1530	PRINT*, "apresilp=", za
1531	zx_t = 0.0
1532	za = 0.0
1533	DO i = 1, klon
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
1538	PRINT*, "Precip=", zx_t
1539	ENDIF
1540	c
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
1553	c-------------------------------------------------------------------
1554	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1555	c-------------------------------------------------------------------
1556
1557	c 1. NUAGES CONVECTIFS
1558	c
1559	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
1560
1561	c Nuages diagnostiques pour Tiedtke
1562	CALL diagcld1(paprs,pplay,
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
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
1578	c facttemps=pdtphys/1.e4
1579	facteur = pdtphys *facttemps
1580	do k=1,klev
1581	do i=1,klon
1582	rnebcon(i,k)=rnebcon(i,k)*facteur
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
1596	ENDIF
1597	c
1598	c 2. NUAGES STARTIFORMES
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
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
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
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
1645	zqsat(i,k)=zx_qs
1646	ENDDO
1647	ENDDO
1648	c
1649	c Calculer les parametres optiques des nuages et quelques
1650	c parametres pour diagnostiques:
1651	c
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
1657	CALL nuage (paprs, pplay,
1658	. t_seri, cldliq, cldfra, cldtau, cldemi,
1659	. cldh, cldl, cldm, cldt, cldq)
1660	endif
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
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)
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)
1674	ENDDO
1675	! if (debut) then
1676	! albsol1 = albsol
1677	! albsollw1 = albsollw
1678	! endif
1679	! albsol = albsol1
1680	! albsollw = albsollw1
1681	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
1682	cIM e (dist, rmu0, fract, co2_ppm, solaire,
1683	e (dist, rmu0, fract,
1684	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
1685	e wo,
1686	e cldfra, cldemi, cldtau,
1687	s heat,heat0,cool,cool0,radsol,albpla,
1688	s topsw,toplw,solsw,sollw,
1689	s sollwdown,
1690	cccIMs topsw0,toplw0,solsw0,sollw0)
1691	s topsw0,toplw0,solsw0,sollw0,
1692	s ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
1693	itaprad = 0
1694	ENDIF
1695	itaprad = itaprad + 1
1696
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
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
1720	c Calculer l'hydrologie de la surface
1721	c
1722	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
1723	c . agesno, ftsol,fqsurf,fsnow, ruis)
1724	c
1725	DO i = 1, klon
1726	zxqsurf(i) = 0.0
1727	zxsnow(i) = 0.0
1728	ENDDO
1729	DO nsrf = 1, nbsrf
1730	DO i = 1, klon
1731	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
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
1738	cXXX DO nsrf = 1, nbsrf
1739	cXXX DO i = 1, klon
1740	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
1741	cXXX fqsurf(i,nsrf) = zxqsurf(i)
1742	cXXX fsnow(i,nsrf) = zxsnow(i)
1743	cXXX ENDIF
1744	cXXX ENDDO
1745	cXXX ENDDO
1746	c
1747	c Calculer le bilan du sol et la derive de temperature (couplage)
1748	c
1749	DO i = 1, klon
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)
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
1772	c igwdim=MAX(1,igwd)
1773	c
1774	CALL drag_noro(klon,klev,dtime,paprs,pplay,
1775	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1776	e igwd,idx,itest,
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
1804	c igwdim=MAX(1,igwd)
1805	c
1806	CALL lift_noro(klon,klev,dtime,paprs,pplay,
1807	e rlat,zmea,zstd,zpic,
1808	e itest,
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
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
1832	cAA
1833	cAA Installation de l'interface online-offline pour traceurs
1834	cAA
1835	c====================================================================
1836	c Calcul des tendances traceurs
1837	c====================================================================
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.
1840	IF (iflag_con.GE.3) THEN
1841	c on ajoute les tendances calculees par KE43
1842	cXXX OM on onhibe la convection sur les traceurs
1843	DO iq=1, nqmax-2 ! Sandrine a -3 ???
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
1850	WRITE(iqn,'(i2.2)') iq
1851	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
1852	ENDDO
1853	CMAF modif pour garder info du nombre de traceurs auxquels
1854	C la physique s'applique
1855	ELSE
1856	CMAF modif pour garder info du nombre de traceurs auxquels
1857	C la physique s'applique
1858	C
1859	call phytrac (rnpb,
1860	I debut,lafin,
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)
1869	ENDIF
1870
1871	IF (offline) THEN
1872
1873	call phystokenc (
1874	I nlon,nlev,pdtphys,rlon,rlat,
1875	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1876	I ycoefh,yu1,yv1,ftsol,pctsrf,
1877	I frac_impa, frac_nucl,
1878	I pphis,paire,dtime,itap)
1879
1880
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
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
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
1924	cccIM cf. FH
1925	c=======================================================================
1926	c SORTIES
1927	c=======================================================================
1928
1929	c Interpollation sur quelques niveaux de pression
1930	c -----------------------------------------------
1931
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)
1939	call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
1940
1941	slp(:) = paprs(:,1)exp(pphis(:)/(289.t_seri(:,1)))
1942	c
1943
1944
1945	c=============================================================
1946	c Ecriture des sorties
1947	c=============================================================
1948
1949	#ifdef histhf
1950	#include "write_histhf.h"
1951	#endif
1952
1953	#include "write_histday.h"
1954	#include "write_histmth.h"
1955	#include "write_histins.h"
1956
1957	c=============================================================
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
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
1990	c====================================================================
1991	c Si c'est la fin, il faut conserver l'etat de redemarrage
1992	c====================================================================
1993	c
1994	IF (lafin) THEN
1995	itau_phy = itau_phy + itap
1996	ccc IF (ok_oasis) CALL quitcpl
1997	CALL phyredem ("restartphy.nc",dtime,radpas,
1998	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
1999	. fsnow, falbe, fevap, rain_fall, snow_fall,
2000	. solsw, sollwdown,dlw,
2001	. radsol,frugs,agesno,
2002	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2003	. t_ancien, q_ancien, rnebcon, ratqs, clwcon)
2004	ENDIF
2005
2006	RETURN
2007	END
2008	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
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)
2017	REAL aire(klon)
2018	REAL qtotal, zx, qcheck
2019	INTEGER i, k
2020	c
2021	zx = 0.0
2022	DO i = 1, klon
2023	zx = zx + aire(i)
2024	ENDDO
2025	qtotal = 0.0
2026	DO k = 1, klev
2027	DO i = 1, klon
2028	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2029	. *(paprs(i,k)-paprs(i,k+1))/RG
2030	ENDDO
2031	ENDDO
2032	c
2033	qcheck = qtotal/zx
2034	c
2035	RETURN
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
2044	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2045	c
2046	INTEGER i, n, ig
2047	c
2048	jjm = jjmp1 - 1
2049	DO n = 1, nfield
2050	DO i=1,iim
2051	ecrit(i,n) = fi(1,n)
2052	ecrit(i+jjm*iim,n) = fi(nlon,n)
2053	ENDDO
2054	DO ig = 1, nlon - 2
2055	ecrit(iim+ig,n) = fi(1+ig,n)
2056	ENDDO
2057	ENDDO
2058	RETURN
2059	END
2060

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

Download in other formats:

Original Format