source: LMDZ.3.3/branches/LF/libf/phylmd/physiq.F @ 467

Last change on this file since 467 was 2, checked in by lmdz, 25 years ago

Initial revision

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