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

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

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