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

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

Interface avec les differentes surface, version de travail.LF

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