source: LMDZ.3.3/trunk/libf/phylmd/physiq.F @ 250

Last change on this file since 250 was 250, checked in by lmdz, 23 years ago

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