source: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F @ 385

Last change on this file since 385 was 385, checked in by lmdzadmin, 22 years ago

Definition d'un temps de relaxation pour le schema de nuages
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 130.6 KB
Line 
1c
2c $Header$
3c
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
15c======================================================================
16c
17c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
18c
19c Objet: Moniteur general de la physique du modele
20cAA      Modifications quant aux traceurs :
21cAA                  -  uniformisation des parametrisations ds phytrac
22cAA                  -  stockage des moyennes des champs necessaires
23cAA                     en mode traceur off-line
24c======================================================================
25c    modif   ( P. Le Van ,  12/10/98 )
26c
27c  Arguments:
28c
29c nlon----input-I-nombre de points horizontaux
30c nlev----input-I-nombre de couches verticales
31c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
32c debut---input-L-variable logique indiquant le premier passage
33c lafin---input-L-variable logique indiquant le dernier passage
34c rjour---input-R-numero du jour de l'experience
35c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
36c pdtphys-input-R-pas d'integration pour la physique (seconde)
37c paprs---input-R-pression pour chaque inter-couche (en Pa)
38c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
39c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
40c pphis---input-R-geopotentiel du sol
41c paire---input-R-aire de chaque maille
42c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
43c u-------input-R-vitesse dans la direction X (de O a E) en m/s
44c v-------input-R-vitesse Y (de S a N) en m/s
45c t-------input-R-temperature (K)
46c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
47c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
48c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
49c omega---input-R-vitesse verticale en Pa/s
50c cufi----input-R-resolution des mailles en x (m)
51c cvfi----input-R-resolution des mailles en y (m)
52c
53c d_u-----output-R-tendance physique de "u" (m/s/s)
54c d_v-----output-R-tendance physique de "v" (m/s/s)
55c d_t-----output-R-tendance physique de "t" (K/s)
56c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
57c d_ps----output-R-tendance physique de la pression au sol
58c======================================================================
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"
69c======================================================================
70      LOGICAL check ! Verifier la conservation du modele en eau
71      PARAMETER (check=.FALSE.)
72      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
73      PARAMETER (ok_stratus=.FALSE.)
74c======================================================================
75c Parametres lies au coupleur OASIS:
76#include "oasis.h"
77      INTEGER,SAVE :: npas, nexca
78      logical rnpb
79      parameter(rnpb=.true.)
80c      PARAMETER (npas=1440)
81c      PARAMETER (nexca=48)
82      EXTERNAL fromcpl, intocpl, inicma
83c      ocean = type de modele ocean a utiliser: force, slab, couple
84      character*6 ocean
85      SAVE ocean
86
87c      parameter (ocean = 'force ')
88c     parameter (ocean = 'couple')
89      logical ok_ocean
90c======================================================================
91c Clef controlant l'activation du cycle diurne:
92ccc      LOGICAL cycle_diurne
93ccc      PARAMETER (cycle_diurne=.FALSE.)
94c======================================================================
95c Modele thermique du sol, a activer pour le cycle diurne:
96ccc      LOGICAL soil_model
97ccc      PARAMETER (soil_model=.FALSE.)
98      logical ok_veget
99      save ok_veget
100c     parameter (ok_veget = .true.)
101c      parameter (ok_veget = .false.)
102c======================================================================
103c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
104c le calcul du rayonnement est celle apres la precipitation des nuages.
105c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
106c la condensation et la precipitation. Cette cle augmente les impacts
107c radiatifs des nuages.
108ccc      LOGICAL new_oliq
109ccc      PARAMETER (new_oliq=.FALSE.)
110c======================================================================
111c Clefs controlant deux parametrisations de l'orographie:
112cc      LOGICAL ok_orodr
113ccc      PARAMETER (ok_orodr=.FALSE.)
114ccc      LOGICAL ok_orolf
115ccc      PARAMETER (ok_orolf=.FALSE.)
116c======================================================================
117      LOGICAL ok_journe ! sortir le fichier journalier
118      save ok_journe
119c      PARAMETER (ok_journe=.true.)
120c
121      LOGICAL ok_mensuel ! sortir le fichier mensuel
122      save ok_mensuel
123c      PARAMETER (ok_mensuel=.true.)
124c
125      LOGICAL ok_instan ! sortir le fichier instantane
126      save ok_instan
127c      PARAMETER (ok_instan=.true.)
128c
129      LOGICAL ok_region ! sortir le fichier regional
130      PARAMETER (ok_region=.FALSE.)
131c======================================================================
132c
133      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
134      PARAMETER (ivap=1)
135      INTEGER iliq          ! indice de traceurs pour eau liquide
136      PARAMETER (iliq=2)
137
138      INTEGER nvm           ! nombre de vegetations
139      PARAMETER (nvm=8)
140      REAL veget(klon,nvm)  ! couverture vegetale
141      SAVE veget
142
143c
144c
145c Variables argument:
146c
147      INTEGER nlon
148      INTEGER nlev
149      INTEGER nqmax
150      REAL rjourvrai, rjour_ecri
151      REAL gmtime
152      REAL pdtphys
153      LOGICAL debut, lafin
154      REAL paprs(klon,klev+1)
155      REAL pplay(klon,klev)
156      REAL pphi(klon,klev)
157      REAL pphis(klon)
158      REAL paire(klon)
159      REAL presnivs(klev)
160      REAL znivsig(klev)
161      REAL zsurf(nbsrf)
162      real cufi(klon), cvfi(klon)
163
164      REAL u(klon,klev)
165      REAL v(klon,klev)
166      REAL t(klon,klev)
167      REAL qx(klon,klev,nqmax)
168
169      REAL t_ancien(klon,klev), q_ancien(klon,klev)
170      SAVE t_ancien, q_ancien
171      LOGICAL ancien_ok
172      SAVE ancien_ok
173
174      REAL d_t_dyn(klon,klev)
175      REAL d_q_dyn(klon,klev)
176
177      REAL omega(klon,klev)
178
179      REAL d_u(klon,klev)
180      REAL d_v(klon,klev)
181      REAL d_t(klon,klev)
182      REAL d_qx(klon,klev,nqmax)
183      REAL d_ps(klon)
184
185      INTEGER        longcles
186      PARAMETER    ( longcles = 20 )
187      REAL clesphy0( longcles      )
188c
189c Variables quasi-arguments
190c
191      REAL xjour
192      SAVE xjour
193c
194c
195c Variables propres a la physique
196c
197      REAL dtime
198      SAVE dtime                  ! pas temporel de la physique
199c
200      INTEGER radpas
201      SAVE radpas                 ! frequence d'appel rayonnement
202c
203      REAL radsol(klon)
204      SAVE radsol                 ! bilan radiatif au sol
205c
206      REAL rlat(klon)
207      SAVE rlat                   ! latitude pour chaque point
208c
209      REAL rlon(klon)
210      SAVE rlon                   ! longitude pour chaque point
211c
212cc      INTEGER iflag_con
213cc      SAVE iflag_con              ! indicateur de la convection
214c
215      INTEGER itap
216      SAVE itap                   ! compteur pour la physique
217c
218      REAL co2_ppm
219      SAVE co2_ppm                ! concentration du CO2
220c
221      REAL solaire
222      SAVE solaire                ! constante solaire
223c
224      REAL ftsol(klon,nbsrf)
225      SAVE ftsol                  ! temperature du sol
226c
227      REAL ftsoil(klon,nsoilmx,nbsrf)
228      SAVE ftsoil                 ! temperature dans le sol
229c
230      REAL fevap(klon,nbsrf)
231      SAVE fevap                 ! evaporation
232      REAL fluxlat(klon,nbsrf)
233      SAVE fluxlat
234c
235      REAL deltat(klon)
236      SAVE deltat                 ! ecart avec la SST de reference
237c
238      REAL fqsol(klon,nbsrf)
239      SAVE fqsol                  ! humidite du sol
240c
241      REAL fsnow(klon,nbsrf)
242      SAVE fsnow                  ! epaisseur neigeuse
243c
244      REAL falbe(klon,nbsrf)
245      SAVE falbe                  ! albedo par type de surface
246      REAL falblw(klon,nbsrf)
247      SAVE falblw                 ! albedo par type de surface
248
249c
250c
251c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
252c
253      REAL zmea(klon)
254      SAVE zmea                   ! orographie moyenne
255c
256      REAL zstd(klon)
257      SAVE zstd                   ! deviation standard de l'OESM
258c
259      REAL zsig(klon)
260      SAVE zsig                   ! pente de l'OESM
261c
262      REAL zgam(klon)
263      save zgam                   ! anisotropie de l'OESM
264c
265      REAL zthe(klon)
266      SAVE zthe                   ! orientation de l'OESM
267c
268      REAL zpic(klon)
269      SAVE zpic                   ! Maximum de l'OESM
270c
271      REAL zval(klon)
272      SAVE zval                   ! Minimum de l'OESM
273c
274      REAL rugoro(klon)
275      SAVE rugoro                 ! longueur de rugosite de l'OESM
276c
277      REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
278c
279      REAL zuthe(klon),zvthe(klon)
280      SAVE zuthe
281      SAVE zvthe
282      INTEGER igwd,idx(klon),itest(klon)
283c
284      REAL agesno(klon,nbsrf)
285      SAVE agesno                 ! age de la neige
286c
287      REAL alb_neig(klon)
288      SAVE alb_neig               ! albedo de la neige
289cKE43
290c Variables liees a la convection de K. Emanuel (sb):
291c
292      REAL ema_workcbmf(klon)   ! cloud base mass flux
293      SAVE ema_workcbmf
294
295      REAL ema_cbmf(klon)       ! cloud base mass flux
296      SAVE ema_cbmf
297
298      REAL ema_pcb(klon)        ! cloud base pressure
299      SAVE ema_pcb
300
301      REAL ema_pct(klon)        ! cloud top pressure
302      SAVE ema_pct
303
304      REAL bas, top             ! cloud base and top levels
305      SAVE bas
306      SAVE top
307
308      REAL Ma(klon,klev)        ! undilute upward mass flux
309      SAVE Ma
310      REAL ema_work1(klon, klev), ema_work2(klon, klev)
311      SAVE ema_work1, ema_work2
312      REAL wdn(klon), tdn(klon), qdn(klon)
313c Variables locales pour la couche limite (al1):
314c
315cAl1      REAL pblh(klon)           ! Hauteur de couche limite
316cAl1      SAVE pblh
317c34EK
318c
319c Variables locales:
320c
321      REAL cdragh(klon) ! drag coefficient pour T and Q
322      REAL cdragm(klon) ! drag coefficient pour vent
323cAA
324cAA  Pour phytrac
325cAA
326      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
327      REAL yu1(klon)            ! vents dans la premiere couche U
328      REAL yv1(klon)            ! vents dans la premiere couche V
329      LOGICAL offline           ! Controle du stockage ds "physique"
330      PARAMETER (offline=.false.)
331      INTEGER physid
332      REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
333      save pfrac_impa
334      REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
335      save pfrac_nucl
336      REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
337      save pfrac_1nucl
338      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
339      REAL frac_nucl(klon,klev) ! idem (nucleation)
340cAA
341      REAL rain_fall(klon) ! pluie
342      REAL snow_fall(klon) ! neige
343      save snow_fall, rain_fall
344      REAL evap(klon), devap(klon) ! evaporation et sa derivee
345      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
346      REAL dlw(klon)    ! derivee infra rouge
347      REAL bils(klon) ! bilan de chaleur au sol
348      REAL fder(klon) ! Derive de flux (sensible et latente)
349      save fder
350      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
351      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
352      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
353      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
354c
355      REAL frugs(klon,nbsrf) ! longueur de rugosite
356      save frugs
357      REAL zxrugs(klon) ! longueur de rugosite
358c
359c Conditions aux limites
360c
361      INTEGER julien
362c
363      INTEGER lmt_pas
364      SAVE lmt_pas                ! frequence de mise a jour
365      REAL pctsrf(klon,nbsrf)
366      SAVE pctsrf                 ! sous-fraction du sol
367      REAL albsol(klon)
368      SAVE albsol                 ! albedo du sol total
369      REAL albsollw(klon)
370      SAVE albsollw                 ! albedo du sol total
371      REAL albsol1(klon)
372      SAVE albsol1                 ! albedo du sol total
373      REAL albsollw1(klon)
374      SAVE albsollw1                 ! albedo du sol total
375
376      REAL wo(klon,klev)
377      SAVE wo                     ! ozone
378c======================================================================
379c
380c Declaration des procedures appelees
381c
382      EXTERNAL angle     ! calculer angle zenithal du soleil
383      EXTERNAL alboc     ! calculer l'albedo sur ocean
384      EXTERNAL albsno    ! calculer albedo sur neige
385      EXTERNAL ajsec     ! ajustement sec
386      EXTERNAL clmain    ! couche limite
387      EXTERNAL condsurf  ! lire les conditions aux limites
388      EXTERNAL conlmd    ! convection (schema LMD)
389cKE43
390      EXTERNAL conema3  ! convect4.3
391      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
392cAA
393      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
394c                          ! stockage des coefficients necessaires au
395c                          ! lessivage OFF-LINE et ON-LINE
396cAA
397      EXTERNAL hgardfou  ! verifier les temperatures
398      EXTERNAL nuage     ! calculer les proprietes radiatives
399      EXTERNAL o3cm      ! initialiser l'ozone
400      EXTERNAL orbite    ! calculer l'orbite terrestre
401      EXTERNAL ozonecm   ! prescrire l'ozone
402      EXTERNAL phyetat0  ! lire l'etat initial de la physique
403      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
404      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
405      EXTERNAL suphec    ! initialiser certaines constantes
406      EXTERNAL transp    ! transport total de l'eau et de l'energie
407      EXTERNAL ecribina  ! ecrire le fichier binaire global
408      EXTERNAL ecribins  ! ecrire le fichier binaire global
409      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
410      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
411c
412c Variables locales
413c
414      real clwcon(klon,klev),rnebcon(klon,klev)
415      real clwcon0(klon,klev),rnebcon0(klon,klev)
416      save rnebcon, clwcon
417
418      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
419      REAL dialiq(klon,klev)  ! eau liquide nuageuse
420      REAL diafra(klon,klev)  ! fraction nuageuse
421      REAL cldliq(klon,klev)  ! eau liquide nuageuse
422      REAL cldfra(klon,klev)  ! fraction nuageuse
423      REAL cldtau(klon,klev)  ! epaisseur optique
424      REAL cldemi(klon,klev)  ! emissivite infrarouge
425c
426C§§§ PB
427      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
428      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
429      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
430      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
431c
432      REAL zxfluxt(klon, klev)
433      REAL zxfluxq(klon, klev)
434      REAL zxfluxu(klon, klev)
435      REAL zxfluxv(klon, klev)
436C§§§
437      REAL heat(klon,klev)    ! chauffage solaire
438      REAL heat0(klon,klev)   ! chauffage solaire ciel clair
439      REAL cool(klon,klev)    ! refroidissement infrarouge
440      REAL cool0(klon,klev)   ! refroidissement infrarouge ciel clair
441      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
442      real sollwdown(klon)    ! downward LW flux at surface
443      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
444      REAL albpla(klon)
445c Le rayonnement n'est pas calcule tous les pas, il faut donc
446c                      sauvegarder les sorties du rayonnement
447      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
448      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
449      INTEGER itaprad
450      SAVE itaprad
451c
452      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
453      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
454c
455      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
456      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
457c
458      REAL zxtsol(klon), zxqsol(klon), zxsnow(klon)
459c
460      REAL dist, rmu0(klon), fract(klon)
461      REAL zdtime, zlongi
462c
463      CHARACTER*2 str2
464      CHARACTER*2 iqn
465c
466      REAL qcheck
467      REAL z_avant(klon), z_apres(klon), z_factor(klon)
468      LOGICAL zx_ajustq
469c
470      REAL za, zb
471      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
472      real zqsat(klon,klev)
473      INTEGER i, k, iq, ig, j, nsrf, ll
474      REAL t_coup
475      PARAMETER (t_coup=234.0)
476c
477      REAL zphi(klon,klev)
478      REAL zx_tmp_x(iim), zx_tmp_yjjmp1
479      REAL zx_relief(iim,jjmp1)
480      REAL zx_aire(iim,jjmp1)
481cKE43
482c Variables locales pour la convection de K. Emanuel (sb):
483c
484      REAL upwd(klon,klev)      ! saturated updraft mass flux
485      REAL dnwd(klon,klev)      ! saturated downdraft mass flux
486      REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
487      REAL tvp(klon,klev)       ! virtual temp of lifted parcel
488      REAL cape(klon)           ! CAPE
489      SAVE cape
490      REAL pbase(klon)          ! cloud base pressure
491      SAVE pbase
492      REAL bbase(klon)          ! cloud base buoyancy
493      SAVE bbase
494      REAL rflag(klon)          ! flag fonctionnement de convect
495      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
496c -- convect43:
497      INTEGER ntra              ! nb traceurs pour convect4.3
498      REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
499      REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
500      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
501      REAL dplcldt(klon), dplcldr(klon)
502c?     .     condm_con(klon,klev),conda_con(klon,klev),
503c?     .     mr_con(klon,klev),ep_con(klon,klev)
504c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
505c --
506c34EK
507c
508c Variables du changement
509c
510c con: convection
511c lsc: condensation a grande echelle (Large-Scale-Condensation)
512c ajs: ajustement sec
513c eva: evaporation de l'eau liquide nuageuse
514c vdf: couche limite (Vertical DiFfusion)
515      REAL d_t_con(klon,klev),d_q_con(klon,klev)
516      REAL d_u_con(klon,klev),d_v_con(klon,klev)
517      REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
518      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
519      REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
520      REAL rneb(klon,klev)
521c
522      REAL pmfu(klon,klev), pmfd(klon,klev)
523      REAL pen_u(klon,klev), pen_d(klon,klev)
524      REAL pde_u(klon,klev), pde_d(klon,klev)
525      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
526      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
527      REAL prfl(klon,klev+1), psfl(klon,klev+1)
528c
529      INTEGER ibas_con(klon), itop_con(klon)
530      REAL rain_con(klon), rain_lsc(klon)
531      REAL snow_con(klon), snow_lsc(klon)
532      REAL d_ts(klon,nbsrf)
533c
534      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
535      REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
536c
537      REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
538      REAL d_t_oro(klon,klev)
539      REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
540      REAL d_t_lif(klon,klev)
541
542      REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
543      real ratqsbas,ratqshaut
544      save ratqsbas,ratqshaut, ratqs
545      real zpt_conv(klon,klev)
546
547c Parametres lies au nouveau schema de nuages (SB, PDF)
548      real fact_cldcon
549      real facttemps
550      logical ok_newmicro
551      save ok_newmicro
552      save fact_cldcon,facttemps
553      real facteur
554
555      integer iflag_cldcon
556      save iflag_cldcon
557
558      logical ptconv(klon,klev)
559
560c
561c Variables liees a l'ecriture de la bande histoire physique
562c
563      INTEGER ecrit_mth
564      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
565c
566      INTEGER ecrit_day
567      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
568c
569      INTEGER ecrit_ins
570      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
571c
572      INTEGER ecrit_reg
573      SAVE ecrit_reg   ! frequence d'ecriture
574c
575      integer itau_w   ! pas de temps ecriture = itap + itau_phy
576c
577c
578c Variables locales pour effectuer les appels en serie
579c
580      REAL t_seri(klon,klev), q_seri(klon,klev)
581      REAL ql_seri(klon,klev),qs_seri(klon,klev)
582      REAL u_seri(klon,klev), v_seri(klon,klev)
583c
584      REAL tr_seri(klon,klev,nbtr)
585      REAL d_tr(klon,klev,nbtr)
586
587      REAL zx_rh(klon,klev)
588
589      INTEGER        length
590      PARAMETER    ( length = 100 )
591      REAL tabcntr0( length       )
592c
593      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
594      REAL zx_tmp_fi2d(klon)
595      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
596      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
597c
598      INTEGER nid_day, nid_mth, nid_ins
599      SAVE nid_day, nid_mth, nid_ins
600c
601      INTEGER nhori, nvert
602      REAL zsto, zout
603      real zjulian
604      save zjulian
605
606      character*20 modname
607      character*80 abort_message
608      logical ok_sync
609      real date0
610      integer idayref
611
612C essai writephys
613      integer fid_day, fid_mth, fid_ins
614      parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
615      integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
616      parameter (prof2d_on = 1, prof3d_on = 2,
617     .           prof2d_av = 3, prof3d_av = 4)
618      character*30 nom_fichier
619      character*10 varname
620      character*40 vartitle
621      character*20 varunits
622C     Variables liees au bilan d'energie et d'enthalpi
623      INTEGER   if_ebil ! level for energy conserv. dignostics
624      SAVE      if_ebil
625      REAL ztsol(klon)
626      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
627     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
628      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
629     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
630      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
631      REAL      d_h_vcol_phy
632      REAL      fs_bound, fq_bound
633      SAVE      d_h_vcol_phy
634      REAL      zero_v(klon)
635      CHARACTER*15 ztit
636      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
637      SAVE      ip_ebil
638      DATA      ip_ebil/2/
639c
640c Declaration des constantes et des fonctions thermodynamiques
641c
642#include "YOMCST.h"
643#include "YOETHF.h"
644#include "FCTTRE.h"
645c======================================================================
646      modname = 'physiq'
647      if_ebil = 2
648      IF (if_ebil.ge.1) THEN
649        DO i=1,klon
650          zero_v(i)=0.
651        END DO
652      END IF
653      ok_sync=.TRUE.
654      IF (nqmax .LT. 2) THEN
655         PRINT*, 'eaux vapeur et liquide sont indispensables'
656         CALL ABORT
657      ENDIF
658      IF (debut) THEN
659         CALL suphec ! initialiser constantes et parametres phys.
660      ENDIF
661
662
663c======================================================================
664      xjour = rjourvrai
665c
666c Si c'est le debut, il faut initialiser plusieurs choses
667c          ********
668c
669       IF (debut) THEN
670C
671         IF (if_ebil.ge.1) d_h_vcol_phy=0.
672c
673c appel a la lecture du run.def physique
674c
675         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
676     .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
677     .                  iflag_cldcon,ratqsbas,ratqshaut)
678
679         DO k = 2, nvm          ! pas de vegetation
680            DO i = 1, klon
681               veget(i,k) = 0.0
682            ENDDO
683         ENDDO
684         DO i = 1, klon
685            veget(i,1) = 1.0    ! il n'y a que du desert
686         ENDDO
687         PRINT*, 'Pas de vegetation; desert partout'
688c
689c
690c Initialiser les compteurs:
691c
692
693         frugs = 0.
694         itap    = 0
695         itaprad = 0
696c
697         CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
698     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
699     .       falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
700     .       dlw,radsol,frugs,agesno,clesphy0,
701     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
702     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
703
704c
705         radpas = NINT( 86400./dtime/nbapp_rad)
706
707c
708         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
709     ,                    ok_instan, ok_region )
710c
711         IF (ABS(dtime-pdtphys).GT.0.001) THEN
712            PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
713            abort_message=' See above '
714            call abort_gcm(modname,abort_message,1)
715         ENDIF
716         IF (nlon .NE. klon) THEN
717            PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
718            abort_message=' See above '
719            call abort_gcm(modname,abort_message,1)
720         ENDIF
721         IF (nlev .NE. klev) THEN
722            PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
723            abort_message=' See above '
724            call abort_gcm(modname,abort_message,1)
725         ENDIF
726c
727         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
728           PRINT*, 'Nbre d appels au rayonnement insuffisant'
729           PRINT*, "Au minimum 4 appels par jour si cycle diurne"
730           abort_message=' See above '
731           call abort_gcm(modname,abort_message,1)
732         ENDIF
733         PRINT*, "Clef pour la convection, iflag_con=", iflag_con
734c
735cKE43
736c Initialisation pour la convection de K.E. (sb):
737         IF (iflag_con.GE.3) THEN
738
739         PRINT*, "*** Convection de Kerry Emanuel 4.3  "
740         PRINT*, "On va utiliser le melange convectif des traceurs qui"
741         PRINT*, "est calcule dans convect4.3"
742         PRINT*, " !!! penser aux logical flags de phytrac"
743
744          DO i = 1, klon
745           ema_cbmf(i) = 0.
746           ema_pcb(i)  = 0.
747           ema_pct(i)  = 0.
748           ema_workcbmf(i) = 0.
749          ENDDO
750         ENDIF
751c34EK
752         IF (ok_orodr) THEN
753         DO i=1,klon
754         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
755         ENDDO
756         CALL SUGWD(klon,klev,paprs,pplay)
757         DO i=1,klon
758         zuthe(i)=0.
759         zvthe(i)=0.
760         if(zstd(i).gt.10.)then
761           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
762           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
763         endif
764         ENDDO
765         ENDIF
766c
767c
768         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
769         PRINT*,'La frequence de lecture surface est de ', lmt_pas
770c
771         ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
772         IF (ok_mensuel) THEN
773         PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
774         ENDIF
775         ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
776         IF (ok_journe) THEN
777         PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
778         ENDIF
779ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
780ccc         ecrit_ins = NINT(86400./dtime *0.25)  ! 4 fois par jour
781         ecrit_ins = NINT(86400./dtime/12.)  ! toutes les deux heures
782         ecrit_ins = NINT(86400./dtime/48.)  ! a chaque pas de temps
783         IF (ok_instan) THEN
784         PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
785         ENDIF
786         ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
787         IF (ok_region) THEN
788         PRINT*, 'La frequence de sortie region est de ', ecrit_reg
789         ENDIF
790
791c
792c Initialiser le couplage si necessaire
793c
794      npas = 0
795      nexca = 0
796      if (ocean == 'couple') then
797        npas = itaufin/ iphysiq
798        nexca = 86400 / dtime
799        write(*,*)' ##### Ocean couple #####'
800        write(*,*)' Valeurs des pas de temps'
801        write(*,*)' npas = ', npas
802        write(*,*)' nexca = ', nexca
803      endif       
804c
805c
806c Gestion calendrier
807
808c
809      IF (ok_journe) THEN
810c
811         idayref = day_ref
812         CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
813c
814         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
815         DO i = 1, iim
816            zx_lon(i,1) = rlon(i+1)
817            zx_lon(i,jjmp1) = rlon(i+1)
818         ENDDO
819         DO ll=1,klev
820            znivsig(ll)=float(ll)
821         ENDDO
822         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
823         write(*,*)'zx_lon = ',zx_lon(:,1)
824         write(*,*)'zx_lat = ',zx_lat(1,:)
825         CALL histbeg("histday", iim,zx_lon(:,1), jjmp1,zx_lat(1,:),
826     .                 1,iim,1,jjmp1, itau_phy, zjulian, dtime,
827     .                 nhori, nid_day)
828         write(*,*)'Journee ', itau_phy, zjulian
829         CALL histvert(nid_day, "presnivs", "Vertical levels", "mb",
830     .                 klev, presnivs, nvert)
831c        call histvert(nid_day, 'sig_s', 'Niveaux sigma','-',
832c    .              klev, znivsig, nvert)
833c
834         zsto = dtime
835         zout = dtime * FLOAT(ecrit_day)
836C Essai writephys
837c        nom_fichier = 'histday1'
838c        call writephy_ini(fid_day,nom_fichier,klon,iim,jjmp1,klev,
839c    .                     rlon,rlat, presnivs,
840c    .                     zjulian, dtime)
841c        call writephy_def(prof2d_on, fid_day, "once", zsto, zout, 0)
842c        call writephy_def(prof3d_on, fid_day, "once", zsto, zout,
843c    .                                                         klev)
844c        call writephy_def(prof2d_av, fid_day, "ave(X)", zsto, zout, 0)
845c        call writephy_def(prof3d_av, fid_day, "ave(X)", zsto, zout,
846c    .                                                         klev)
847 
848c
849         CALL histdef(nid_day, "phis", "Surface geop. height", "-",
850     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
851     .                "once", zsto,zout)
852c
853         CALL histdef(nid_day, "aire", "Grid area", "-",
854     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
855     .                "once", zsto,zout)
856c
857c Champs 2D:
858c
859         CALL histdef(nid_day, "tsol", "Surface Temperature", "K",
860     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
861     .                "ave(X)", zsto,zout)
862c
863         CALL histdef(nid_day, "tter", "Surface Temperature", "K",
864     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
865     .                "ave(X)", zsto,zout)
866c
867         CALL histdef(nid_day, "tlic", "Surface Temperature", "K",
868     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
869     .                "ave(X)", zsto,zout)
870c
871         CALL histdef(nid_day, "toce", "Surface Temperature", "K",
872     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
873     .                "ave(X)", zsto,zout)
874c
875         CALL histdef(nid_day, "tsic", "Surface Temperature", "K",
876     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
877     .                "ave(X)", zsto,zout)
878c
879         CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
880     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
881     .                "ave(X)", zsto,zout)
882c
883         CALL histdef(nid_day, "rain", "Precipitation", "mm/day",
884     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
885     .                "ave(X)", zsto,zout)
886c
887         CALL histdef(nid_day, "snow", "Snow fall", "mm/day",
888     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
889     .                "ave(X)", zsto,zout)
890c
891         CALL histdef(nid_day, "snow_cov", "Snow cover", "mm",
892     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
893     .                "ave(X)", zsto,zout)
894c
895         CALL histdef(nid_day, "evap", "Evaporation", "mm/day",
896     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
897     .                "ave(X)", zsto,zout)
898c
899         CALL histdef(nid_day, "tops", "Solar rad. at TOA", "W/m2",
900     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
901     .                "ave(X)", zsto,zout)
902c
903         CALL histdef(nid_day, "topl", "IR rad. at TOA", "W/m2",
904     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
905     .                "ave(X)", zsto,zout)
906c
907         CALL histdef(nid_day, "sols", "Solar rad. at surf.", "W/m2",
908     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
909     .                "ave(X)", zsto,zout)
910c
911         CALL histdef(nid_day, "soll", "IR rad. at surface", "W/m2",
912     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
913     .                "ave(X)", zsto,zout)
914c
915         CALL histdef(nid_day, "solldown", "Down. IR rad. at surface",
916     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32,
917     .                "ave(X)", zsto,zout)
918c
919         CALL histdef(nid_day, "bils", "Surf. total heat flux", "W/m2",
920     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
921     .                "ave(X)", zsto,zout)
922c
923         CALL histdef(nid_day, "sens", "Sensible heat flux", "W/m2",
924     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
925     .                "ave(X)", zsto,zout)
926c
927         CALL histdef(nid_day, "fder", "Heat flux derivation", "W/m2",
928     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
929     .                "ave(X)", zsto,zout)
930c
931         CALL histdef(nid_day, "frtu", "Zonal wind stress", "Pa",
932     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
933     .                "ave(X)", zsto,zout)
934c
935         CALL histdef(nid_day, "frtv", "Meridional wind stress", "Pa",
936     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
937     .                "ave(X)", zsto,zout)
938c
939C §§§ PB flux pour chauqe sous surface
940C
941         DO nsrf = 1, nbsrf
942C
943           call histdef(nid_day, "pourc_"//clnsurf(nsrf),
944     $         "Fraction"//clnsurf(nsrf), "W/m2", 
945     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
946     $         "ave(X)", zsto,zout)
947C
948           call histdef(nid_day, "tsol_"//clnsurf(nsrf),
949     $         "Fraction"//clnsurf(nsrf), "W/m2", 
950     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
951     $         "ave(X)", zsto,zout)
952C
953           call histdef(nid_day, "sens_"//clnsurf(nsrf),
954     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
955     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
956     $         "ave(X)", zsto,zout)
957c
958           call histdef(nid_day, "lat_"//clnsurf(nsrf),
959     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
960     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
961     $         "ave(X)", zsto,zout)
962C
963           call histdef(nid_day, "taux_"//clnsurf(nsrf),
964     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
965     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
966     $         "ave(X)", zsto,zout)
967
968           call histdef(nid_day, "tauy_"//clnsurf(nsrf),
969     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
970     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
971     $         "ave(X)", zsto,zout)
972C
973           call histdef(nid_day, "albe_"//clnsurf(nsrf),
974     $         "Albedo surf. "//clnsurf(nsrf), "W/m2", 
975     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
976     $         "ave(X)", zsto,zout)
977C
978           call histdef(nid_day, "rugs_"//clnsurf(nsrf),
979     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
980     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
981     $         "ave(X)", zsto,zout)
982
983C§§§
984         END DO
985           
986         CALL histdef(nid_day, "sicf", "Sea-ice fraction", "-",
987     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
988     .                "ave(X)", zsto,zout)
989c
990         CALL histdef(nid_day, "cldl", "Low-level cloudiness", "-",
991     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
992     .                "ave(X)", zsto,zout)
993c
994         CALL histdef(nid_day, "cldm", "Mid-level cloudiness", "-",
995     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
996     .                "ave(X)", zsto,zout)
997c
998         CALL histdef(nid_day, "cldh", "High-level cloudiness", "-",
999     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1000     .                "ave(X)", zsto,zout)
1001c
1002         CALL histdef(nid_day, "cldt", "Total cloudiness", "-",
1003     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1004     .                "ave(X)", zsto,zout)
1005c
1006         CALL histdef(nid_day, "cldq", "Cloud liquid water path", "-",
1007     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1008     .                "ave(X)", zsto,zout)
1009c
1010c Champs 3D:
1011c
1012         CALL histdef(nid_day, "temp", "Air temperature", "K",
1013     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1014     .                "ave(X)", zsto,zout)
1015c
1016         CALL histdef(nid_day, "ovap", "Specific humidity", "Kg/Kg",
1017     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1018     .                "ave(X)", zsto,zout)
1019c
1020         CALL histdef(nid_day, "geop", "Geopotential height", "m",
1021     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1022     .                "ave(X)", zsto,zout)
1023c
1024         CALL histdef(nid_day, "vitu", "Zonal wind", "m/s",
1025     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1026     .                "ave(X)", zsto,zout)
1027c
1028         CALL histdef(nid_day, "vitv", "Meridional wind", "m/s",
1029     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1030     .                "ave(X)", zsto,zout)
1031c
1032         CALL histdef(nid_day, "vitw", "Vertical wind", "m/s",
1033     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1034     .                "ave(X)", zsto,zout)
1035c
1036         CALL histdef(nid_day, "pres", "Air pressure", "Pa",
1037     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1038     .                "ave(X)", zsto,zout)
1039c
1040         CALL histend(nid_day)
1041c
1042         ndex2d = 0
1043         ndex3d = 0
1044c
1045      ENDIF ! fin de test sur ok_journe
1046c
1047      IF (ok_mensuel) THEN
1048c
1049         idayref = day_ref
1050         CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
1051c
1052         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
1053         DO i = 1, iim
1054            zx_lon(i,1) = rlon(i+1)
1055            zx_lon(i,jjmp1) = rlon(i+1)
1056         ENDDO
1057         DO ll=1,klev
1058            znivsig(ll)=float(ll)
1059         ENDDO
1060         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
1061         CALL histbeg("histmth.nc", iim,zx_lon(:,1), jjmp1,zx_lat(1,:),
1062     .                 1,iim,1,jjmp1, itau_phy, zjulian, dtime,
1063     .                 nhori, nid_mth)
1064         write(*,*)'Mensuel ', itau_phy, zjulian
1065         CALL histvert(nid_mth, "presnivs", "Vertical levels", "mb",
1066     .                 klev, presnivs, nvert)
1067c        call histvert(nid_mth, 'sig_s', 'Niveaux sigma','-',
1068c    .              klev, znivsig, nvert)
1069c
1070         zsto = dtime
1071         zout = dtime * ecrit_mth
1072c
1073         CALL histdef(nid_mth, "phis", "Surface geop. height", "-",
1074     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1075     .                "once",  zsto,zout)
1076c
1077         CALL histdef(nid_mth, "aire", "Grid area", "-",
1078     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1079     .                "once",  zsto,zout)
1080c
1081c Champs 2D:
1082c
1083         CALL histdef(nid_mth, "tsol", "Surface Temperature", "K",
1084     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1085     .                "ave(X)", zsto,zout)
1086c
1087         CALL histdef(nid_mth, "psol", "Surface Pressure", "Pa",
1088     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1089     .                "ave(X)", zsto,zout)
1090c
1091         CALL histdef(nid_mth, "qsol", "Surface humidity", "mm",
1092     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1093     .                "ave(X)", zsto,zout)
1094c
1095         CALL histdef(nid_mth, "rain", "Precipitation", "mm/day",
1096     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1097     .                "ave(X)", zsto,zout)
1098c
1099         CALL histdef(nid_mth, "plul", "Large-scale Precip.", "mm/day",
1100     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1101     .                "ave(X)", zsto,zout)
1102c
1103         CALL histdef(nid_mth, "pluc", "Convective Precip.", "mm/day",
1104     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1105     .                "ave(X)", zsto,zout)
1106c
1107         CALL histdef(nid_mth, "snow", "Snow fall", "mm/day",
1108     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1109     .                "ave(X)", zsto,zout)
1110c
1111         CALL histdef(nid_mth, "snow_cov", "Snow cover", "mm",
1112     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1113     .                "ave(X)", zsto,zout)
1114c
1115         CALL histdef(nid_mth, "evap", "Evaporation", "mm/day",
1116     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1117     .                "ave(X)", zsto,zout)
1118c
1119         CALL histdef(nid_mth, "tops", "Solar rad. at TOA", "W/m2",
1120     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1121     .                "ave(X)", zsto,zout)
1122c
1123         CALL histdef(nid_mth, "topl", "IR rad. at TOA", "W/m2",
1124     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1125     .                "ave(X)", zsto,zout)
1126c
1127         CALL histdef(nid_mth, "sols", "Solar rad. at surf.", "W/m2",
1128     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1129     .                "ave(X)", zsto,zout)
1130c
1131         CALL histdef(nid_mth, "soll", "IR rad. at surface", "W/m2",
1132     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1133     .                "ave(X)", zsto,zout)
1134c
1135         CALL histdef(nid_mth, "solldown", "Down. IR rad. at surface",
1136     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32,
1137     .                "ave(X)", zsto,zout)
1138c
1139         CALL histdef(nid_mth, "tops0", "Solar rad. at TOA", "W/m2",
1140     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1141     .                "ave(X)", zsto,zout)
1142c
1143         CALL histdef(nid_mth, "topl0", "IR rad. at TOA", "W/m2",
1144     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1145     .                "ave(X)", zsto,zout)
1146c
1147         CALL histdef(nid_mth, "sols0", "Solar rad. at surf.", "W/m2",
1148     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1149     .                "ave(X)", zsto,zout)
1150c
1151         CALL histdef(nid_mth, "soll0", "IR rad. at surface", "W/m2",
1152     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1153     .                "ave(X)", zsto,zout)
1154c
1155         CALL histdef(nid_mth, "bils", "Surf. total heat flux", "W/m2",
1156     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1157     .                "ave(X)", zsto,zout)
1158c
1159         CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
1160     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1161     .                "ave(X)", zsto,zout)
1162c
1163         CALL histdef(nid_mth, "fder", "Heat flux derivation", "W/m2",
1164     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1165     .                "ave(X)", zsto,zout)
1166c
1167         CALL histdef(nid_mth, "frtu", "Zonal wind stress", "Pa",
1168     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1169     .                "ave(X)", zsto,zout)
1170c
1171         CALL histdef(nid_mth, "frtv", "Meridional wind stress", "Pa",
1172     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1173     .                "ave(X)", zsto,zout)
1174c
1175         DO nsrf = 1, nbsrf
1176C
1177           call histdef(nid_mth, "pourc_"//clnsurf(nsrf),
1178     $         "Fraction "//clnsurf(nsrf), "W/m2", 
1179     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1180     $         "ave(X)", zsto,zout)
1181C
1182           call histdef(nid_mth, "tsol_"//clnsurf(nsrf),
1183     $         "Fraction "//clnsurf(nsrf), "W/m2", 
1184     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1185     $         "ave(X)", zsto,zout)
1186C
1187           call histdef(nid_mth, "sens_"//clnsurf(nsrf),
1188     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
1189     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1190     $         "ave(X)", zsto,zout)
1191c
1192           call histdef(nid_mth, "lat_"//clnsurf(nsrf),
1193     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
1194     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1195     $         "ave(X)", zsto,zout)
1196C
1197           call histdef(nid_mth, "taux_"//clnsurf(nsrf),
1198     $         "Zonal wind stress"//clnsurf(nsrf), "Pa", 
1199     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1200     $         "ave(X)", zsto,zout)
1201
1202           call histdef(nid_mth, "tauy_"//clnsurf(nsrf),
1203     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
1204     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1205     $         "ave(X)", zsto,zout)
1206c
1207           call histdef(nid_mth, "albe_"//clnsurf(nsrf),
1208     $         "Albedo surf. "//clnsurf(nsrf), "W/m2", 
1209     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1210     $         "ave(X)", zsto,zout)
1211c
1212           call histdef(nid_mth, "rugs_"//clnsurf(nsrf),
1213     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
1214     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1215     $         "ave(X)", zsto,zout)
1216c
1217         CALL histdef(nid_mth, "ages_"//clnsurf(nsrf), "Snow age","day",
1218     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1219     .                "ave(X)", zsto,zout)
1220
1221         END DO
1222C
1223         CALL histdef(nid_mth, "sicf", "Sea-ice fraction", "-",
1224     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1225     .                "ave(X)", zsto,zout)
1226c
1227         CALL histdef(nid_mth, "albs", "Surface albedo", "-",
1228     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1229     .                "ave(X)", zsto,zout)
1230         CALL histdef(nid_mth, "albslw", "Surface albedo LW", "-",
1231     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1232     .                "ave(X)", zsto,zout)
1233c
1234         CALL histdef(nid_mth, "cdrm", "Momentum drag coef.", "-",
1235     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1236     .                "ave(X)", zsto,zout)
1237c
1238         CALL histdef(nid_mth, "cdrh", "Heat drag coef.", "-",
1239     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1240     .                "ave(X)", zsto,zout)
1241c
1242         CALL histdef(nid_mth, "cldl", "Low-level cloudiness", "-",
1243     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1244     .                "ave(X)", zsto,zout)
1245c
1246         CALL histdef(nid_mth, "cldm", "Mid-level cloudiness", "-",
1247     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1248     .                "ave(X)", zsto,zout)
1249c
1250         CALL histdef(nid_mth, "cldh", "High-level cloudiness", "-",
1251     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1252     .                "ave(X)", zsto,zout)
1253c
1254         CALL histdef(nid_mth, "cldt", "Total cloudiness", "-",
1255     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1256     .                "ave(X)", zsto,zout)
1257c
1258         CALL histdef(nid_mth, "cldq", "Cloud liquid water path", "-",
1259     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1260     .                "ave(X)", zsto,zout)
1261c
1262         CALL histdef(nid_mth, "ue", "Zonal energy transport", "-",
1263     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1264     .                "ave(X)", zsto,zout)
1265c
1266         CALL histdef(nid_mth, "ve", "Merid energy transport", "-",
1267     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1268     .                "ave(X)", zsto,zout)
1269c
1270         CALL histdef(nid_mth, "uq", "Zonal humidity transport", "-",
1271     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1272     .                "ave(X)", zsto,zout)
1273c
1274         CALL histdef(nid_mth, "vq", "Merid humidity transport", "-",
1275     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1276     .                "ave(X)", zsto,zout)
1277cKE43
1278      IF (iflag_con .GE. 3) THEN ! sb
1279c
1280         CALL histdef(nid_mth, "cape", "Conv avlbl pot ener", "J/Kg",
1281     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1282     .                "ave(X)", zsto,zout)
1283c
1284         CALL histdef(nid_mth, "pbase", "Cld base pressure", "hPa",
1285     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1286     .                "ave(X)", zsto,zout)
1287c
1288         CALL histdef(nid_mth, "ptop", "Cld top pressure", "hPa",
1289     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1290     .                "ave(X)", zsto,zout)
1291c
1292         CALL histdef(nid_mth, "fbase", "Cld base mass flux", "Kg/m2/s",
1293     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1294     .                "ave(X)", zsto,zout)
1295c
1296c
1297      ENDIF
1298c34EK
1299c
1300c Champs 3D:
1301c
1302         CALL histdef(nid_mth, "temp", "Air temperature", "K",
1303     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1304     .                "ave(X)", zsto,zout)
1305c
1306         CALL histdef(nid_mth, "ovap", "Specific humidity", "Kg/Kg",
1307     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1308     .                "ave(X)", zsto,zout)
1309c
1310         CALL histdef(nid_mth, "geop", "Geopotential height", "m",
1311     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1312     .                "ave(X)", zsto,zout)
1313c
1314         CALL histdef(nid_mth, "vitu", "Zonal wind", "m/s",
1315     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1316     .                "ave(X)", zsto,zout)
1317c
1318         CALL histdef(nid_mth, "vitv", "Meridional wind", "m/s",
1319     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1320     .                "ave(X)", zsto,zout)
1321c
1322         CALL histdef(nid_mth, "vitw", "Vertical wind", "m/s",
1323     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1324     .                "ave(X)", zsto,zout)
1325c
1326         CALL histdef(nid_mth, "pres", "Air pressure", "Pa",
1327     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1328     .                "ave(X)", zsto,zout)
1329c
1330         CALL histdef(nid_mth, "rneb", "Cloud fraction", "-",
1331     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1332     .                "ave(X)", zsto,zout)
1333c
1334         CALL histdef(nid_mth, "rhum", "Relative humidity", "-",
1335     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1336     .                "ave(X)", zsto,zout)
1337c
1338         CALL histdef(nid_mth, "oliq", "Liquid water content", "kg/kg",
1339     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1340     .                "ave(X)", zsto,zout)
1341c
1342         CALL histdef(nid_mth, "dtdyn", "Dynamics dT", "K/s",
1343     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1344     .                "ave(X)", zsto,zout)
1345c
1346         CALL histdef(nid_mth, "dqdyn", "Dynamics dQ", "Kg/Kg/s",
1347     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1348     .                "ave(X)", zsto,zout)
1349c
1350         CALL histdef(nid_mth, "dtcon", "Convection dT", "K/s",
1351     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1352     .                "ave(X)", zsto,zout)
1353c
1354         CALL histdef(nid_mth, "ducon", "Convection du", "m/s2",
1355     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1356     .                "ave(X)", zsto,zout)
1357c
1358         CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
1359     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1360     .                "ave(X)", zsto,zout)
1361c
1362         CALL histdef(nid_mth, "dtlsc", "Condensation dT", "K/s",
1363     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1364     .                "ave(X)", zsto,zout)
1365c
1366         CALL histdef(nid_mth, "dqlsc", "Condensation dQ", "Kg/Kg/s",
1367     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1368     .                "ave(X)", zsto,zout)
1369c
1370         CALL histdef(nid_mth, "dtvdf", "Boundary-layer dT", "K/s",
1371     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1372     .                "ave(X)", zsto,zout)
1373c
1374         CALL histdef(nid_mth, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
1375     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1376     .                "ave(X)", zsto,zout)
1377c
1378         CALL histdef(nid_mth, "dteva", "Reevaporation dT", "K/s",
1379     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1380     .                "ave(X)", zsto,zout)
1381c
1382         CALL histdef(nid_mth, "dqeva", "Reevaporation dQ", "Kg/Kg/s",
1383     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1384     .                "ave(X)", zsto,zout)
1385
1386         CALL histdef(nid_mth, "ptconv", "POINTS CONVECTIFS"," ",
1387     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1388     .                "ave(X)", zsto,zout)
1389
1390         CALL histdef(nid_mth, "ratqs", "RATQS"," ",
1391     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1392     .                "ave(X)", zsto,zout)
1393
1394c
1395         CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
1396     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1397     .                "ave(X)", zsto,zout)
1398
1399         CALL histdef(nid_mth, "dqajs", "Dry adjust. dQ", "Kg/Kg/s",
1400     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1401     .                "ave(X)", zsto,zout)
1402c
1403         CALL histdef(nid_mth, "dtswr", "SW radiation dT", "K/s",
1404     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1405     .                "ave(X)", zsto,zout)
1406c
1407         CALL histdef(nid_mth, "dtsw0", "SW radiation dT", "K/s",
1408     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1409     .                "ave(X)", zsto,zout)
1410c
1411         CALL histdef(nid_mth, "dtlwr", "LW radiation dT", "K/s",
1412     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1413     .                "ave(X)", zsto,zout)
1414c
1415         CALL histdef(nid_mth, "dtlw0", "LW radiation dT", "K/s",
1416     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1417     .                "ave(X)", zsto,zout)
1418c
1419         CALL histdef(nid_mth, "duvdf", "Boundary-layer dU", "m/s2",
1420     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1421     .                "ave(X)", zsto,zout)
1422c
1423         CALL histdef(nid_mth, "dvvdf", "Boundary-layer dV", "m/s2",
1424     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1425     .                "ave(X)", zsto,zout)
1426c
1427         IF (ok_orodr) THEN
1428         CALL histdef(nid_mth, "duoro", "Orography dU", "m/s2",
1429     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1430     .                "ave(X)", zsto,zout)
1431c
1432         CALL histdef(nid_mth, "dvoro", "Orography dV", "m/s2",
1433     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1434     .                "ave(X)", zsto,zout)
1435c
1436         ENDIF
1437C
1438         IF (ok_orolf) THEN
1439         CALL histdef(nid_mth, "dulif", "Orography dU", "m/s2",
1440     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1441     .                "ave(X)", zsto,zout)
1442c
1443         CALL histdef(nid_mth, "dvlif", "Orography dV", "m/s2",
1444     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1445     .                "ave(X)", zsto,zout)
1446         ENDIF
1447C
1448         CALL histdef(nid_mth, "ozone", "Ozone concentration", "-",
1449     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1450     .                "ave(X)", zsto,zout)
1451c
1452         if (nqmax.GE.3) THEN
1453         DO iq=1,nqmax-2
1454         IF (iq.LE.99) THEN
1455         WRITE(str2,'(i2.2)') iq
1456         CALL histdef(nid_mth, "trac"//str2, "Tracer No."//str2, "-",
1457     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1458     .                "ave(X)", zsto,zout)
1459         ELSE
1460         PRINT*, "Trop de traceurs"
1461         CALL abort
1462         ENDIF
1463         ENDDO
1464         ENDIF
1465c
1466cKE43
1467      IF (iflag_con.GE.3) THEN ! (sb)
1468c
1469         CALL histdef(nid_mth, "upwd", "saturated updraft", "Kg/m2/s",
1470     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1471     .                "ave(X)", zsto,zout)
1472c
1473         CALL histdef(nid_mth, "dnwd", "saturated downdraft","Kg/m2/s",
1474     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1475     .                "ave(X)", zsto,zout)
1476c
1477         CALL histdef(nid_mth, "dnwd0", "unsat. downdraft", "Kg/m2/s",
1478     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1479     .                "ave(X)", zsto,zout)
1480c
1481         CALL histdef(nid_mth,"Ma","undilute adiab updraft","Kg/m2/s",
1482     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1483     .                "ave(X)", zsto,zout)
1484c
1485c
1486      ENDIF
1487c34EK
1488         CALL histend(nid_mth)
1489c
1490         ndex2d = 0
1491         ndex3d = 0
1492c
1493      ENDIF ! fin de test sur ok_mensuel
1494c
1495c
1496      IF (ok_instan) THEN
1497c
1498         idayref = day_ref
1499         CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
1500c
1501         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
1502         DO i = 1, iim
1503            zx_lon(i,1) = rlon(i+1)
1504            zx_lon(i,jjmp1) = rlon(i+1)
1505         ENDDO
1506         DO ll=1,klev
1507            znivsig(ll)=float(ll)
1508         ENDDO
1509         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
1510         CALL histbeg("histins", iim,zx_lon(:,1), jjmp1,zx_lat(1,:),
1511     .                 1,iim,1,jjmp1, itau_phy, zjulian, dtime,
1512     .                 nhori, nid_ins)
1513         write(*,*)'Inst ', itau_phy, zjulian
1514         CALL histvert(nid_ins, "presnivs", "Vertical levels", "mb",
1515     .                 klev, presnivs, nvert)
1516c        call histvert(nid_ins, 'sig_s', 'Niveaux sigma','-',
1517c    .              klev, znivsig, nvert)
1518c
1519c
1520         zsto = dtime * ecrit_ins
1521         zout = dtime * ecrit_ins
1522C
1523         CALL histdef(nid_ins, "phis", "Surface geop. height", "-",
1524     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1525     .                "once", zsto,zout)
1526c
1527         CALL histdef(nid_ins, "aire", "Grid area", "-",
1528     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1529     .                "once", zsto,zout)
1530c
1531c Champs 2D:
1532c
1533        CALL histdef(nid_ins, "tsol", "Surface Temperature", "K",
1534     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1535     .                "inst(X)", zsto,zout)
1536c
1537        CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
1538     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1539     .                "inst(X)", zsto,zout)
1540c
1541         CALL histdef(nid_ins, "plul", "Large-scale Precip.", "mm/day",
1542     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1543     .                "inst(X)", zsto,zout)
1544c
1545         CALL histdef(nid_ins, "pluc", "Convective Precip.", "mm/day",
1546     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1547     .                "inst(X)", zsto,zout)
1548
1549        CALL histdef(nid_ins, "qsol", "Surface humidity", "mm",
1550     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1551     .                "inst(X)", zsto,zout)
1552c
1553         CALL histdef(nid_ins, "cdrm", "Momentum drag coef.", "-",
1554     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1555     .                "inst(X)", zsto,zout)
1556c
1557         CALL histdef(nid_ins, "cdrh", "Heat drag coef.", "-",
1558     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1559     .                "inst(X)", zsto,zout)
1560c
1561         CALL histdef(nid_ins, "rain", "Precipitation", "mm/day",
1562     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1563     .                "inst(X)", zsto,zout)
1564c
1565         CALL histdef(nid_ins, "snow", "Snow fall", "mm/day",
1566     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1567     .                "inst(X)", zsto,zout)
1568c
1569         CALL histdef(nid_ins, "snow_cov", "Snow cover", "mm",
1570     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1571     .                "inst(X)", zsto,zout)
1572c
1573         CALL histdef(nid_ins, "topl", "OLR", "W/m2",
1574     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1575     .                "inst(X)", zsto,zout)
1576c
1577         CALL histdef(nid_ins, "evap", "Evaporation", "mm/day",
1578     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1579     .                "inst(X)", zsto,zout)
1580c
1581         CALL histdef(nid_ins, "sols", "Solar rad. at surf.", "W/m2",
1582     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1583     .                "inst(X)", zsto,zout)
1584c
1585         CALL histdef(nid_ins, "soll", "IR rad. at surface", "W/m2",
1586     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1587     .                "inst(X)", zsto,zout)
1588c
1589         CALL histdef(nid_ins, "solldown", "Down. IR rad. at surface",
1590     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32,
1591     .                "inst(X)", zsto,zout)
1592c
1593         CALL histdef(nid_ins, "bils", "Surf. total heat flux", "W/m2",
1594     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1595     .                "inst(X)", zsto,zout)
1596c
1597         CALL histdef(nid_ins, "sens", "Sensible heat flux", "W/m2",
1598     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1599     .                "inst(X)", zsto,zout)
1600c
1601         CALL histdef(nid_ins, "fder", "Heat flux derivation", "W/m2",
1602     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1603     .                "inst(X)", zsto,zout)
1604c
1605      CALL histdef(nid_ins, "dtsvdfo", "Boundary-layer dTs(o)", "K/s",
1606     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1607     .                "inst(X)", zsto,zout)
1608c
1609      CALL histdef(nid_ins, "dtsvdft", "Boundary-layer dTs(t)", "K/s",
1610     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1611     .                "inst(X)", zsto,zout)
1612c
1613      CALL histdef(nid_ins, "dtsvdfg", "Boundary-layer dTs(g)", "K/s",
1614     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1615     .                "inst(X)", zsto,zout)
1616c
1617      CALL histdef(nid_ins, "dtsvdfi", "Boundary-layer dTs(g)", "K/s",
1618     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1619     .                "inst(X)", zsto,zout)
1620
1621         DO nsrf = 1, nbsrf
1622C
1623           call histdef(nid_ins, "pourc_"//clnsurf(nsrf),
1624     $         "Fraction"//clnsurf(nsrf), "W/m2", 
1625     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1626     $         "inst(X)", zsto,zout)
1627
1628           call histdef(nid_ins, "sens_"//clnsurf(nsrf),
1629     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
1630     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1631     $         "inst(X)", zsto,zout)
1632c
1633           call histdef(nid_ins, "tsol_"//clnsurf(nsrf),
1634     $         "Surface Temperature"//clnsurf(nsrf), "W/m2", 
1635     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1636     $         "inst(X)", zsto,zout)
1637c
1638           call histdef(nid_ins, "lat_"//clnsurf(nsrf),
1639     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
1640     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1641     $         "inst(X)", zsto,zout)
1642C
1643           call histdef(nid_ins, "taux_"//clnsurf(nsrf),
1644     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
1645     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1646     $         "inst(X)", zsto,zout)
1647
1648           call histdef(nid_ins, "tauy_"//clnsurf(nsrf),
1649     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
1650     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1651     $         "inst(X)", zsto,zout)
1652c
1653           call histdef(nid_ins, "albe_"//clnsurf(nsrf),
1654     $         "Albedo "//clnsurf(nsrf), "-", 
1655     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1656     $         "inst(X)", zsto,zout)
1657c
1658           call histdef(nid_ins, "rugs_"//clnsurf(nsrf),
1659     $         "rugosite "//clnsurf(nsrf), "-", 
1660     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1661     $         "inst(X)", zsto,zout)
1662C§§§
1663         END DO
1664         CALL histdef(nid_ins, "rugs", "rugosity", "-",
1665     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1666     .                "inst(X)", zsto,zout)
1667
1668c
1669         CALL histdef(nid_ins, "albs", "Surface albedo", "-",
1670     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1671     .                "inst(X)", zsto,zout)
1672         CALL histdef(nid_ins, "albslw", "Surface albedo LW", "-",
1673     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1674     .                "inst(X)", zsto,zout)
1675c
1676c
1677c Champs 3D:
1678c
1679         CALL histdef(nid_ins, "temp", "Temperature", "K",
1680     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1681     .                "inst(X)", zsto,zout)
1682c
1683         CALL histdef(nid_ins, "vitu", "Zonal wind", "m/s",
1684     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1685     .                "inst(X)", zsto,zout)
1686c
1687         CALL histdef(nid_ins, "vitv", "Merid wind", "m/s",
1688     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1689     .                "inst(X)", zsto,zout)
1690c
1691         CALL histdef(nid_ins, "geop", "Geopotential height", "m",
1692     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1693     .                "inst(X)", zsto,zout)
1694c
1695         CALL histdef(nid_ins, "pres", "Air pressure", "Pa",
1696     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1697     .                "inst(X)", zsto,zout)
1698c
1699         CALL histdef(nid_ins, "dtvdf", "Boundary-layer dT", "K/s",
1700     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1701     .                "inst(X)", zsto,zout)
1702c
1703         CALL histdef(nid_ins, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
1704     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1705     .                "inst(X)", zsto,zout)
1706c
1707
1708         CALL histend(nid_ins)
1709c
1710         ndex2d = 0
1711         ndex3d = 0
1712c
1713      ENDIF
1714
1715c$$$PB Positionner date0 pour initialisation de ORCHIDEE
1716      date0 = zjulian
1717C      date0 = day_ini
1718      WRITE(*,*) 'physiq date0 : ',date0
1719c
1720c
1721c
1722c Prescrire l'ozone dans l'atmosphere
1723c
1724c
1725cc         DO i = 1, klon
1726cc         DO k = 1, klev
1727cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1728cc         ENDDO
1729cc         ENDDO
1730c
1731c
1732      ENDIF
1733c
1734c   ****************     Fin  de   IF ( debut  )   ***************
1735c
1736c
1737c Mettre a zero des variables de sortie (pour securite)
1738c
1739      DO i = 1, klon
1740         d_ps(i) = 0.0
1741      ENDDO
1742      DO k = 1, klev
1743      DO i = 1, klon
1744         d_t(i,k) = 0.0
1745         d_u(i,k) = 0.0
1746         d_v(i,k) = 0.0
1747      ENDDO
1748      ENDDO
1749      DO iq = 1, nqmax
1750      DO k = 1, klev
1751      DO i = 1, klon
1752         d_qx(i,k,iq) = 0.0
1753      ENDDO
1754      ENDDO
1755      ENDDO
1756c
1757c Ne pas affecter les valeurs entrees de u, v, h, et q
1758c
1759      DO k = 1, klev
1760      DO i = 1, klon
1761         t_seri(i,k)  = t(i,k)
1762         u_seri(i,k)  = u(i,k)
1763         v_seri(i,k)  = v(i,k)
1764         q_seri(i,k)  = qx(i,k,ivap)
1765         ql_seri(i,k) = qx(i,k,iliq)
1766         qs_seri(i,k) = 0.
1767      ENDDO
1768      ENDDO
1769      IF (nqmax.GE.3) THEN
1770      DO iq = 3, nqmax
1771      DO  k = 1, klev
1772      DO  i = 1, klon
1773         tr_seri(i,k,iq-2) = qx(i,k,iq)
1774      ENDDO
1775      ENDDO
1776      ENDDO
1777      ELSE
1778      DO k = 1, klev
1779      DO i = 1, klon
1780         tr_seri(i,k,1) = 0.0
1781      ENDDO
1782      ENDDO
1783      ENDIF
1784C
1785      IF (if_ebil.ge.1) THEN
1786        DO i = 1, klon
1787          ztsol(i) = 0.
1788        ENDDO
1789        DO i = 1, klon
1790          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1791        ENDDO
1792        ztit='after dynamic'
1793        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
1794     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1795     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1796C     Comme les tendances de la physique sont ajoute dans la dynamique,
1797C     on devrait avoir que la variation d'entalpie par la dynamique
1798C     est egale a la variation de la physique au pas de temps precedent.
1799C     Donc la somme de ces 2 variations devrait etre nulle.
1800        call diagphy(paire,ztit,ip_ebil
1801     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1802     e      , zero_v, zero_v, zero_v, ztsol
1803     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1804     s      , fs_bound, fq_bound )
1805      END IF
1806
1807c Diagnostiquer la tendance dynamique
1808c
1809      IF (ancien_ok) THEN
1810         DO k = 1, klev
1811         DO i = 1, klon
1812            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1813            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1814         ENDDO
1815         ENDDO
1816      ELSE
1817         DO k = 1, klev
1818         DO i = 1, klon
1819            d_t_dyn(i,k) = 0.0
1820            d_q_dyn(i,k) = 0.0
1821         ENDDO
1822         ENDDO
1823         ancien_ok = .TRUE.
1824      ENDIF
1825c
1826c Ajouter le geopotentiel du sol:
1827c
1828      DO k = 1, klev
1829      DO i = 1, klon
1830         zphi(i,k) = pphi(i,k) + pphis(i)
1831      ENDDO
1832      ENDDO
1833c
1834c Verifier les temperatures
1835c
1836      CALL hgardfou(t_seri,ftsol,'debutphy')
1837c
1838c Incrementer le compteur de la physique
1839c
1840      itap   = itap + 1
1841      julien = MOD(NINT(xjour),360)
1842c
1843c Mettre en action les conditions aux limites (albedo, sst, etc.).
1844c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1845c
1846      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1847         PRINT *,' PHYS cond  julien ',julien
1848         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1849      ENDIF
1850c
1851c Re-evaporer l'eau liquide nuageuse
1852c
1853      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1854      DO i = 1, klon
1855         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1856c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1857         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1858         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1859         zb = MAX(0.0,ql_seri(i,k))
1860         za = - MAX(0.0,ql_seri(i,k))
1861     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1862         t_seri(i,k) = t_seri(i,k) + za
1863         q_seri(i,k) = q_seri(i,k) + zb
1864         ql_seri(i,k) = 0.0
1865         d_t_eva(i,k) = za
1866         d_q_eva(i,k) = zb
1867      ENDDO
1868      ENDDO
1869c
1870      IF (if_ebil.ge.2) THEN
1871        ztit='after reevap'
1872        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1873     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1874     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1875         call diagphy(paire,ztit,ip_ebil
1876     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1877     e      , zero_v, zero_v, zero_v, ztsol
1878     e      , d_h_vcol, d_qt, d_ec
1879     s      , fs_bound, fq_bound )
1880C
1881      END IF
1882C
1883c
1884c Appeler la diffusion verticale (programme de couche limite)
1885c
1886      DO i = 1, klon
1887c       if (.not. ok_veget) then
1888c          frugs(i,is_ter) = SQRT(frugs(i,is_ter)**2+rugoro(i)**2)
1889c       endif
1890c         frugs(i,is_lic) = rugoro(i)
1891c         frugs(i,is_oce) = rugmer(i)
1892c         frugs(i,is_sic) = 0.001
1893         zxrugs(i) = 0.0
1894      ENDDO
1895      DO nsrf = 1, nbsrf
1896      DO i = 1, klon
1897c$$$         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001
1898        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1899      ENDDO
1900      ENDDO
1901      DO nsrf = 1, nbsrf
1902      DO i = 1, klon
1903            zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1904      ENDDO
1905      ENDDO
1906c
1907C calculs necessaires au calcul de l'albedo dans l'interface
1908c
1909      CALL orbite(FLOAT(julien),zlongi,dist)
1910      IF (cycle_diurne) THEN
1911        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1912        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1913      ELSE
1914        rmu0 = -999.999
1915      ENDIF
1916
1917      fder = dlw
1918
1919      CALL clmain(dtime,itap,date0,pctsrf,
1920     e            t_seri,q_seri,u_seri,v_seri,
1921     e            julien, rmu0,
1922     e            ok_veget, ocean, npas, nexca, ftsol,
1923     $            soil_model,ftsoil,
1924     $            paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,falblw,
1925     $            fluxlat,
1926     e            rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1927     e            rlon, rlat, cufi, cvfi, frugs,
1928     e            debut, lafin, agesno,rugoro ,
1929     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1930     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1931     s            dsens, devap,
1932     s            ycoefh,yu1,yv1)
1933
1934c
1935C§§§ PB
1936C§§§ Incrementation des flux
1937C§§
1938      zxfluxt=0.
1939      zxfluxq=0.
1940      zxfluxu=0.
1941      zxfluxv=0.
1942      DO nsrf = 1, nbsrf
1943        DO k = 1, klev
1944          DO i = 1, klon
1945            zxfluxt(i,k) = zxfluxt(i,k) +
1946     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1947            zxfluxq(i,k) = zxfluxq(i,k) +
1948     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1949            zxfluxu(i,k) = zxfluxu(i,k) +
1950     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1951            zxfluxv(i,k) = zxfluxv(i,k) +
1952     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1953          END DO
1954        END DO
1955      END DO
1956      DO i = 1, klon
1957         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1958c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1959         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1960         fder(i) = dlw(i) + dsens(i) + devap(i)
1961      ENDDO
1962
1963
1964      DO k = 1, klev
1965      DO i = 1, klon
1966         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1967         q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1968         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1969         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1970      ENDDO
1971      ENDDO
1972c
1973      IF (if_ebil.ge.2) THEN
1974        ztit='after clmain'
1975        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1976     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1977     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1978         call diagphy(paire,ztit,ip_ebil
1979     e      , zero_v, zero_v, zero_v, zero_v, sens
1980     e      , evap  , zero_v, zero_v, ztsol
1981     e      , d_h_vcol, d_qt, d_ec
1982     s      , fs_bound, fq_bound )
1983      END IF
1984C
1985c
1986c Incrementer la temperature du sol
1987c
1988      DO i = 1, klon
1989         zxtsol(i) = 0.0
1990         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1991     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
1992     $       THEN
1993             WRITE(*,*) 'physiq : pb sous surface au point ', i,
1994     $           pctsrf(i, 1 : nbsrf)
1995         ENDIF
1996      ENDDO
1997      DO nsrf = 1, nbsrf
1998      DO i = 1, klon
1999c$$$        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
2000            ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
2001            zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2002c$$$        ENDIF
2003      ENDDO
2004      ENDDO
2005
2006c
2007c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
2008c
2009      DO nsrf = 1, nbsrf
2010        DO i = 1, klon
2011          IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
2012        ENDDO
2013      ENDDO
2014
2015c
2016c Calculer la derive du flux infrarouge
2017c
2018c$$$      DO nsrf = 1, nbsrf
2019      DO i = 1, klon
2020c$$$        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
2021            dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3
2022c$$$     .          *(ftsol(i,nsrf)-zxtsol(i))
2023c$$$     .          *pctsrf(i,nsrf)
2024c$$$        ENDIF
2025c$$$      ENDDO
2026      ENDDO
2027c
2028c Appeler la convection (au choix)
2029c
2030      DO k = 1, klev
2031      DO i = 1, klon
2032         conv_q(i,k) = d_q_dyn(i,k)
2033     .               + d_q_vdf(i,k)/dtime
2034         conv_t(i,k) = d_t_dyn(i,k)
2035     .               + d_t_vdf(i,k)/dtime
2036      ENDDO
2037      ENDDO
2038      IF (check) THEN
2039         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
2040         PRINT*, "avantcon=", za
2041      ENDIF
2042      zx_ajustq = .FALSE.
2043      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2044      IF (zx_ajustq) THEN
2045         DO i = 1, klon
2046            z_avant(i) = 0.0
2047         ENDDO
2048         DO k = 1, klev
2049         DO i = 1, klon
2050            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2051     .                        *(paprs(i,k)-paprs(i,k+1))/RG
2052         ENDDO
2053         ENDDO
2054      ENDIF
2055      IF (iflag_con.EQ.1) THEN
2056          stop'reactiver le call conlmd dans physiq.F'
2057c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2058c    .             d_t_con, d_q_con,
2059c    .             rain_con, snow_con, ibas_con, itop_con)
2060      ELSE IF (iflag_con.EQ.2) THEN
2061      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2062     e            conv_t, conv_q, zxfluxq(1,1), omega,
2063     s            d_t_con, d_q_con, rain_con, snow_con,
2064     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2065     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
2066      WHERE (rain_con < 0.) rain_con = 0.
2067      WHERE (snow_con < 0.) snow_con = 0.
2068      DO i = 1, klon
2069         ibas_con(i) = klev+1 - kcbot(i)
2070         itop_con(i) = klev+1 - kctop(i)
2071      ENDDO
2072      ELSE IF (iflag_con.GE.3) THEN
2073c nb of tracers for the KE convection:
2074          if (nqmax .GE. 4) then
2075              ntra = nbtr
2076          else
2077              ntra = 1
2078          endif
2079          if (iflag_con.eq.4) then ! vectorise
2080          CALL conemav (dtime,paprs,pplay,t_seri,q_seri,
2081     .        u_seri,v_seri,tr_seri,nbtr,
2082     .        ema_work1,ema_work2,
2083     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2084     .        rain_con, snow_con, ibas_con, itop_con,
2085     .        upwd,dnwd,dnwd0,
2086c    .        Ma,cape,tvp,(/(nint(rflag(i)),i=1,size(rflag))/),
2087     .        Ma,cape,tvp,iflagctrl,
2088     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
2089
2090          else
2091
2092c          print*,'Avant conema OUI'
2093          CALL conema3 (dtime,
2094     .        paprs,pplay,t_seri,q_seri,
2095     .        u_seri,v_seri,tr_seri,nbtr,
2096     .        ema_work1,ema_work2,
2097     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2098     .        rain_con, snow_con, ibas_con, itop_con,
2099     .        upwd,dnwd,dnwd0,bas,top,
2100     .        Ma,cape,tvp,rflag,
2101     .        pbase
2102     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2103     .        ,clwcon0)
2104c          print*,'Apres conema3 '
2105
2106c Calculer l'humidite relative pour diagnostique
2107c
2108      DO k = 1, klev
2109      DO i = 1, klon
2110         zx_t = t_seri(i,k)
2111         IF (thermcep) THEN
2112            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2113            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2114            zx_qs  = MIN(0.5,zx_qs)
2115            zcor   = 1./(1.-retv*zx_qs)
2116            zx_qs  = zx_qs*zcor
2117         ELSE
2118           IF (zx_t.LT.t_coup) THEN
2119              zx_qs = qsats(zx_t)/pplay(i,k)
2120           ELSE
2121              zx_qs = qsatl(zx_t)/pplay(i,k)
2122           ENDIF
2123         ENDIF
2124         zqsat(i,k)=zx_qs
2125      ENDDO
2126      ENDDO
2127
2128c   calcul des propriétés des nuages convectifs
2129             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2130             call clouds_gno
2131     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2132
2133          endif
2134          DO i = 1, klon
2135            ema_pcb(i)  = pbase(i)
2136          ENDDO
2137          DO i = 1, klon
2138            ema_pct(i)  = paprs(i,itop_con(i))
2139          ENDDO
2140          DO i = 1, klon
2141            ema_cbmf(i) = ema_workcbmf(i)
2142          ENDDO     
2143      ELSE
2144          PRINT*, "iflag_con non-prevu", iflag_con
2145          CALL abort
2146      ENDIF
2147
2148c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2149c    .              d_u_con, d_v_con)
2150      DO k = 1, klev
2151        DO i = 1, klon
2152         t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2153         q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2154         u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2155         v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2156        ENDDO
2157      ENDDO
2158c
2159      IF (if_ebil.ge.2) THEN
2160        ztit='after convect'
2161        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2162     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2163     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2164         call diagphy(paire,ztit,ip_ebil
2165     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2166     e      , zero_v, rain_con, snow_con, ztsol
2167     e      , d_h_vcol, d_qt, d_ec
2168     s      , fs_bound, fq_bound )
2169      END IF
2170C
2171      IF (check) THEN
2172          za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
2173          PRINT*, "aprescon=", za
2174          zx_t = 0.0
2175          za = 0.0
2176          DO i = 1, klon
2177            za = za + paire(i)/FLOAT(klon)
2178            zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
2179          ENDDO
2180          zx_t = zx_t/za*dtime
2181          PRINT*, "Precip=", zx_t
2182      ENDIF
2183      IF (zx_ajustq) THEN
2184          DO i = 1, klon
2185            z_apres(i) = 0.0
2186          ENDDO
2187          DO k = 1, klev
2188            DO i = 1, klon
2189              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2190     .            *(paprs(i,k)-paprs(i,k+1))/RG
2191            ENDDO
2192          ENDDO
2193          DO i = 1, klon
2194            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2195     .          /z_apres(i)
2196          ENDDO
2197          DO k = 1, klev
2198            DO i = 1, klon
2199              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2200     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
2201                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
2202              ENDIF
2203            ENDDO
2204          ENDDO
2205      ENDIF
2206      zx_ajustq=.FALSE.
2207c
2208      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
2209c
2210          IF (iflag_con .NE. 2 .AND. debut) THEN
2211              PRINT*, 'Pour l instant, seul conflx fonctionne ',
2212     $            'avec traceurs', iflag_con
2213              PRINT*,' Mettre iflag_con',
2214     $            ' = 2 dans run.def et repasser'
2215c              CALL abort
2216              ENDIF
2217c
2218      ENDIF !--nqmax.GT.2
2219c
2220c Appeler l'ajustement sec
2221c
2222      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
2223      DO k = 1, klev
2224      DO i = 1, klon
2225         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
2226         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
2227      ENDDO
2228      ENDDO
2229c
2230      IF (if_ebil.ge.2) THEN
2231        ztit='after dry_adjust'
2232        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2233     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2234     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2235      END IF
2236
2237
2238c-------------------------------------------------------------------------
2239c  Caclul des ratqs
2240c-------------------------------------------------------------------------
2241
2242c      print*,'calcul des ratqs'
2243c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2244c   ----------------
2245c   on ecrase le tableau ratqsc calcule par clouds_gno
2246      if (iflag_cldcon.eq.1) then
2247         do k=1,klev
2248         do i=1,klon
2249            if(ptconv(i,k)) then
2250              ratqsc(i,k)=ratqsbas
2251     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2252            else
2253               ratqsc(i,k)=0.
2254            endif
2255         enddo
2256         enddo
2257      endif
2258
2259c   ratqs stables
2260c   -------------
2261      do k=1,klev
2262         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
2263     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)       
2264      enddo
2265
2266
2267c  ratqs final
2268c  -----------
2269      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
2270c   les ratqs sont une conbinaison de ratqss et ratqsc
2271c   ratqs final
2272c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
2273c   relaxation des ratqs
2274c         facttemps=exp(-pdtphys/1.e4)
2275         facteur=exp(-pdtphys*facttemps)
2276         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
2277         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
2278c         print*,'calcul des ratqs fini'
2279      else
2280c   on ne prend que le ratqs stable pour fisrtilp
2281         ratqs(:,:)=ratqss(:,:)
2282      endif
2283
2284
2285c
2286c Appeler le processus de condensation a grande echelle
2287c et le processus de precipitation
2288c-------------------------------------------------------------------------
2289      CALL fisrtilp(dtime,paprs,pplay,
2290     .           t_seri, q_seri,ptconv,ratqs,
2291     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2292     .           rain_lsc, snow_lsc,
2293     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
2294     .           frac_impa, frac_nucl,
2295     .           prfl, psfl, rhcl)
2296
2297      WHERE (rain_lsc < 0) rain_lsc = 0.
2298      WHERE (snow_lsc < 0) snow_lsc = 0.
2299      DO k = 1, klev
2300      DO i = 1, klon
2301         t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2302         q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2303         ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2304         cldfra(i,k) = rneb(i,k)
2305         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2306      ENDDO
2307      ENDDO
2308      IF (check) THEN
2309         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
2310         PRINT*, "apresilp=", za
2311         zx_t = 0.0
2312         za = 0.0
2313         DO i = 1, klon
2314            za = za + paire(i)/FLOAT(klon)
2315            zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
2316        ENDDO
2317         zx_t = zx_t/za*dtime
2318         PRINT*, "Precip=", zx_t
2319      ENDIF
2320c
2321      IF (if_ebil.ge.2) THEN
2322        ztit='after fisrt'
2323        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2324     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2325     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2326        call diagphy(paire,ztit,ip_ebil
2327     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2328     e      , zero_v, rain_lsc, snow_lsc, ztsol
2329     e      , d_h_vcol, d_qt, d_ec
2330     s      , fs_bound, fq_bound )
2331      END IF
2332c
2333c-------------------------------------------------------------------
2334c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2335c-------------------------------------------------------------------
2336
2337c 1. NUAGES CONVECTIFS
2338c
2339      IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2340
2341c Nuages diagnostiques pour Tiedtke
2342      CALL diagcld1(paprs,pplay,
2343     .             rain_con,snow_con,ibas_con,itop_con,
2344     .             diafra,dialiq)
2345      DO k = 1, klev
2346      DO i = 1, klon
2347      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2348         cldliq(i,k) = dialiq(i,k)
2349         cldfra(i,k) = diafra(i,k)
2350      ENDIF
2351      ENDDO
2352      ENDDO
2353
2354      ELSE IF (iflag_cldcon.eq.3) THEN
2355c  On prend pour les nuages convectifs le max du calcul de la
2356c  convection et du calcul du pas de temps précédent diminué d'un facteur
2357c  facttemps
2358c      facttemps=pdtphys/1.e4
2359      facteur = pdtphys *facttemps
2360      do k=1,klev
2361         do i=1,klon
2362            rnebcon(i,k)=rnebcon(i,k)*facteur
2363            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
2364     s      then
2365                rnebcon(i,k)=rnebcon0(i,k)
2366                clwcon(i,k)=clwcon0(i,k)
2367            endif
2368         enddo
2369      enddo
2370
2371c   On prend la somme des fractions nuageuses et des contenus en eau
2372      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2373      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2374
2375
2376      ENDIF
2377c
2378c 2. NUAGES STARTIFORMES
2379c
2380      IF (ok_stratus) THEN
2381      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2382      DO k = 1, klev
2383      DO i = 1, klon
2384      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2385         cldliq(i,k) = dialiq(i,k)
2386         cldfra(i,k) = diafra(i,k)
2387      ENDIF
2388      ENDDO
2389      ENDDO
2390      ENDIF
2391c
2392c Precipitation totale
2393c
2394      DO i = 1, klon
2395         rain_fall(i) = rain_con(i) + rain_lsc(i)
2396         snow_fall(i) = snow_con(i) + snow_lsc(i)
2397      ENDDO
2398c
2399      IF (if_ebil.ge.2) THEN
2400        ztit="after diagcld"
2401        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2402     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2403     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2404      END IF
2405c
2406c Calculer l'humidite relative pour diagnostique
2407c
2408      DO k = 1, klev
2409      DO i = 1, klon
2410         zx_t = t_seri(i,k)
2411         IF (thermcep) THEN
2412            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2413            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2414            zx_qs  = MIN(0.5,zx_qs)
2415            zcor   = 1./(1.-retv*zx_qs)
2416            zx_qs  = zx_qs*zcor
2417         ELSE
2418           IF (zx_t.LT.t_coup) THEN
2419              zx_qs = qsats(zx_t)/pplay(i,k)
2420           ELSE
2421              zx_qs = qsatl(zx_t)/pplay(i,k)
2422           ENDIF
2423         ENDIF
2424         zx_rh(i,k) = q_seri(i,k)/zx_qs
2425         zqsat(i,k)=zx_qs
2426      ENDDO
2427      ENDDO
2428c
2429c Calculer les parametres optiques des nuages et quelques
2430c parametres pour diagnostiques:
2431c
2432      if (ok_newmicro) then
2433      CALL newmicro (paprs, pplay,ok_newmicro,
2434     .            t_seri, cldliq, cldfra, cldtau, cldemi,
2435     .            cldh, cldl, cldm, cldt, cldq)
2436      else
2437      CALL nuage (paprs, pplay,
2438     .            t_seri, cldliq, cldfra, cldtau, cldemi,
2439     .            cldh, cldl, cldm, cldt, cldq)
2440      endif
2441c
2442c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2443c
2444      IF (MOD(itaprad,radpas).EQ.0) THEN
2445      DO i = 1, klon
2446         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2447     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
2448     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
2449     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
2450         albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2451     .               + falblw(i,is_lic) * pctsrf(i,is_lic)
2452     .               + falblw(i,is_ter) * pctsrf(i,is_ter)
2453     .               + falblw(i,is_sic) * pctsrf(i,is_sic)
2454      ENDDO
2455!      if (debut) then
2456!        albsol1 = albsol
2457!        albsollw1 = albsollw
2458!      endif
2459!      albsol = albsol1
2460!      albsollw = albsollw1
2461      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2462     e            (dist, rmu0, fract, co2_ppm, solaire,
2463     e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2464     e             wo,
2465     e             cldfra, cldemi, cldtau,
2466     s             heat,heat0,cool,cool0,radsol,albpla,
2467     s             topsw,toplw,solsw,sollw,
2468     s             sollwdown,
2469     s             topsw0,toplw0,solsw0,sollw0)
2470      itaprad = 0
2471      ENDIF
2472      itaprad = itaprad + 1
2473c
2474c Ajouter la tendance des rayonnements (tous les pas)
2475c
2476      DO k = 1, klev
2477      DO i = 1, klon
2478         t_seri(i,k) = t_seri(i,k)
2479     .               + (heat(i,k)-cool(i,k)) * dtime/86400.
2480      ENDDO
2481      ENDDO
2482c
2483      IF (if_ebil.ge.2) THEN
2484        ztit='after rad'
2485        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2486     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2487     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2488        call diagphy(paire,ztit,ip_ebil
2489     e      , topsw, toplw, solsw, sollw, zero_v
2490     e      , zero_v, zero_v, zero_v, ztsol
2491     e      , d_h_vcol, d_qt, d_ec
2492     s      , fs_bound, fq_bound )
2493      END IF
2494c
2495c
2496c Calculer l'hydrologie de la surface
2497c
2498c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2499c     .            agesno, ftsol,fqsol,fsnow, ruis)
2500c
2501      DO i = 1, klon
2502         zxqsol(i) = 0.0
2503         zxsnow(i) = 0.0
2504      ENDDO
2505      DO nsrf = 1, nbsrf
2506      DO i = 1, klon
2507         zxqsol(i) = zxqsol(i) + fqsol(i,nsrf)*pctsrf(i,nsrf)
2508         zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2509      ENDDO
2510      ENDDO
2511c
2512c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2513c
2514c$$$      DO nsrf = 1, nbsrf
2515c$$$      DO i = 1, klon
2516c$$$         IF (pctsrf(i,nsrf).LT.epsfra) THEN
2517c$$$            fqsol(i,nsrf) = zxqsol(i)
2518c$$$            fsnow(i,nsrf) = zxsnow(i)
2519c$$$         ENDIF
2520c$$$      ENDDO
2521c$$$      ENDDO
2522c
2523c Calculer le bilan du sol et la derive de temperature (couplage)
2524c
2525      DO i = 1, klon
2526         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2527      ENDDO
2528c
2529cmoddeblott(jan95)
2530c Appeler le programme de parametrisation de l'orographie
2531c a l'echelle sous-maille:
2532c
2533      IF (ok_orodr) THEN
2534c
2535c  selection des points pour lesquels le shema est actif:
2536        igwd=0
2537        DO i=1,klon
2538        itest(i)=0
2539c        IF ((zstd(i).gt.10.0)) THEN
2540        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2541          itest(i)=1
2542          igwd=igwd+1
2543          idx(igwd)=i
2544        ENDIF
2545        ENDDO
2546c        igwdim=MAX(1,igwd)
2547c
2548        CALL drag_noro(klon,klev,dtime,paprs,pplay,
2549     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
2550     e                   igwd,idx,itest,
2551     e                   t_seri, u_seri, v_seri,
2552     s                   zulow, zvlow, zustr, zvstr,
2553     s                   d_t_oro, d_u_oro, d_v_oro)
2554c
2555c  ajout des tendances
2556        DO k = 1, klev
2557        DO i = 1, klon
2558           t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2559           u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2560           v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2561        ENDDO
2562        ENDDO
2563c
2564      ENDIF ! fin de test sur ok_orodr
2565c
2566      IF (ok_orolf) THEN
2567c
2568c  selection des points pour lesquels le shema est actif:
2569        igwd=0
2570        DO i=1,klon
2571        itest(i)=0
2572        IF ((zpic(i)-zmea(i)).GT.100.) THEN
2573          itest(i)=1
2574          igwd=igwd+1
2575          idx(igwd)=i
2576        ENDIF
2577        ENDDO
2578c        igwdim=MAX(1,igwd)
2579c
2580        CALL lift_noro(klon,klev,dtime,paprs,pplay,
2581     e                   rlat,zmea,zstd,zpic,
2582     e                   itest,
2583     e                   t_seri, u_seri, v_seri,
2584     s                   zulow, zvlow, zustr, zvstr,
2585     s                   d_t_lif, d_u_lif, d_v_lif)
2586c
2587c  ajout des tendances
2588        DO k = 1, klev
2589        DO i = 1, klon
2590           t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2591           u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2592           v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2593        ENDDO
2594        ENDDO
2595c
2596      ENDIF ! fin de test sur ok_orolf
2597c
2598      IF (if_ebil.ge.2) THEN
2599        ztit='after orography'
2600        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2601     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2602     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2603      END IF
2604c
2605c
2606cAA
2607cAA Installation de l'interface online-offline pour traceurs
2608cAA
2609c====================================================================
2610c   Calcul  des tendances traceurs
2611c====================================================================
2612C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2613cKE43       des traceurs => il faut donc mettre des flags a .false.
2614      IF (iflag_con.GE.3) THEN
2615c           on ajoute les tendances calculees par KE43
2616c$$$ OM on onhibe la convection sur les traceurs
2617        DO iq=1, nqmax-2 ! Sandrine a -3 ???
2618c$$$ OM on inhibe la convection sur les traceur
2619c$$$        DO k = 1, nlev
2620c$$$        DO i = 1, klon
2621c$$$          tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2622c$$$        ENDDO
2623c$$$        ENDDO
2624        WRITE(iqn,'(i2.2)') iq
2625        CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2626        ENDDO
2627CMAF modif pour garder info du nombre de traceurs auxquels
2628C la physique s'applique
2629      ELSE
2630CMAF modif pour garder info du nombre de traceurs auxquels
2631C la physique s'applique
2632C
2633      call phytrac (rnpb,
2634     I                   debut,lafin,
2635     I                   nqmax-2,
2636     I                   nlon,nlev,dtime,
2637     I                   t,paprs,pplay,
2638     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2639     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2640     I                   frac_impa, frac_nucl,
2641     I                   rlon,presnivs,paire,pphis,
2642     O                   tr_seri)
2643      ENDIF
2644
2645      IF (offline) THEN
2646
2647         call phystokenc (
2648     I                   nlon,nlev,pdtphys,rlon,rlat,
2649     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2650     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
2651     I                   frac_impa, frac_nucl,
2652     I                   pphis,paire,dtime,itap)
2653
2654
2655      ENDIF
2656
2657c
2658c Calculer le transport de l'eau et de l'energie (diagnostique)
2659c
2660      CALL transp (paprs,zxtsol,
2661     e                   t_seri, q_seri, u_seri, v_seri, zphi,
2662     s                   ve, vq, ue, uq)
2663c
2664c Accumuler les variables a stocker dans les fichiers histoire:
2665c
2666c
2667c
2668      IF (if_ebil.ge.1) THEN
2669        ztit='after physic'
2670        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
2671     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2672     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2673C     Comme les tendances de la physique sont ajoute dans la dynamique,
2674C     on devrait avoir que la variation d'entalpie par la dynamique
2675C     est egale a la variation de la physique au pas de temps precedent.
2676C     Donc la somme de ces 2 variations devrait etre nulle.
2677        call diagphy(paire,ztit,ip_ebil
2678     e      , topsw, toplw, solsw, sollw, sens
2679     e      , evap, rain_fall, snow_fall, ztsol
2680     e      , d_h_vcol, d_qt, d_ec
2681     s      , fs_bound, fq_bound )
2682C
2683      d_h_vcol_phy=d_h_vcol
2684C
2685      END IF
2686C
2687      IF (ok_journe) THEN
2688c
2689      ndex2d = 0
2690      ndex3d = 0
2691c
2692c Champs 2D:
2693c
2694         zsto = dtime
2695         zout = dtime * FLOAT(ecrit_day)
2696         itau_w = itau_phy + itap
2697
2698         i = NINT(zout/zsto)
2699         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2700       CALL histwrite(nid_day,"phis",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2701         varname = 'phis'
2702         vartitle= 'Surface geop. height'
2703         varunits= '-'
2704c        call writephy(fid_day,prof2d_on,varname,pphis,vartitle,
2705c    .                                                    varunits)
2706c
2707         i = NINT(zout/zsto)
2708         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2709       CALL histwrite(nid_day,"aire",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2710         varname = 'aire'
2711         vartitle= 'Grid area'
2712         varunits= '-'
2713c        call writephy(fid_day,prof2d_on,varname,paire,vartitle,
2714c    .                                                    varunits)
2715C
2716      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2717      CALL histwrite(nid_day,"tsol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2718c     call writephy(fid_day,prof2d_av,'tsol',zxtsol,
2719c    .              'Surface Temperature','K')
2720c
2721C
2722      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_ter)
2723      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d ,zx_tmp_2d)
2724      CALL histwrite(nid_day,"tter",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2725c     call writephy(fid_day,prof2d_av,'tter',ftsol(1 : klon, is_ter),
2726c    .              'Surface Temperature','K')
2727C
2728      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_lic)
2729      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2730      CALL histwrite(nid_day,"tlic",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2731c     call writephy(fid_day,prof2d_av,'tlic',ftsol(1 : klon, is_lic),
2732c    .              'Surface Temperature','K')
2733C
2734      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_oce)
2735      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2736      CALL histwrite(nid_day,"toce",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2737c     call writephy(fid_day,prof2d_av,'toce',ftsol(1 : klon, is_oce),
2738c    .              'Surface Temperature','K')
2739C
2740      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_sic)
2741      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2742      CALL histwrite(nid_day,"tsic",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2743c     call writephy(fid_day,prof2d_av,'tsic',ftsol(1 : klon, is_sic),
2744c    .              'Surface Temperature','K')
2745C
2746      DO i = 1, klon
2747         zx_tmp_fi2d(i) = paprs(i,1)
2748      ENDDO
2749      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2750      CALL histwrite(nid_day,"psol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2751c Essai writephys
2752      varname = 'psol'
2753      vartitle= 'pression au sol'
2754      varunits= 'hPa'
2755c     call writephy(fid_day,prof2d_av,varname,zx_tmp_fi2d,vartitle,
2756c    .                                                    varunits)
2757c
2758      DO i = 1, klon
2759         zx_tmp_fi2d(i) = (rain_fall(i) + snow_fall(i))* 86400.
2760      ENDDO
2761      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2762      CALL histwrite(nid_day,"rain",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2763c     call writephy(fid_day,prof2d_av,'rain',zx_tmp_fi2d,
2764c    .              'Precipitation','mm/day')
2765
2766
2767c
2768      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
2769      CALL histwrite(nid_day,"snow",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2770c     call writephy(fid_day,prof2d_av,'snow',snow_fall,
2771c    .              'Snow','mm/day')
2772c
2773      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxsnow,zx_tmp_2d)
2774      CALL histwrite(nid_day,"snow_cov",itau_w,zx_tmp_2d,iim*jjmp1,
2775     .               ndex2d)
2776c     call writephy(fid_day,prof2d_av,'snow_cov',zxsnow,
2777c    .              'Snow cover','mm')
2778c
2779      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2780      CALL histwrite(nid_day,"evap",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2781c     call writephy(fid_day,prof2d_av,'evap',evap,
2782c    .              'Evaporation','mm/day')
2783c
2784      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
2785      CALL histwrite(nid_day,"tops",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2786c     call writephy(fid_day,prof2d_av,'tops',topsw,
2787c    .              'Solar rad. at TOA','W/m2')
2788c
2789      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2790      CALL histwrite(nid_day,"topl",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2791c     call writephy(fid_day,prof2d_av,'topl',toplw,
2792c    .              'IR rad. at TOA','W/m2')
2793c
2794      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2795      CALL histwrite(nid_day,"sols",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2796c     call writephy(fid_day,prof2d_av,'sols',solsw,
2797c    .              'Solar rad. at surf.','W/m2')
2798c
2799      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2800      CALL histwrite(nid_day,"soll",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2801c     call writephy(fid_day,prof2d_av,'soll',sollw,
2802c    .              'IR rad. at surface','W/m2')
2803c
2804      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
2805      CALL histwrite(nid_day,"solldown",itau_w,zx_tmp_2d,iim*jjmp1,
2806     .               ndex2d)
2807c     call writephy(fid_day,prof2d_av,'solldown',sollwdown,
2808c    .              'Down. IR rad. at surface','W/m2')
2809c
2810      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2811      CALL histwrite(nid_day,"bils",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2812c     call writephy(fid_day,prof2d_av,'bils',bils,
2813c    .              'Surf. total heat flux','W/m2')
2814c
2815      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2816      CALL histwrite(nid_day,"sens",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2817c     call writephy(fid_day,prof2d_av,'sens',sens,
2818c    .              'Sensible heat flux','W/m2')
2819c
2820      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2821      CALL histwrite(nid_day,"fder",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2822c     call writephy(fid_day,prof2d_av,'fder',fder,
2823c    .              'Heat flux derivation','W/m2')
2824c
2825c
2826      DO nsrf = 1, nbsrf
2827C§§§
2828        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2829        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2830        CALL histwrite(nid_day,"pourc_"//clnsurf(nsrf),itau_w,
2831     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2832c       call writephy(fid_day,prof2d_av,'pourc_'//clnsurf(nsrf),
2833c    .                pctsrf( 1 : klon, nsrf),
2834c    .                'Fraction'//clnsurf(nsrf),'-')
2835C
2836        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
2837        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2838        CALL histwrite(nid_day,"tsol_"//clnsurf(nsrf),itau_w,
2839     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2840c       call writephy(fid_day,prof2d_av,'tsol_'//clnsurf(nsrf),
2841c    .                ftsol( 1 : klon, nsrf),
2842c    .                'Surf. Temp'//clnsurf(nsrf),'K')
2843C
2844        zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
2845        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2846        CALL histwrite(nid_day,"sens_"//clnsurf(nsrf),itau_w,
2847     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2848c       call writephy(fid_day,prof2d_av,'sens_'//clnsurf(nsrf),
2849c    .                fluxt( 1 : klon, 1, nsrf),
2850c    .                'Sensible heat flux '//clnsurf(nsrf),'W/m2')
2851
2852        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
2853        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2854        CALL histwrite(nid_day,"lat_"//clnsurf(nsrf),itau_w,
2855     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2856c       call writephy(fid_day,prof2d_av,'lat_'//clnsurf(nsrf),
2857c    .                fluxlat( 1 : klon, nsrf),
2858c    .                'Latent heat flux '//clnsurf(nsrf),'W/m2')
2859C
2860        zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)
2861        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2862        CALL histwrite(nid_day,"taux_"//clnsurf(nsrf),itau_w,
2863     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2864c       call writephy(fid_day,prof2d_av,'taux_'//clnsurf(nsrf),
2865c    .                fluxu( 1 : klon, 1, nsrf),
2866c    .                'Zonal wind stress '//clnsurf(nsrf),'Pa')
2867C     
2868        zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)
2869        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2870        CALL histwrite(nid_day,"tauy_"//clnsurf(nsrf),itau_w,
2871     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2872c       call writephy(fid_day,prof2d_av,'tauy_'//clnsurf(nsrf),
2873c    .                fluxv( 1 : klon, 1, nsrf),
2874c    .                'Meridional wind stress '//clnsurf(nsrf),'Pa')
2875C
2876        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
2877        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2878        CALL histwrite(nid_day,"albe_"//clnsurf(nsrf),itau_w,
2879     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2880c       call writephy(fid_day,prof2d_av,'albe_'//clnsurf(nsrf),
2881c    .                falbe( 1 : klon, nsrf),
2882c    .                'Albedo surf. SW'//clnsurf(nsrf),'-')
2883c       call writephy(fid_day,prof2d_av,'alblw_'//clnsurf(nsrf),
2884c    .                falblw( 1 : klon, nsrf),
2885c    .                'Albedo surf. LW'//clnsurf(nsrf),'-')
2886C
2887        zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
2888        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2889        CALL histwrite(nid_day,"rugs_"//clnsurf(nsrf),itau_w,
2890     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2891c       call writephy(fid_day,prof2d_av,'rugs_'//clnsurf(nsrf),
2892c    .                frugs( 1 : klon, nsrf),
2893c    .                'Rugosity '//clnsurf(nsrf),' - ')
2894C
2895      END DO 
2896C
2897c$$$      DO i = 1, klon
2898c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
2899c$$$      ENDDO
2900c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2901c$$$      CALL histwrite(nid_day,"sicf",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2902c
2903      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
2904      CALL histwrite(nid_day,"cldl",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2905c     call writephy(fid_day,prof2d_av,'cldl',cldl,
2906c    .              'Low-level cloudiness','-')
2907c
2908      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
2909      CALL histwrite(nid_day,"cldm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2910c     call writephy(fid_day,prof2d_av,'cldm',cldm,
2911c    .              'Mid-level cloudiness','-')
2912c
2913      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
2914      CALL histwrite(nid_day,"cldh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2915c     call writephy(fid_day,prof2d_av,'cldh',cldh,
2916c    .              'High-level cloudiness','-')
2917c
2918      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
2919      CALL histwrite(nid_day,"cldt",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2920c     call writephy(fid_day,prof2d_av,'cldt',cldt,
2921c    .              'Total cloudiness','-')
2922c
2923      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
2924      CALL histwrite(nid_day,"cldq",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2925c     call writephy(fid_day,prof2d_av,'cldq',cldq,
2926c    .              'Cloud liquid water path','-')
2927c
2928c Champs 3D:
2929c
2930      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2931      CALL histwrite(nid_day,"temp",itau_w,zx_tmp_3d,
2932     .                                   iim*jjmp1*klev,ndex3d)
2933c Essai writephys
2934      varname = 'temp'
2935      vartitle= 'temperature 3D'
2936      varunits= 'K'
2937c     call writephy(fid_day,prof3d_av,varname,t_seri,vartitle,varunits)
2938c
2939      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
2940      CALL histwrite(nid_day,"ovap",itau_w,zx_tmp_3d,
2941     .                                   iim*jjmp1*klev,ndex3d)
2942c     call writephy(fid_day,prof3d_av,'ovap',qx(1,1,ivap),
2943c    .              'Specific humidity','Kg/Kg')
2944c
2945      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2946      CALL histwrite(nid_day,"geop",itau_w,zx_tmp_3d,
2947     .                                   iim*jjmp1*klev,ndex3d)
2948c     call writephy(fid_day,prof3d_av,'geop',zphi,
2949c    .              'Geopotential height','m')
2950c
2951      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2952      CALL histwrite(nid_day,"vitu",itau_w,zx_tmp_3d,
2953     .                                   iim*jjmp1*klev,ndex3d)
2954c     call writephy(fid_day,prof3d_av,'vitu',u_seri,
2955c    .              'Zonal wind','m/s')
2956c
2957      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2958      CALL histwrite(nid_day,"vitv",itau_w,zx_tmp_3d,
2959     .                                   iim*jjmp1*klev,ndex3d)
2960c     call writephy(fid_day,prof3d_av,'vitv',v_seri,
2961c    .              'Meridional wind','m/s')
2962c
2963      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
2964      CALL histwrite(nid_day,"vitw",itau_w,zx_tmp_3d,
2965     .                                   iim*jjmp1*klev,ndex3d)
2966c     call writephy(fid_day,prof3d_av,'vitw',omega,
2967c    .              'Vertical wind','m/s')
2968c
2969      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2970      CALL histwrite(nid_day,"pres",itau_w,zx_tmp_3d,
2971     .                                   iim*jjmp1*klev,ndex3d)
2972c     call writephy(fid_day,prof3d_av,'pres',pplay,
2973c    .              'Air pressure','Pa')
2974
2975c
2976      if (ok_sync) then
2977c       call writephy_sync(fid_day)
2978        call histsync(nid_day)
2979      endif
2980      ENDIF
2981C
2982      IF (ok_mensuel) THEN
2983c
2984      ndex2d = 0
2985      ndex3d = 0
2986c
2987c Champs 2D:
2988c
2989         zsto = dtime
2990         zout = dtime * ecrit_mth
2991         itau_w = itau_phy + itap
2992
2993      i = NINT(zout/zsto)
2994      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2995      CALL histwrite(nid_mth,"phis",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
2996C
2997      i = NINT(zout/zsto)
2998      CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2999      CALL histwrite(nid_mth,"aire",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3000
3001      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
3002      CALL histwrite(nid_mth,"tsol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3003c
3004      DO i = 1, klon
3005         zx_tmp_fi2d(i) = paprs(i,1)
3006      ENDDO
3007      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3008      CALL histwrite(nid_mth,"psol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3009c
3010      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxqsol,zx_tmp_2d)
3011      CALL histwrite(nid_mth,"qsol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3012c
3013      DO i = 1, klon
3014         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
3015      ENDDO
3016      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3017      CALL histwrite(nid_mth,"rain",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3018c
3019      DO i = 1, klon
3020         zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
3021      ENDDO
3022      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3023      CALL histwrite(nid_mth,"plul",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3024c
3025      DO i = 1, klon
3026         zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
3027      ENDDO
3028      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3029      CALL histwrite(nid_mth,"pluc",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3030c
3031      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
3032      CALL histwrite(nid_mth,"snow",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3033c
3034      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxsnow,zx_tmp_2d)
3035      CALL histwrite(nid_mth,"snow_cov",itau_w,zx_tmp_2d,iim*jjmp1,
3036     .               ndex2d)
3037c
3038      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
3039      CALL histwrite(nid_mth,"evap",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3040c
3041      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
3042      CALL histwrite(nid_mth,"tops",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3043c
3044      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
3045      CALL histwrite(nid_mth,"topl",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3046c
3047      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
3048      CALL histwrite(nid_mth,"sols",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3049c
3050      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
3051      CALL histwrite(nid_mth,"soll",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3052c
3053      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
3054      CALL histwrite(nid_mth,"solldown",itau_w,zx_tmp_2d,iim*jjmp1,
3055     .               ndex2d)
3056c
3057      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
3058      CALL histwrite(nid_mth,"tops0",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3059c
3060      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw0,zx_tmp_2d)
3061      CALL histwrite(nid_mth,"topl0",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3062c
3063      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw0,zx_tmp_2d)
3064      CALL histwrite(nid_mth,"sols0",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3065c
3066      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw0,zx_tmp_2d)
3067      CALL histwrite(nid_mth,"soll0",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3068c
3069      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
3070      CALL histwrite(nid_mth,"bils",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3071c
3072      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
3073      CALL histwrite(nid_mth,"sens",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3074c
3075      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
3076      CALL histwrite(nid_mth,"fder",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3077c
3078c
3079c      DO i = 1, klon
3080c         zx_tmp_fi2d(i) = fluxu(i,1)
3081c      ENDDO
3082c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3083c      CALL histwrite(nid_mth,"frtu",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3084c
3085c      DO i = 1, klon
3086c         zx_tmp_fi2d(i) = fluxv(i,1)
3087c      ENDDO
3088c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3089c      CALL histwrite(nid_mth,"frtv",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3090c
3091      DO nsrf = 1, nbsrf
3092C§§§
3093        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
3094        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3095        CALL histwrite(nid_mth,"pourc_"//clnsurf(nsrf),itau_w,
3096     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3097C
3098        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
3099        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3100        CALL histwrite(nid_mth,"tsol_"//clnsurf(nsrf),itau_w,
3101     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3102C
3103        zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
3104        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3105        CALL histwrite(nid_mth,"sens_"//clnsurf(nsrf),itau_w,
3106     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3107C
3108        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
3109        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3110        CALL histwrite(nid_mth,"lat_"//clnsurf(nsrf),itau_w,
3111     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3112C
3113        zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)
3114        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3115        CALL histwrite(nid_mth,"taux_"//clnsurf(nsrf),itau_w,
3116     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3117C     
3118        zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)
3119        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3120        CALL histwrite(nid_mth,"tauy_"//clnsurf(nsrf),itau_w,
3121     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3122C
3123        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
3124        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3125        CALL histwrite(nid_mth,"albe_"//clnsurf(nsrf),itau_w,
3126     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3127C
3128        zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
3129        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3130        CALL histwrite(nid_mth,"rugs_"//clnsurf(nsrf),itau_w,
3131     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3132c
3133      zx_tmp_fi2d(1 : klon) = agesno( 1 : klon, nsrf)
3134      CALL gr_fi_ecrit(1, klon,iim,jjmp1, agesno,zx_tmp_2d)
3135      CALL histwrite(nid_mth,"ages_"//clnsurf(nsrf),itau_w
3136     $    ,zx_tmp_2d,iim*jjmp1,ndex2d)
3137
3138      END DO 
3139c$$$      DO i = 1, klon
3140c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
3141c$$$      ENDDO
3142c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3143c$$$      CALL histwrite(nid_mth,"sicf",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3144c
3145      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
3146      CALL histwrite(nid_mth,"albs",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3147      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsollw,zx_tmp_2d)
3148      CALL histwrite(nid_mth,"albslw",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3149c
3150      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
3151      CALL histwrite(nid_mth,"cdrm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3152c
3153      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
3154      CALL histwrite(nid_mth,"cdrh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3155c
3156      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
3157      CALL histwrite(nid_mth,"cldl",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3158c
3159      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
3160      CALL histwrite(nid_mth,"cldm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3161c
3162      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
3163      CALL histwrite(nid_mth,"cldh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3164c
3165      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
3166      CALL histwrite(nid_mth,"cldt",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3167c
3168      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
3169      CALL histwrite(nid_mth,"cldq",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3170c
3171      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ue,zx_tmp_2d)
3172      CALL histwrite(nid_mth,"ue",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3173c
3174      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ve,zx_tmp_2d)
3175      CALL histwrite(nid_mth,"ve",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3176c
3177      CALL gr_fi_ecrit(1, klon,iim,jjmp1, uq,zx_tmp_2d)
3178      CALL histwrite(nid_mth,"uq",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3179c
3180      CALL gr_fi_ecrit(1, klon,iim,jjmp1, vq,zx_tmp_2d)
3181      CALL histwrite(nid_mth,"vq",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3182cKE43
3183      IF (iflag_con .GE. 3) THEN ! sb
3184c
3185      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cape,zx_tmp_2d)
3186      CALL histwrite(nid_mth,"cape",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3187c
3188      CALL gr_fi_ecrit(1, klon,iim,jjmp1,pbase,zx_tmp_2d)
3189      CALL histwrite(nid_mth,"pbase",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3190c
3191      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_pct,zx_tmp_2d)
3192      CALL histwrite(nid_mth,"ptop",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3193c
3194      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_cbmf,zx_tmp_2d)
3195      CALL histwrite(nid_mth,"fbase",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3196c
3197c
3198      ENDIF
3199c34EK
3200c
3201c Champs 3D:
3202C
3203      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
3204      CALL histwrite(nid_mth,"temp",itau_w,zx_tmp_3d,
3205     .                                   iim*jjmp1*klev,ndex3d)
3206c
3207      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
3208      CALL histwrite(nid_mth,"ovap",itau_w,zx_tmp_3d,
3209     .                                   iim*jjmp1*klev,ndex3d)
3210c
3211      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
3212      CALL histwrite(nid_mth,"geop",itau_w,zx_tmp_3d,
3213     .                                   iim*jjmp1*klev,ndex3d)
3214c
3215      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
3216      CALL histwrite(nid_mth,"vitu",itau_w,zx_tmp_3d,
3217     .                                   iim*jjmp1*klev,ndex3d)
3218c
3219      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
3220      CALL histwrite(nid_mth,"vitv",itau_w,zx_tmp_3d,
3221     .                                   iim*jjmp1*klev,ndex3d)
3222c
3223      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
3224      CALL histwrite(nid_mth,"vitw",itau_w,zx_tmp_3d,
3225     .                                   iim*jjmp1*klev,ndex3d)
3226c
3227      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
3228      CALL histwrite(nid_mth,"pres",itau_w,zx_tmp_3d,
3229     .                                   iim*jjmp1*klev,ndex3d)
3230c
3231      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldfra, zx_tmp_3d)
3232      CALL histwrite(nid_mth,"rneb",itau_w,zx_tmp_3d,
3233     .                                   iim*jjmp1*klev,ndex3d)
3234c
3235      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_rh, zx_tmp_3d)
3236      CALL histwrite(nid_mth,"rhum",itau_w,zx_tmp_3d,
3237     .                                   iim*jjmp1*klev,ndex3d)
3238c
3239      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldliq, zx_tmp_3d)
3240      CALL histwrite(nid_mth,"oliq",itau_w,zx_tmp_3d,
3241     .                                   iim*jjmp1*klev,ndex3d)
3242c
3243      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_dyn, zx_tmp_3d)
3244      CALL histwrite(nid_mth,"dtdyn",itau_w,zx_tmp_3d,
3245     .                                   iim*jjmp1*klev,ndex3d)
3246c
3247      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_dyn, zx_tmp_3d)
3248      CALL histwrite(nid_mth,"dqdyn",itau_w,zx_tmp_3d,
3249     .                                   iim*jjmp1*klev,ndex3d)
3250c
3251      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_con, zx_tmp_3d)
3252      CALL histwrite(nid_mth,"dtcon",itau_w,zx_tmp_3d,
3253     .                                   iim*jjmp1*klev,ndex3d)
3254c
3255      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_con, zx_tmp_3d)
3256      CALL histwrite(nid_mth,"dqcon",itau_w,zx_tmp_3d,
3257     .                                   iim*jjmp1*klev,ndex3d)
3258c
3259      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_lsc, zx_tmp_3d)
3260      CALL histwrite(nid_mth,"dtlsc",itau_w,zx_tmp_3d,
3261     .                                   iim*jjmp1*klev,ndex3d)
3262c
3263      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_lsc, zx_tmp_3d)
3264      CALL histwrite(nid_mth,"dqlsc",itau_w,zx_tmp_3d,
3265     .                                   iim*jjmp1*klev,ndex3d)
3266c
3267      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
3268      CALL histwrite(nid_mth,"dtvdf",itau_w,zx_tmp_3d,
3269     .                                   iim*jjmp1*klev,ndex3d)
3270c
3271      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
3272      CALL histwrite(nid_mth,"dqvdf",itau_w,zx_tmp_3d,
3273     .                                   iim*jjmp1*klev,ndex3d)
3274c
3275      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_eva, zx_tmp_3d)
3276      CALL histwrite(nid_mth,"dteva",itau_w,zx_tmp_3d,
3277     .                                   iim*jjmp1*klev,ndex3d)
3278c
3279      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_eva, zx_tmp_3d)
3280      CALL histwrite(nid_mth,"dqeva",itau_w,zx_tmp_3d,
3281     .                                   iim*jjmp1*klev,ndex3d)
3282c
3283      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zpt_conv, zx_tmp_3d)
3284      CALL histwrite(nid_mth,"ptconv",itau_w,zx_tmp_3d,
3285     .                                   iim*(jjmp1)*klev,ndex3d)
3286c
3287      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, ratqs, zx_tmp_3d)
3288      CALL histwrite(nid_mth,"ratqs",itau_w,zx_tmp_3d,
3289     .                                   iim*(jjmp1)*klev,ndex3d)
3290c
3291      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
3292      CALL histwrite(nid_mth,"dtajs",itau_w,zx_tmp_3d,
3293     .                                   iim*jjmp1*klev,ndex3d)
3294c
3295      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_ajs, zx_tmp_3d)
3296      CALL histwrite(nid_mth,"dqajs",itau_w,zx_tmp_3d,
3297     .                                   iim*jjmp1*klev,ndex3d)
3298c
3299      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat, zx_tmp_3d)
3300      CALL histwrite(nid_mth,"dtswr",itau_w,zx_tmp_3d,
3301     .                                   iim*jjmp1*klev,ndex3d)
3302c
3303      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat0, zx_tmp_3d)
3304      CALL histwrite(nid_mth,"dtsw0",itau_w,zx_tmp_3d,
3305     .                                   iim*jjmp1*klev,ndex3d)
3306c
3307      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool, zx_tmp_3d)
3308      CALL histwrite(nid_mth,"dtlwr",itau_w,zx_tmp_3d,
3309     .                                   iim*jjmp1*klev,ndex3d)
3310c
3311      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool0, zx_tmp_3d)
3312      CALL histwrite(nid_mth,"dtlw0",itau_w,zx_tmp_3d,
3313     .                                   iim*jjmp1*klev,ndex3d)
3314c
3315      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_vdf, zx_tmp_3d)
3316      CALL histwrite(nid_mth,"duvdf",itau_w,zx_tmp_3d,
3317     .                                   iim*jjmp1*klev,ndex3d)
3318c
3319      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_vdf, zx_tmp_3d)
3320      CALL histwrite(nid_mth,"dvvdf",itau_w,zx_tmp_3d,
3321     .                                   iim*jjmp1*klev,ndex3d)
3322c
3323      IF (ok_orodr) THEN
3324      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_oro, zx_tmp_3d)
3325      CALL histwrite(nid_mth,"duoro",itau_w,zx_tmp_3d,
3326     .                                   iim*jjmp1*klev,ndex3d)
3327c
3328      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_oro, zx_tmp_3d)
3329      CALL histwrite(nid_mth,"dvoro",itau_w,zx_tmp_3d,
3330     .                                   iim*jjmp1*klev,ndex3d)
3331c
3332      ENDIF
3333C
3334      IF (ok_orolf) THEN
3335      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_lif, zx_tmp_3d)
3336      CALL histwrite(nid_mth,"dulif",itau_w,zx_tmp_3d,
3337     .                                   iim*jjmp1*klev,ndex3d)
3338c
3339      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_lif, zx_tmp_3d)
3340      CALL histwrite(nid_mth,"dvlif",itau_w,zx_tmp_3d,
3341     .                                   iim*jjmp1*klev,ndex3d)
3342      ENDIF
3343C
3344      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, wo, zx_tmp_3d)
3345      CALL histwrite(nid_mth,"ozone",itau_w,zx_tmp_3d,
3346     .                                   iim*jjmp1*klev,ndex3d)
3347c
3348      IF (nqmax.GE.3) THEN
3349      DO iq=1,nqmax-2
3350      IF (iq.LE.99) THEN
3351         CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,iq+2), zx_tmp_3d)
3352         WRITE(str2,'(i2.2)') iq
3353         CALL histwrite(nid_mth,"trac"//str2,itau_w,zx_tmp_3d,
3354     .                                   iim*jjmp1*klev,ndex3d)
3355      ELSE
3356         PRINT*, "Trop de traceurs"
3357         CALL abort
3358      ENDIF
3359      ENDDO
3360      ENDIF
3361cKE43
3362      IF (iflag_con.GE.3) THEN ! (sb)
3363c
3364      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, upwd, zx_tmp_3d)
3365      CALL histwrite(nid_mth,"upwd",itau_w,zx_tmp_3d,
3366     .                                   iim*jjmp1*klev,ndex3d)
3367c
3368      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd, zx_tmp_3d)
3369      CALL histwrite(nid_mth,"dnwd",itau_w,zx_tmp_3d,
3370     .                                   iim*jjmp1*klev,ndex3d)
3371c
3372      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd0, zx_tmp_3d)
3373      CALL histwrite(nid_mth,"dnwd0",itau_w,zx_tmp_3d,
3374     .                                   iim*jjmp1*klev,ndex3d)
3375c
3376      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, Ma, zx_tmp_3d)
3377      CALL histwrite(nid_mth,"Ma",itau_w,zx_tmp_3d,
3378     .                                   iim*jjmp1*klev,ndex3d)
3379c
3380c
3381      ENDIF
3382c34EK
3383c
3384      if (ok_sync) then
3385        call histsync(nid_mth)
3386      endif
3387      ENDIF
3388c
3389      IF (ok_instan) THEN
3390c
3391      ndex2d = 0
3392      ndex3d = 0
3393c
3394c Champs 2D:
3395c
3396         zsto = dtime * ecrit_ins
3397         zout = dtime * ecrit_ins
3398         itau_w = itau_phy + itap
3399
3400         i = NINT(zout/zsto)
3401      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
3402      CALL histwrite(nid_ins,"phis",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3403c
3404         i = NINT(zout/zsto)
3405      CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
3406      CALL histwrite(nid_ins,"aire",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3407
3408      DO i = 1, klon
3409         zx_tmp_fi2d(i) = paprs(i,1)
3410      ENDDO
3411      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3412      CALL histwrite(nid_ins,"psol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3413c
3414      DO i = 1, klon
3415         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
3416      ENDDO
3417      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3418      CALL histwrite(nid_ins,"rain",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3419c
3420      DO i = 1, klon
3421         zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
3422      ENDDO
3423      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3424      CALL histwrite(nid_ins,"plul",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3425c
3426      DO i = 1, klon
3427         zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
3428      ENDDO
3429      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
3430      CALL histwrite(nid_ins,"pluc",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3431
3432      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
3433      CALL histwrite(nid_ins,"tsol",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3434c
3435      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
3436      CALL histwrite(nid_ins,"snow",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3437
3438      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
3439      CALL histwrite(nid_ins,"cdrm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3440c
3441      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
3442      CALL histwrite(nid_ins,"cdrh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3443c
3444      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
3445      CALL histwrite(nid_ins,"topl",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3446c
3447      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
3448      CALL histwrite(nid_ins,"evap",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3449c
3450      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
3451      CALL histwrite(nid_ins,"sols",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3452c
3453      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
3454      CALL histwrite(nid_ins,"soll",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3455c
3456      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
3457      CALL histwrite(nid_ins,"solldown",itau_w,zx_tmp_2d,iim*jjmp1,
3458     .                ndex2d)
3459c
3460      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
3461      CALL histwrite(nid_ins,"bils",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3462c
3463      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
3464      CALL histwrite(nid_ins,"sens",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3465c
3466      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
3467      CALL histwrite(nid_ins,"fder",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3468c
3469      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_oce),zx_tmp_2d)
3470      CALL histwrite(nid_ins,"dtsvdfo",itau_w,zx_tmp_2d,iim*jjmp1,
3471     .               ndex2d)
3472c
3473      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_ter),zx_tmp_2d)
3474      CALL histwrite(nid_ins,"dtsvdft",itau_w,zx_tmp_2d,iim*jjmp1,
3475     .               ndex2d)
3476c
3477      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_lic),zx_tmp_2d)
3478      CALL histwrite(nid_ins,"dtsvdfg",itau_w,zx_tmp_2d,iim*jjmp1,
3479     .               ndex2d)
3480c
3481      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_sic),zx_tmp_2d)
3482      CALL histwrite(nid_ins,"dtsvdfi",itau_w,zx_tmp_2d,iim*jjmp1,
3483     .               ndex2d)
3484
3485      DO nsrf = 1, nbsrf
3486C§§§
3487        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
3488        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3489        CALL histwrite(nid_ins,"pourc_"//clnsurf(nsrf),itau_w,
3490     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3491C
3492        zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
3493        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3494        CALL histwrite(nid_ins,"sens_"//clnsurf(nsrf),itau_w,
3495     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3496C
3497        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
3498        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3499        CALL histwrite(nid_ins,"lat_"//clnsurf(nsrf),itau_w,
3500     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3501C
3502        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
3503        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3504        CALL histwrite(nid_ins,"tsol_"//clnsurf(nsrf),itau_w,
3505     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3506C
3507        zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)
3508        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3509        CALL histwrite(nid_ins,"taux_"//clnsurf(nsrf),itau_w,
3510     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3511C     
3512        zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)
3513        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3514        CALL histwrite(nid_ins,"tauy_"//clnsurf(nsrf),itau_w,
3515     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3516C
3517        zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
3518        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3519        CALL histwrite(nid_ins,"rugs_"//clnsurf(nsrf),itau_w,
3520     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3521C
3522        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
3523        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
3524        CALL histwrite(nid_ins,"albe_"//clnsurf(nsrf),itau_w,
3525     $      zx_tmp_2d,iim*jjmp1,ndex2d)
3526C
3527      END DO 
3528      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
3529      CALL histwrite(nid_ins,"albs",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3530      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsollw,zx_tmp_2d)
3531      CALL histwrite(nid_ins,"albslw",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3532c
3533      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxsnow,zx_tmp_2d)
3534      CALL histwrite(nid_ins,"snow_cov",itau_w,zx_tmp_2d,iim*jjmp1,
3535     .               ndex2d)
3536c
3537      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxrugs,zx_tmp_2d)
3538      CALL histwrite(nid_ins,"rugs",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
3539c
3540c Champs 3D:
3541c
3542      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
3543      CALL histwrite(nid_ins,"temp",itau_w,zx_tmp_3d,
3544     .                                   iim*jjmp1*klev,ndex3d)
3545c
3546      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
3547      CALL histwrite(nid_ins,"vitu",itau_w,zx_tmp_3d,
3548     .                                   iim*jjmp1*klev,ndex3d)
3549c
3550      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
3551      CALL histwrite(nid_ins,"vitv",itau_w,zx_tmp_3d,
3552     .                                   iim*jjmp1*klev,ndex3d)
3553c
3554      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
3555      CALL histwrite(nid_ins,"geop",itau_w,zx_tmp_3d,
3556     .                                   iim*jjmp1*klev,ndex3d)
3557c
3558      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
3559      CALL histwrite(nid_ins,"pres",itau_w,zx_tmp_3d,
3560     .                                   iim*jjmp1*klev,ndex3d)
3561c
3562      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
3563      CALL histwrite(nid_ins,"dtvdf",itau_w,zx_tmp_3d,
3564     .                                   iim*jjmp1*klev,ndex3d)
3565c
3566      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
3567      CALL histwrite(nid_ins,"dqvdf",itau_w,zx_tmp_3d,
3568     .                                   iim*jjmp1*klev,ndex3d)
3569
3570c
3571      if (ok_sync) then
3572        call histsync(nid_ins)
3573      endif
3574      ENDIF
3575c
3576c
3577c Ecrire la bande regionale (binaire grads)
3578      IF (ok_region .AND. mod(itap,ecrit_reg).eq.0) THEN
3579         CALL ecriregs(84,zxtsol)
3580         CALL ecriregs(84,paprs(1,1))
3581         CALL ecriregs(84,topsw)
3582         CALL ecriregs(84,toplw)
3583         CALL ecriregs(84,solsw)
3584         CALL ecriregs(84,sollw)
3585         CALL ecriregs(84,rain_fall)
3586         CALL ecriregs(84,snow_fall)
3587         CALL ecriregs(84,evap)
3588         CALL ecriregs(84,sens)
3589         CALL ecriregs(84,bils)
3590         CALL ecriregs(84,pctsrf(1,is_sic))
3591         CALL ecriregs(84,zxfluxu(1,1))
3592         CALL ecriregs(84,zxfluxv(1,1))
3593         CALL ecriregs(84,ue)
3594         CALL ecriregs(84,ve)
3595         CALL ecriregs(84,uq)
3596         CALL ecriregs(84,vq)
3597c
3598         CALL ecrirega(84,u_seri)
3599         CALL ecrirega(84,v_seri)
3600         CALL ecrirega(84,omega)
3601         CALL ecrirega(84,t_seri)
3602         CALL ecrirega(84,zphi)
3603         CALL ecrirega(84,q_seri)
3604         CALL ecrirega(84,cldfra)
3605         CALL ecrirega(84,cldliq)
3606         CALL ecrirega(84,pplay)
3607
3608
3609cc         CALL ecrirega(84,d_t_dyn)
3610cc         CALL ecrirega(84,d_q_dyn)
3611cc         CALL ecrirega(84,heat)
3612cc         CALL ecrirega(84,cool)
3613cc         CALL ecrirega(84,d_t_con)
3614cc         CALL ecrirega(84,d_q_con)
3615cc         CALL ecrirega(84,d_t_lsc)
3616cc         CALL ecrirega(84,d_q_lsc)
3617      ENDIF
3618c
3619c Convertir les incrementations en tendances
3620c
3621      DO k = 1, klev
3622      DO i = 1, klon
3623         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3624         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3625         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3626         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3627         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3628      ENDDO
3629      ENDDO
3630c
3631      IF (nqmax.GE.3) THEN
3632      DO iq = 3, nqmax
3633      DO  k = 1, klev
3634      DO  i = 1, klon
3635         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3636      ENDDO
3637      ENDDO
3638      ENDDO
3639      ENDIF
3640c
3641c Sauvegarder les valeurs de t et q a la fin de la physique:
3642c
3643      DO k = 1, klev
3644      DO i = 1, klon
3645         t_ancien(i,k) = t_seri(i,k)
3646         q_ancien(i,k) = q_seri(i,k)
3647      ENDDO
3648      ENDDO
3649c
3650c====================================================================
3651c Si c'est la fin, il faut conserver l'etat de redemarrage
3652c====================================================================
3653c
3654      IF (lafin) THEN
3655         itau_phy = itau_phy + itap
3656ccc         IF (ok_oasis) CALL quitcpl
3657         CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
3658     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
3659     .      falbe, fevap, rain_fall, snow_fall,
3660     .      solsw, sollwdown,dlw,
3661     .      radsol,frugs,agesno,
3662     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3663     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
3664      ENDIF
3665
3666      RETURN
3667      END
3668      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3669      IMPLICIT none
3670c
3671c Calculer et imprimer l'eau totale. A utiliser pour verifier
3672c la conservation de l'eau
3673c
3674#include "YOMCST.h"
3675      INTEGER klon,klev
3676      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3677      REAL aire(klon)
3678      REAL qtotal, zx, qcheck
3679      INTEGER i, k
3680c
3681      zx = 0.0
3682      DO i = 1, klon
3683         zx = zx + aire(i)
3684      ENDDO
3685      qtotal = 0.0
3686      DO k = 1, klev
3687      DO i = 1, klon
3688         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3689     .                     *(paprs(i,k)-paprs(i,k+1))/RG
3690      ENDDO
3691      ENDDO
3692c
3693      qcheck = qtotal/zx
3694c
3695      RETURN
3696      END
3697      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3698      IMPLICIT none
3699c
3700c Tranformer une variable de la grille physique a
3701c la grille d'ecriture
3702c
3703      INTEGER nfield,nlon,iim,jjmp1, jjm
3704      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3705c
3706      INTEGER i, n, ig
3707c
3708      jjm = jjmp1 - 1
3709      DO n = 1, nfield
3710         DO i=1,iim
3711            ecrit(i,n) = fi(1,n)
3712            ecrit(i+jjm*iim,n) = fi(nlon,n)
3713         ENDDO
3714         DO ig = 1, nlon - 2
3715           ecrit(iim+ig,n) = fi(1+ig,n)
3716         ENDDO
3717      ENDDO
3718      RETURN
3719      END
3720
Note: See TracBrowser for help on using the repository browser.