source: LMDZ.3.3/branches/rel-1-0-patch/libf/phylmd/physiq.F @ 349

Last change on this file since 349 was 349, checked in by lmdz, 22 years ago

Regle le probleme de decalage de 1 jour au debut de chaque simulation
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      real 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      integer idayref
594
595      character*20 modname
596      character*80 abort_message
597      logical ok_sync
598
599c
600c Declaration des constantes et des fonctions thermodynamiques
601c
602#include "YOMCST.h"
603#include "YOETHF.h"
604#include "FCTTRE.h"
605c======================================================================
606      modname = 'physiq'
607      ok_sync=.TRUE.
608      IF (nqmax .LT. 2) THEN
609         PRINT*, 'eaux vapeur et liquide sont indispensables'
610         CALL ABORT
611      ENDIF
612      IF (debut) THEN
613         CALL suphec ! initialiser constantes et parametres phys.
614      ENDIF
615c======================================================================
616      xjour = rjourvrai
617c
618c Si c'est le debut, il faut initialiser plusieurs choses
619c          ********
620c
621       IF (debut) THEN
622c
623 
624         IF (ok_oasis) THEN
625            PRINT*, "Attentions! les parametres suivants sont fixes:"
626            PRINT *,'***********************************************'
627            PRINT*, "npas, nexca, itimestep=", npas, nexca, itimestep
628            PRINT*, "Changer-les manuellement s il le faut"
629            PRINT *,'***********************************************'
630            CALL inicma( npas, nexca, itimestep)
631         ENDIF
632c
633         IF (ok_ocean) THEN
634            PRINT*, '************************'
635            PRINT*, 'SLAB OCEAN est active, prenez precautions !'
636            PRINT*, '************************'
637         ENDIF
638c
639         DO k = 2, nvm          ! pas de vegetation
640            DO i = 1, klon
641               veget(i,k) = 0.0
642            ENDDO
643         ENDDO
644         DO i = 1, klon
645            veget(i,1) = 1.0    ! il n'y a que du desert
646         ENDDO
647         PRINT*, 'Pas de vegetation; desert partout'
648c
649c Initialiser les compteurs:
650c
651
652         itap    = 0
653         itaprad = 0
654c
655         CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
656     .         rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
657     .        radsol,rugmer,agesno,clesphy0,
658     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
659     .       t_ancien, q_ancien, ancien_ok )
660
661c
662         radpas = NINT( 86400./dtime/nbapp_rad)
663
664c
665         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
666     ,                    ok_instan, ok_region )
667c
668         IF (ABS(dtime-pdtphys).GT.0.001) THEN
669            PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
670            abort_message=' See above '
671            call abort_gcm(modname,abort_message,1)
672         ENDIF
673         IF (nlon .NE. klon) THEN
674            PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
675            abort_message=' See above '
676            call abort_gcm(modname,abort_message,1)
677         ENDIF
678         IF (nlev .NE. klev) THEN
679            PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
680            abort_message=' See above '
681            call abort_gcm(modname,abort_message,1)
682         ENDIF
683c
684         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
685           PRINT*, 'Nbre d appels au rayonnement insuffisant'
686           PRINT*, "Au minimum 4 appels par jour si cycle diurne"
687           abort_message=' See above '
688           call abort_gcm(modname,abort_message,1)
689         ENDIF
690         PRINT*, "Clef pour la convection, iflag_con=", iflag_con
691c
692cKE43
693c Initialisation pour la convection de K.E. (sb):
694         IF (iflag_con.EQ.4) THEN
695
696         PRINT*, "*** Convection de Kerry Emanuel 4.3  "
697         PRINT*, "On va utiliser le melange convectif des traceurs qui"
698         PRINT*, "est calcule dans convect4.3"
699         PRINT*, " !!! penser aux logical flags de phytrac"
700
701          DO i = 1, klon
702           ema_cbmf(i) = 0.
703           ema_pcb(i)  = 0.
704           ema_pct(i)  = 0.
705           ema_workcbmf(i) = 0.
706          ENDDO
707         ENDIF
708c34EK
709         IF (ok_orodr) THEN
710         DO i=1,klon
711         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
712         ENDDO
713         CALL SUGWD(klon,klev,paprs,pplay)
714         DO i=1,klon
715         zuthe(i)=0.
716         zvthe(i)=0.
717         if(zstd(i).gt.10.)then
718           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
719           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
720         endif
721         ENDDO
722         ENDIF
723c
724         IF (soil_model) THEN
725            DO nsrf = 1, nbsrf
726            CALL soil(dtime, nsrf, fsnow(1,nsrf),
727     .                ftsol(1,nsrf), ftsoil(1,1,nsrf),
728     .                soilcap(1,nsrf), soilflux(1,nsrf))
729            ENDDO
730         ENDIF
731c
732         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
733         PRINT*,'La frequence de lecture surface est de ', lmt_pas
734c
735         ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
736         IF (ok_mensuel) THEN
737         PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
738         ENDIF
739         ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
740         IF (ok_journe) THEN
741         PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
742         ENDIF
743ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
744         ecrit_ins = NINT(86400./dtime *0.25)  ! tous les jours
745         IF (ok_instan) THEN
746         PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
747         ENDIF
748         ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
749         IF (ok_region) THEN
750         PRINT*, 'La frequence de sortie region est de ', ecrit_reg
751         ENDIF
752c
753c
754      IF (ok_journe) THEN
755c
756         idayref = day_ini
757         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
758c
759         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
760         DO i = 1, iim
761            zx_lon(i,1) = rlon(i+1)
762            zx_lon(i,jjmp1) = rlon(i+1)
763         ENDDO
764         DO ll=1,klev
765            znivsig(ll)=float(ll)
766         ENDDO
767         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
768         CALL histbeg("histday", iim,zx_lon, jjmp1,zx_lat,
769     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
770     .                 nhori, nid_day)
771         CALL histvert(nid_day, "presnivs", "Vertical levels", "mb",
772     .                 klev, presnivs, nvert)
773c        call histvert(nid_day, 'sig_s', 'Niveaux sigma','-',
774c    .              klev, znivsig, nvert)
775c
776         zsto = dtime
777         zout = dtime * FLOAT(ecrit_day)
778c
779         CALL histdef(nid_day, "phis", "Surface geop. height", "-",
780     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
781     .                "once", zsto,zout)
782c
783         CALL histdef(nid_day, "aire", "Grid area", "-",
784     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
785     .                "once", zsto,zout)
786c
787c Champs 2D:
788c
789         CALL histdef(nid_day, "tsol", "Surface Temperature", "K",
790     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
791     .                "ave(X)", zsto,zout)
792c
793         CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
794     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
795     .                "ave(X)", zsto,zout)
796c
797         CALL histdef(nid_day, "rain", "Precipitation", "mm/day",
798     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
799     .                "ave(X)", zsto,zout)
800c
801         CALL histdef(nid_day, "snow", "Snow fall", "mm/day",
802     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
803     .                "ave(X)", zsto,zout)
804c
805         CALL histdef(nid_day, "evap", "Evaporation", "mm/day",
806     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
807     .                "ave(X)", zsto,zout)
808c
809         CALL histdef(nid_day, "tops", "Solar rad. at TOA", "W/m2",
810     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
811     .                "ave(X)", zsto,zout)
812c
813         CALL histdef(nid_day, "topl", "IR rad. at TOA", "W/m2",
814     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
815     .                "ave(X)", zsto,zout)
816c
817         CALL histdef(nid_day, "sols", "Solar rad. at surf.", "W/m2",
818     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
819     .                "ave(X)", zsto,zout)
820c
821         CALL histdef(nid_day, "soll", "IR rad. at surface", "W/m2",
822     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
823     .                "ave(X)", zsto,zout)
824c
825         CALL histdef(nid_day, "bils", "Surf. total heat flux", "W/m2",
826     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
827     .                "ave(X)", zsto,zout)
828c
829         CALL histdef(nid_day, "sens", "Sensible heat flux", "W/m2",
830     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
831     .                "ave(X)", zsto,zout)
832c
833         CALL histdef(nid_day, "fder", "Heat flux derivation", "W/m2",
834     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
835     .                "ave(X)", zsto,zout)
836c
837         CALL histdef(nid_day, "frtu", "Zonal wind stress", "Pa",
838     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
839     .                "ave(X)", zsto,zout)
840c
841         CALL histdef(nid_day, "frtv", "Meridional wind stress", "Pa",
842     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
843     .                "ave(X)", zsto,zout)
844c
845         CALL histdef(nid_day, "ruis", "Runoff", "mm/day",
846     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
847     .                "ave(X)", zsto,zout)
848c
849         CALL histdef(nid_day, "sicf", "Sea-ice fraction", "-",
850     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
851     .                "ave(X)", zsto,zout)
852c
853         CALL histdef(nid_day, "cldl", "Low-level cloudiness", "-",
854     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
855     .                "ave(X)", zsto,zout)
856c
857         CALL histdef(nid_day, "cldm", "Mid-level cloudiness", "-",
858     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
859     .                "ave(X)", zsto,zout)
860c
861         CALL histdef(nid_day, "cldh", "High-level cloudiness", "-",
862     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
863     .                "ave(X)", zsto,zout)
864c
865         CALL histdef(nid_day, "cldt", "Total cloudiness", "-",
866     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
867     .                "ave(X)", zsto,zout)
868c
869         CALL histdef(nid_day, "cldq", "Cloud liquid water path", "-",
870     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
871     .                "ave(X)", zsto,zout)
872c
873c Champs 3D:
874c
875         CALL histdef(nid_day, "temp", "Air temperature", "K",
876     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
877     .                "ave(X)", zsto,zout)
878c
879         CALL histdef(nid_day, "ovap", "Specific humidity", "Kg/Kg",
880     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
881     .                "ave(X)", zsto,zout)
882c
883         CALL histdef(nid_day, "geop", "Geopotential height", "m",
884     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
885     .                "ave(X)", zsto,zout)
886c
887         CALL histdef(nid_day, "vitu", "Zonal wind", "m/s",
888     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
889     .                "ave(X)", zsto,zout)
890c
891         CALL histdef(nid_day, "vitv", "Meridional wind", "m/s",
892     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
893     .                "ave(X)", zsto,zout)
894c
895         CALL histdef(nid_day, "vitw", "Vertical wind", "m/s",
896     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
897     .                "ave(X)", zsto,zout)
898c
899         CALL histdef(nid_day, "pres", "Air pressure", "Pa",
900     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
901     .                "ave(X)", zsto,zout)
902c
903         CALL histend(nid_day)
904c
905         ndex2d = 0
906         ndex3d = 0
907c
908      ENDIF ! fin de test sur ok_journe
909c
910      IF (ok_mensuel) THEN
911c
912         idayref = day_ini
913         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
914c
915         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
916         DO i = 1, iim
917            zx_lon(i,1) = rlon(i+1)
918            zx_lon(i,jjmp1) = rlon(i+1)
919         ENDDO
920         DO ll=1,klev
921            znivsig(ll)=float(ll)
922         ENDDO
923         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
924         CALL histbeg("histmth", iim,zx_lon, jjmp1,zx_lat,
925     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
926     .                 nhori, nid_mth)
927         CALL histvert(nid_mth, "presnivs", "Vertical levels", "mb",
928     .                 klev, presnivs, nvert)
929c        call histvert(nid_mth, 'sig_s', 'Niveaux sigma','-',
930c    .              klev, znivsig, nvert)
931c
932         zsto = dtime
933         zout = dtime * FLOAT(ecrit_mth)
934c
935         CALL histdef(nid_mth, "phis", "Surface geop. height", "-",
936     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
937     .                "once",  zsto,zout)
938c
939         CALL histdef(nid_mth, "aire", "Grid area", "-",
940     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
941     .                "once",  zsto,zout)
942c
943c Champs 2D:
944c
945         CALL histdef(nid_mth, "tsol", "Surface Temperature", "K",
946     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
947     .                "ave(X)", zsto,zout)
948c
949         CALL histdef(nid_mth, "psol", "Surface Pressure", "Pa",
950     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
951     .                "ave(X)", zsto,zout)
952c
953         CALL histdef(nid_mth, "qsol", "Surface humidity", "mm",
954     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
955     .                "ave(X)", zsto,zout)
956c
957         CALL histdef(nid_mth, "rain", "Precipitation", "mm/day",
958     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
959     .                "ave(X)", zsto,zout)
960c
961         CALL histdef(nid_mth, "plul", "Large-scale Precip.", "mm/day",
962     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
963     .                "ave(X)", zsto,zout)
964c
965         CALL histdef(nid_mth, "pluc", "Convective Precip.", "mm/day",
966     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
967     .                "ave(X)", zsto,zout)
968c
969         CALL histdef(nid_mth, "snow", "Snow fall", "mm/day",
970     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
971     .                "ave(X)", zsto,zout)
972c
973         CALL histdef(nid_mth, "ages", "Snow age", "day",
974     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
975     .                "ave(X)", zsto,zout)
976c
977         CALL histdef(nid_mth, "evap", "Evaporation", "mm/day",
978     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
979     .                "ave(X)", zsto,zout)
980c
981         CALL histdef(nid_mth, "tops", "Solar rad. at TOA", "W/m2",
982     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
983     .                "ave(X)", zsto,zout)
984c
985         CALL histdef(nid_mth, "topl", "IR rad. at TOA", "W/m2",
986     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
987     .                "ave(X)", zsto,zout)
988c
989         CALL histdef(nid_mth, "sols", "Solar rad. at surf.", "W/m2",
990     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
991     .                "ave(X)", zsto,zout)
992c
993         CALL histdef(nid_mth, "soll", "IR rad. at surface", "W/m2",
994     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
995     .                "ave(X)", zsto,zout)
996c
997         CALL histdef(nid_mth, "tops0", "Solar rad. at TOA", "W/m2",
998     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
999     .                "ave(X)", zsto,zout)
1000c
1001         CALL histdef(nid_mth, "topl0", "IR rad. at TOA", "W/m2",
1002     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1003     .                "ave(X)", zsto,zout)
1004c
1005         CALL histdef(nid_mth, "sols0", "Solar rad. at surf.", "W/m2",
1006     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1007     .                "ave(X)", zsto,zout)
1008c
1009         CALL histdef(nid_mth, "soll0", "IR rad. at surface", "W/m2",
1010     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1011     .                "ave(X)", zsto,zout)
1012c
1013         CALL histdef(nid_mth, "bils", "Surf. total heat flux", "W/m2",
1014     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1015     .                "ave(X)", zsto,zout)
1016c
1017         CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
1018     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1019     .                "ave(X)", zsto,zout)
1020c
1021         CALL histdef(nid_mth, "fder", "Heat flux derivation", "W/m2",
1022     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1023     .                "ave(X)", zsto,zout)
1024c
1025         CALL histdef(nid_mth, "frtu", "Zonal wind stress", "Pa",
1026     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1027     .                "ave(X)", zsto,zout)
1028c
1029         CALL histdef(nid_mth, "frtv", "Meridional wind stress", "Pa",
1030     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1031     .                "ave(X)", zsto,zout)
1032c
1033         CALL histdef(nid_mth, "ruis", "Runoff", "mm/day",
1034     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1035     .                "ave(X)", zsto,zout)
1036c
1037         CALL histdef(nid_mth, "sicf", "Sea-ice fraction", "-",
1038     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1039     .                "ave(X)", zsto,zout)
1040c
1041         CALL histdef(nid_mth, "albs", "Surface albedo", "-",
1042     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1043     .                "ave(X)", zsto,zout)
1044c
1045         CALL histdef(nid_mth, "cdrm", "Momentum drag coef.", "-",
1046     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1047     .                "ave(X)", zsto,zout)
1048c
1049         CALL histdef(nid_mth, "cdrh", "Heat drag coef.", "-",
1050     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1051     .                "ave(X)", zsto,zout)
1052c
1053         CALL histdef(nid_mth, "cldl", "Low-level cloudiness", "-",
1054     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1055     .                "ave(X)", zsto,zout)
1056c
1057         CALL histdef(nid_mth, "cldm", "Mid-level cloudiness", "-",
1058     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1059     .                "ave(X)", zsto,zout)
1060c
1061         CALL histdef(nid_mth, "cldh", "High-level cloudiness", "-",
1062     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1063     .                "ave(X)", zsto,zout)
1064c
1065         CALL histdef(nid_mth, "cldt", "Total cloudiness", "-",
1066     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1067     .                "ave(X)", zsto,zout)
1068c
1069         CALL histdef(nid_mth, "cldq", "Cloud liquid water path", "-",
1070     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1071     .                "ave(X)", zsto,zout)
1072c
1073         CALL histdef(nid_mth, "ue", "Zonal energy transport", "-",
1074     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1075     .                "ave(X)", zsto,zout)
1076c
1077         CALL histdef(nid_mth, "ve", "Merid energy transport", "-",
1078     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1079     .                "ave(X)", zsto,zout)
1080c
1081         CALL histdef(nid_mth, "uq", "Zonal humidity transport", "-",
1082     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1083     .                "ave(X)", zsto,zout)
1084c
1085         CALL histdef(nid_mth, "vq", "Merid humidity transport", "-",
1086     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1087     .                "ave(X)", zsto,zout)
1088cKE43
1089      IF (iflag_con .EQ. 4) THEN ! sb
1090c
1091         CALL histdef(nid_mth, "cape", "Conv avlbl pot ener", "J/Kg",
1092     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1093     .                "ave(X)", zsto,zout)
1094c
1095         CALL histdef(nid_mth, "pbase", "Cld base pressure", "hPa",
1096     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1097     .                "ave(X)", zsto,zout)
1098c
1099         CALL histdef(nid_mth, "ptop", "Cld top pressure", "hPa",
1100     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1101     .                "ave(X)", zsto,zout)
1102c
1103         CALL histdef(nid_mth, "fbase", "Cld base mass flux", "Kg/m2/s",
1104     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1105     .                "ave(X)", zsto,zout)
1106c
1107c
1108      ENDIF
1109c34EK
1110c
1111c Champs 3D:
1112c
1113         CALL histdef(nid_mth, "temp", "Air temperature", "K",
1114     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1115     .                "ave(X)", zsto,zout)
1116c
1117         CALL histdef(nid_mth, "ovap", "Specific humidity", "Kg/Kg",
1118     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1119     .                "ave(X)", zsto,zout)
1120c
1121         CALL histdef(nid_mth, "geop", "Geopotential height", "m",
1122     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1123     .                "ave(X)", zsto,zout)
1124c
1125         CALL histdef(nid_mth, "vitu", "Zonal wind", "m/s",
1126     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1127     .                "ave(X)", zsto,zout)
1128c
1129         CALL histdef(nid_mth, "vitv", "Meridional wind", "m/s",
1130     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1131     .                "ave(X)", zsto,zout)
1132c
1133         CALL histdef(nid_mth, "vitw", "Vertical wind", "m/s",
1134     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1135     .                "ave(X)", zsto,zout)
1136c
1137         CALL histdef(nid_mth, "pres", "Air pressure", "Pa",
1138     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1139     .                "ave(X)", zsto,zout)
1140c
1141         CALL histdef(nid_mth, "rneb", "Cloud fraction", "-",
1142     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1143     .                "ave(X)", zsto,zout)
1144c
1145         CALL histdef(nid_mth, "rhum", "Relative humidity", "-",
1146     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1147     .                "ave(X)", zsto,zout)
1148c
1149         CALL histdef(nid_mth, "oliq", "Liquid water content", "kg/kg",
1150     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1151     .                "ave(X)", zsto,zout)
1152c
1153         CALL histdef(nid_mth, "dtdyn", "Dynamics dT", "K/s",
1154     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1155     .                "ave(X)", zsto,zout)
1156c
1157         CALL histdef(nid_mth, "dqdyn", "Dynamics dQ", "Kg/Kg/s",
1158     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1159     .                "ave(X)", zsto,zout)
1160c
1161         CALL histdef(nid_mth, "dtcon", "Convection dT", "K/s",
1162     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1163     .                "ave(X)", zsto,zout)
1164c
1165         CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
1166     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1167     .                "ave(X)", zsto,zout)
1168c
1169         CALL histdef(nid_mth, "dtlsc", "Condensation dT", "K/s",
1170     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1171     .                "ave(X)", zsto,zout)
1172c
1173         CALL histdef(nid_mth, "dqlsc", "Condensation dQ", "Kg/Kg/s",
1174     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1175     .                "ave(X)", zsto,zout)
1176c
1177         CALL histdef(nid_mth, "dtvdf", "Boundary-layer dT", "K/s",
1178     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1179     .                "ave(X)", zsto,zout)
1180c
1181         CALL histdef(nid_mth, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
1182     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1183     .                "ave(X)", zsto,zout)
1184c
1185         CALL histdef(nid_mth, "dteva", "Reevaporation dT", "K/s",
1186     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1187     .                "ave(X)", zsto,zout)
1188c
1189         CALL histdef(nid_mth, "dqeva", "Reevaporation dQ", "Kg/Kg/s",
1190     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1191     .                "ave(X)", zsto,zout)
1192
1193         CALL histdef(nid_mth, "ptconv", "POINTS CONVECTIFS"," ",
1194     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1195     .                "ave(X)", zsto,zout)
1196
1197         CALL histdef(nid_mth, "ratqs", "RATQS"," ",
1198     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1199     .                "ave(X)", zsto,zout)
1200
1201c
1202         CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
1203     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1204     .                "ave(X)", zsto,zout)
1205
1206         CALL histdef(nid_mth, "dqajs", "Dry adjust. dQ", "Kg/Kg/s",
1207     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1208     .                "ave(X)", zsto,zout)
1209c
1210         CALL histdef(nid_mth, "dtswr", "SW radiation dT", "K/s",
1211     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1212     .                "ave(X)", zsto,zout)
1213c
1214         CALL histdef(nid_mth, "dtsw0", "SW radiation dT", "K/s",
1215     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1216     .                "ave(X)", zsto,zout)
1217c
1218         CALL histdef(nid_mth, "dtlwr", "LW radiation dT", "K/s",
1219     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1220     .                "ave(X)", zsto,zout)
1221c
1222         CALL histdef(nid_mth, "dtlw0", "LW radiation dT", "K/s",
1223     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1224     .                "ave(X)", zsto,zout)
1225c
1226         CALL histdef(nid_mth, "duvdf", "Boundary-layer dU", "m/s2",
1227     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1228     .                "ave(X)", zsto,zout)
1229c
1230         CALL histdef(nid_mth, "dvvdf", "Boundary-layer dV", "m/s2",
1231     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1232     .                "ave(X)", zsto,zout)
1233c
1234         IF (ok_orodr) THEN
1235         CALL histdef(nid_mth, "duoro", "Orography dU", "m/s2",
1236     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1237     .                "ave(X)", zsto,zout)
1238c
1239         CALL histdef(nid_mth, "dvoro", "Orography dV", "m/s2",
1240     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1241     .                "ave(X)", zsto,zout)
1242c
1243         ENDIF
1244C
1245         IF (ok_orolf) THEN
1246         CALL histdef(nid_mth, "dulif", "Orography dU", "m/s2",
1247     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1248     .                "ave(X)", zsto,zout)
1249c
1250         CALL histdef(nid_mth, "dvlif", "Orography dV", "m/s2",
1251     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1252     .                "ave(X)", zsto,zout)
1253         ENDIF
1254C
1255         CALL histdef(nid_mth, "ozone", "Ozone concentration", "-",
1256     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1257     .                "ave(X)", zsto,zout)
1258c
1259         if (nqmax.GE.3) THEN
1260         DO iq=1,nqmax-2
1261         IF (iq.LE.99) THEN
1262         WRITE(str2,'(i2.2)') iq
1263         CALL histdef(nid_mth, "trac"//str2, "Tracer No."//str2, "-",
1264     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1265     .                "ave(X)", zsto,zout)
1266         ELSE
1267         PRINT*, "Trop de traceurs"
1268         CALL abort
1269         ENDIF
1270         ENDDO
1271         ENDIF
1272c
1273cKE43
1274      IF (iflag_con.EQ.4) THEN ! (sb)
1275c
1276         CALL histdef(nid_mth, "upwd", "saturated updraft", "Kg/m2/s",
1277     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1278     .                "ave(X)", zsto,zout)
1279c
1280         CALL histdef(nid_mth, "dnwd", "saturated downdraft","Kg/m2/s",
1281     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1282     .                "ave(X)", zsto,zout)
1283c
1284         CALL histdef(nid_mth, "dnwd0", "unsat. downdraft", "Kg/m2/s",
1285     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1286     .                "ave(X)", zsto,zout)
1287c
1288         CALL histdef(nid_mth,"Ma","undilute adiab updraft","Kg/m2/s",
1289     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1290     .                "ave(X)", zsto,zout)
1291c
1292c
1293      ENDIF
1294c34EK
1295         CALL histend(nid_mth)
1296c
1297         ndex2d = 0
1298         ndex3d = 0
1299c
1300      ENDIF ! fin de test sur ok_mensuel
1301c
1302c
1303      IF (ok_instan) THEN
1304c
1305         idayref = day_ini
1306         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
1307c
1308         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
1309         DO i = 1, iim
1310            zx_lon(i,1) = rlon(i+1)
1311            zx_lon(i,jjmp1) = rlon(i+1)
1312         ENDDO
1313         DO ll=1,klev
1314            znivsig(ll)=float(ll)
1315         ENDDO
1316         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
1317         CALL histbeg("histins", iim,zx_lon, jjmp1,zx_lat,
1318     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
1319     .                 nhori, nid_ins)
1320         CALL histvert(nid_ins, "presnivs", "Vertical levels", "mb",
1321     .                 klev, presnivs, nvert)
1322c        call histvert(nid_ins, 'sig_s', 'Niveaux sigma','-',
1323c    .              klev, znivsig, nvert)
1324c
1325         zsto = dtime * ecrit_ins
1326         zout = dtime * ecrit_ins
1327C
1328         CALL histdef(nid_ins, "phis", "Surface geop. height", "-",
1329     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1330     .                "once", zsto,zout)
1331c
1332         CALL histdef(nid_ins, "aire", "Grid area", "-",
1333     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1334     .                "once", zsto,zout)
1335c
1336c Champs 2D:
1337c
1338         CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
1339     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1340     .                "inst(X)", zsto,zout)
1341c
1342         CALL histdef(nid_ins, "topl", "OLR", "W/m2",
1343     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
1344     .                "inst(X)", zsto,zout)
1345c
1346c Champs 3D:
1347c
1348         CALL histdef(nid_ins, "temp", "Temperature", "K",
1349     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1350     .                "inst(X)", zsto,zout)
1351c
1352         CALL histdef(nid_ins, "vitu", "Zonal wind", "m/s",
1353     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1354     .                "inst(X)", zsto,zout)
1355c
1356         CALL histdef(nid_ins, "vitv", "Merid wind", "m/s",
1357     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1358     .                "inst(X)", zsto,zout)
1359c
1360         CALL histdef(nid_ins, "geop", "Geopotential height", "m",
1361     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1362     .                "inst(X)", zsto,zout)
1363c
1364         CALL histdef(nid_ins, "pres", "Air pressure", "Pa",
1365     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
1366     .                "inst(X)", zsto,zout)
1367c
1368         CALL histend(nid_ins)
1369c
1370         ndex2d = 0
1371         ndex3d = 0
1372c
1373      ENDIF
1374c
1375c
1376c
1377c Prescrire l'ozone dans l'atmosphere
1378c
1379c
1380cc         DO i = 1, klon
1381cc         DO k = 1, klev
1382cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1383cc         ENDDO
1384cc         ENDDO
1385c
1386         IF (ok_oasis) THEN
1387         DO i = 1, klon
1388           oas_sols(i) = 0.0
1389           oas_nsol(i) = 0.0
1390           oas_rain(i) = 0.0
1391           oas_snow(i) = 0.0
1392           oas_evap(i) = 0.0
1393           oas_ruis(i) = 0.0
1394           oas_tsol(i) = 0.0
1395           oas_fder(i) = 0.0
1396           oas_albe(i) = 0.0
1397           oas_taux(i) = 0.0
1398           oas_tauy(i) = 0.0
1399         ENDDO
1400         ENDIF
1401c
1402      ENDIF
1403c
1404c   ****************     Fin  de   IF ( debut  )   ***************
1405c
1406c
1407c Mettre a zero des variables de sortie (pour securite)
1408c
1409      DO i = 1, klon
1410         d_ps(i) = 0.0
1411      ENDDO
1412      DO k = 1, klev
1413      DO i = 1, klon
1414         d_t(i,k) = 0.0
1415         d_u(i,k) = 0.0
1416         d_v(i,k) = 0.0
1417      ENDDO
1418      ENDDO
1419      DO iq = 1, nqmax
1420      DO k = 1, klev
1421      DO i = 1, klon
1422         d_qx(i,k,iq) = 0.0
1423      ENDDO
1424      ENDDO
1425      ENDDO
1426c
1427c Ne pas affecter les valeurs entrees de u, v, h, et q
1428c
1429      DO k = 1, klev
1430      DO i = 1, klon
1431         t_seri(i,k)  = t(i,k)
1432         u_seri(i,k)  = u(i,k)
1433         v_seri(i,k)  = v(i,k)
1434         q_seri(i,k)  = qx(i,k,ivap)
1435         ql_seri(i,k) = qx(i,k,iliq)
1436      ENDDO
1437      ENDDO
1438      IF (nqmax.GE.3) THEN
1439      DO iq = 3, nqmax
1440      DO  k = 1, klev
1441      DO  i = 1, klon
1442         tr_seri(i,k,iq-2) = qx(i,k,iq)
1443      ENDDO
1444      ENDDO
1445      ENDDO
1446      ELSE
1447      DO k = 1, klev
1448      DO i = 1, klon
1449         tr_seri(i,k,1) = 0.0
1450      ENDDO
1451      ENDDO
1452      ENDIF
1453c
1454c Diagnostiquer la tendance dynamique
1455c
1456      IF (ancien_ok) THEN
1457         DO k = 1, klev
1458         DO i = 1, klon
1459            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1460            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1461         ENDDO
1462         ENDDO
1463      ELSE
1464         DO k = 1, klev
1465         DO i = 1, klon
1466            d_t_dyn(i,k) = 0.0
1467            d_q_dyn(i,k) = 0.0
1468         ENDDO
1469         ENDDO
1470         ancien_ok = .TRUE.
1471      ENDIF
1472c
1473c Ajouter le geopotentiel du sol:
1474c
1475      DO k = 1, klev
1476      DO i = 1, klon
1477         zphi(i,k) = pphi(i,k) + pphis(i)
1478      ENDDO
1479      ENDDO
1480c
1481c Verifier les temperatures
1482c
1483
1484      CALL hgardfou(t_seri,ftsol,'debutphy')
1485c
1486c Incrementer le compteur de la physique
1487c
1488      itap   = itap + 1
1489      julien = MOD(NINT(xjour),360)
1490c
1491c Mettre en action les conditions aux limites (albedo, sst, etc.).
1492c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1493c
1494      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1495         idayvrai = NINT(xjour)
1496         PRINT *,' PHYS cond  julien ',julien,idayvrai,ok_limitvrai
1497         CALL condsurf(julien,idayvrai, pctsrf ,
1498     .                  lmt_sst,lmt_alb,lmt_rug,lmt_bils  )
1499         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1500      ENDIF
1501cccccccccc
1502      IF (ok_oasis .AND. MOD(itap-1,nexca).EQ.0) THEN
1503C
1504         CALL fromcpl(itap,jjmp1*iim,
1505     .        cpl_sst,cpl_sic,cpl_alb_sst,cpl_alb_sic)
1506         DO i = 1, iim-1 ! un seul point pour le pole nord
1507            cpl_sst(i,1) = cpl_sst(iim,1)
1508            cpl_sic(i,1) = cpl_sic(iim,1)
1509            cpl_alb_sst(i,1) = cpl_alb_sst(iim,1)
1510            cpl_alb_sic(i,1) = cpl_alb_sic(iim,1)
1511         ENDDO
1512         DO i = 2, iim ! un seul point pour le pole sud
1513            cpl_sst(i,jjmp1) = cpl_sst(1,jjmp1)
1514            cpl_sic(i,jjmp1) = cpl_sic(1,jjmp1)
1515            cpl_alb_sst(i,jjmp1) = cpl_alb_sst(1,jjmp1)
1516            cpl_alb_sic(i,jjmp1) = cpl_alb_sic(1,jjmp1)
1517         ENDDO
1518c
1519         ig = 1
1520         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
1521     .       pctsrf(ig,is_sic).GT.epsfra) THEN
1522            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
1523     .                        - (cpl_sic(1,1)-pctsrf(ig,is_sic))
1524            pctsrf(ig,is_sic) = cpl_sic(1,1)
1525            lmt_sst(ig) = cpl_sst(1,1)
1526         ENDIF
1527         DO j = 2, jjm
1528         DO i = 1, iim
1529         ig = ig + 1
1530         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
1531     .       pctsrf(ig,is_sic).GT.epsfra) THEN
1532           pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
1533     .                       - (cpl_sic(i,j)-pctsrf(ig,is_sic))
1534           pctsrf(ig,is_sic) = cpl_sic(i,j)
1535           lmt_sst(ig) = cpl_sst(i,j)
1536         ENDIF
1537         ENDDO
1538         ENDDO
1539         ig = ig + 1
1540         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
1541     .       pctsrf(ig,is_sic).GT.epsfra) THEN
1542            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
1543     .                        - (cpl_sic(1,jjmp1)-pctsrf(ig,is_sic))
1544            pctsrf(ig,is_sic) = cpl_sic(1,jjmp1)
1545            lmt_sst(ig) = cpl_sst(1,jjmp1)
1546         ENDIF
1547c
1548      ENDIF ! ok_oasis
1549cccccccccc
1550c
1551 
1552      IF (ok_ocean) THEN
1553         DO i = 1, klon
1554            ftsol(i,is_oce) = lmt_sst(i) + deltat(i)
1555         ENDDO
1556
1557      ELSE
1558         DO i = 1, klon
1559            ftsol(i,is_oce) = lmt_sst(i)
1560         ENDDO
1561
1562      ENDIF
1563c
1564c Re-evaporer l'eau liquide nuageuse
1565c
1566      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1567      DO i = 1, klon
1568         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1569         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1570         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1571         zb = MAX(0.0,ql_seri(i,k))
1572         za = - MAX(0.0,ql_seri(i,k))
1573     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1574         t_seri(i,k) = t_seri(i,k) + za
1575         q_seri(i,k) = q_seri(i,k) + zb
1576         ql_seri(i,k) = 0.0
1577         d_t_eva(i,k) = za
1578         d_q_eva(i,k) = zb
1579      ENDDO
1580      ENDDO
1581c
1582c Appeler la diffusion verticale (programme de couche limite)
1583c
1584      DO i = 1, klon
1585         frugs(i,is_ter) = SQRT(lmt_rug(i)**2+rugoro(i)**2)
1586         frugs(i,is_lic) = rugoro(i)
1587         frugs(i,is_oce) = rugmer(i)
1588         frugs(i,is_sic) = 0.001
1589         zxrugs(i) = 0.0
1590      ENDDO
1591      DO nsrf = 1, nbsrf
1592      DO i = 1, klon
1593         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1594      ENDDO
1595      ENDDO
1596      DO nsrf = 1, nbsrf
1597      DO i = 1, klon
1598         zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1599      ENDDO
1600      ENDDO
1601c
1602      CALL clmain(dtime,pctsrf,
1603     e            t_seri,q_seri,u_seri,v_seri,soil_model,
1604     e            ftsol,soilcap,soilflux,paprs,pplay,radsol,
1605     e            fsnow,fqsol,
1606     e            rlat, frugs,
1607     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1608     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,rugmer,
1609     s            dsens, devap,
1610     s            ycoefh,yu1,yv1)
1611c
1612      DO i = 1, klon
1613         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1614         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1615         fder(i) = dsens(i) + devap(i)
1616      ENDDO
1617      DO k = 1, klev
1618      DO i = 1, klon
1619         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1620         q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1621         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1622         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1623      ENDDO
1624      ENDDO
1625c
1626c Incrementer la temperature du sol
1627c
1628      DO i = 1, klon
1629         zxtsol(i) = 0.0
1630         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1631     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
1632     $       THEN
1633             WRITE(*,*) 'physiq : pb sous surface au point ', i,
1634     $           pctsrf(i, 1 : nbsrf)
1635         ENDIF
1636      ENDDO
1637      DO nsrf = 1, nbsrf
1638      DO i = 1, klon
1639         ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1640         zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1641      ENDDO
1642      ENDDO
1643
1644c
1645c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1646c
1647      DO nsrf = 1, nbsrf
1648        DO i = 1, klon
1649          IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1650        ENDDO
1651      ENDDO
1652
1653c
1654c Appeler le modele du sol
1655c
1656      IF (soil_model) THEN
1657         DO nsrf = 1, nbsrf
1658         CALL soil(dtime, nsrf, fsnow(1,nsrf),
1659     .             ftsol(1,nsrf), ftsoil(1,1,nsrf),
1660     .             soilcap(1,nsrf), soilflux(1,nsrf))
1661         ENDDO
1662      ENDIF
1663c
1664c Calculer la derive du flux infrarouge
1665c
1666      DO nsrf = 1, nbsrf
1667      DO i = 1, klon
1668         fder(i) = fder(i) - 4.0*RSIGMA*zxtsol(i)**3 *
1669     .                       (ftsol(i,nsrf)-zxtsol(i))
1670     .                      *pctsrf(i,nsrf)
1671      ENDDO
1672      ENDDO
1673c
1674c Appeler la convection (au choix)
1675c
1676      DO k = 1, klev
1677      DO i = 1, klon
1678         conv_q(i,k) = d_q_dyn(i,k)
1679     .               + d_q_vdf(i,k)/dtime
1680         conv_t(i,k) = d_t_dyn(i,k)
1681     .               + d_t_vdf(i,k)/dtime
1682      ENDDO
1683      ENDDO
1684      IF (check) THEN
1685         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1686         PRINT*, "avantcon=", za
1687      ENDIF
1688      zx_ajustq = .FALSE.
1689      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1690      IF (zx_ajustq) THEN
1691         DO i = 1, klon
1692            z_avant(i) = 0.0
1693         ENDDO
1694         DO k = 1, klev
1695         DO i = 1, klon
1696            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1697     .                        *(paprs(i,k)-paprs(i,k+1))/RG
1698         ENDDO
1699         ENDDO
1700      ENDIF
1701      IF (iflag_con.EQ.1) THEN
1702          stop'reactiver le call conlmd dans physiq.F'
1703c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1704c    .             d_t_con, d_q_con,
1705c    .             rain_con, snow_con, ibas_con, itop_con)
1706      ELSE IF (iflag_con.EQ.2) THEN
1707      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1708     e            conv_t, conv_q, fluxq(1,1), omega,
1709     s            d_t_con, d_q_con, rain_con, snow_con,
1710     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1711     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
1712      DO i = 1, klon
1713         ibas_con(i) = klev+1 - kcbot(i)
1714         itop_con(i) = klev+1 - kctop(i)
1715      ENDDO
1716      ELSE IF (iflag_con.EQ.3) THEN
1717          stop'reactiver le call conlmd dans physiq.F'
1718c     CALL conccm (dtime,paprs,pplay,t_seri,q_seri,conv_q,
1719c    s             d_t_con, d_q_con,
1720c    s             rain_con, snow_con, ibas_con, itop_con)
1721cKE43
1722      ELSE IF (iflag_con.EQ.4) THEN
1723c nb of tracers for the KE convection:
1724          if (nqmax .GE. 4) then
1725              ntra = nbtr
1726          else
1727              ntra = 1
1728          endif
1729cke43 (arguments inutiles enleves => des SAVE dans conema43?)
1730c$$$          CALL conema43(dtime,paprs,pplay,t_seri,q_seri,
1731c$$$     $        u_seri,v_seri,tr_seri,nbtr,
1732c$$$     .        ema_workcbmf,
1733c$$$     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1734c$$$     .        wdn, tdn, qdn,
1735c$$$     .        rain_con, snow_con, ibas_con, itop_con,
1736c$$$     .        upwd,dnwd,dnwd0,bas,top,Ma,cape,tvp,rflag,
1737c$$$     .        pbase
1738c$$$     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
1739c$$$     .        pori_con,plcl_con,dtma_con,dtlcl_con)
1740          CALL conema (dtime,paprs,pplay,t_seri,q_seri,
1741     $        u_seri,v_seri,tr_seri,nbtr,
1742     .        ema_work1,ema_work2,
1743     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1744c$$$     .        wdn, tdn, qdn,
1745     .        rain_con, snow_con, ibas_con, itop_con,
1746     .        upwd,dnwd,dnwd0,bas,top,Ma,cape,tvp,rflag,
1747     .        pbase
1748     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
1749c$$$     .        pori_con,plcl_con,dtma_con,dtlcl_con)
1750          DO i = 1, klon
1751            ema_pcb(i)  = pbase(i)
1752          ENDDO
1753          DO i = 1, klon
1754            ema_pct(i)  = paprs(i,itop_con(i))
1755          ENDDO
1756          DO i = 1, klon
1757            ema_cbmf(i) = ema_workcbmf(i)
1758          ENDDO     
1759      ELSE
1760          PRINT*, "iflag_con non-prevu", iflag_con
1761          CALL abort
1762      ENDIF
1763      CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1764     .              d_u_con, d_v_con)
1765      DO k = 1, klev
1766        DO i = 1, klon
1767          t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1768          q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1769          u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1770          v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1771        ENDDO
1772      ENDDO
1773      IF (check) THEN
1774          za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1775          PRINT*, "aprescon=", za
1776          zx_t = 0.0
1777          za = 0.0
1778          DO i = 1, klon
1779            za = za + paire(i)/FLOAT(klon)
1780            zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1781          ENDDO
1782          zx_t = zx_t/za*dtime
1783          PRINT*, "Precip=", zx_t
1784      ENDIF
1785      IF (zx_ajustq) THEN
1786          DO i = 1, klon
1787            z_apres(i) = 0.0
1788          ENDDO
1789          DO k = 1, klev
1790            DO i = 1, klon
1791              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1792     .            *(paprs(i,k)-paprs(i,k+1))/RG
1793            ENDDO
1794          ENDDO
1795          DO i = 1, klon
1796            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1797     .          /z_apres(i)
1798          ENDDO
1799          DO k = 1, klev
1800            DO i = 1, klon
1801              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1802     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
1803                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
1804              ENDIF
1805            ENDDO
1806          ENDDO
1807      ENDIF
1808      zx_ajustq=.FALSE.
1809c
1810      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1811c
1812          IF (iflag_con .NE. 2 .AND.  iflag_con .NE. 4 ) THEN
1813              PRINT*, 'Pour l instant, seul conflx fonctionne ',
1814     $            'avec traceurs', iflag_con
1815              PRINT*,' Mettre iflag_con',
1816     $            ' = 2  ou 4 dans run.def et repasser'
1817              CALL abort
1818              ENDIF
1819c
1820      ENDIF !--nqmax.GT.2
1821c
1822c Appeler l'ajustement sec
1823c
1824      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1825      DO k = 1, klev
1826      DO i = 1, klon
1827         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1828         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1829      ENDDO
1830      ENDDO
1831
1832c   RATQS
1833      call calcratqs (
1834     I            paprs,pplay,q_seri,d_t_con,d_t_ajs
1835     O           ,ratqs,zpt_conv)
1836c
1837c Appeler le processus de condensation a grande echelle
1838c et le processus de precipitation
1839c
1840      CALL fisrtilp_tr(dtime,paprs,pplay,
1841     .           t_seri, q_seri,ratqs,
1842     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1843     .           rain_lsc, snow_lsc,
1844     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
1845     .           frac_impa, frac_nucl,
1846     .           prfl, psfl)
1847      DO k = 1, klev
1848      DO i = 1, klon
1849         t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1850         q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1851         ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1852         cldfra(i,k) = rneb(i,k)
1853         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1854      ENDDO
1855      ENDDO
1856      IF (check) THEN
1857         za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1858         PRINT*, "apresilp=", za
1859         zx_t = 0.0
1860         za = 0.0
1861         DO i = 1, klon
1862            za = za + paire(i)/FLOAT(klon)
1863            zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
1864        ENDDO
1865         zx_t = zx_t/za*dtime
1866         PRINT*, "Precip=", zx_t
1867      ENDIF
1868c
1869c Nuages diagnostiques:
1870c
1871      IF (iflag_con.EQ.2) THEN ! seulement pour Tiedtke
1872      CALL diagcld1(paprs,pplay,
1873     .             rain_con,snow_con,ibas_con,itop_con,
1874     .             diafra,dialiq)
1875      DO k = 1, klev
1876      DO i = 1, klon
1877      IF (diafra(i,k).GT.cldfra(i,k)) THEN
1878         cldliq(i,k) = dialiq(i,k)
1879         cldfra(i,k) = diafra(i,k)
1880      ENDIF
1881      ENDDO
1882      ENDDO
1883      ENDIF
1884c
1885c Nuages stratus artificiels:
1886c
1887      IF (ok_stratus) THEN
1888      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
1889      DO k = 1, klev
1890      DO i = 1, klon
1891      IF (diafra(i,k).GT.cldfra(i,k)) THEN
1892         cldliq(i,k) = dialiq(i,k)
1893         cldfra(i,k) = diafra(i,k)
1894      ENDIF
1895      ENDDO
1896      ENDDO
1897      ENDIF
1898c
1899c Precipitation totale
1900c
1901      DO i = 1, klon
1902         rain_fall(i) = rain_con(i) + rain_lsc(i)
1903         snow_fall(i) = snow_con(i) + snow_lsc(i)
1904      ENDDO
1905c
1906c Calculer l'humidite relative pour diagnostique
1907c
1908      DO k = 1, klev
1909      DO i = 1, klon
1910         zx_t = t_seri(i,k)
1911         IF (thermcep) THEN
1912            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1913            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1914            zx_qs  = MIN(0.5,zx_qs)
1915            zcor   = 1./(1.-retv*zx_qs)
1916            zx_qs  = zx_qs*zcor
1917         ELSE
1918           IF (zx_t.LT.t_coup) THEN
1919              zx_qs = qsats(zx_t)/pplay(i,k)
1920           ELSE
1921              zx_qs = qsatl(zx_t)/pplay(i,k)
1922           ENDIF
1923         ENDIF
1924         zx_rh(i,k) = q_seri(i,k)/zx_qs
1925      ENDDO
1926      ENDDO
1927c
1928c Calculer les parametres optiques des nuages et quelques
1929c parametres pour diagnostiques:
1930c
1931      CALL nuage (paprs, pplay,
1932     .            t_seri, cldliq, cldfra, cldtau, cldemi,
1933     .            cldh, cldl, cldm, cldt, cldq)
1934c
1935c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1936c
1937      IF (MOD(itaprad,radpas).EQ.0) THEN
1938      CALL orbite(FLOAT(julien),zlongi,dist)
1939      IF (cycle_diurne) THEN
1940        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1941        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1942c        CALL zenith(zlongi,gmtime,rlat,rlon,rmu0,fract) !va disparaitre
1943        CALL alboc_cd(rmu0,alb_eau)
1944      ELSE
1945        CALL angle(zlongi,rlat,fract,rmu0)
1946        CALL alboc(FLOAT(julien),rlat,alb_eau)
1947      ENDIF
1948      CALL albsno(veget,agesno,alb_neig)
1949      DO i = 1, klon
1950         zx_alb_oce = alb_eau(i)
1951         IF (pctsrf(i,is_oce).GT.epsfra .AND. ftsol(i,is_oce).LT.271.35)
1952     .   zx_alb_oce = 0.6 ! pour slab_ocean
1953         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_lic)/(fsnow(i,is_lic)+10.0)))
1954         zx_alb_lic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
1955         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_ter)/(fsnow(i,is_ter)+10.0)))
1956         zx_alb_ter = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
1957         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_sic)/(fsnow(i,is_sic)+10.0)))
1958         zx_alb_sic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
1959         albsol(i) = zx_alb_oce * pctsrf(i,is_oce)
1960     .             + zx_alb_lic * pctsrf(i,is_lic)
1961     .             + zx_alb_ter * pctsrf(i,is_ter)
1962     .             + zx_alb_sic * pctsrf(i,is_sic)
1963      ENDDO
1964      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
1965     e            (dist, rmu0, fract, co2_ppm, solaire,
1966     e             paprs, pplay,zxtsol,albsol, t_seri,q_seri,wo,
1967     e             cldfra, cldemi, cldtau,
1968     s             heat,heat0,cool,cool0,radsol,albpla,
1969     s             topsw,toplw,solsw,sollw,
1970     s             topsw0,toplw0,solsw0,sollw0)
1971      itaprad = 0
1972      ENDIF
1973      itaprad = itaprad + 1
1974c
1975c Ajouter la tendance des rayonnements (tous les pas)
1976c
1977      DO k = 1, klev
1978      DO i = 1, klon
1979         t_seri(i,k) = t_seri(i,k)
1980     .               + (heat(i,k)-cool(i,k)) * dtime/86400.
1981      ENDDO
1982      ENDDO
1983c
1984c Calculer l'hydrologie de la surface
1985c
1986      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, evap,
1987     .            agesno, ftsol,fqsol,fsnow, ruis)
1988c
1989      DO i = 1, klon
1990         zxqsol(i) = 0.0
1991         zxsnow(i) = 0.0
1992      ENDDO
1993      DO nsrf = 1, nbsrf
1994      DO i = 1, klon
1995         zxqsol(i) = zxqsol(i) + fqsol(i,nsrf)*pctsrf(i,nsrf)
1996         zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
1997      ENDDO
1998      ENDDO
1999c
2000c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2001c
2002      DO nsrf = 1, nbsrf
2003      DO i = 1, klon
2004         IF (pctsrf(i,nsrf).LT.epsfra) THEN
2005            fqsol(i,nsrf) = zxqsol(i)
2006            fsnow(i,nsrf) = zxsnow(i)
2007         ENDIF
2008      ENDDO
2009      ENDDO
2010c
2011c Calculer le bilan du sol et la derive de temperature (couplage)
2012c
2013      DO i = 1, klon
2014         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2015      ENDDO
2016      IF (ok_ocean) THEN
2017         DO i = 1, klon
2018            cthermiq = cyang
2019            IF (ftsol(i,is_oce).LT. 271.35) cthermiq = cbing
2020            IF (pctsrf(i,is_oce).GT.epsfra) deltat(i) = deltat(i) +
2021     .                          (bils(i)-lmt_bils(i))/cthermiq * dtime
2022            IF (deltat(i).GT.15.0 ) deltat(i) = 15.0
2023            IF (deltat(i).LT.-15.0) deltat(i) = -15.0
2024         ENDDO
2025      ENDIF
2026c
2027cmoddeblott(jan95)
2028c Appeler le programme de parametrisation de l'orographie
2029c a l'echelle sous-maille:
2030c
2031      IF (ok_orodr) THEN
2032c
2033c  selection des points pour lesquels le shema est actif:
2034        igwd=0
2035        DO i=1,klon
2036        itest(i)=0
2037c        IF ((zstd(i).gt.10.0)) THEN
2038        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2039          itest(i)=1
2040          igwd=igwd+1
2041          idx(igwd)=i
2042        ENDIF
2043        ENDDO
2044        igwdim=MAX(1,igwd)
2045c
2046        CALL drag_noro(klon,klev,dtime,paprs,pplay,
2047     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
2048     e                   igwd,igwdim,idx,itest,
2049     e                   t_seri, u_seri, v_seri,
2050     s                   zulow, zvlow, zustr, zvstr,
2051     s                   d_t_oro, d_u_oro, d_v_oro)
2052c
2053c  ajout des tendances
2054        DO k = 1, klev
2055        DO i = 1, klon
2056           t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2057           u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2058           v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2059        ENDDO
2060        ENDDO
2061c
2062      ENDIF ! fin de test sur ok_orodr
2063c
2064      IF (ok_orolf) THEN
2065c
2066c  selection des points pour lesquels le shema est actif:
2067        igwd=0
2068        DO i=1,klon
2069        itest(i)=0
2070        IF ((zpic(i)-zmea(i)).GT.100.) THEN
2071          itest(i)=1
2072          igwd=igwd+1
2073          idx(igwd)=i
2074        ENDIF
2075        ENDDO
2076        igwdim=MAX(1,igwd)
2077c
2078        CALL lift_noro(klon,klev,dtime,paprs,pplay,
2079     e                   rlat,zmea,zstd, zsig, zgam, zthe,zpic,zval,
2080     e                   igwd,igwdim,idx,itest,
2081     e                   t_seri, u_seri, v_seri,
2082     s                   zulow, zvlow, zustr, zvstr,
2083     s                   d_t_lif, d_u_lif, d_v_lif)
2084c
2085c  ajout des tendances
2086        DO k = 1, klev
2087        DO i = 1, klon
2088           t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2089           u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2090           v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2091        ENDDO
2092        ENDDO
2093c
2094      ENDIF ! fin de test sur ok_orolf
2095c
2096cAA
2097cAA Installation de l'interface online-offline pour traceurs
2098cAA
2099c====================================================================
2100c   Calcul  des tendances traceurs
2101c====================================================================
2102C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2103cKE43       des traceurs => il faut donc mettre des flags a .false.
2104      IF (iflag_con.EQ.4) THEN
2105c           on ajoute les tendances calculees par KE43
2106        DO iq=1, nqmax-2 ! Sandrine a -3 ???
2107        DO k = 1, nlev
2108        DO i = 1, klon
2109          tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2110        ENDDO
2111        ENDDO
2112        WRITE(iqn,'(i2.2)') iq
2113        CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2114        ENDDO
2115CMAF modif pour garder info du nombre de traceurs auxquels
2116C la physique s'applique
2117      ELSE
2118CMAF modif pour garder info du nombre de traceurs auxquels
2119C la physique s'applique
2120C
2121      call phytrac (rnpb,
2122     I                   debut,lafin,
2123     I                   nqmax-2,
2124     I                   nlon,nlev,dtime,
2125     I                   t,paprs,pplay,
2126     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2127     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2128     I                   frac_impa, frac_nucl,
2129     I                   rlon,presnivs,paire,pphis,
2130     O                   tr_seri)
2131      ENDIF
2132
2133      IF (offline) THEN
2134
2135         call phystokenc (
2136     I                   nlon,nlev,pdtphys,rlon,rlat,
2137     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2138     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
2139     I                   frac_impa, frac_nucl,
2140     I                   pphis,paire,dtime,itap)
2141
2142
2143      ENDIF
2144
2145c
2146c Calculer le transport de l'eau et de l'energie (diagnostique)
2147c
2148      CALL transp (paprs,zxtsol,
2149     e                   t_seri, q_seri, u_seri, v_seri, zphi,
2150     s                   ve, vq, ue, uq)
2151c
2152c Accumuler les variables a stocker dans les fichiers histoire:
2153c
2154      IF (ok_oasis) THEN ! couplage oasis
2155      DO i = 1, klon
2156        oas_sols(i) = oas_sols(i) + solsw(i)          / FLOAT(nexca)
2157        oas_nsol(i) = oas_nsol(i) + (bils(i)-solsw(i))/ FLOAT(nexca)
2158        oas_rain(i) = oas_rain(i) + rain_fall(i)      / FLOAT(nexca)
2159        oas_snow(i) = oas_snow(i) + snow_fall(i)      / FLOAT(nexca)
2160        oas_evap(i) = oas_evap(i) + evap(i)           / FLOAT(nexca)
2161        oas_tsol(i) = oas_tsol(i) + zxtsol(i)         / FLOAT(nexca)
2162        oas_fder(i) = oas_fder(i) + fder(i)           / FLOAT(nexca)
2163        oas_albe(i) = oas_albe(i) + albsol(i)         / FLOAT(nexca)
2164        oas_taux(i) = oas_taux(i) + fluxu(i,1)        / FLOAT(nexca)
2165        oas_tauy(i) = oas_tauy(i) + fluxv(i,1)        / FLOAT(nexca)
2166        oas_ruis(i) = oas_ruis(i) + ruis(i)       /FLOAT(nexca)/dtime
2167      ENDDO
2168      ENDIF
2169c
2170c
2171      IF (ok_journe) THEN
2172c
2173      ndex2d = 0
2174      ndex3d = 0
2175c
2176c Champs 2D:
2177c
2178         i = NINT(zout/zsto)
2179         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2180         CALL histwrite(nid_day,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2181c
2182         i = NINT(zout/zsto)
2183         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2184         CALL histwrite(nid_day,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2185C
2186      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2187      CALL histwrite(nid_day,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2188c
2189      DO i = 1, klon
2190         zx_tmp_fi2d(i) = paprs(i,1)
2191      ENDDO
2192      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2193      CALL histwrite(nid_day,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2194c
2195      DO i = 1, klon
2196         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
2197      ENDDO
2198      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2199      CALL histwrite(nid_day,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2200c
2201      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
2202      CALL histwrite(nid_day,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2203c
2204      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2205      CALL histwrite(nid_day,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2206c
2207      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
2208      CALL histwrite(nid_day,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2209c
2210      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2211      CALL histwrite(nid_day,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2212c
2213      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2214      CALL histwrite(nid_day,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2215c
2216      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2217      CALL histwrite(nid_day,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2218c
2219      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2220      CALL histwrite(nid_day,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2221c
2222      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2223      CALL histwrite(nid_day,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2224c
2225      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2226      CALL histwrite(nid_day,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2227c
2228      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
2229      CALL histwrite(nid_day,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2230c
2231      DO i = 1, klon
2232         zx_tmp_fi2d(i) = fluxu(i,1)
2233      ENDDO
2234      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2235      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2236c
2237      DO i = 1, klon
2238         zx_tmp_fi2d(i) = fluxv(i,1)
2239      ENDDO
2240      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2241      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2242c
2243      DO i = 1, klon
2244         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
2245      ENDDO
2246      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2247      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2248c
2249      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
2250      CALL histwrite(nid_day,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2251c
2252      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
2253      CALL histwrite(nid_day,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2254c
2255      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
2256      CALL histwrite(nid_day,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2257c
2258      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
2259      CALL histwrite(nid_day,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2260c
2261      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
2262      CALL histwrite(nid_day,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2263c
2264c Champs 3D:
2265c
2266      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2267      CALL histwrite(nid_day,"temp",itap,zx_tmp_3d,
2268     .                                   iim*jjmp1*klev,ndex3d)
2269c
2270      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
2271      CALL histwrite(nid_day,"ovap",itap,zx_tmp_3d,
2272     .                                   iim*jjmp1*klev,ndex3d)
2273c
2274      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2275      CALL histwrite(nid_day,"geop",itap,zx_tmp_3d,
2276     .                                   iim*jjmp1*klev,ndex3d)
2277c
2278      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2279      CALL histwrite(nid_day,"vitu",itap,zx_tmp_3d,
2280     .                                   iim*jjmp1*klev,ndex3d)
2281c
2282      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2283      CALL histwrite(nid_day,"vitv",itap,zx_tmp_3d,
2284     .                                   iim*jjmp1*klev,ndex3d)
2285c
2286      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
2287      CALL histwrite(nid_day,"vitw",itap,zx_tmp_3d,
2288     .                                   iim*jjmp1*klev,ndex3d)
2289c
2290      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2291      CALL histwrite(nid_day,"pres",itap,zx_tmp_3d,
2292     .                                   iim*jjmp1*klev,ndex3d)
2293c
2294      if (ok_sync) then
2295        call histsync(nid_day)
2296      endif
2297      ENDIF
2298C
2299      IF (ok_mensuel) THEN
2300c
2301      ndex2d = 0
2302      ndex3d = 0
2303c
2304c Champs 2D:
2305c
2306         i = NINT(zout/zsto)
2307         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2308         CALL histwrite(nid_mth,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2309C
2310         i = NINT(zout/zsto)
2311         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2312         CALL histwrite(nid_mth,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2313
2314      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
2315      CALL histwrite(nid_mth,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2316c
2317      DO i = 1, klon
2318         zx_tmp_fi2d(i) = paprs(i,1)
2319      ENDDO
2320      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2321      CALL histwrite(nid_mth,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2322c
2323      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxqsol,zx_tmp_2d)
2324      CALL histwrite(nid_mth,"qsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2325c
2326      DO i = 1, klon
2327         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
2328      ENDDO
2329      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2330      CALL histwrite(nid_mth,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2331c
2332      DO i = 1, klon
2333         zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
2334      ENDDO
2335      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2336      CALL histwrite(nid_mth,"plul",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2337c
2338      DO i = 1, klon
2339         zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
2340      ENDDO
2341      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2342      CALL histwrite(nid_mth,"pluc",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2343c
2344      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
2345      CALL histwrite(nid_mth,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2346c
2347      CALL gr_fi_ecrit(1, klon,iim,jjmp1, agesno,zx_tmp_2d)
2348      CALL histwrite(nid_mth,"ages",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2349c
2350      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
2351      CALL histwrite(nid_mth,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2352c
2353      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
2354      CALL histwrite(nid_mth,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2355c
2356      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2357      CALL histwrite(nid_mth,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2358c
2359      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
2360      CALL histwrite(nid_mth,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2361c
2362      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
2363      CALL histwrite(nid_mth,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2364c
2365      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
2366      CALL histwrite(nid_mth,"tops0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2367c
2368      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw0,zx_tmp_2d)
2369      CALL histwrite(nid_mth,"topl0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2370c
2371      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw0,zx_tmp_2d)
2372      CALL histwrite(nid_mth,"sols0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2373c
2374      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw0,zx_tmp_2d)
2375      CALL histwrite(nid_mth,"soll0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2376c
2377      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
2378      CALL histwrite(nid_mth,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2379c
2380      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
2381      CALL histwrite(nid_mth,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2382c
2383      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
2384      CALL histwrite(nid_mth,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2385c
2386      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
2387      CALL histwrite(nid_mth,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2388c
2389      DO i = 1, klon
2390         zx_tmp_fi2d(i) = fluxu(i,1)
2391      ENDDO
2392      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2393      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2394c
2395      DO i = 1, klon
2396         zx_tmp_fi2d(i) = fluxv(i,1)
2397      ENDDO
2398      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2399      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2400c
2401      DO i = 1, klon
2402         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
2403      ENDDO
2404      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2405      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2406c
2407      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
2408      CALL histwrite(nid_mth,"albs",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2409c
2410      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
2411      CALL histwrite(nid_mth,"cdrm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2412c
2413      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
2414      CALL histwrite(nid_mth,"cdrh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2415c
2416      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
2417      CALL histwrite(nid_mth,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2418c
2419      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
2420      CALL histwrite(nid_mth,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2421c
2422      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
2423      CALL histwrite(nid_mth,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2424c
2425      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
2426      CALL histwrite(nid_mth,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2427c
2428      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
2429      CALL histwrite(nid_mth,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2430c
2431      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ue,zx_tmp_2d)
2432      CALL histwrite(nid_mth,"ue",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2433c
2434      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ve,zx_tmp_2d)
2435      CALL histwrite(nid_mth,"ve",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2436c
2437      CALL gr_fi_ecrit(1, klon,iim,jjmp1, uq,zx_tmp_2d)
2438      CALL histwrite(nid_mth,"uq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2439c
2440      CALL gr_fi_ecrit(1, klon,iim,jjmp1, vq,zx_tmp_2d)
2441      CALL histwrite(nid_mth,"vq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2442cKE43
2443      IF (iflag_con .EQ. 4) THEN ! sb
2444c
2445      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cape,zx_tmp_2d)
2446      CALL histwrite(nid_mth,"cape",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2447c
2448      CALL gr_fi_ecrit(1, klon,iim,jjmp1,pbase,zx_tmp_2d)
2449      CALL histwrite(nid_mth,"pbase",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2450c
2451      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_pct,zx_tmp_2d)
2452      CALL histwrite(nid_mth,"ptop",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2453c
2454      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_cbmf,zx_tmp_2d)
2455      CALL histwrite(nid_mth,"fbase",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2456c
2457c
2458      ENDIF
2459c34EK
2460c
2461c Champs 3D:
2462C
2463      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2464      CALL histwrite(nid_mth,"temp",itap,zx_tmp_3d,
2465     .                                   iim*jjmp1*klev,ndex3d)
2466c
2467      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
2468      CALL histwrite(nid_mth,"ovap",itap,zx_tmp_3d,
2469     .                                   iim*jjmp1*klev,ndex3d)
2470c
2471      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2472      CALL histwrite(nid_mth,"geop",itap,zx_tmp_3d,
2473     .                                   iim*jjmp1*klev,ndex3d)
2474c
2475      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2476      CALL histwrite(nid_mth,"vitu",itap,zx_tmp_3d,
2477     .                                   iim*jjmp1*klev,ndex3d)
2478c
2479      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2480      CALL histwrite(nid_mth,"vitv",itap,zx_tmp_3d,
2481     .                                   iim*jjmp1*klev,ndex3d)
2482c
2483      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
2484      CALL histwrite(nid_mth,"vitw",itap,zx_tmp_3d,
2485     .                                   iim*jjmp1*klev,ndex3d)
2486c
2487      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2488      CALL histwrite(nid_mth,"pres",itap,zx_tmp_3d,
2489     .                                   iim*jjmp1*klev,ndex3d)
2490c
2491      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldfra, zx_tmp_3d)
2492      CALL histwrite(nid_mth,"rneb",itap,zx_tmp_3d,
2493     .                                   iim*jjmp1*klev,ndex3d)
2494c
2495      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_rh, zx_tmp_3d)
2496      CALL histwrite(nid_mth,"rhum",itap,zx_tmp_3d,
2497     .                                   iim*jjmp1*klev,ndex3d)
2498c
2499      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldliq, zx_tmp_3d)
2500      CALL histwrite(nid_mth,"oliq",itap,zx_tmp_3d,
2501     .                                   iim*jjmp1*klev,ndex3d)
2502c
2503      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_dyn, zx_tmp_3d)
2504      CALL histwrite(nid_mth,"dtdyn",itap,zx_tmp_3d,
2505     .                                   iim*jjmp1*klev,ndex3d)
2506c
2507      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_dyn, zx_tmp_3d)
2508      CALL histwrite(nid_mth,"dqdyn",itap,zx_tmp_3d,
2509     .                                   iim*jjmp1*klev,ndex3d)
2510c
2511      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_con, zx_tmp_3d)
2512      CALL histwrite(nid_mth,"dtcon",itap,zx_tmp_3d,
2513     .                                   iim*jjmp1*klev,ndex3d)
2514c
2515      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_con, zx_tmp_3d)
2516      CALL histwrite(nid_mth,"dqcon",itap,zx_tmp_3d,
2517     .                                   iim*jjmp1*klev,ndex3d)
2518c
2519      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_lsc, zx_tmp_3d)
2520      CALL histwrite(nid_mth,"dtlsc",itap,zx_tmp_3d,
2521     .                                   iim*jjmp1*klev,ndex3d)
2522c
2523      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_lsc, zx_tmp_3d)
2524      CALL histwrite(nid_mth,"dqlsc",itap,zx_tmp_3d,
2525     .                                   iim*jjmp1*klev,ndex3d)
2526c
2527      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
2528      CALL histwrite(nid_mth,"dtvdf",itap,zx_tmp_3d,
2529     .                                   iim*jjmp1*klev,ndex3d)
2530c
2531      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
2532      CALL histwrite(nid_mth,"dqvdf",itap,zx_tmp_3d,
2533     .                                   iim*jjmp1*klev,ndex3d)
2534c
2535      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_eva, zx_tmp_3d)
2536      CALL histwrite(nid_mth,"dteva",itap,zx_tmp_3d,
2537     .                                   iim*jjmp1*klev,ndex3d)
2538c
2539      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_eva, zx_tmp_3d)
2540      CALL histwrite(nid_mth,"dqeva",itap,zx_tmp_3d,
2541     .                                   iim*jjmp1*klev,ndex3d)
2542c
2543      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zpt_conv, zx_tmp_3d)
2544      CALL histwrite(nid_mth,"ptconv",itap,zx_tmp_3d,
2545     .                                   iim*jjmp1*klev,ndex3d)
2546c
2547      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, ratqs, zx_tmp_3d)
2548      CALL histwrite(nid_mth,"ratqs",itap,zx_tmp_3d,
2549     .                                   iim*jjmp1*klev,ndex3d)
2550c
2551      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
2552      CALL histwrite(nid_mth,"dtajs",itap,zx_tmp_3d,
2553     .                                   iim*jjmp1*klev,ndex3d)
2554c
2555      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_ajs, zx_tmp_3d)
2556      CALL histwrite(nid_mth,"dqajs",itap,zx_tmp_3d,
2557     .                                   iim*jjmp1*klev,ndex3d)
2558c
2559      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat, zx_tmp_3d)
2560      CALL histwrite(nid_mth,"dtswr",itap,zx_tmp_3d,
2561     .                                   iim*jjmp1*klev,ndex3d)
2562c
2563      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat0, zx_tmp_3d)
2564      CALL histwrite(nid_mth,"dtsw0",itap,zx_tmp_3d,
2565     .                                   iim*jjmp1*klev,ndex3d)
2566c
2567      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool, zx_tmp_3d)
2568      CALL histwrite(nid_mth,"dtlwr",itap,zx_tmp_3d,
2569     .                                   iim*jjmp1*klev,ndex3d)
2570c
2571      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool0, zx_tmp_3d)
2572      CALL histwrite(nid_mth,"dtlw0",itap,zx_tmp_3d,
2573     .                                   iim*jjmp1*klev,ndex3d)
2574c
2575      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_vdf, zx_tmp_3d)
2576      CALL histwrite(nid_mth,"duvdf",itap,zx_tmp_3d,
2577     .                                   iim*jjmp1*klev,ndex3d)
2578c
2579      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_vdf, zx_tmp_3d)
2580      CALL histwrite(nid_mth,"dvvdf",itap,zx_tmp_3d,
2581     .                                   iim*jjmp1*klev,ndex3d)
2582c
2583      IF (ok_orodr) THEN
2584      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_oro, zx_tmp_3d)
2585      CALL histwrite(nid_mth,"duoro",itap,zx_tmp_3d,
2586     .                                   iim*jjmp1*klev,ndex3d)
2587c
2588      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_oro, zx_tmp_3d)
2589      CALL histwrite(nid_mth,"dvoro",itap,zx_tmp_3d,
2590     .                                   iim*jjmp1*klev,ndex3d)
2591c
2592      ENDIF
2593C
2594      IF (ok_orolf) THEN
2595      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_lif, zx_tmp_3d)
2596      CALL histwrite(nid_mth,"dulif",itap,zx_tmp_3d,
2597     .                                   iim*jjmp1*klev,ndex3d)
2598c
2599      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_lif, zx_tmp_3d)
2600      CALL histwrite(nid_mth,"dvlif",itap,zx_tmp_3d,
2601     .                                   iim*jjmp1*klev,ndex3d)
2602      ENDIF
2603C
2604      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, wo, zx_tmp_3d)
2605      CALL histwrite(nid_mth,"ozone",itap,zx_tmp_3d,
2606     .                                   iim*jjmp1*klev,ndex3d)
2607c
2608      IF (nqmax.GE.3) THEN
2609      DO iq=1,nqmax-2
2610      IF (iq.LE.99) THEN
2611         CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,iq+2), zx_tmp_3d)
2612         WRITE(str2,'(i2.2)') iq
2613         CALL histwrite(nid_mth,"trac"//str2,itap,zx_tmp_3d,
2614     .                                   iim*jjmp1*klev,ndex3d)
2615      ELSE
2616         PRINT*, "Trop de traceurs"
2617         CALL abort
2618      ENDIF
2619      ENDDO
2620      ENDIF
2621cKE43
2622      IF (iflag_con.EQ.4) THEN ! (sb)
2623c
2624      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, upwd, zx_tmp_3d)
2625      CALL histwrite(nid_mth,"upwd",itap,zx_tmp_3d,
2626     .                                   iim*jjmp1*klev,ndex3d)
2627c
2628      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd, zx_tmp_3d)
2629      CALL histwrite(nid_mth,"dnwd",itap,zx_tmp_3d,
2630     .                                   iim*jjmp1*klev,ndex3d)
2631c
2632      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd0, zx_tmp_3d)
2633      CALL histwrite(nid_mth,"dnwd0",itap,zx_tmp_3d,
2634     .                                   iim*jjmp1*klev,ndex3d)
2635c
2636      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, Ma, zx_tmp_3d)
2637      CALL histwrite(nid_mth,"Ma",itap,zx_tmp_3d,
2638     .                                   iim*jjmp1*klev,ndex3d)
2639c
2640c
2641      ENDIF
2642c34EK
2643c
2644      if (ok_sync) then
2645        call histsync(nid_mth)
2646      endif
2647      ENDIF
2648c
2649      IF (ok_instan) THEN
2650c
2651      ndex2d = 0
2652      ndex3d = 0
2653c
2654c Champs 2D:
2655c
2656         i = NINT(zout/zsto)
2657         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
2658         CALL histwrite(nid_ins,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2659c
2660         i = NINT(zout/zsto)
2661         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
2662         CALL histwrite(nid_ins,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
2663
2664      DO i = 1, klon
2665         zx_tmp_fi2d(i) = paprs(i,1)
2666      ENDDO
2667      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
2668      CALL histwrite(nid_ins,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2669c
2670      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
2671      CALL histwrite(nid_ins,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
2672c
2673c Champs 3D:
2674c
2675      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
2676      CALL histwrite(nid_ins,"temp",itap,zx_tmp_3d,
2677     .                                   iim*jjmp1*klev,ndex3d)
2678c
2679      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
2680      CALL histwrite(nid_ins,"vitu",itap,zx_tmp_3d,
2681     .                                   iim*jjmp1*klev,ndex3d)
2682c
2683      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
2684      CALL histwrite(nid_ins,"vitv",itap,zx_tmp_3d,
2685     .                                   iim*jjmp1*klev,ndex3d)
2686c
2687      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
2688      CALL histwrite(nid_ins,"geop",itap,zx_tmp_3d,
2689     .                                   iim*jjmp1*klev,ndex3d)
2690c
2691      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
2692      CALL histwrite(nid_ins,"pres",itap,zx_tmp_3d,
2693     .                                   iim*jjmp1*klev,ndex3d)
2694c
2695      if (ok_sync) then
2696        call histsync(nid_ins)
2697      endif
2698      ENDIF
2699c
2700      IF (ok_oasis .AND. mod(itap,nexca).EQ.0) THEN
2701c
2702c Je ne traite pas le ruissellement, pour l'instant (qui m'aidera ?)
2703         DO i = 1, klon
2704            oas_ruisoce(i) = 0.0
2705            oas_ruisriv(i) = 0.0
2706         ENDDO
2707c
2708         ig = 1
2709         DO i = 1, iim
2710            z_sols(i,1) = oas_sols(ig)
2711            z_nsol(i,1) = oas_nsol(ig)
2712            z_rain(i,1) = oas_rain(ig)
2713            z_snow(i,1) = oas_snow(ig)
2714            z_evap(i,1) = oas_evap(ig)
2715            z_ruisoce(i,1) = oas_ruisoce(ig)
2716            z_ruisriv(i,1) = oas_ruisriv(ig)
2717            z_tsol(i,1) = oas_tsol(ig)
2718            z_fder(i,1) = oas_fder(ig)
2719            z_albe(i,1) = oas_albe(ig)
2720            z_taux(i,1) = oas_taux(ig)
2721            z_tauy(i,1) = oas_tauy(ig)
2722         ENDDO
2723         DO j = 2, jjm
2724         DO i = 1, iim
2725            ig = ig + 1
2726            z_sols(i,j) = oas_sols(ig)
2727            z_nsol(i,j) = oas_nsol(ig)
2728            z_rain(i,j) = oas_rain(ig)
2729            z_snow(i,j) = oas_snow(ig)
2730            z_evap(i,j) = oas_evap(ig)
2731            z_ruisoce(i,j) = oas_ruisoce(ig)
2732            z_ruisriv(i,j) = oas_ruisriv(ig)
2733            z_tsol(i,j) = oas_tsol(ig)
2734            z_fder(i,j) = oas_fder(ig)
2735            z_albe(i,j) = oas_albe(ig)
2736            z_taux(i,j) = oas_taux(ig)
2737            z_tauy(i,j) = oas_tauy(ig)
2738         ENDDO
2739         ENDDO
2740         ig = ig + 1
2741         DO i = 1, iim
2742            z_sols(i,jjmp1)    = oas_sols(ig)
2743            z_nsol(i,jjmp1)    = oas_nsol(ig)
2744            z_rain(i,jjmp1)    = oas_rain(ig)
2745            z_snow(i,jjmp1)    = oas_snow(ig)
2746            z_evap(i,jjmp1)    = oas_evap(ig)
2747            z_ruisoce(i,jjmp1) = oas_ruisoce(ig)
2748            z_ruisriv(i,jjmp1) = oas_ruisriv(ig)
2749            z_tsol(i,jjmp1)    = oas_tsol(ig)
2750            z_fder(i,jjmp1)    = oas_fder(ig)
2751            z_albe(i,jjmp1)    = oas_albe(ig)
2752            z_taux(i,jjmp1)    = oas_taux(ig)
2753            z_tauy(i,jjmp1)    = oas_tauy(ig)
2754         ENDDO
2755c
2756c Passer les champs au coupleur:
2757c
2758         CALL intocpl(itap,jjmp1*iim,
2759     .                   z_sols, z_nsol,
2760     .                   z_rain, z_snow, z_evap,
2761     .                   z_ruisoce, z_ruisriv,
2762     .                   z_tsol, z_fder, z_albe,
2763     .                   z_taux, z_tauy)
2764         DO i = 1, klon
2765           oas_sols(i) = 0.0
2766           oas_nsol(i) = 0.0
2767           oas_rain(i) = 0.0
2768           oas_snow(i) = 0.0
2769           oas_evap(i) = 0.0
2770           oas_ruis(i) = 0.0
2771           oas_tsol(i) = 0.0
2772           oas_fder(i) = 0.0
2773           oas_albe(i) = 0.0
2774           oas_taux(i) = 0.0
2775           oas_tauy(i) = 0.0
2776         ENDDO
2777      ENDIF
2778c
2779c Ecrire la bande regionale (binaire grads)
2780      IF (ok_region .AND. mod(itap,ecrit_reg).eq.0) THEN
2781         CALL ecriregs(84,zxtsol)
2782         CALL ecriregs(84,paprs(1,1))
2783         CALL ecriregs(84,topsw)
2784         CALL ecriregs(84,toplw)
2785         CALL ecriregs(84,solsw)
2786         CALL ecriregs(84,sollw)
2787         CALL ecriregs(84,rain_fall)
2788         CALL ecriregs(84,snow_fall)
2789         CALL ecriregs(84,evap)
2790         CALL ecriregs(84,sens)
2791         CALL ecriregs(84,bils)
2792         CALL ecriregs(84,pctsrf(1,is_sic))
2793         CALL ecriregs(84,fluxu(1,1))
2794         CALL ecriregs(84,fluxv(1,1))
2795         CALL ecriregs(84,ue)
2796         CALL ecriregs(84,ve)
2797         CALL ecriregs(84,uq)
2798         CALL ecriregs(84,vq)
2799c
2800         CALL ecrirega(84,u_seri)
2801         CALL ecrirega(84,v_seri)
2802         CALL ecrirega(84,omega)
2803         CALL ecrirega(84,t_seri)
2804         CALL ecrirega(84,zphi)
2805         CALL ecrirega(84,q_seri)
2806         CALL ecrirega(84,cldfra)
2807         CALL ecrirega(84,cldliq)
2808         CALL ecrirega(84,pplay)
2809
2810
2811cc         CALL ecrirega(84,d_t_dyn)
2812cc         CALL ecrirega(84,d_q_dyn)
2813cc         CALL ecrirega(84,heat)
2814cc         CALL ecrirega(84,cool)
2815cc         CALL ecrirega(84,d_t_con)
2816cc         CALL ecrirega(84,d_q_con)
2817cc         CALL ecrirega(84,d_t_lsc)
2818cc         CALL ecrirega(84,d_q_lsc)
2819      ENDIF
2820c
2821c Convertir les incrementations en tendances
2822c
2823      DO k = 1, klev
2824      DO i = 1, klon
2825         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2826         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2827         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2828         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2829         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2830      ENDDO
2831      ENDDO
2832c
2833      IF (nqmax.GE.3) THEN
2834      DO iq = 3, nqmax
2835      DO  k = 1, klev
2836      DO  i = 1, klon
2837         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2838      ENDDO
2839      ENDDO
2840      ENDDO
2841      ENDIF
2842c
2843c Sauvegarder les valeurs de t et q a la fin de la physique:
2844c
2845      DO k = 1, klev
2846      DO i = 1, klon
2847         t_ancien(i,k) = t_seri(i,k)
2848         q_ancien(i,k) = q_seri(i,k)
2849      ENDDO
2850      ENDDO
2851c
2852c====================================================================
2853c Si c'est la fin, il faut conserver l'etat de redemarrage
2854c====================================================================
2855c
2856      IF (lafin) THEN
2857ccc         IF (ok_oasis) CALL quitcpl
2858         CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
2859     .        rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
2860     .        radsol,rugmer,agesno,
2861     .        zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2862     .        t_ancien, q_ancien)
2863      ENDIF
2864
2865      RETURN
2866      END
2867      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2868      IMPLICIT none
2869c
2870c Calculer et imprimer l'eau totale. A utiliser pour verifier
2871c la conservation de l'eau
2872c
2873#include "YOMCST.h"
2874      INTEGER klon,klev
2875      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2876      REAL aire(klon)
2877      REAL qtotal, zx, qcheck
2878      INTEGER i, k
2879c
2880      zx = 0.0
2881      DO i = 1, klon
2882         zx = zx + aire(i)
2883      ENDDO
2884      qtotal = 0.0
2885      DO k = 1, klev
2886      DO i = 1, klon
2887         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2888     .                     *(paprs(i,k)-paprs(i,k+1))/RG
2889      ENDDO
2890      ENDDO
2891c
2892      qcheck = qtotal/zx
2893c
2894      RETURN
2895      END
2896      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2897      IMPLICIT none
2898c
2899c Tranformer une variable de la grille physique a
2900c la grille d'ecriture
2901c
2902      INTEGER nfield,nlon,iim,jjmp1, jjm
2903      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2904c
2905      INTEGER i, j, n, ig
2906c
2907      jjm = jjmp1 - 1
2908      DO n = 1, nfield
2909         DO i=1,iim
2910            ecrit(i,n) = fi(1,n)
2911            ecrit(i+jjm*iim,n) = fi(nlon,n)
2912         ENDDO
2913         DO ig = 1, nlon - 2
2914           ecrit(iim+ig,n) = fi(1+ig,n)
2915         ENDDO
2916      ENDDO
2917      RETURN
2918      END
2919
Note: See TracBrowser for help on using the repository browser.