Context Navigation

physiq.F @ 1962

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

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

Download in other formats:

Original Format