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

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

Un peu de menage
Elimination du runoff
LF

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