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

Last change on this file since 120 was 120, checked in by lmdzadmin, 24 years ago

* empty log message *

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 99.7 KB
Line 
1      SUBROUTINE physiq (nlon,nlev,nqmax  ,
2     .            debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
3     .            paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
4     .            u,v,t,qx,
5     .            omega,
6     .            d_u, d_v, d_t, d_qx, d_ps)
7      USE ioipsl
8      IMPLICIT none
9c======================================================================
10c
11c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
12c
13c Objet: Moniteur general de la physique du modele
14cAA      Modifications quant aux traceurs :
15cAA                  -  uniformisation des parametrisations ds phytrac
16cAA                  -  stockage des moyennes des champs necessaires
17cAA                     en mode traceur off-line
18c======================================================================
19c    modif   ( P. Le Van ,  12/10/98 )
20c
21c  Arguments:
22c
23c nlon----input-I-nombre de points horizontaux
24c nlev----input-I-nombre de couches verticales
25c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
26c debut---input-L-variable logique indiquant le premier passage
27c lafin---input-L-variable logique indiquant le dernier passage
28c rjour---input-R-numero du jour de l'experience
29c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
30c pdtphys-input-R-pas d'integration pour la physique (seconde)
31c paprs---input-R-pression pour chaque inter-couche (en Pa)
32c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
33c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
34c pphis---input-R-geopotentiel du sol
35c paire---input-R-aire de chaque maille
36c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
37c u-------input-R-vitesse dans la direction X (de O a E) en m/s
38c v-------input-R-vitesse Y (de S a N) en m/s
39c t-------input-R-temperature (K)
40c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
41c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
42c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
43c omega---input-R-vitesse verticale en Pa/s
44c
45c d_u-----output-R-tendance physique de "u" (m/s/s)
46c d_v-----output-R-tendance physique de "v" (m/s/s)
47c d_t-----output-R-tendance physique de "t" (K/s)
48c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
49c d_ps----output-R-tendance physique de la pression au sol
50c======================================================================
51#include "dimensions.h"
52      integer jjmp1
53      parameter (jjmp1=jjm+1-1/jjm)
54#include "dimphy.h"
55#include "regdim.h"
56#include "indicesol.h"
57#include "dimsoil.h"
58#include "clesphys.h"
59#include "control.h"
60#include "temps.h"
61c======================================================================
62      LOGICAL check ! Verifier la conservation du modele en eau
63      PARAMETER (check=.FALSE.)
64      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
65      PARAMETER (ok_stratus=.FALSE.)
66c======================================================================
67c Parametres lies au coupleur OASIS:
68#include "oasis.h"
69      INTEGER npas, nexca, itimestep
70      logical rnpb
71      parameter(rnpb=.true.)
72c      PARAMETER (npas=1440)
73c      PARAMETER (nexca=48)
74c      PARAMETER (itimestep=1800)
75      EXTERNAL fromcpl, intocpl, inicma
76      REAL cpl_sst(iim,jjmp1), cpl_sic(iim,jjmp1)
77      REAL cpl_alb_sst(iim,jjmp1), cpl_alb_sic(iim,jjmp1)
78      character *6 ocean
79      parameter (ocean = 'force ')
80c======================================================================
81c ok_ocean indique l'utilisation du modele oceanique "slab ocean",
82c il faut bien sur s'assurer que le bilan energetique de reference
83c a la surface de l'ocean est bien present dans le fichier des
84c conditions aux limites, ainsi que l'indicateur du sol ne contient
85c pas de glace oceanique (pas de valeur 3).
86c
87      LOGICAL ok_ocean
88      PARAMETER (ok_ocean=.FALSE.)
89      REAL cyang  ! capacite thermique de l'ocean superficiel
90      PARAMETER (cyang=30.0 * 4.228e+06)
91      REAL cbing  ! capacite thermique de la glace oceanique
92      PARAMETER (cbing=1.0 * 4.228e+06)
93      REAL cthermiq
94c======================================================================
95c Clef controlant l'activation du cycle diurne:
96ccc      LOGICAL cycle_diurne
97ccc      PARAMETER (cycle_diurne=.FALSE.)
98c======================================================================
99c Modele thermique du sol, a activer pour le cycle diurne:
100ccc      LOGICAL soil_model
101ccc      PARAMETER (soil_model=.FALSE.)
102      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
103      SAVE soilcap, soilflux
104      logical ok_veget
105      parameter (ok_veget = .false.)
106c======================================================================
107c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
108c le calcul du rayonnement est celle apres la precipitation des nuages.
109c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
110c la condensation et la precipitation. Cette cle augmente les impacts
111c radiatifs des nuages.
112ccc      LOGICAL new_oliq
113ccc      PARAMETER (new_oliq=.FALSE.)
114c======================================================================
115c Clefs controlant deux parametrisations de l'orographie:
116cc      LOGICAL ok_orodr
117ccc      PARAMETER (ok_orodr=.FALSE.)
118ccc      LOGICAL ok_orolf
119ccc      PARAMETER (ok_orolf=.FALSE.)
120c======================================================================
121      LOGICAL ok_journe ! sortir le fichier journalier
122      PARAMETER (ok_journe=.FALSE.)
123c
124      LOGICAL ok_mensuel ! sortir le fichier mensuel
125      PARAMETER (ok_mensuel=.TRUE.)
126c
127      LOGICAL ok_instan ! sortir le fichier instantane
128      PARAMETER (ok_instan=.false.)
129c
130      LOGICAL ok_region ! sortir le fichier regional
131      PARAMETER (ok_region=.FALSE.)
132c======================================================================
133c
134      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
135      PARAMETER (ivap=1)
136      INTEGER iliq          ! indice de traceurs pour eau liquide
137      PARAMETER (iliq=2)
138
139      INTEGER nvm           ! nombre de vegetations
140      PARAMETER (nvm=8)
141      REAL veget(klon,nvm)  ! couverture vegetale
142      SAVE veget
143
144c
145c
146c Variables argument:
147c
148      INTEGER nlon
149      INTEGER nlev
150      INTEGER nqmax
151      REAL rjourvrai, rjour_ecri
152      REAL gmtime
153      REAL pdtphys
154      LOGICAL debut, lafin
155      REAL paprs(klon,klev+1)
156      REAL pplay(klon,klev)
157      REAL pphi(klon,klev)
158      REAL pphis(klon)
159      REAL paire(klon)
160      REAL presnivs(klev)
161      REAL znivsig(klev)
162      REAL zsurf(nbsrf)
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_u_dyn(klon,klev)
175      REAL d_v_dyn(klon,klev)
176      REAL d_t_dyn(klon,klev)
177      REAL d_q_dyn(klon,klev)
178
179      REAL omega(klon,klev)
180
181      REAL d_u(klon,klev)
182      REAL d_v(klon,klev)
183      REAL d_t(klon,klev)
184      REAL d_qx(klon,klev,nqmax)
185      REAL d_ps(klon)
186
187      INTEGER        longcles
188      PARAMETER    ( longcles = 20 )
189      REAL clesphy0( longcles      )
190c
191c Variables quasi-arguments
192c
193      REAL xjour
194      SAVE xjour
195c
196c
197c Variables propres a la physique
198c
199      REAL dtime
200      SAVE dtime                  ! pas temporel de la physique
201c
202      INTEGER radpas
203      SAVE radpas                 ! frequence d'appel rayonnement
204c
205      REAL radsol(klon)
206      SAVE radsol                 ! bilan radiatif au sol
207c
208      REAL rlat(klon)
209      SAVE rlat                   ! latitude pour chaque point
210c
211      REAL rlon(klon)
212      SAVE rlon                   ! longitude pour chaque point
213c
214cc      INTEGER iflag_con
215cc      SAVE iflag_con              ! indicateur de la convection
216c
217      INTEGER itap
218      SAVE itap                   ! compteur pour la physique
219c
220      REAL co2_ppm
221      SAVE co2_ppm                ! concentration du CO2
222c
223      REAL solaire
224      SAVE solaire                ! constante solaire
225c
226      REAL ftsol(klon,nbsrf)
227      SAVE ftsol                  ! temperature du sol
228c
229      REAL ftsoil(klon,nsoilmx,nbsrf)
230      SAVE ftsoil                 ! temperature dans le sol
231c
232      REAL fevap(klon,nbsrf)
233      SAVE fevap                 ! evaporation
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
246c
247      REAL rugmer(klon)
248      SAVE rugmer                 ! longeur de rugosite sur mer (m)
249c
250c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
251c
252      REAL zmea(klon)
253      SAVE zmea                   ! orographie moyenne
254c
255      REAL zstd(klon)
256      SAVE zstd                   ! deviation standard de l'OESM
257c
258      REAL zsig(klon)
259      SAVE zsig                   ! pente de l'OESM
260c
261      REAL zgam(klon)
262      save zgam                   ! anisotropie de l'OESM
263c
264      REAL zthe(klon)
265      SAVE zthe                   ! orientation de l'OESM
266c
267      REAL zpic(klon)
268      SAVE zpic                   ! Maximum de l'OESM
269c
270      REAL zval(klon)
271      SAVE zval                   ! Minimum de l'OESM
272c
273      REAL rugoro(klon)
274      SAVE rugoro                 ! longueur de rugosite de l'OESM
275c
276      REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
277c
278      REAL zuthe(klon),zvthe(klon)
279      SAVE zuthe
280      SAVE zvthe
281      INTEGER igwd,igwdim,idx(klon),itest(klon)
282c
283      REAL agesno(klon)
284      SAVE agesno                 ! age de la neige
285c
286      REAL alb_neig(klon)
287      SAVE alb_neig               ! albedo de la neige
288c
289c Variables locales:
290c
291      REAL cdragh(klon) ! drag coefficient pour T and Q
292      REAL cdragm(klon) ! drag coefficient pour vent
293cAA
294cAA  Pour phytrac
295cAA
296      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
297      REAL yu1(klon)            ! vents dans la premiere couche U
298      REAL yv1(klon)            ! vents dans la premiere couche V
299      LOGICAL offline           ! Controle du stockage ds "physique"
300      PARAMETER (offline=.false.)
301      INTEGER physid
302      REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
303      save pfrac_impa
304      REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
305      save pfrac_nucl
306      REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
307      save pfrac_1nucl
308      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
309      REAL frac_nucl(klon,klev) ! idem (nucleation)
310cAA
311      REAL rain_fall(klon) ! pluie
312      REAL snow_fall(klon) ! neige
313      REAL evap(klon), devap(klon) ! evaporation et sa derivee
314      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
315      REAL bils(klon) ! bilan de chaleur au sol
316      REAL fder(klon) ! Derive de flux (sensible et latente)
317      REAL ruis(klon) ! ruissellement
318      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
319      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
320      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
321      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
322c
323      REAL frugs(klon,nbsrf) ! longueur de rugosite
324      REAL zxrugs(klon) ! longueur de rugosite
325c
326c Conditions aux limites
327c
328      INTEGER julien
329      INTEGER idayvrai
330      SAVE idayvrai
331c
332      INTEGER lmt_pas
333      SAVE lmt_pas                ! frequence de mise a jour
334      REAL pctsrf(klon,nbsrf)
335      SAVE pctsrf                 ! sous-fraction du sol
336      REAL lmt_sst(klon)
337      SAVE lmt_sst                ! temperature de la surface ocean
338      REAL lmt_bils(klon)
339      SAVE lmt_bils               ! bilan de chaleur au sol
340      REAL lmt_alb(klon)
341      SAVE lmt_alb                ! temperature de la surface ocean
342      REAL lmt_rug(klon)
343      SAVE lmt_rug                ! longueur de rugosite
344      REAL alb_eau(klon)
345      SAVE alb_eau                ! albedo sur l'ocean
346      REAL albsol(klon)
347      SAVE albsol                 ! albedo du sol total
348      REAL wo(klon,klev)
349      SAVE wo                     ! ozone
350c======================================================================
351c
352c Declaration des procedures appelees
353c
354      EXTERNAL angle     ! calculer angle zenithal du soleil
355      EXTERNAL alboc     ! calculer l'albedo sur ocean
356      EXTERNAL albsno    ! calculer albedo sur neige
357      EXTERNAL ajsec     ! ajustement sec
358      EXTERNAL clmain    ! couche limite
359      EXTERNAL condsurf  ! lire les conditions aux limites
360      EXTERNAL conlmd    ! convection (schema LMD)
361      EXTERNAL diagcld   ! nuages diagnostiques
362      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
363cAA
364      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
365c                          ! stockage des coefficients necessaires au
366c                          ! lessivage OFF-LINE et ON-LINE
367cAA
368      EXTERNAL hgardfou  ! verifier les temperatures
369      EXTERNAL hydrol    ! hydrologie du sol
370      EXTERNAL nuage     ! calculer les proprietes radiatives
371      EXTERNAL o3cm      ! initialiser l'ozone
372      EXTERNAL orbite    ! calculer l'orbite terrestre
373      EXTERNAL ozonecm   ! prescrire l'ozone
374      EXTERNAL phyetat0  ! lire l'etat initial de la physique
375      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
376      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
377      EXTERNAL suphec    ! initialiser certaines constantes
378      EXTERNAL transp    ! transport total de l'eau et de l'energie
379      EXTERNAL ecribina  ! ecrire le fichier binaire global
380      EXTERNAL ecribins  ! ecrire le fichier binaire global
381      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
382      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
383c
384c Variables locales
385c
386      REAL dialiq(klon,klev)  ! eau liquide nuageuse
387      REAL diafra(klon,klev)  ! fraction nuageuse
388      REAL cldliq(klon,klev)  ! eau liquide nuageuse
389      REAL cldfra(klon,klev)  ! fraction nuageuse
390      REAL cldtau(klon,klev)  ! epaisseur optique
391      REAL cldemi(klon,klev)  ! emissivite infrarouge
392c
393C§§§ PB
394      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
395      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
396      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
397      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
398c
399      REAL zxfluxt(klon, klev)
400      REAL zxfluxq(klon, klev)
401      REAL zxfluxu(klon, klev)
402      REAL zxfluxv(klon, klev)
403C§§§
404      REAL heat(klon,klev)    ! chauffage solaire
405      REAL heat0(klon,klev)   ! chauffage solaire ciel clair
406      REAL cool(klon,klev)    ! refroidissement infrarouge
407      REAL cool0(klon,klev)   ! refroidissement infrarouge ciel clair
408      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
409      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
410      REAL albpla(klon)
411c Le rayonnement n'est pas calcule tous les pas, il faut donc
412c                      sauvegarder les sorties du rayonnement
413      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw
414      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
415      INTEGER itaprad
416      SAVE itaprad
417c
418      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
419      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
420c
421      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
422      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
423c
424      REAL zx_alb_lic, zx_alb_oce, zx_alb_ter, zx_alb_sic
425      REAL zxtsol(klon), zxqsol(klon), zxsnow(klon)
426c
427      REAL dist, rmu0(klon), fract(klon)
428      REAL zdtime, zlongi
429c
430      CHARACTER*2 str2
431      CHARACTER*2 iqn
432c
433      REAL qcheck
434      REAL z_avant(klon), z_apres(klon), z_factor(klon)
435      LOGICAL zx_ajustq
436c
437      REAL za, zb
438      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
439      INTEGER i, k, iq, ig, j, nsrf, ll
440      REAL t_coup
441      PARAMETER (t_coup=234.0)
442c
443      REAL zphi(klon,klev)
444      REAL zx_relief(iim,jjmp1)
445      REAL zx_aire(iim,jjmp1)
446c
447c Variables du changement
448c
449c con: convection
450c lsc: condensation a grande echelle (Large-Scale-Condensation)
451c ajs: ajustement sec
452c eva: evaporation de l'eau liquide nuageuse
453c vdf: couche limite (Vertical DiFfusion)
454      REAL d_t_con(klon,klev),d_q_con(klon,klev)
455      REAL d_u_con(klon,klev),d_v_con(klon,klev)
456      REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
457      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
458      REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
459      REAL rneb(klon,klev)
460c
461      REAL pmfu(klon,klev), pmfd(klon,klev)
462      REAL pen_u(klon,klev), pen_d(klon,klev)
463      REAL pde_u(klon,klev), pde_d(klon,klev)
464      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
465      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
466      REAL prfl(klon,klev+1), psfl(klon,klev+1)
467c
468      INTEGER ibas_con(klon), itop_con(klon)
469      REAL rain_con(klon), rain_lsc(klon)
470      REAL snow_con(klon), snow_lsc(klon)
471      REAL d_ts(klon,nbsrf)
472c
473      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
474      REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
475c
476      REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
477      REAL d_t_oro(klon,klev)
478      REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
479      REAL d_t_lif(klon,klev)
480c
481c Variables liees a l'ecriture de la bande histoire physique
482c
483      INTEGER ecrit_mth
484      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
485c
486      INTEGER ecrit_day
487      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
488c
489      INTEGER ecrit_ins
490      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
491c
492      INTEGER ecrit_reg
493      SAVE ecrit_reg   ! frequence d'ecriture
494c
495      REAL oas_sols(klon), z_sols(iim,jjmp1)
496      SAVE oas_sols
497      REAL oas_nsol(klon), z_nsol(iim,jjmp1)
498      SAVE oas_nsol
499      REAL oas_rain(klon), z_rain(iim,jjmp1)
500      SAVE oas_rain
501      REAL oas_snow(klon), z_snow(iim,jjmp1)
502      SAVE oas_snow
503      REAL oas_evap(klon), z_evap(iim,jjmp1)
504      SAVE oas_evap
505      REAL oas_ruis(klon), z_ruis(iim,jjmp1)
506      SAVE oas_ruis
507      REAL oas_tsol(klon), z_tsol(iim,jjmp1)
508      SAVE oas_tsol
509      REAL oas_fder(klon), z_fder(iim,jjmp1)
510      SAVE oas_fder
511      REAL oas_albe(klon), z_albe(iim,jjmp1)
512      SAVE oas_albe
513      REAL oas_taux(klon), z_taux(iim,jjmp1)
514      SAVE oas_taux
515      REAL oas_tauy(klon), z_tauy(iim,jjmp1)
516      SAVE oas_tauy
517      REAL oas_ruisoce(klon), z_ruisoce(iim,jjmp1)
518      SAVE oas_ruisoce
519      REAL oas_ruisriv(klon), z_ruisriv(iim,jjmp1)
520      SAVE oas_ruisriv
521c
522c
523c Variables locales pour effectuer les appels en serie
524c
525      REAL t_seri(klon,klev), q_seri(klon,klev)
526      REAL ql_seri(klon,klev)
527      REAL u_seri(klon,klev), v_seri(klon,klev)
528c
529      REAL tr_seri(klon,klev,nbtr)
530      REAL d_tr(klon,klev,nbtr)
531      REAL source_tr(klon,nbtr)
532
533      REAL zx_rh(klon,klev)
534      REAL dtimeday,dtimecri,dtimexp9,fecri_pas,fecri86400,fecritday
535
536      INTEGER        length
537      PARAMETER    ( length = 100 )
538      REAL tabcntr0( length       )
539c
540      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
541      REAL zx_tmp_fi2d(klon)
542      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
543      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
544c
545      INTEGER nid_day, nid_mth, nid_ins
546      SAVE nid_day, nid_mth, nid_ins
547c
548      INTEGER nhori, nvert,nvertsf
549      REAL zsto, zout, zjulian
550
551      character*20 modname
552      character*80 abort_message
553      logical ok_sync
554
555c
556c Declaration des constantes et des fonctions thermodynamiques
557c
558#include "YOMCST.h"
559#include "YOETHF.h"
560#include "FCTTRE.h"
561c======================================================================
562      modname = 'physiq'
563      ok_sync=.TRUE.
564      IF (nqmax .LT. 2) THEN
565         PRINT*, 'eaux vapeur et liquide sont indispensables'
566         CALL ABORT
567      ENDIF
568      IF (debut) THEN
569         CALL suphec ! initialiser constantes et parametres phys.
570      ENDIF
571
572
573c======================================================================
574      xjour = rjourvrai
575c
576c Si c'est le debut, il faut initialiser plusieurs choses
577c          ********
578c
579       IF (debut) THEN
580
581         DO k = 2, nvm          ! pas de vegetation
582            DO i = 1, klon
583               veget(i,k) = 0.0
584            ENDDO
585         ENDDO
586         DO i = 1, klon
587            veget(i,1) = 1.0    ! il n'y a que du desert
588         ENDDO
589         PRINT*, 'Pas de vegetation; desert partout'
590c
591c
592c Initialiser les compteurs:
593c
594
595         itap    = 0
596         itaprad = 0
597c
598         CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
599     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
600     .       falbe, fevap, rain_fall,snow_fall,sollw, solsw,
601     .       radsol,rugmer,agesno,clesphy0,
602     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
603     .       t_ancien, q_ancien, ancien_ok )
604
605c
606         radpas = NINT( 86400./dtime/nbapp_rad)
607
608c
609         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
610     ,                    ok_instan, ok_region )
611c
612         IF (ABS(dtime-pdtphys).GT.0.001) THEN
613            PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
614            abort_message=' See above '
615            call abort_gcm(modname,abort_message,1)
616         ENDIF
617         IF (nlon .NE. klon) THEN
618            PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
619            abort_message=' See above '
620            call abort_gcm(modname,abort_message,1)
621         ENDIF
622         IF (nlev .NE. klev) THEN
623            PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
624            abort_message=' See above '
625            call abort_gcm(modname,abort_message,1)
626         ENDIF
627c
628         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
629           PRINT*, 'Nbre d appels au rayonnement insuffisant'
630           PRINT*, "Au minimum 4 appels par jour si cycle diurne"
631           abort_message=' See above '
632           call abort_gcm(modname,abort_message,1)
633         ENDIF
634         PRINT*, "Clef pour la convection, iflag_con=", iflag_con
635c
636         IF (ok_orodr) THEN
637         DO i=1,klon
638         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
639         ENDDO
640         CALL SUGWD(klon,klev,paprs,pplay)
641         DO i=1,klon
642         zuthe(i)=0.
643         zvthe(i)=0.
644         if(zstd(i).gt.10.)then
645           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
646           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
647         endif
648         ENDDO
649         ENDIF
650c
651c
652         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
653         PRINT*,'La frequence de lecture surface est de ', lmt_pas
654c
655         ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
656         IF (ok_mensuel) THEN
657         PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
658         ENDIF
659         ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
660         IF (ok_journe) THEN
661         PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
662         ENDIF
663ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
664ccc         ecrit_ins = NINT(86400./dtime *0.25)  ! 4 fois par jour
665         ecrit_ins = NINT(86400./dtime/48.)  ! a chaque pas de temps
666         IF (ok_instan) THEN
667         PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
668         ENDIF
669         ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
670         IF (ok_region) THEN
671         PRINT*, 'La frequence de sortie region est de ', ecrit_reg
672         ENDIF
673
674c
675c Initialiser le couplage si necessaire
676c
677      npas = 0
678      nexca = 0
679      if (ocean == 'couple') then
680        npas = itaufin/ iphysiq
681        nexca = 86400 / dtime
682        write(*,*)' ##### Ocean couple #####'
683        write(*,*)' Valeurs des pas de temps'
684        write(*,*)' npas = ', npas
685        write(*,*)' nexca = ', nexca
686      endif       
687c
688c
689      IF (ok_journe) THEN
690c
691         CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
692         zjulian = zjulian + day_ini
693c
694         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
695         DO i = 1, iim
696            zx_lon(i,1) = rlon(i+1)
697            zx_lon(i,jjmp1) = rlon(i+1)
698         ENDDO
699         DO ll=1,klev
700            znivsig(ll)=float(ll)
701         ENDDO
702         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
703         CALL histbeg("histday", iim,zx_lon, jjmp1,zx_lat,
704     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
705     .                 nhori, nid_day)
706c         CALL histvert(nid_day, "presnivs", "Vertical levels", "mb",
707c     .                 klev, presnivs, nvert)
708         call histvert(nid_day, 'sig_s', 'Niveaux sigma','-',
709     .              klev, znivsig, nvert)
710c
711         zsto = dtime
712         zout = dtime * FLOAT(ecrit_day)
713c
714         CALL histdef(nid_day, "phis", "Surface geop. height", "-",
715     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
716     .                "once", zsto,zout)
717c
718         CALL histdef(nid_day, "aire", "Grid area", "-",
719     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
720     .                "once", zsto,zout)
721c
722c Champs 2D:
723c
724         CALL histdef(nid_day, "tsol", "Surface Temperature", "K",
725     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
726     .                "ave(X)", zsto,zout)
727c
728         CALL histdef(nid_day, "tter", "Surface Temperature", "K",
729     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
730     .                "ave(X)", zsto,zout)
731c
732         CALL histdef(nid_day, "tlic", "Surface Temperature", "K",
733     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
734     .                "ave(X)", zsto,zout)
735c
736         CALL histdef(nid_day, "toce", "Surface Temperature", "K",
737     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
738     .                "ave(X)", zsto,zout)
739c
740         CALL histdef(nid_day, "tsic", "Surface Temperature", "K",
741     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
742     .                "ave(X)", zsto,zout)
743c
744         CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
745     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
746     .                "ave(X)", zsto,zout)
747c
748         CALL histdef(nid_day, "rain", "Precipitation", "mm/day",
749     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
750     .                "ave(X)", zsto,zout)
751c
752         CALL histdef(nid_day, "snow", "Snow fall", "mm/day",
753     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
754     .                "ave(X)", zsto,zout)
755c
756         CALL histdef(nid_day, "evap", "Evaporation", "mm/day",
757     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
758     .                "ave(X)", zsto,zout)
759c
760         CALL histdef(nid_day, "tops", "Solar rad. at TOA", "W/m2",
761     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
762     .                "ave(X)", zsto,zout)
763c
764         CALL histdef(nid_day, "topl", "IR rad. at TOA", "W/m2",
765     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
766     .                "ave(X)", zsto,zout)
767c
768         CALL histdef(nid_day, "sols", "Solar rad. at surf.", "W/m2",
769     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
770     .                "ave(X)", zsto,zout)
771c
772         CALL histdef(nid_day, "soll", "IR rad. at surface", "W/m2",
773     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
774     .                "ave(X)", zsto,zout)
775c
776         CALL histdef(nid_day, "bils", "Surf. total heat flux", "W/m2",
777     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
778     .                "ave(X)", zsto,zout)
779c
780         CALL histdef(nid_day, "sens", "Sensible heat flux", "W/m2",
781     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
782     .                "ave(X)", zsto,zout)
783c
784         CALL histdef(nid_day, "fder", "Heat flux derivation", "W/m2",
785     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
786     .                "ave(X)", zsto,zout)
787c
788         CALL histdef(nid_day, "frtu", "Zonal wind stress", "Pa",
789     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
790     .                "ave(X)", zsto,zout)
791c
792         CALL histdef(nid_day, "frtv", "Meridional wind stress", "Pa",
793     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
794     .                "ave(X)", zsto,zout)
795c
796C §§§ PB flux pour chauqe sous surface
797C
798         DO nsrf = 1, nbsrf
799C
800           call histdef(nid_day, "pourc_"//clnsurf(nsrf),
801     $         "Fraction"//clnsurf(nsrf), "W/m2", 
802     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
803     $         "ave(X)", zsto,zout)
804
805           call histdef(nid_day, "sens_"//clnsurf(nsrf),
806     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
807     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
808     $         "ave(X)", zsto,zout)
809c
810           call histdef(nid_day, "lat_"//clnsurf(nsrf),
811     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
812     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
813     $         "ave(X)", zsto,zout)
814C
815           call histdef(nid_day, "taux_"//clnsurf(nsrf),
816     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
817     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
818     $         "ave(X)", zsto,zout)
819
820           call histdef(nid_day, "tauy_"//clnsurf(nsrf),
821     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
822     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
823     $         "ave(X)", zsto,zout)
824C§§§
825         END DO
826           
827         CALL histdef(nid_day, "ruis", "Runoff", "mm/day",
828     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
829     .                "ave(X)", zsto,zout)
830c
831         CALL histdef(nid_day, "sicf", "Sea-ice fraction", "-",
832     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
833     .                "ave(X)", zsto,zout)
834c
835         CALL histdef(nid_day, "cldl", "Low-level cloudiness", "-",
836     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
837     .                "ave(X)", zsto,zout)
838c
839         CALL histdef(nid_day, "cldm", "Mid-level cloudiness", "-",
840     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
841     .                "ave(X)", zsto,zout)
842c
843         CALL histdef(nid_day, "cldh", "High-level cloudiness", "-",
844     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
845     .                "ave(X)", zsto,zout)
846c
847         CALL histdef(nid_day, "cldt", "Total cloudiness", "-",
848     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
849     .                "ave(X)", zsto,zout)
850c
851         CALL histdef(nid_day, "cldq", "Cloud liquid water path", "-",
852     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
853     .                "ave(X)", zsto,zout)
854c
855c Champs 3D:
856c
857         CALL histdef(nid_day, "temp", "Air temperature", "K",
858     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
859     .                "ave(X)", zsto,zout)
860c
861         CALL histdef(nid_day, "ovap", "Specific humidity", "Kg/Kg",
862     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
863     .                "ave(X)", zsto,zout)
864c
865         CALL histdef(nid_day, "geop", "Geopotential height", "m",
866     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
867     .                "ave(X)", zsto,zout)
868c
869         CALL histdef(nid_day, "vitu", "Zonal wind", "m/s",
870     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
871     .                "ave(X)", zsto,zout)
872c
873         CALL histdef(nid_day, "vitv", "Meridional wind", "m/s",
874     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
875     .                "ave(X)", zsto,zout)
876c
877         CALL histdef(nid_day, "vitw", "Vertical wind", "m/s",
878     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
879     .                "ave(X)", zsto,zout)
880c
881         CALL histdef(nid_day, "pres", "Air pressure", "Pa",
882     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
883     .                "ave(X)", zsto,zout)
884c
885         CALL histend(nid_day)
886c
887         ndex2d = 0
888         ndex3d = 0
889c
890      ENDIF ! fin de test sur ok_journe
891c
892      IF (ok_mensuel) THEN
893c
894         CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
895         zjulian = zjulian + day_ini
896c
897         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
898         DO i = 1, iim
899            zx_lon(i,1) = rlon(i+1)
900            zx_lon(i,jjmp1) = rlon(i+1)
901         ENDDO
902         DO ll=1,klev
903            znivsig(ll)=float(ll)
904         ENDDO
905         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
906         CALL histbeg("histmth", iim,zx_lon, jjmp1,zx_lat,
907     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
908     .                 nhori, nid_mth)
909c         CALL histvert(nid_mth, "presnivs", "Vertical levels", "mb",
910c     .                 klev, presnivs, nvert)
911         call histvert(nid_mth, 'sig_s', 'Niveaux sigma','-',
912     .              klev, znivsig, nvert)
913c
914         zsto = dtime
915         zout = dtime * ecrit_mth
916c
917         CALL histdef(nid_mth, "phis", "Surface geop. height", "-",
918     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
919     .                "once",  zsto,zout)
920c
921         CALL histdef(nid_mth, "aire", "Grid area", "-",
922     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
923     .                "once",  zsto,zout)
924c
925c Champs 2D:
926c
927         CALL histdef(nid_mth, "tsol", "Surface Temperature", "K",
928     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
929     .                "ave(X)", zsto,zout)
930c
931         CALL histdef(nid_mth, "psol", "Surface Pressure", "Pa",
932     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
933     .                "ave(X)", zsto,zout)
934c
935         CALL histdef(nid_mth, "qsol", "Surface humidity", "mm",
936     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
937     .                "ave(X)", zsto,zout)
938c
939         CALL histdef(nid_mth, "rain", "Precipitation", "mm/day",
940     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
941     .                "ave(X)", zsto,zout)
942c
943         CALL histdef(nid_mth, "plul", "Large-scale Precip.", "mm/day",
944     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
945     .                "ave(X)", zsto,zout)
946c
947         CALL histdef(nid_mth, "pluc", "Convective Precip.", "mm/day",
948     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
949     .                "ave(X)", zsto,zout)
950c
951         CALL histdef(nid_mth, "snow", "Snow fall", "mm/day",
952     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
953     .                "ave(X)", zsto,zout)
954c
955         CALL histdef(nid_mth, "ages", "Snow age", "day",
956     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
957     .                "ave(X)", zsto,zout)
958c
959         CALL histdef(nid_mth, "evap", "Evaporation", "mm/day",
960     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
961     .                "ave(X)", zsto,zout)
962c
963         CALL histdef(nid_mth, "tops", "Solar rad. at TOA", "W/m2",
964     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
965     .                "ave(X)", zsto,zout)
966c
967         CALL histdef(nid_mth, "topl", "IR rad. at TOA", "W/m2",
968     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
969     .                "ave(X)", zsto,zout)
970c
971         CALL histdef(nid_mth, "sols", "Solar rad. at surf.", "W/m2",
972     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
973     .                "ave(X)", zsto,zout)
974c
975         CALL histdef(nid_mth, "soll", "IR rad. at surface", "W/m2",
976     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
977     .                "ave(X)", zsto,zout)
978c
979         CALL histdef(nid_mth, "tops0", "Solar rad. at TOA", "W/m2",
980     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
981     .                "ave(X)", zsto,zout)
982c
983         CALL histdef(nid_mth, "topl0", "IR rad. at TOA", "W/m2",
984     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
985     .                "ave(X)", zsto,zout)
986c
987         CALL histdef(nid_mth, "sols0", "Solar rad. at surf.", "W/m2",
988     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
989     .                "ave(X)", zsto,zout)
990c
991         CALL histdef(nid_mth, "soll0", "IR rad. at surface", "W/m2",
992     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
993     .                "ave(X)", zsto,zout)
994c
995         CALL histdef(nid_mth, "bils", "Surf. total heat flux", "W/m2",
996     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
997     .                "ave(X)", zsto,zout)
998c
999         CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
1000     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1001     .                "ave(X)", zsto,zout)
1002c
1003         CALL histdef(nid_mth, "fder", "Heat flux derivation", "W/m2",
1004     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1005     .                "ave(X)", zsto,zout)
1006c
1007         CALL histdef(nid_mth, "frtu", "Zonal wind stress", "Pa",
1008     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1009     .                "ave(X)", zsto,zout)
1010c
1011         CALL histdef(nid_mth, "frtv", "Meridional wind stress", "Pa",
1012     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1013     .                "ave(X)", zsto,zout)
1014c
1015         DO nsrf = 1, nbsrf
1016C
1017           call histdef(nid_mth, "pourc_"//clnsurf(nsrf),
1018     $         "Fraction "//clnsurf(nsrf), "W/m2", 
1019     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1020     $         "ave(X)", zsto,zout)
1021C
1022           call histdef(nid_mth, "sens_"//clnsurf(nsrf),
1023     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
1024     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1025     $         "ave(X)", zsto,zout)
1026c
1027           call histdef(nid_mth, "lat_"//clnsurf(nsrf),
1028     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
1029     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1030     $         "ave(X)", zsto,zout)
1031C
1032           call histdef(nid_mth, "taux_"//clnsurf(nsrf),
1033     $         "Zonal wind stress"//clnsurf(nsrf), "Pa", 
1034     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1035     $         "ave(X)", zsto,zout)
1036
1037           call histdef(nid_mth, "tauy_"//clnsurf(nsrf),
1038     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
1039     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1040     $         "ave(X)", zsto,zout)
1041         END DO
1042C
1043         CALL histdef(nid_mth, "ruis", "Runoff", "mm/day",
1044     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1045     .                "ave(X)", zsto,zout)
1046c
1047         CALL histdef(nid_mth, "sicf", "Sea-ice fraction", "-",
1048     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1049     .                "ave(X)", zsto,zout)
1050c
1051         CALL histdef(nid_mth, "albs", "Surface albedo", "-",
1052     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1053     .                "ave(X)", zsto,zout)
1054c
1055         CALL histdef(nid_mth, "cdrm", "Momentum drag coef.", "-",
1056     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1057     .                "ave(X)", zsto,zout)
1058c
1059         CALL histdef(nid_mth, "cdrh", "Heat drag coef.", "-",
1060     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1061     .                "ave(X)", zsto,zout)
1062c
1063         CALL histdef(nid_mth, "cldl", "Low-level cloudiness", "-",
1064     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1065     .                "ave(X)", zsto,zout)
1066c
1067         CALL histdef(nid_mth, "cldm", "Mid-level cloudiness", "-",
1068     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1069     .                "ave(X)", zsto,zout)
1070c
1071         CALL histdef(nid_mth, "cldh", "High-level cloudiness", "-",
1072     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1073     .                "ave(X)", zsto,zout)
1074c
1075         CALL histdef(nid_mth, "cldt", "Total cloudiness", "-",
1076     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1077     .                "ave(X)", zsto,zout)
1078c
1079         CALL histdef(nid_mth, "cldq", "Cloud liquid water path", "-",
1080     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1081     .                "ave(X)", zsto,zout)
1082c
1083         CALL histdef(nid_mth, "ue", "Zonal energy transport", "-",
1084     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1085     .                "ave(X)", zsto,zout)
1086c
1087         CALL histdef(nid_mth, "ve", "Merid energy transport", "-",
1088     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1089     .                "ave(X)", zsto,zout)
1090c
1091         CALL histdef(nid_mth, "uq", "Zonal humidity transport", "-",
1092     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1093     .                "ave(X)", zsto,zout)
1094c
1095         CALL histdef(nid_mth, "vq", "Merid humidity transport", "-",
1096     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1097     .                "ave(X)", zsto,zout)
1098c
1099c Champs 3D:
1100c
1101         CALL histdef(nid_mth, "temp", "Air temperature", "K",
1102     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1103     .                "ave(X)", zsto,zout)
1104c
1105         CALL histdef(nid_mth, "ovap", "Specific humidity", "Kg/Kg",
1106     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1107     .                "ave(X)", zsto,zout)
1108c
1109         CALL histdef(nid_mth, "geop", "Geopotential height", "m",
1110     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1111     .                "ave(X)", zsto,zout)
1112c
1113         CALL histdef(nid_mth, "vitu", "Zonal wind", "m/s",
1114     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1115     .                "ave(X)", zsto,zout)
1116c
1117         CALL histdef(nid_mth, "vitv", "Meridional wind", "m/s",
1118     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1119     .                "ave(X)", zsto,zout)
1120c
1121         CALL histdef(nid_mth, "vitw", "Vertical wind", "m/s",
1122     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1123     .                "ave(X)", zsto,zout)
1124c
1125         CALL histdef(nid_mth, "pres", "Air pressure", "Pa",
1126     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1127     .                "ave(X)", zsto,zout)
1128c
1129         CALL histdef(nid_mth, "rneb", "Cloud fraction", "-",
1130     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1131     .                "ave(X)", zsto,zout)
1132c
1133         CALL histdef(nid_mth, "rhum", "Relative humidity", "-",
1134     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1135     .                "ave(X)", zsto,zout)
1136c
1137         CALL histdef(nid_mth, "oliq", "Liquid water content", "kg/kg",
1138     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1139     .                "ave(X)", zsto,zout)
1140c
1141         CALL histdef(nid_mth, "dtdyn", "Dynamics dT", "K/s",
1142     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1143     .                "ave(X)", zsto,zout)
1144c
1145         CALL histdef(nid_mth, "dqdyn", "Dynamics dQ", "Kg/Kg/s",
1146     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1147     .                "ave(X)", zsto,zout)
1148c
1149         CALL histdef(nid_mth, "dtcon", "Convection dT", "K/s",
1150     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1151     .                "ave(X)", zsto,zout)
1152c
1153         CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
1154     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1155     .                "ave(X)", zsto,zout)
1156c
1157         CALL histdef(nid_mth, "dtlsc", "Condensation dT", "K/s",
1158     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1159     .                "ave(X)", zsto,zout)
1160c
1161         CALL histdef(nid_mth, "dqlsc", "Condensation dQ", "Kg/Kg/s",
1162     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1163     .                "ave(X)", zsto,zout)
1164c
1165         CALL histdef(nid_mth, "dtvdf", "Boundary-layer dT", "K/s",
1166     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1167     .                "ave(X)", zsto,zout)
1168c
1169         CALL histdef(nid_mth, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
1170     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1171     .                "ave(X)", zsto,zout)
1172c
1173         CALL histdef(nid_mth, "dteva", "Reevaporation dT", "K/s",
1174     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1175     .                "ave(X)", zsto,zout)
1176c
1177         CALL histdef(nid_mth, "dqeva", "Reevaporation dQ", "Kg/Kg/s",
1178     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1179     .                "ave(X)", zsto,zout)
1180c
1181         CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
1182     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1183     .                "ave(X)", zsto,zout)
1184
1185         CALL histdef(nid_mth, "dqajs", "Dry adjust. dQ", "Kg/Kg/s",
1186     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1187     .                "ave(X)", zsto,zout)
1188c
1189         CALL histdef(nid_mth, "dtswr", "SW radiation dT", "K/s",
1190     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1191     .                "ave(X)", zsto,zout)
1192c
1193         CALL histdef(nid_mth, "dtsw0", "SW radiation dT", "K/s",
1194     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1195     .                "ave(X)", zsto,zout)
1196c
1197         CALL histdef(nid_mth, "dtlwr", "LW radiation dT", "K/s",
1198     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1199     .                "ave(X)", zsto,zout)
1200c
1201         CALL histdef(nid_mth, "dtlw0", "LW radiation dT", "K/s",
1202     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1203     .                "ave(X)", zsto,zout)
1204c
1205         CALL histdef(nid_mth, "duvdf", "Boundary-layer dU", "m/s2",
1206     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1207     .                "ave(X)", zsto,zout)
1208c
1209         CALL histdef(nid_mth, "dvvdf", "Boundary-layer dV", "m/s2",
1210     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1211     .                "ave(X)", zsto,zout)
1212c
1213         IF (ok_orodr) THEN
1214         CALL histdef(nid_mth, "duoro", "Orography dU", "m/s2",
1215     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1216     .                "ave(X)", zsto,zout)
1217c
1218         CALL histdef(nid_mth, "dvoro", "Orography dV", "m/s2",
1219     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1220     .                "ave(X)", zsto,zout)
1221c
1222         ENDIF
1223C
1224         IF (ok_orolf) THEN
1225         CALL histdef(nid_mth, "dulif", "Orography dU", "m/s2",
1226     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1227     .                "ave(X)", zsto,zout)
1228c
1229         CALL histdef(nid_mth, "dvlif", "Orography dV", "m/s2",
1230     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1231     .                "ave(X)", zsto,zout)
1232         ENDIF
1233C
1234         CALL histdef(nid_mth, "ozone", "Ozone concentration", "-",
1235     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1236     .                "ave(X)", zsto,zout)
1237c
1238         if (nqmax.GE.3) THEN
1239         DO iq=1,nqmax-2
1240         IF (iq.LE.99) THEN
1241         WRITE(str2,'(i2.2)') iq
1242         CALL histdef(nid_mth, "trac"//str2, "Tracer No."//str2, "-",
1243     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1244     .                "ave(X)", zsto,zout)
1245         ELSE
1246         PRINT*, "Trop de traceurs"
1247         CALL abort
1248         ENDIF
1249         ENDDO
1250         ENDIF
1251c
1252         CALL histend(nid_mth)
1253c
1254         ndex2d = 0
1255         ndex3d = 0
1256c
1257      ENDIF ! fin de test sur ok_mensuel
1258c
1259c
1260      IF (ok_instan) THEN
1261c
1262         CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
1263         zjulian = zjulian + day_ini
1264c
1265         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
1266         DO i = 1, iim
1267            zx_lon(i,1) = rlon(i+1)
1268            zx_lon(i,jjmp1) = rlon(i+1)
1269         ENDDO
1270         DO ll=1,klev
1271            znivsig(ll)=float(ll)
1272         ENDDO
1273         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
1274         CALL histbeg("histins", iim,zx_lon, jjmp1,zx_lat,
1275     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
1276     .                 nhori, nid_ins)
1277c         CALL histvert(nid_ins, "presnivs", "Vertical levels", "mb",
1278c     .                 klev, presnivs, nvert)
1279         call histvert(nid_ins, 'sig_s', 'Niveaux sigma','-',
1280     .              klev, znivsig, nvert)
1281c
1282c
1283         zsto = dtime * ecrit_ins
1284         zout = dtime * ecrit_ins
1285C
1286         CALL histdef(nid_ins, "phis", "Surface geop. height", "-",
1287     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1288     .                "once", zsto,zout)
1289c
1290         CALL histdef(nid_ins, "aire", "Grid area", "-",
1291     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1292     .                "once", zsto,zout)
1293c
1294c Champs 2D:
1295c
1296        CALL histdef(nid_ins, "tsol", "Surface Temperature", "K",
1297     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1298     .                "inst(X)", zsto,zout)
1299c
1300        CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
1301     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1302     .                "inst(X)", zsto,zout)
1303c
1304         CALL histdef(nid_ins, "topl", "OLR", "W/m2",
1305     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1306     .                "inst(X)", zsto,zout)
1307c
1308         CALL histdef(nid_ins, "evap", "Evaporation", "mm/day",
1309     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1310     .                "inst(X)", zsto,zout)
1311c
1312         CALL histdef(nid_ins, "sols", "Solar rad. at surf.", "W/m2",
1313     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1314     .                "inst(X)", zsto,zout)
1315c
1316         CALL histdef(nid_ins, "soll", "IR rad. at surface", "W/m2",
1317     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1318     .                "inst(X)", zsto,zout)
1319c
1320         CALL histdef(nid_ins, "bils", "Surf. total heat flux", "W/m2",
1321     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1322     .                "inst(X)", zsto,zout)
1323c
1324         CALL histdef(nid_ins, "sens", "Sensible heat flux", "W/m2",
1325     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1326     .                "inst(X)", zsto,zout)
1327c
1328         CALL histdef(nid_ins, "fder", "Heat flux derivation", "W/m2",
1329     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1330     .                "inst(X)", zsto,zout)
1331c
1332      CALL histdef(nid_ins, "dtsvdfo", "Boundary-layer dTs(o)", "K/s",
1333     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1334     .                "inst(X)", zsto,zout)
1335c
1336      CALL histdef(nid_ins, "dtsvdft", "Boundary-layer dTs(t)", "K/s",
1337     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1338     .                "inst(X)", zsto,zout)
1339c
1340      CALL histdef(nid_ins, "dtsvdfg", "Boundary-layer dTs(g)", "K/s",
1341     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1342     .                "inst(X)", zsto,zout)
1343c
1344      CALL histdef(nid_ins, "dtsvdfi", "Boundary-layer dTs(g)", "K/s",
1345     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1346     .                "inst(X)", zsto,zout)
1347
1348         DO nsrf = 1, nbsrf
1349C
1350           call histdef(nid_ins, "pourc_"//clnsurf(nsrf),
1351     $         "Fraction"//clnsurf(nsrf), "W/m2", 
1352     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1353     $         "inst(X)", zsto,zout)
1354
1355           call histdef(nid_ins, "sens_"//clnsurf(nsrf),
1356     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2", 
1357     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1358     $         "inst(X)", zsto,zout)
1359c
1360           call histdef(nid_ins, "tsol_"//clnsurf(nsrf),
1361     $         "Surface Temperature"//clnsurf(nsrf), "W/m2", 
1362     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1363     $         "inst(X)", zsto,zout)
1364c
1365           call histdef(nid_ins, "lat_"//clnsurf(nsrf),
1366     $         "Latent heat flux "//clnsurf(nsrf), "W/m2", 
1367     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1368     $         "inst(X)", zsto,zout)
1369C
1370           call histdef(nid_ins, "taux_"//clnsurf(nsrf),
1371     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
1372     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1373     $         "inst(X)", zsto,zout)
1374
1375           call histdef(nid_ins, "tauy_"//clnsurf(nsrf),
1376     $         "Meridional xind stress "//clnsurf(nsrf), "Pa", 
1377     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
1378     $         "inst(X)", zsto,zout)
1379C§§§
1380         END DO
1381c
1382c Champs 3D:
1383c
1384         CALL histdef(nid_ins, "temp", "Temperature", "K",
1385     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1386     .                "inst(X)", zsto,zout)
1387c
1388         CALL histdef(nid_ins, "vitu", "Zonal wind", "m/s",
1389     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1390     .                "inst(X)", zsto,zout)
1391c
1392         CALL histdef(nid_ins, "vitv", "Merid wind", "m/s",
1393     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1394     .                "inst(X)", zsto,zout)
1395c
1396         CALL histdef(nid_ins, "geop", "Geopotential height", "m",
1397     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1398     .                "inst(X)", zsto,zout)
1399c
1400         CALL histdef(nid_ins, "pres", "Air pressure", "Pa",
1401     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1402     .                "inst(X)", zsto,zout)
1403c
1404         CALL histdef(nid_ins, "dtvdf", "Boundary-layer dT", "K/s",
1405     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1406     .                "inst(X)", zsto,zout)
1407c
1408         CALL histdef(nid_ins, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
1409     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1410     .                "inst(X)", zsto,zout)
1411c
1412
1413         CALL histend(nid_ins)
1414c
1415         ndex2d = 0
1416         ndex3d = 0
1417c
1418      ENDIF
1419c
1420c
1421c
1422c Prescrire l'ozone dans l'atmosphere
1423c
1424c
1425cc         DO i = 1, klon
1426cc         DO k = 1, klev
1427cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1428cc         ENDDO
1429cc         ENDDO
1430c
1431c
1432      ENDIF
1433c
1434c   ****************     Fin  de   IF ( debut  )   ***************
1435c
1436c
1437c Mettre a zero des variables de sortie (pour securite)
1438c
1439      DO i = 1, klon
1440         d_ps(i) = 0.0
1441      ENDDO
1442      DO k = 1, klev
1443      DO i = 1, klon
1444         d_t(i,k) = 0.0
1445         d_u(i,k) = 0.0
1446         d_v(i,k) = 0.0
1447      ENDDO
1448      ENDDO
1449      DO iq = 1, nqmax
1450      DO k = 1, klev
1451      DO i = 1, klon
1452         d_qx(i,k,iq) = 0.0
1453      ENDDO
1454      ENDDO
1455      ENDDO
1456c
1457c Ne pas affecter les valeurs entrees de u, v, h, et q
1458c
1459      DO k = 1, klev
1460      DO i = 1, klon
1461         t_seri(i,k)  = t(i,k)
1462         u_seri(i,k)  = u(i,k)
1463         v_seri(i,k)  = v(i,k)
1464         q_seri(i,k)  = qx(i,k,ivap)
1465         ql_seri(i,k) = qx(i,k,iliq)
1466      ENDDO
1467      ENDDO
1468      IF (nqmax.GE.3) THEN
1469      DO iq = 3, nqmax
1470      DO  k = 1, klev
1471      DO  i = 1, klon
1472         tr_seri(i,k,iq-2) = qx(i,k,iq)
1473      ENDDO
1474      ENDDO
1475      ENDDO
1476      ELSE
1477      DO k = 1, klev
1478      DO i = 1, klon
1479         tr_seri(i,k,1) = 0.0
1480      ENDDO
1481      ENDDO
1482      ENDIF
1483c
1484c Diagnostiquer la tendance dynamique
1485c
1486      IF (ancien_ok) THEN
1487         DO k = 1, klev
1488         DO i = 1, klon
1489            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1490            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1491         ENDDO
1492         ENDDO
1493      ELSE
1494         DO k = 1, klev
1495         DO i = 1, klon
1496            d_t_dyn(i,k) = 0.0
1497            d_q_dyn(i,k) = 0.0
1498         ENDDO
1499         ENDDO
1500         ancien_ok = .TRUE.
1501      ENDIF
1502c
1503c Ajouter le geopotentiel du sol:
1504c
1505      DO k = 1, klev
1506      DO i = 1, klon
1507         zphi(i,k) = pphi(i,k) + pphis(i)
1508      ENDDO
1509      ENDDO
1510c
1511c Verifier les temperatures
1512c
1513
1514      CALL hgardfou(t_seri,ftsol,'debutphy')
1515c
1516c Incrementer le compteur de la physique
1517c
1518      itap   = itap + 1
1519      julien = MOD(NINT(xjour),360)
1520c
1521c Mettre en action les conditions aux limites (albedo, sst, etc.).
1522c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1523c
1524      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1525         idayvrai = NINT(xjour)
1526         PRINT *,' PHYS cond  julien ',julien,idayvrai
1527c         CALL condsurf(julien,idayvrai, pctsrf ,
1528c     .                  lmt_sst,lmt_alb,lmt_rug,lmt_bils  )
1529         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1530      ENDIF
1531c
1532c Re-evaporer l'eau liquide nuageuse
1533c
1534      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1535      DO i = 1, klon
1536         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1537         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1538         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1539         zb = MAX(0.0,ql_seri(i,k))
1540         za = - MAX(0.0,ql_seri(i,k))
1541     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1542         t_seri(i,k) = t_seri(i,k) + za
1543         q_seri(i,k) = q_seri(i,k) + zb
1544         ql_seri(i,k) = 0.0
1545         d_t_eva(i,k) = za
1546         d_q_eva(i,k) = zb
1547      ENDDO
1548      ENDDO
1549c
1550c Appeler la diffusion verticale (programme de couche limite)
1551c
1552      DO i = 1, klon
1553        if (.not. ok_veget) then
1554          frugs(i,is_ter) = SQRT(frugs(i,is_ter)**2+rugoro(i)**2)
1555        endif
1556         frugs(i,is_lic) = rugoro(i)
1557         frugs(i,is_oce) = rugmer(i)
1558         frugs(i,is_sic) = 0.001
1559         zxrugs(i) = 0.0
1560      ENDDO
1561      DO nsrf = 1, nbsrf
1562      DO i = 1, klon
1563         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1564      ENDDO
1565      ENDDO
1566      DO nsrf = 1, nbsrf
1567      DO i = 1, klon
1568         zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1569      ENDDO
1570      ENDDO
1571c
1572C calculs necessaires au calcul de l'albedo dans l'interface
1573c
1574      CALL orbite(FLOAT(julien),zlongi,dist)
1575      IF (cycle_diurne) THEN
1576        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1577        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1578      ELSE
1579        rmu0 = -999.999
1580      ENDIF
1581
1582      CALL clmain(dtime,itap,pctsrf,
1583     e            t_seri,q_seri,u_seri,v_seri,
1584     e            julien, rmu0,
1585     e            ok_veget, ocean, npas, nexca, ftsol,
1586     e            paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,
1587     e            rain_fall, snow_fall, solsw, sollw, fder,
1588     e            rlon, rlat, frugs,
1589     e            debut, lafin, agesno,
1590     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1591     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,rugmer,
1592     s            dsens, devap,
1593     s            ycoefh,yu1,yv1)
1594c
1595C§§§ PB
1596C§§§ Incrementation des flux
1597C§§
1598      zxfluxt=0.
1599      zxfluxq=0.
1600      zxfluxu=0.
1601      zxfluxv=0.
1602      DO nsrf = 1, nbsrf
1603        DO k = 1, klev
1604          DO i = 1, klon
1605            zxfluxt(i,k) = zxfluxt(i,k) +
1606     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1607            zxfluxq(i,k) = zxfluxq(i,k) +
1608     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1609            zxfluxu(i,k) = zxfluxu(i,k) +
1610     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1611            zxfluxv(i,k) = zxfluxv(i,k) +
1612     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1613          END DO
1614        END DO
1615      END DO
1616      DO i = 1, klon
1617         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1618c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1619         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1620         fder(i) = dsens(i) + devap(i)
1621      ENDDO
1622
1623      DO k = 1, klev
1624      DO i = 1, klon
1625         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1626         q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1627         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1628         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1629      ENDDO
1630      ENDDO
1631c
1632c Incrementer la temperature du sol
1633c
1634      DO i = 1, klon
1635         zxtsol(i) = 0.0
1636         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1637     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
1638     $       THEN
1639             WRITE(*,*) 'physiq : pb sous surface au point ', i,
1640     $           pctsrf(i, 1 : nbsrf)
1641         ENDIF
1642      ENDDO
1643      DO nsrf = 1, nbsrf
1644      DO i = 1, klon
1645         ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1646         zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1647      ENDDO
1648      ENDDO
1649
1650c
1651c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1652c
1653      DO nsrf = 1, nbsrf
1654      DO i = 1, klon
1655         IF (pctsrf(i,nsrf).LT.epsfra) ftsol(i,nsrf) = zxtsol(i)
1656      ENDDO
1657      ENDDO
1658
1659c
1660c Calculer la derive du flux infrarouge
1661c
1662      DO nsrf = 1, nbsrf
1663      DO i = 1, klon
1664         fder(i) = fder(i) - 4.0*RSIGMA*zxtsol(i)**3 *
1665     .                       (ftsol(i,nsrf)-zxtsol(i))
1666     .                      *pctsrf(i,nsrf)
1667      ENDDO
1668      ENDDO
1669c
1670c Appeler la convection (au choix)
1671c
1672      DO k = 1, klev
1673      DO i = 1, klon
1674         conv_q(i,k) = d_q_dyn(i,k)
1675     .               + d_q_vdf(i,k)/dtime
1676         conv_t(i,k) = d_t_dyn(i,k)
1677     .               + d_t_vdf(i,k)/dtime
1678      ENDDO
1679      ENDDO
1680      IF (check) THEN
1681         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1682         PRINT*, "avantcon=", za
1683      ENDIF
1684      zx_ajustq = .FALSE.
1685      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1686      IF (zx_ajustq) THEN
1687         DO i = 1, klon
1688            z_avant(i) = 0.0
1689         ENDDO
1690         DO k = 1, klev
1691         DO i = 1, klon
1692            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1693     .                        *(paprs(i,k)-paprs(i,k+1))/RG
1694         ENDDO
1695         ENDDO
1696      ENDIF
1697      IF (iflag_con.EQ.1) THEN
1698          stop'reactiver le call conlmd dans physiq.F'
1699c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1700c    .             d_t_con, d_q_con,
1701c    .             rain_con, snow_con, ibas_con, itop_con)
1702      ELSE IF (iflag_con.EQ.2) THEN
1703      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1704     e            conv_t, conv_q, zxfluxq(1,1), omega,
1705     s            d_t_con, d_q_con, rain_con, snow_con,
1706     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1707     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
1708      DO i = 1, klon
1709         ibas_con(i) = klev+1 - kcbot(i)
1710         itop_con(i) = klev+1 - kctop(i)
1711      ENDDO
1712      ELSE IF (iflag_con.EQ.3) THEN
1713          stop'reactiver le call conlmd dans physiq.F'
1714c     CALL conccm (dtime,paprs,pplay,t_seri,q_seri,conv_q,
1715c    s             d_t_con, d_q_con,
1716c    s             rain_con, snow_con, ibas_con, itop_con)
1717      ELSE
1718      PRINT*, "iflag_con non-prevu", iflag_con
1719      CALL abort
1720      ENDIF
1721      CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1722     .              d_u_con, d_v_con)
1723      DO k = 1, klev
1724      DO i = 1, klon
1725         t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1726         q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1727         u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1728         v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1729      ENDDO
1730      ENDDO
1731      IF (check) THEN
1732         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1733         PRINT*, "aprescon=", za
1734         zx_t = 0.0
1735         za = 0.0
1736         DO i = 1, klon
1737            za = za + paire(i)/FLOAT(klon)
1738            zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1739        ENDDO
1740         zx_t = zx_t/za*dtime
1741         PRINT*, "Precip=", zx_t
1742      ENDIF
1743      IF (zx_ajustq) THEN
1744         DO i = 1, klon
1745            z_apres(i) = 0.0
1746         ENDDO
1747         DO k = 1, klev
1748         DO i = 1, klon
1749            z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1750     .                        *(paprs(i,k)-paprs(i,k+1))/RG
1751         ENDDO
1752         ENDDO
1753         DO i = 1, klon
1754         z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1755     .                /z_apres(i)
1756         ENDDO
1757         DO k = 1, klev
1758         DO i = 1, klon
1759         IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1760     .       z_factor(i).LT.(1.0-1.0E-08)) THEN
1761               q_seri(i,k) = q_seri(i,k) * z_factor(i)
1762         ENDIF
1763         ENDDO
1764         ENDDO
1765      ENDIF
1766      zx_ajustq=.FALSE.
1767c
1768      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1769c
1770      IF (iflag_con.NE.2) THEN
1771         PRINT*, "Pour l instant, seul conflx fonctionne avec traceurs"
1772         PRINT*,' Mettre iflag_con = 2  dans  run.def et repasser  !'
1773         CALL abort
1774      ENDIF
1775c
1776      ENDIF !--nqmax.GT.2
1777c
1778c Appeler l'ajustement sec
1779c
1780      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1781      DO k = 1, klev
1782      DO i = 1, klon
1783         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1784         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1785      ENDDO
1786      ENDDO
1787c
1788c Appeler le processus de condensation a grande echelle
1789c et le processus de precipitation
1790c
1791      CALL fisrtilp_tr(dtime,paprs,pplay,
1792     .           t_seri, q_seri,
1793     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1794     .           rain_lsc, snow_lsc,
1795     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
1796     .           frac_impa, frac_nucl,
1797     .           prfl, psfl)
1798      DO k = 1, klev
1799      DO i = 1, klon
1800         t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1801         q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1802         ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1803         cldfra(i,k) = rneb(i,k)
1804         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1805      ENDDO
1806      ENDDO
1807      IF (check) THEN
1808         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1809         PRINT*, "apresilp=", za
1810         zx_t = 0.0
1811         za = 0.0
1812         DO i = 1, klon
1813            za = za + paire(i)/FLOAT(klon)
1814            zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
1815        ENDDO
1816         zx_t = zx_t/za*dtime
1817         PRINT*, "Precip=", zx_t
1818      ENDIF
1819c
1820c Nuages diagnostiques:
1821c
1822      IF (iflag_con.EQ.2) THEN ! seulement pour Tiedtke
1823      CALL diagcld1(paprs,pplay,
1824     .             rain_con,snow_con,ibas_con,itop_con,
1825     .             diafra,dialiq)
1826      DO k = 1, klev
1827      DO i = 1, klon
1828      IF (diafra(i,k).GT.cldfra(i,k)) THEN
1829         cldliq(i,k) = dialiq(i,k)
1830         cldfra(i,k) = diafra(i,k)
1831      ENDIF
1832      ENDDO
1833      ENDDO
1834      ENDIF
1835c
1836c Nuages stratus artificiels:
1837c
1838      IF (ok_stratus) THEN
1839      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
1840      DO k = 1, klev
1841      DO i = 1, klon
1842      IF (diafra(i,k).GT.cldfra(i,k)) THEN
1843         cldliq(i,k) = dialiq(i,k)
1844         cldfra(i,k) = diafra(i,k)
1845      ENDIF
1846      ENDDO
1847      ENDDO
1848      ENDIF
1849c
1850c Precipitation totale
1851c
1852      DO i = 1, klon
1853         rain_fall(i) = rain_con(i) + rain_lsc(i)
1854         snow_fall(i) = snow_con(i) + snow_lsc(i)
1855      ENDDO
1856c
1857c Calculer l'humidite relative pour diagnostique
1858c
1859      DO k = 1, klev
1860      DO i = 1, klon
1861         zx_t = t_seri(i,k)
1862         IF (thermcep) THEN
1863            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1864            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1865            zx_qs  = MIN(0.5,zx_qs)
1866            zcor   = 1./(1.-retv*zx_qs)
1867            zx_qs  = zx_qs*zcor
1868         ELSE
1869           IF (zx_t.LT.t_coup) THEN
1870              zx_qs = qsats(zx_t)/pplay(i,k)
1871           ELSE
1872              zx_qs = qsatl(zx_t)/pplay(i,k)
1873           ENDIF
1874         ENDIF
1875         zx_rh(i,k) = q_seri(i,k)/zx_qs
1876      ENDDO
1877      ENDDO
1878c
1879c Calculer les parametres optiques des nuages et quelques
1880c parametres pour diagnostiques:
1881c
1882      CALL nuage (paprs, pplay,
1883     .            t_seri, cldliq, cldfra, cldtau, cldemi,
1884     .            cldh, cldl, cldm, cldt, cldq)
1885c
1886c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1887c
1888      IF (MOD(itaprad,radpas).EQ.0) THEN
1889c      CALL orbite(FLOAT(julien),zlongi,dist)
1890c      IF (cycle_diurne) THEN
1891c        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1892c        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1893c        CALL zenith(zlongi,gmtime,rlat,rlon,rmu0,fract) !va disparaitre
1894c        CALL alboc_cd(rmu0,alb_eau)
1895c      ELSE
1896c        CALL angle(zlongi,rlat,fract,rmu0)
1897c        CALL alboc(FLOAT(julien),rlat,alb_eau)
1898c      ENDIF
1899c      CALL albsno(veget,agesno,alb_neig)
1900      DO i = 1, klon
1901c         falbe(i,is_oce) = alb_eau(i)
1902c         IF (pctsrf(i,is_oce).GT.epsfra .AND. ftsol(i,is_oce).LT.271.35)
1903c     .   falbe(i,is_oce) = 0.6 ! pour slab_ocean
1904c         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_lic)/(fsnow(i,is_lic)+10.0)))
1905c         falbe(i,is_lic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
1906c         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_ter)/(fsnow(i,is_ter)+10.0)))
1907c         falbe(i,is_ter) = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
1908c         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_sic)/(fsnow(i,is_sic)+10.0)))
1909         falbe(i,is_sic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
1910         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
1911     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
1912     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
1913     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
1914      ENDDO
1915c      DO nsrf = 1, nbsrf
1916c        DO i = 1, klon
1917c           albsol(i) = albsol(i) + falbe(i,nsrf)*pctsrf(i,nsrf)
1918c        ENDDO
1919c      ENDDO
1920      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
1921     e            (dist, rmu0, fract, co2_ppm, solaire,
1922     e             paprs, pplay,zxtsol,albsol, t_seri,q_seri,wo,
1923     e             cldfra, cldemi, cldtau,
1924     s             heat,heat0,cool,cool0,radsol,albpla,
1925     s             topsw,toplw,solsw,sollw,
1926     s             topsw0,toplw0,solsw0,sollw0)
1927      itaprad = 0
1928      ENDIF
1929      itaprad = itaprad + 1
1930c
1931c Ajouter la tendance des rayonnements (tous les pas)
1932c
1933      DO k = 1, klev
1934      DO i = 1, klon
1935         t_seri(i,k) = t_seri(i,k)
1936     .               + (heat(i,k)-cool(i,k)) * dtime/86400.
1937      ENDDO
1938      ENDDO
1939c
1940c Calculer l'hydrologie de la surface
1941c
1942c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
1943c     .            agesno, ftsol,fqsol,fsnow, ruis)
1944c
1945      DO i = 1, klon
1946         zxqsol(i) = 0.0
1947         zxsnow(i) = 0.0
1948      ENDDO
1949      DO nsrf = 1, nbsrf
1950      DO i = 1, klon
1951         zxqsol(i) = zxqsol(i) + fqsol(i,nsrf)*pctsrf(i,nsrf)
1952         zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
1953      ENDDO
1954      ENDDO
1955c
1956c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
1957c
1958      DO nsrf = 1, nbsrf
1959      DO i = 1, klon
1960         IF (pctsrf(i,nsrf).LT.epsfra) THEN
1961            fqsol(i,nsrf) = zxqsol(i)
1962            fsnow(i,nsrf) = zxsnow(i)
1963         ENDIF
1964      ENDDO
1965      ENDDO
1966c
1967c Calculer le bilan du sol et la derive de temperature (couplage)
1968c
1969      DO i = 1, klon
1970         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
1971      ENDDO
1972      IF (ok_ocean) THEN
1973         DO i = 1, klon
1974            cthermiq = cyang
1975            IF (ftsol(i,is_oce).LT. 271.35) cthermiq = cbing
1976            IF (pctsrf(i,is_oce).GT.epsfra) deltat(i) = deltat(i) +
1977     .                          (bils(i)-lmt_bils(i))/cthermiq * dtime
1978            IF (deltat(i).GT.15.0 ) deltat(i) = 15.0
1979            IF (deltat(i).LT.-15.0) deltat(i) = -15.0
1980         ENDDO
1981      ENDIF
1982c
1983cmoddeblott(jan95)
1984c Appeler le programme de parametrisation de l'orographie
1985c a l'echelle sous-maille:
1986c
1987      IF (ok_orodr) THEN
1988c
1989c  selection des points pour lesquels le shema est actif:
1990        igwd=0
1991        DO i=1,klon
1992        itest(i)=0
1993c        IF ((zstd(i).gt.10.0)) THEN
1994        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1995          itest(i)=1
1996          igwd=igwd+1
1997          idx(igwd)=i
1998        ENDIF
1999        ENDDO
2000        igwdim=MAX(1,igwd)
2001c
2002        CALL drag_noro(klon,klev,dtime,paprs,pplay,
2003     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
2004     e                   igwd,igwdim,idx,itest,
2005     e                   t_seri, u_seri, v_seri,
2006     s                   zulow, zvlow, zustr, zvstr,
2007     s                   d_t_oro, d_u_oro, d_v_oro)
2008c
2009c  ajout des tendances
2010        DO k = 1, klev
2011        DO i = 1, klon
2012           t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2013           u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2014           v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2015        ENDDO
2016        ENDDO
2017c
2018      ENDIF ! fin de test sur ok_orodr
2019c
2020      IF (ok_orolf) THEN
2021c
2022c  selection des points pour lesquels le shema est actif:
2023        igwd=0
2024        DO i=1,klon
2025        itest(i)=0
2026        IF ((zpic(i)-zmea(i)).GT.100.) THEN
2027          itest(i)=1
2028          igwd=igwd+1
2029          idx(igwd)=i
2030        ENDIF
2031        ENDDO
2032        igwdim=MAX(1,igwd)
2033c
2034        CALL lift_noro(klon,klev,dtime,paprs,pplay,
2035     e                   rlat,zmea,zstd, zsig, zgam, zthe,zpic,zval,
2036     e                   igwd,igwdim,idx,itest,
2037     e                   t_seri, u_seri, v_seri,
2038     s                   zulow, zvlow, zustr, zvstr,
2039     s                   d_t_lif, d_u_lif, d_v_lif)
2040c
2041c  ajout des tendances
2042        DO k = 1, klev
2043        DO i = 1, klon
2044           t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2045           u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2046           v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2047        ENDDO
2048        ENDDO
2049c
2050      ENDIF ! fin de test sur ok_orolf
2051c
2052cAA
2053cAA Installation de l'interface online-offline pour traceurs
2054cAA
2055c====================================================================
2056c   Calcul  des tendances traceurs
2057c====================================================================
2058CMAF modif pour garder info du nombre de traceurs auxquels
2059C la physique s'applique
2060C
2061      call phytrac (rnpb,
2062     I                   debut,
2063     I                   nqmax-2,
2064     I                   nlon,nlev,dtime,
2065     I                   t,paprs,pplay,
2066     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2067     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2068     I                   frac_impa, frac_nucl,
2069     I                   rlon,presnivs,paire,pphis,
2070     O                   tr_seri)
2071
2072      IF (offline) THEN
2073
2074         call phystokenc (
2075     I                   nlon,nlev,pdtphys,rlon,rlat,
2076     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2077     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
2078     I                   frac_impa, frac_nucl,
2079     I                   pphis,paire,dtime,itap,
2080     O                   physid)
2081
2082      ENDIF
2083
2084c
2085c Calculer le transport de l'eau et de l'energie (diagnostique)
2086c
2087      CALL transp (paprs,zxtsol,
2088     e                   t_seri, q_seri, u_seri, v_seri, zphi,
2089     s                   ve, vq, ue, uq)
2090c
2091c Accumuler les variables a stocker dans les fichiers histoire:
2092c
2093c
2094c
2095
2096      IF (ok_journe) THEN
2097c
2098      ndex2d = 0
2099      ndex3d = 0
2100c
2101c Champs 2D:
2102c
2103         i = NINT(zout/zsto)
2104         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2105         CALL histwrite(nid_day,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2106c
2107         i = NINT(zout/zsto)
2108         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2109         CALL histwrite(nid_day,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2110C
2111      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2112      CALL histwrite(nid_day,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2113c
2114C
2115      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_ter)
2116      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d ,zx_tmp_2d)
2117      CALL histwrite(nid_day,"tter",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2118C
2119      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_lic)
2120      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2121      CALL histwrite(nid_day,"tlic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2122C
2123      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_oce)
2124      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2125      CALL histwrite(nid_day,"toce",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2126C
2127      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_sic)
2128      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2129      CALL histwrite(nid_day,"tsic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2130C
2131      DO i = 1, klon
2132         zx_tmp_fi2d(i) = paprs(i,1)
2133      ENDDO
2134      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2135      CALL histwrite(nid_day,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2136c
2137      DO i = 1, klon
2138         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
2139      ENDDO
2140      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2141      CALL histwrite(nid_day,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2142c
2143      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
2144      CALL histwrite(nid_day,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2145c
2146      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2147      CALL histwrite(nid_day,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2148c
2149      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
2150      CALL histwrite(nid_day,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2151c
2152      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2153      CALL histwrite(nid_day,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2154c
2155      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2156      CALL histwrite(nid_day,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2157c
2158      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2159      CALL histwrite(nid_day,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2160c
2161      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2162      CALL histwrite(nid_day,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2163c
2164      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2165      CALL histwrite(nid_day,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2166c
2167      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2168      CALL histwrite(nid_day,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2169c
2170      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
2171      CALL histwrite(nid_day,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2172c
2173c      DO i = 1, klon
2174c         zx_tmp_fi2d(i) = fluxu(i,1)
2175c      ENDDO
2176c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2177c      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2178c
2179c      DO i = 1, klon
2180c         zx_tmp_fi2d(i) = fluxv(i,1)
2181c      ENDDO
2182c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2183c      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2184c
2185      DO nsrf = 1, nbsrf
2186C§§§
2187        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2188        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2189        CALL histwrite(nid_day,"pourc_"//clnsurf(nsrf),itap,
2190     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2191C
2192        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
2193        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2194        CALL histwrite(nid_day,"sens_"//clnsurf(nsrf),itap,
2195     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2196C
2197        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
2198        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2199        CALL histwrite(nid_day,"lat_"//clnsurf(nsrf),itap,
2200     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2201C
2202        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
2203        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2204        CALL histwrite(nid_day,"taux_"//clnsurf(nsrf),itap,
2205     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2206C     
2207        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
2208        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2209        CALL histwrite(nid_day,"tauy_"//clnsurf(nsrf),itap,
2210     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2211C
2212      END DO 
2213C
2214c$$$      DO i = 1, klon
2215c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
2216c$$$      ENDDO
2217c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2218c$$$      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2219c
2220      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
2221      CALL histwrite(nid_day,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2222c
2223      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
2224      CALL histwrite(nid_day,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2225c
2226      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
2227      CALL histwrite(nid_day,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2228c
2229      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
2230      CALL histwrite(nid_day,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2231c
2232      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
2233      CALL histwrite(nid_day,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2234c
2235c Champs 3D:
2236c
2237      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2238      CALL histwrite(nid_day,"temp",itap,zx_tmp_3d,
2239     .                                   iim*jjmp1*klev,ndex3d)
2240c
2241      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
2242      CALL histwrite(nid_day,"ovap",itap,zx_tmp_3d,
2243     .                                   iim*jjmp1*klev,ndex3d)
2244c
2245      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2246      CALL histwrite(nid_day,"geop",itap,zx_tmp_3d,
2247     .                                   iim*jjmp1*klev,ndex3d)
2248c
2249      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2250      CALL histwrite(nid_day,"vitu",itap,zx_tmp_3d,
2251     .                                   iim*jjmp1*klev,ndex3d)
2252c
2253      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2254      CALL histwrite(nid_day,"vitv",itap,zx_tmp_3d,
2255     .                                   iim*jjmp1*klev,ndex3d)
2256c
2257      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
2258      CALL histwrite(nid_day,"vitw",itap,zx_tmp_3d,
2259     .                                   iim*jjmp1*klev,ndex3d)
2260c
2261      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2262      CALL histwrite(nid_day,"pres",itap,zx_tmp_3d,
2263     .                                   iim*jjmp1*klev,ndex3d)
2264c
2265      if (ok_sync) then
2266        call histsync(nid_day)
2267      endif
2268      ENDIF
2269C
2270      IF (ok_mensuel) THEN
2271c
2272      ndex2d = 0
2273      ndex3d = 0
2274c
2275c Champs 2D:
2276c
2277         i = NINT(zout/zsto)
2278         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2279         CALL histwrite(nid_mth,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2280C
2281         i = NINT(zout/zsto)
2282         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2283         CALL histwrite(nid_mth,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2284
2285      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2286      CALL histwrite(nid_mth,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2287c
2288      DO i = 1, klon
2289         zx_tmp_fi2d(i) = paprs(i,1)
2290      ENDDO
2291      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2292      CALL histwrite(nid_mth,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2293c
2294      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxqsol,zx_tmp_2d)
2295      CALL histwrite(nid_mth,"qsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2296c
2297      DO i = 1, klon
2298         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
2299      ENDDO
2300      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2301      CALL histwrite(nid_mth,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2302c
2303      DO i = 1, klon
2304         zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
2305      ENDDO
2306      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2307      CALL histwrite(nid_mth,"plul",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2308c
2309      DO i = 1, klon
2310         zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
2311      ENDDO
2312      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2313      CALL histwrite(nid_mth,"pluc",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2314c
2315      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
2316      CALL histwrite(nid_mth,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2317c
2318      CALL gr_fi_ecrit(1, klon,iim,jjmp1, agesno,zx_tmp_2d)
2319      CALL histwrite(nid_mth,"ages",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2320c
2321      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2322      CALL histwrite(nid_mth,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2323c
2324      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
2325      CALL histwrite(nid_mth,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2326c
2327      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2328      CALL histwrite(nid_mth,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2329c
2330      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2331      CALL histwrite(nid_mth,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2332c
2333      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2334      CALL histwrite(nid_mth,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2335c
2336      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
2337      CALL histwrite(nid_mth,"tops0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2338c
2339      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw0,zx_tmp_2d)
2340      CALL histwrite(nid_mth,"topl0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2341c
2342      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw0,zx_tmp_2d)
2343      CALL histwrite(nid_mth,"sols0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2344c
2345      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw0,zx_tmp_2d)
2346      CALL histwrite(nid_mth,"soll0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2347c
2348      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2349      CALL histwrite(nid_mth,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2350c
2351      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2352      CALL histwrite(nid_mth,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2353c
2354      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2355      CALL histwrite(nid_mth,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2356c
2357      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
2358      CALL histwrite(nid_mth,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2359c
2360c      DO i = 1, klon
2361c         zx_tmp_fi2d(i) = fluxu(i,1)
2362c      ENDDO
2363c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2364c      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2365c
2366c      DO i = 1, klon
2367c         zx_tmp_fi2d(i) = fluxv(i,1)
2368c      ENDDO
2369c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2370c      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2371c
2372      DO nsrf = 1, nbsrf
2373C§§§
2374        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2375        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2376        CALL histwrite(nid_mth,"pourc_"//clnsurf(nsrf),itap,
2377     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2378C
2379        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
2380        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2381        CALL histwrite(nid_mth,"sens_"//clnsurf(nsrf),itap,
2382     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2383C
2384        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
2385        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2386        CALL histwrite(nid_mth,"lat_"//clnsurf(nsrf),itap,
2387     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2388C
2389        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
2390        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2391        CALL histwrite(nid_mth,"taux_"//clnsurf(nsrf),itap,
2392     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2393C     
2394        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
2395        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2396        CALL histwrite(nid_mth,"tauy_"//clnsurf(nsrf),itap,
2397     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2398C
2399      END DO 
2400c$$$      DO i = 1, klon
2401c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
2402c$$$      ENDDO
2403c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2404c$$$      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2405c
2406      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
2407      CALL histwrite(nid_mth,"albs",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2408c
2409      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
2410      CALL histwrite(nid_mth,"cdrm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2411c
2412      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
2413      CALL histwrite(nid_mth,"cdrh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2414c
2415      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
2416      CALL histwrite(nid_mth,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2417c
2418      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
2419      CALL histwrite(nid_mth,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2420c
2421      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
2422      CALL histwrite(nid_mth,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2423c
2424      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
2425      CALL histwrite(nid_mth,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2426c
2427      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
2428      CALL histwrite(nid_mth,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2429c
2430      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ue,zx_tmp_2d)
2431      CALL histwrite(nid_mth,"ue",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2432c
2433      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ve,zx_tmp_2d)
2434      CALL histwrite(nid_mth,"ve",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2435c
2436      CALL gr_fi_ecrit(1, klon,iim,jjmp1, uq,zx_tmp_2d)
2437      CALL histwrite(nid_mth,"uq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2438c
2439      CALL gr_fi_ecrit(1, klon,iim,jjmp1, vq,zx_tmp_2d)
2440      CALL histwrite(nid_mth,"vq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2441c
2442c Champs 3D:
2443C
2444      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2445      CALL histwrite(nid_mth,"temp",itap,zx_tmp_3d,
2446     .                                   iim*jjmp1*klev,ndex3d)
2447c
2448      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
2449      CALL histwrite(nid_mth,"ovap",itap,zx_tmp_3d,
2450     .                                   iim*jjmp1*klev,ndex3d)
2451c
2452      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2453      CALL histwrite(nid_mth,"geop",itap,zx_tmp_3d,
2454     .                                   iim*jjmp1*klev,ndex3d)
2455c
2456      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2457      CALL histwrite(nid_mth,"vitu",itap,zx_tmp_3d,
2458     .                                   iim*jjmp1*klev,ndex3d)
2459c
2460      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2461      CALL histwrite(nid_mth,"vitv",itap,zx_tmp_3d,
2462     .                                   iim*jjmp1*klev,ndex3d)
2463c
2464      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
2465      CALL histwrite(nid_mth,"vitw",itap,zx_tmp_3d,
2466     .                                   iim*jjmp1*klev,ndex3d)
2467c
2468      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2469      CALL histwrite(nid_mth,"pres",itap,zx_tmp_3d,
2470     .                                   iim*jjmp1*klev,ndex3d)
2471c
2472      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldfra, zx_tmp_3d)
2473      CALL histwrite(nid_mth,"rneb",itap,zx_tmp_3d,
2474     .                                   iim*jjmp1*klev,ndex3d)
2475c
2476      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_rh, zx_tmp_3d)
2477      CALL histwrite(nid_mth,"rhum",itap,zx_tmp_3d,
2478     .                                   iim*jjmp1*klev,ndex3d)
2479c
2480      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldliq, zx_tmp_3d)
2481      CALL histwrite(nid_mth,"oliq",itap,zx_tmp_3d,
2482     .                                   iim*jjmp1*klev,ndex3d)
2483c
2484      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_dyn, zx_tmp_3d)
2485      CALL histwrite(nid_mth,"dtdyn",itap,zx_tmp_3d,
2486     .                                   iim*jjmp1*klev,ndex3d)
2487c
2488      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_dyn, zx_tmp_3d)
2489      CALL histwrite(nid_mth,"dqdyn",itap,zx_tmp_3d,
2490     .                                   iim*jjmp1*klev,ndex3d)
2491c
2492      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_con, zx_tmp_3d)
2493      CALL histwrite(nid_mth,"dtcon",itap,zx_tmp_3d,
2494     .                                   iim*jjmp1*klev,ndex3d)
2495c
2496      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_con, zx_tmp_3d)
2497      CALL histwrite(nid_mth,"dqcon",itap,zx_tmp_3d,
2498     .                                   iim*jjmp1*klev,ndex3d)
2499c
2500      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_lsc, zx_tmp_3d)
2501      CALL histwrite(nid_mth,"dtlsc",itap,zx_tmp_3d,
2502     .                                   iim*jjmp1*klev,ndex3d)
2503c
2504      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_lsc, zx_tmp_3d)
2505      CALL histwrite(nid_mth,"dqlsc",itap,zx_tmp_3d,
2506     .                                   iim*jjmp1*klev,ndex3d)
2507c
2508      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
2509      CALL histwrite(nid_mth,"dtvdf",itap,zx_tmp_3d,
2510     .                                   iim*jjmp1*klev,ndex3d)
2511c
2512      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
2513      CALL histwrite(nid_mth,"dqvdf",itap,zx_tmp_3d,
2514     .                                   iim*jjmp1*klev,ndex3d)
2515c
2516      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_eva, zx_tmp_3d)
2517      CALL histwrite(nid_mth,"dteva",itap,zx_tmp_3d,
2518     .                                   iim*jjmp1*klev,ndex3d)
2519c
2520      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_eva, zx_tmp_3d)
2521      CALL histwrite(nid_mth,"dqeva",itap,zx_tmp_3d,
2522     .                                   iim*jjmp1*klev,ndex3d)
2523c
2524      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
2525      CALL histwrite(nid_mth,"dtajs",itap,zx_tmp_3d,
2526     .                                   iim*jjmp1*klev,ndex3d)
2527c
2528      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_ajs, zx_tmp_3d)
2529      CALL histwrite(nid_mth,"dqajs",itap,zx_tmp_3d,
2530     .                                   iim*jjmp1*klev,ndex3d)
2531c
2532      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat, zx_tmp_3d)
2533      CALL histwrite(nid_mth,"dtswr",itap,zx_tmp_3d,
2534     .                                   iim*jjmp1*klev,ndex3d)
2535c
2536      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat0, zx_tmp_3d)
2537      CALL histwrite(nid_mth,"dtsw0",itap,zx_tmp_3d,
2538     .                                   iim*jjmp1*klev,ndex3d)
2539c
2540      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool, zx_tmp_3d)
2541      CALL histwrite(nid_mth,"dtlwr",itap,zx_tmp_3d,
2542     .                                   iim*jjmp1*klev,ndex3d)
2543c
2544      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool0, zx_tmp_3d)
2545      CALL histwrite(nid_mth,"dtlw0",itap,zx_tmp_3d,
2546     .                                   iim*jjmp1*klev,ndex3d)
2547c
2548      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_vdf, zx_tmp_3d)
2549      CALL histwrite(nid_mth,"duvdf",itap,zx_tmp_3d,
2550     .                                   iim*jjmp1*klev,ndex3d)
2551c
2552      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_vdf, zx_tmp_3d)
2553      CALL histwrite(nid_mth,"dvvdf",itap,zx_tmp_3d,
2554     .                                   iim*jjmp1*klev,ndex3d)
2555c
2556      IF (ok_orodr) THEN
2557      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_oro, zx_tmp_3d)
2558      CALL histwrite(nid_mth,"duoro",itap,zx_tmp_3d,
2559     .                                   iim*jjmp1*klev,ndex3d)
2560c
2561      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_oro, zx_tmp_3d)
2562      CALL histwrite(nid_mth,"dvoro",itap,zx_tmp_3d,
2563     .                                   iim*jjmp1*klev,ndex3d)
2564c
2565      ENDIF
2566C
2567      IF (ok_orolf) THEN
2568      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_lif, zx_tmp_3d)
2569      CALL histwrite(nid_mth,"dulif",itap,zx_tmp_3d,
2570     .                                   iim*jjmp1*klev,ndex3d)
2571c
2572      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_lif, zx_tmp_3d)
2573      CALL histwrite(nid_mth,"dvlif",itap,zx_tmp_3d,
2574     .                                   iim*jjmp1*klev,ndex3d)
2575      ENDIF
2576C
2577      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, wo, zx_tmp_3d)
2578      CALL histwrite(nid_mth,"ozone",itap,zx_tmp_3d,
2579     .                                   iim*jjmp1*klev,ndex3d)
2580c
2581      IF (nqmax.GE.3) THEN
2582      DO iq=1,nqmax-2
2583      IF (iq.LE.99) THEN
2584         CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,iq+2), zx_tmp_3d)
2585         WRITE(str2,'(i2.2)') iq
2586         CALL histwrite(nid_mth,"trac"//str2,itap,zx_tmp_3d,
2587     .                                   iim*jjmp1*klev,ndex3d)
2588      ELSE
2589         PRINT*, "Trop de traceurs"
2590         CALL abort
2591      ENDIF
2592      ENDDO
2593      ENDIF
2594c
2595      if (ok_sync) then
2596        call histsync(nid_mth)
2597      endif
2598      ENDIF
2599c
2600      IF (ok_instan) THEN
2601c
2602      ndex2d = 0
2603      ndex3d = 0
2604c
2605c Champs 2D:
2606c
2607         i = NINT(zout/zsto)
2608         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2609         CALL histwrite(nid_ins,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2610c
2611         i = NINT(zout/zsto)
2612         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2613         CALL histwrite(nid_ins,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2614
2615      DO i = 1, klon
2616         zx_tmp_fi2d(i) = paprs(i,1)
2617      ENDDO
2618      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2619      CALL histwrite(nid_ins,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2620c
2621      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2622      CALL histwrite(nid_ins,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2623c
2624      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2625      CALL histwrite(nid_ins,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2626c
2627      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2628      CALL histwrite(nid_ins,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2629c
2630      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2631      CALL histwrite(nid_ins,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2632c
2633      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2634      CALL histwrite(nid_ins,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2635c
2636      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2637      CALL histwrite(nid_ins,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2638c
2639      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2640      CALL histwrite(nid_ins,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2641c
2642      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2643      CALL histwrite(nid_ins,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2644c
2645      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_oce),zx_tmp_2d)
2646      CALL histwrite(nid_ins,"dtsvdfo",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2647c
2648      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_ter),zx_tmp_2d)
2649      CALL histwrite(nid_ins,"dtsvdft",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2650c
2651      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_lic),zx_tmp_2d)
2652      CALL histwrite(nid_ins,"dtsvdfg",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2653c
2654      CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_sic),zx_tmp_2d)
2655      CALL histwrite(nid_ins,"dtsvdfi",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2656
2657      DO nsrf = 1, nbsrf
2658C§§§
2659        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2660        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2661        CALL histwrite(nid_ins,"pourc_"//clnsurf(nsrf),itap,
2662     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2663C
2664        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
2665        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2666        CALL histwrite(nid_ins,"sens_"//clnsurf(nsrf),itap,
2667     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2668C
2669        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
2670        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2671        CALL histwrite(nid_ins,"lat_"//clnsurf(nsrf),itap,
2672     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2673C
2674        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
2675        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2676        CALL histwrite(nid_ins,"tsol_"//clnsurf(nsrf),itap,
2677     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2678C
2679        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
2680        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2681        CALL histwrite(nid_ins,"taux_"//clnsurf(nsrf),itap,
2682     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2683C     
2684        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
2685        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
2686        CALL histwrite(nid_ins,"tauy_"//clnsurf(nsrf),itap,
2687     $      zx_tmp_2d,iim*jjmp1,ndex2d)
2688C
2689      END DO 
2690
2691c
2692c Champs 3D:
2693c
2694      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2695      CALL histwrite(nid_ins,"temp",itap,zx_tmp_3d,
2696     .                                   iim*jjmp1*klev,ndex3d)
2697c
2698      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2699      CALL histwrite(nid_ins,"vitu",itap,zx_tmp_3d,
2700     .                                   iim*jjmp1*klev,ndex3d)
2701c
2702      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2703      CALL histwrite(nid_ins,"vitv",itap,zx_tmp_3d,
2704     .                                   iim*jjmp1*klev,ndex3d)
2705c
2706      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2707      CALL histwrite(nid_ins,"geop",itap,zx_tmp_3d,
2708     .                                   iim*jjmp1*klev,ndex3d)
2709c
2710      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2711      CALL histwrite(nid_ins,"pres",itap,zx_tmp_3d,
2712     .                                   iim*jjmp1*klev,ndex3d)
2713c
2714      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
2715      CALL histwrite(nid_ins,"dtvdf",itap,zx_tmp_3d,
2716     .                                   iim*jjmp1*klev,ndex3d)
2717c
2718      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
2719      CALL histwrite(nid_ins,"dqvdf",itap,zx_tmp_3d,
2720     .                                   iim*jjmp1*klev,ndex3d)
2721
2722c
2723      if (ok_sync) then
2724        call histsync(nid_ins)
2725      endif
2726      ENDIF
2727c
2728c
2729c Ecrire la bande regionale (binaire grads)
2730      IF (ok_region .AND. mod(itap,ecrit_reg).eq.0) THEN
2731         CALL ecriregs(84,zxtsol)
2732         CALL ecriregs(84,paprs(1,1))
2733         CALL ecriregs(84,topsw)
2734         CALL ecriregs(84,toplw)
2735         CALL ecriregs(84,solsw)
2736         CALL ecriregs(84,sollw)
2737         CALL ecriregs(84,rain_fall)
2738         CALL ecriregs(84,snow_fall)
2739         CALL ecriregs(84,evap)
2740         CALL ecriregs(84,sens)
2741         CALL ecriregs(84,bils)
2742         CALL ecriregs(84,pctsrf(1,is_sic))
2743         CALL ecriregs(84,zxfluxu(1,1))
2744         CALL ecriregs(84,zxfluxv(1,1))
2745         CALL ecriregs(84,ue)
2746         CALL ecriregs(84,ve)
2747         CALL ecriregs(84,uq)
2748         CALL ecriregs(84,vq)
2749c
2750         CALL ecrirega(84,u_seri)
2751         CALL ecrirega(84,v_seri)
2752         CALL ecrirega(84,omega)
2753         CALL ecrirega(84,t_seri)
2754         CALL ecrirega(84,zphi)
2755         CALL ecrirega(84,q_seri)
2756         CALL ecrirega(84,cldfra)
2757         CALL ecrirega(84,cldliq)
2758         CALL ecrirega(84,pplay)
2759
2760
2761cc         CALL ecrirega(84,d_t_dyn)
2762cc         CALL ecrirega(84,d_q_dyn)
2763cc         CALL ecrirega(84,heat)
2764cc         CALL ecrirega(84,cool)
2765cc         CALL ecrirega(84,d_t_con)
2766cc         CALL ecrirega(84,d_q_con)
2767cc         CALL ecrirega(84,d_t_lsc)
2768cc         CALL ecrirega(84,d_q_lsc)
2769      ENDIF
2770c
2771c Convertir les incrementations en tendances
2772c
2773      DO k = 1, klev
2774      DO i = 1, klon
2775         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2776         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2777         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2778         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2779         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2780      ENDDO
2781      ENDDO
2782c
2783      IF (nqmax.GE.3) THEN
2784      DO iq = 3, nqmax
2785      DO  k = 1, klev
2786      DO  i = 1, klon
2787         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2788      ENDDO
2789      ENDDO
2790      ENDDO
2791      ENDIF
2792c
2793c Sauvegarder les valeurs de t et q a la fin de la physique:
2794c
2795      DO k = 1, klev
2796      DO i = 1, klon
2797         t_ancien(i,k) = t_seri(i,k)
2798         q_ancien(i,k) = q_seri(i,k)
2799      ENDDO
2800      ENDDO
2801c
2802c====================================================================
2803c Si c'est la fin, il faut conserver l'etat de redemarrage
2804c====================================================================
2805c
2806      IF (lafin) THEN
2807ccc         IF (ok_oasis) CALL quitcpl
2808         CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
2809     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
2810     .      falbe, fevap, rain_fall, snow_fall,
2811     .      solsw, sollw,
2812     .      radsol,frugs,agesno,
2813     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2814     .      t_ancien, q_ancien)
2815      ENDIF
2816
2817      RETURN
2818      END
2819      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2820      IMPLICIT none
2821c
2822c Calculer et imprimer l'eau totale. A utiliser pour verifier
2823c la conservation de l'eau
2824c
2825#include "YOMCST.h"
2826      INTEGER klon,klev
2827      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2828      REAL aire(klon)
2829      REAL qtotal, zx, qcheck
2830      INTEGER i, k
2831c
2832      zx = 0.0
2833      DO i = 1, klon
2834         zx = zx + aire(i)
2835      ENDDO
2836      qtotal = 0.0
2837      DO k = 1, klev
2838      DO i = 1, klon
2839         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2840     .                     *(paprs(i,k)-paprs(i,k+1))/RG
2841      ENDDO
2842      ENDDO
2843c
2844      qcheck = qtotal/zx
2845c
2846      RETURN
2847      END
2848      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2849      IMPLICIT none
2850c
2851c Tranformer une variable de la grille physique a
2852c la grille d'ecriture
2853c
2854      INTEGER nfield,nlon,iim,jjmp1, jjm
2855      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2856c
2857      INTEGER i, j, n, ig
2858c
2859      jjm = jjmp1 - 1
2860      DO n = 1, nfield
2861         DO i=1,iim
2862            ecrit(i,n) = fi(1,n)
2863            ecrit(i+jjm*iim,n) = fi(nlon,n)
2864         ENDDO
2865         DO ig = 1, nlon - 2
2866           ecrit(iim+ig,n) = fi(1+ig,n)
2867         ENDDO
2868      ENDDO
2869      RETURN
2870      END
2871
Note: See TracBrowser for help on using the repository browser.