source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/physiq.F @ 4640

Last change on this file since 4640 was 710, checked in by Laurent Fairhead, 18 years ago

Integration des modifications d'OM et JLD pour une meilleure prise en compte
de la conservation du flux d'eau
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 89.0 KB
Line 
1!
2! $Header$
3!
4C
5      SUBROUTINE physiq (nlon,nlev,nqmax,
6     .            debut,lafin,rjourvrai,gmtime,pdtphys,
7     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
8     .            u,v,t,qx,
9     .            omega,
10#ifdef INCA_CH4
11     .            flxmass_w,
12#endif
13     .            d_u, d_v, d_t, d_qx, d_ps)
14
15      USE ioipsl
16      USE histcom
17#ifdef INCA
18      USE chemshut
19#ifdef INCA_CH4
20!      USE obs_pos
21#endif
22#endif
23      IMPLICIT none
24c======================================================================
25c
26c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
27c
28c Objet: Moniteur general de la physique du modele
29cAA      Modifications quant aux traceurs :
30cAA                  -  uniformisation des parametrisations ds phytrac
31cAA                  -  stockage des moyennes des champs necessaires
32cAA                     en mode traceur off-line
33c======================================================================
34c   CLEFS CPP POUR LES IO
35c   =====================
36c #define histhf
37#define histday
38#define histmth
39#define histins
40c #define histISCCP
41c #define histREGDYN
42c #define histmthNMC
43c======================================================================
44c    modif   ( P. Le Van ,  12/10/98 )
45c
46c  Arguments:
47c
48c nlon----input-I-nombre de points horizontaux
49c nlev----input-I-nombre de couches verticales
50c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
51c debut---input-L-variable logique indiquant le premier passage
52c lafin---input-L-variable logique indiquant le dernier passage
53c rjour---input-R-numero du jour de l'experience
54c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
55c pdtphys-input-R-pas d'integration pour la physique (seconde)
56c paprs---input-R-pression pour chaque inter-couche (en Pa)
57c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
58c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
59c pphis---input-R-geopotentiel du sol
60c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
61c u-------input-R-vitesse dans la direction X (de O a E) en m/s
62c v-------input-R-vitesse Y (de S a N) en m/s
63c t-------input-R-temperature (K)
64c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
65c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
66c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
67c omega---input-R-vitesse verticale en Pa/s
68c
69c d_u-----output-R-tendance physique de "u" (m/s/s)
70c d_v-----output-R-tendance physique de "v" (m/s/s)
71c d_t-----output-R-tendance physique de "t" (K/s)
72c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
73c d_ps----output-R-tendance physique de la pression au sol
74c======================================================================
75#include "dimensions.h"
76      integer jjmp1
77      parameter (jjmp1=jjm+1-1/jjm)
78#include "dimphy.h"
79#include "regdim.h"
80#include "indicesol.h"
81#include "dimsoil.h"
82#include "clesphys.h"
83#include "control.h"
84#include "temps.h"
85#include "comgeomphy.h"
86#include "advtrac.h"
87#include "iniprint.h"
88c======================================================================
89      LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
90      PARAMETER (ok_cvl=.TRUE.)
91      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
92      PARAMETER (ok_gust=.FALSE.)
93c======================================================================
94      LOGICAL check ! Verifier la conservation du modele en eau
95      PARAMETER (check=.FALSE.)
96      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
97      PARAMETER (ok_stratus=.FALSE.)
98c======================================================================
99c Parametres lies au coupleur OASIS:
100#include "oasis.h"
101      INTEGER,SAVE :: npas, nexca
102      logical rnpb
103#ifdef INCA
104      parameter(rnpb=.false.)
105#else
106      parameter(rnpb=.true.)
107#endif
108c      ocean = type de modele ocean a utiliser: force, slab, couple
109      character*6 ocean
110      SAVE ocean
111
112c      parameter (ocean = 'force ')
113c     parameter (ocean = 'couple')
114      logical ok_ocean
115c======================================================================
116c Clef controlant l'activation du cycle diurne:
117ccc      LOGICAL cycle_diurne
118ccc      PARAMETER (cycle_diurne=.FALSE.)
119c======================================================================
120c Modele thermique du sol, a activer pour le cycle diurne:
121ccc      LOGICAL soil_model
122ccc      PARAMETER (soil_model=.FALSE.)
123      logical ok_veget
124      save ok_veget
125c     parameter (ok_veget = .true.)
126c      parameter (ok_veget = .false.)
127c======================================================================
128c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
129c le calcul du rayonnement est celle apres la precipitation des nuages.
130c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
131c la condensation et la precipitation. Cette cle augmente les impacts
132c radiatifs des nuages.
133ccc      LOGICAL new_oliq
134ccc      PARAMETER (new_oliq=.FALSE.)
135c======================================================================
136c Clefs controlant deux parametrisations de l'orographie:
137cc      LOGICAL ok_orodr
138ccc      PARAMETER (ok_orodr=.FALSE.)
139ccc      LOGICAL ok_orolf
140ccc      PARAMETER (ok_orolf=.FALSE.)
141c======================================================================
142      LOGICAL ok_journe ! sortir le fichier journalier
143      save ok_journe
144c      PARAMETER (ok_journe=.true.)
145c
146      LOGICAL ok_mensuel ! sortir le fichier mensuel
147      save ok_mensuel
148c      PARAMETER (ok_mensuel=.true.)
149c
150      LOGICAL ok_instan ! sortir le fichier instantane
151      save ok_instan
152c      PARAMETER (ok_instan=.true.)
153c
154      LOGICAL ok_region ! sortir le fichier regional
155      PARAMETER (ok_region=.FALSE.)
156c======================================================================
157c
158      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
159      PARAMETER (ivap=1)
160      INTEGER iliq          ! indice de traceurs pour eau liquide
161      PARAMETER (iliq=2)
162
163c
164c
165c Variables argument:
166c
167      INTEGER nlon
168      INTEGER nlev
169      INTEGER nqmax
170      REAL rjourvrai
171      REAL gmtime
172      REAL pdtphys
173      LOGICAL debut, lafin
174      REAL paprs(klon,klev+1)
175      REAL pplay(klon,klev)
176      REAL pphi(klon,klev)
177      REAL pphis(klon)
178      REAL presnivs(klev)
179      REAL znivsig(klev)
180      REAL zsurf(nbsrf)
181
182      REAL u(klon,klev)
183      REAL v(klon,klev)
184      REAL t(klon,klev)
185      REAL qx(klon,klev,nqmax)
186
187      REAL t_ancien(klon,klev), q_ancien(klon,klev)
188      SAVE t_ancien, q_ancien
189      LOGICAL ancien_ok
190      SAVE ancien_ok
191
192      REAL d_t_dyn(klon,klev)
193      REAL d_q_dyn(klon,klev)
194
195      REAL omega(klon,klev)
196
197#ifdef INCA_CH4
198      REAL flxmass_w(klon,klev)
199#endif
200      REAL d_u(klon,klev)
201      REAL d_v(klon,klev)
202      REAL d_t(klon,klev)
203      REAL d_qx(klon,klev,nqmax)
204      REAL d_ps(klon)
205
206      INTEGER klevp1, klevm1
207      PARAMETER(klevp1=klev+1,klevm1=klev-1)
208#include "raddim.h"
209c
210cIM 080304   REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
211      REAL swdn0(klon,klevp1), swdn(klon,klevp1)
212      REAL swup0(klon,klevp1), swup(klon,klevp1)
213      SAVE swdn0 , swdn, swup0, swup
214c
215      REAL SWdn200clr(klon), SWdn200(klon)
216      REAL SWup200clr(klon), SWup200(klon)
217      SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
218c
219      REAL lwdn0(klon,klevp1), lwdn(klon,klevp1)
220      REAL lwup0(klon,klevp1), lwup(klon,klevp1)
221      SAVE lwdn0 , lwdn, lwup0, lwup
222c
223      REAL LWdn200clr(klon), LWdn200(klon)
224      REAL LWup200clr(klon), LWup200(klon)
225      SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
226c
227      REAL LWdnTOA(klon), LWdnTOAclr(klon)
228      SAVE LWdnTOA, LWdnTOAclr
229c
230c vents meridien et zonal a un niveau de pression
231c
232      integer nlevSTD
233      PARAMETER(nlevSTD=17)
234      real rlevSTD(nlevSTD)
235      DATA rlevSTD/100000., 92500., 85000., 70000.,
236     .60000., 50000., 40000., 30000., 25000., 20000.,
237     .15000., 10000., 7000., 5000., 3000., 2000., 1000./
238      CHARACTER*5 clevSTD(nlevSTD), aa, bb
239      DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
240     .'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
241     .'70  ','50  ','30  ','20  ','10  '/
242c
243      real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
244      real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
245      real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
246c
247cIM ENSEMBLES BEG
248c
249      integer nlevENS
250      PARAMETER(nlevENS=4)
251      integer indENS(nlevENS)
252      save indENS
253      real rlevENS(nlevENS)
254      DATA rlevENS/85000., 70000., 50000., 20000./
255      CHARACTER*3 clev(nlevENS)
256      DATA clev/'850','700','500','200'/
257 
258      real tlev(klon,nlevENS), qlev(klon,nlevENS), rhlev(klon,nlevENS)
259      real ulev(klon,nlevENS), vlev(klon,nlevENS), philev(klon,nlevENS)
260      real wlev(klon,nlevENS)
261cIM ENSEMBLES END
262c
263c prw: precipitable water
264      real prw(klon)
265
266      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
267      REAL convfra(klon,klev)  ! fraction nuageuse convective
268
269      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
270      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
271      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
272      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
273
274      INTEGER linv, kp1
275c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
276c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
277      REAL flwp(klon), fiwp(klon)
278      REAL flwc(klon,klev), fiwc(klon,klev)
279      REAL flwp_c(klon), fiwp_c(klon)
280      REAL flwc_c(klon,klev), fiwc_c(klon,klev)
281      REAL flwp_s(klon), fiwp_s(klon)
282      REAL flwc_s(klon,klev), fiwc_s(klon,klev)
283
284c ISCCP simulator v3.4
285c dans clesphys.h top_height, overlap
286cv3.4
287      INTEGER debug, debugcol
288      INTEGER npoints
289      PARAMETER(npoints=klon)
290c
291      INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
292      INTEGER nregISCtot
293      PARAMETER(nregISCtot=1)
294c
295c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
296c y compris pour 1 point
297c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
298c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
299      INTEGER imin_debut, nbpti
300      INTEGER jmin_debut, nbptj
301c
302      REAL nbsunlit(nregISCtot,klon)  !nbsunlit : moyenne de sunlit
303      INTEGER ncol, seed(klon)
304
305c ncol = nb. de sous-colonnes pour chaque maille du GCM
306c     PARAMETER(ncol=100)
307c     PARAMETER(ncol=625)
308c     PARAMETER(ncol=10)
309      PARAMETER(ncol=25)
310      REAL tautab(0:255)
311      INTEGER invtau(-20:45000)
312      REAL emsfc_lw
313      PARAMETER(emsfc_lw=0.99)
314      REAL    ran0                      ! type for random number fuction
315c
316      REAL cldtot(klon,klev)
317c variables de haut en bas pour le simulateur ISCCP
318      REAL dtau_s(klon,klev) !tau nuages startiformes
319      REAL dtau_c(klon,klev) !tau nuages convectifs
320      REAL dem_s(klon,klev)  !emissivite nuages startiformes
321      REAL dem_c(klon,klev)  !emissivite nuages convectifs
322c
323c variables de haut en bas pour le simulateur ISCCP
324      REAL pfull(klon,klev)
325      REAL phalf(klon,klev+1)
326      REAL qv(klon,klev)
327      REAL cc(klon,klev)
328      REAL conv(klon,klev)
329      REAL dtau_sH2B(klon,klev)
330      REAL dtau_cH2B(klon,klev)
331      REAL at(klon,klev)
332      REAL dem_sH2B(klon,klev)
333      REAL dem_cH2B(klon,klev)
334
335c output from ISCCP simulator
336      REAL fq_isccp(klon,7,7)
337      REAL totalcldarea(klon)
338      REAL meanptop(klon)
339      REAL meantaucld(klon)
340      REAL boxtau(klon,ncol)
341      REAL boxptop(klon,ncol)
342c
343      INTEGER l, ni, nj, kmax, lmax
344      PARAMETER(kmax=8, lmax=8)
345      INTEGER kmaxm1, lmaxm1
346      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
347      INTEGER iimx7, jjmx7, jjmp1x7
348      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1,
349     .jjmp1x7=jjmp1*lmaxm1)
350      REAL fq4d(iim,jjmp1,kmaxm1,lmaxm1)
351      REAL fq3d(iimx7, jjmp1x7)
352c
353      INTEGER iw, iwmax
354      REAL wmin, pas_w
355c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
356      PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
357      REAL o500(klon)
358c
359cIM: nbregdyn = nbre regions pour calculs statistiques sur output du ISCCP
360cIM: dynamiques 
361      INTEGER nreg, nbregdyn
362      PARAMETER(nbregdyn=5)
363      REAL histoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
364      REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
365      REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbregdyn)
366      SAVE nhistoWt
367
368      INTEGER linv
369      INTEGER pct_ocean(klon,nbregdyn)
370      REAL rlonPOS(klon)
371 
372c sorties ISCCP
373
374      logical ok_isccp
375      real ecrit_isccp
376      integer nid_isccp
377      save ok_isccp, ecrit_isccp, nid_isccp       
378
379#ifdef histISCCP
380      data ok_isccp/.true./     
381#else
382      data ok_isccp/.false./
383#endif
384
385c sorties statistiques regime dynamique
386      logical ok_regdyn
387      real ecrit_regdyn
388      integer nid_regdyn
389      save ok_regdyn, ecrit_regdyn, nid_regdyn
390
391#ifdef histREGDYN
392c     data ok_regdyn,ecrit_regdyn/.true.,0.125/
393c     data ok_regdyn,ecrit_regdyn/.true.,1./
394       data ok_regdyn/.true./
395#else
396      data ok_regdyn/.false./
397#endif
398
399      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
400      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
401      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
402
403c cldtopres pression au sommet des nuages
404      REAL cldtopres(lmaxm1)
405      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
406
407      INTEGER komega, nhoriRD
408
409c taulev: numero du niveau de tau dans les sorties ISCCP
410      CHARACTER *4 taulev(kmaxm1)
411      DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
412
413      REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
414      INTEGER nhorix7
415cIM: region='3d' <==> sorties en global
416      CHARACTER*3 region
417      PARAMETER(region='3d')
418c
419      logical ok_hf
420      real ecrit_hf
421      integer nid_hf, nid_hf3d
422      save ok_hf, ecrit_hf, nid_hf, nid_hf3d
423
424c  QUESTION : noms de variables ?
425
426#ifdef histhf
427      data ok_hf,ecrit_hf/.true.,0.25/
428#else
429      data ok_hf/.false./
430#endif
431
432      INTEGER        longcles
433      PARAMETER    ( longcles = 20 )
434      REAL clesphy0( longcles      )
435c
436c Variables quasi-arguments
437c
438      REAL xjour
439      SAVE xjour
440c
441c
442c Variables propres a la physique
443c
444      REAL dtime
445      SAVE dtime                  ! pas temporel de la physique
446c
447      INTEGER radpas
448      SAVE radpas                 ! frequence d'appel rayonnement
449c
450      REAL radsol(klon)
451      SAVE radsol               ! bilan radiatif au sol calcule par code radiatif
452c
453      REAL rlat(klon)
454      SAVE rlat                   ! latitude pour chaque point
455c
456      REAL rlon(klon)
457      SAVE rlon                   ! longitude pour chaque point
458c
459cc      INTEGER iflag_con
460cc      SAVE iflag_con              ! indicateur de la convection
461c
462      INTEGER itap
463      SAVE itap                   ! compteur pour la physique
464c
465      REAL co2_ppm_etat0
466c
467      REAL solaire_etat0
468c
469      real slp(klon) ! sea level pressure
470
471      REAL ftsol(klon,nbsrf)
472      SAVE ftsol                  ! temperature du sol
473c
474      REAL ftsoil(klon,nsoilmx,nbsrf)
475      SAVE ftsoil                 ! temperature dans le sol
476c
477      REAL fevap(klon,nbsrf)
478      SAVE fevap                 ! evaporation
479      REAL fluxlat(klon,nbsrf)
480      SAVE fluxlat
481c
482      REAL deltat(klon)
483      SAVE deltat                 ! ecart avec la SST de reference
484c
485      REAL fqsurf(klon,nbsrf)
486      SAVE fqsurf                 ! humidite de l'air au contact de la surface
487c
488      REAL qsol(klon)
489      SAVE qsol                  ! hauteur d'eau dans le sol
490c
491      REAL fsnow(klon,nbsrf)
492      SAVE fsnow                  ! epaisseur neigeuse
493c
494      REAL falbe(klon,nbsrf)
495      SAVE falbe                  ! albedo par type de surface
496      REAL falblw(klon,nbsrf)
497      SAVE falblw                 ! albedo par type de surface
498
499c
500c
501c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
502c
503      REAL zmea(klon)
504      SAVE zmea                   ! orographie moyenne
505c
506      REAL zstd(klon)
507      SAVE zstd                   ! deviation standard de l'OESM
508c
509      REAL zsig(klon)
510      SAVE zsig                   ! pente de l'OESM
511c
512      REAL zgam(klon)
513      save zgam                   ! anisotropie de l'OESM
514c
515      REAL zthe(klon)
516      SAVE zthe                   ! orientation de l'OESM
517c
518      REAL zpic(klon)
519      SAVE zpic                   ! Maximum de l'OESM
520c
521      REAL zval(klon)
522      SAVE zval                   ! Minimum de l'OESM
523c
524      REAL rugoro(klon)
525      SAVE rugoro                 ! longueur de rugosite de l'OESM
526c
527      REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
528c
529      REAL zuthe(klon),zvthe(klon)
530      SAVE zuthe
531      SAVE zvthe
532      INTEGER igwd,idx(klon),itest(klon)
533c
534      REAL agesno(klon,nbsrf)
535      SAVE agesno                 ! age de la neige
536c
537      REAL alb_neig(klon)
538      SAVE alb_neig               ! albedo de la neige
539c
540      REAL run_off_lic_0(klon)
541      SAVE run_off_lic_0
542cKE43
543c Variables liees a la convection de K. Emanuel (sb):
544c
545      REAL ema_workcbmf(klon)   ! cloud base mass flux
546      SAVE ema_workcbmf
547
548      REAL ema_cbmf(klon)       ! cloud base mass flux
549      SAVE ema_cbmf
550
551      REAL ema_pcb(klon)        ! cloud base pressure
552      SAVE ema_pcb
553
554      REAL ema_pct(klon)        ! cloud top pressure
555      SAVE ema_pct
556
557      REAL bas, top             ! cloud base and top levels
558      SAVE bas
559      SAVE top
560
561      REAL Ma(klon,klev)        ! undilute upward mass flux
562      SAVE Ma
563      REAL qcondc(klon,klev)    ! in-cld water content from convect
564      SAVE qcondc
565      REAL ema_work1(klon, klev), ema_work2(klon, klev)
566      SAVE ema_work1, ema_work2
567      REAL wdn(klon), tdn(klon), qdn(klon)
568
569      REAL wd(klon) ! sb
570      SAVE wd       ! sb
571
572c Variables locales pour la couche limite (al1):
573c
574cAl1      REAL pblh(klon)           ! Hauteur de couche limite
575cAl1      SAVE pblh
576c34EK
577c
578c Variables locales:
579c
580      REAL cdragh(klon) ! drag coefficient pour T and Q
581      REAL cdragm(klon) ! drag coefficient pour vent
582cAA
583cAA  Pour phytrac
584cAA
585      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
586      REAL yu1(klon)            ! vents dans la premiere couche U
587      REAL yv1(klon)            ! vents dans la premiere couche V
588      REAL ffonte(klon,nbsrf)    !Flux thermique utilise pour fondre la neige
589      REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface
590c                               !et necessaire pour limiter la
591c                               !hauteur de neige, en kg/m2/s
592      REAL fqfonte(klon,nbsrf)  !Quantite d'eau de fonte des glaciers
593      REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
594
595      LOGICAL offline           ! Controle du stockage ds "physique"
596      PARAMETER (offline=.false.)
597      INTEGER physid
598      REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
599      save pfrac_impa
600      REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
601      save pfrac_nucl
602      REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
603      save pfrac_1nucl
604      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
605      REAL frac_nucl(klon,klev) ! idem (nucleation)
606#ifdef INCA
607      REAL          :: calday
608#endif
609
610cAA
611      REAL rain_fall(klon) ! pluie
612      REAL snow_fall(klon) ! neige
613      save snow_fall, rain_fall
614cIM 050204 BEG
615      REAL total_rain(klon), nday_rain(klon)
616      save total_rain, nday_rain
617cIM 050204 END
618      REAL evap(klon), devap(klon) ! evaporation et sa derivee
619      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
620      REAL dlw(klon)    ! derivee infra rouge
621      REAL bils(klon) ! bilan de chaleur au sol
622      REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
623C                             ! type de sous-surface et pondere par la fraction
624      REAL fder(klon) ! Derive de flux (sensible et latente)
625      save fder
626      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
627      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
628      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
629      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
630c
631      REAL frugs(klon,nbsrf) ! longueur de rugosite
632      save frugs
633      REAL zxrugs(klon) ! longueur de rugosite
634c
635c Conditions aux limites
636c
637      INTEGER julien
638c
639      INTEGER lmt_pas
640      SAVE lmt_pas                ! frequence de mise a jour
641      REAL pctsrf(klon,nbsrf)
642cIM
643      REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
644      REAL paire_ter(klon)        !surfaces terre
645cIM
646      SAVE pctsrf                 ! sous-fraction du sol
647      REAL albsol(klon)
648      SAVE albsol                 ! albedo du sol total
649      REAL albsollw(klon)
650      SAVE albsollw                 ! albedo du sol total
651
652      REAL wo(klon,klev)
653      SAVE wo                     ! ozone
654c======================================================================
655c
656c Declaration des procedures appelees
657c
658      EXTERNAL angle     ! calculer angle zenithal du soleil
659      EXTERNAL alboc     ! calculer l'albedo sur ocean
660      EXTERNAL ajsec     ! ajustement sec
661      EXTERNAL clmain    ! couche limite
662      EXTERNAL conlmd    ! convection (schema LMD)
663cKE43
664      EXTERNAL conema3  ! convect4.3
665      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
666cAA
667      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
668c                          ! stockage des coefficients necessaires au
669c                          ! lessivage OFF-LINE et ON-LINE
670      EXTERNAL hgardfou  ! verifier les temperatures
671      EXTERNAL nuage     ! calculer les proprietes radiatives
672      EXTERNAL o3cm      ! initialiser l'ozone
673      EXTERNAL orbite    ! calculer l'orbite terrestre
674      EXTERNAL ozonecm   ! prescrire l'ozone
675      EXTERNAL phyetat0  ! lire l'etat initial de la physique
676      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
677      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
678      EXTERNAL suphec    ! initialiser certaines constantes
679      EXTERNAL transp    ! transport total de l'eau et de l'energie
680      EXTERNAL ecribina  ! ecrire le fichier binaire global
681      EXTERNAL ecribins  ! ecrire le fichier binaire global
682      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
683      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
684cIM
685      EXTERNAL haut2bas  !variables de haut en bas
686      INTEGER lnblnk1
687      EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
688                         !caracter
689c
690c Variables locales
691c
692      real clwcon(klon,klev),rnebcon(klon,klev)
693      real clwcon0(klon,klev),rnebcon0(klon,klev)
694      save rnebcon, clwcon
695
696      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
697      REAL dialiq(klon,klev)  ! eau liquide nuageuse
698      REAL diafra(klon,klev)  ! fraction nuageuse
699      REAL cldliq(klon,klev)  ! eau liquide nuageuse
700      REAL cldfra(klon,klev)  ! fraction nuageuse
701      REAL cldtau(klon,klev)  ! epaisseur optique
702      REAL cldemi(klon,klev)  ! emissivite infrarouge
703c
704CXXX PB
705      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
706      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
707      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
708      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
709c
710      REAL zxfluxt(klon, klev)
711      REAL zxfluxq(klon, klev)
712      REAL zxfluxu(klon, klev)
713      REAL zxfluxv(klon, klev)
714CXXX
715      REAL heat(klon,klev)    ! chauffage solaire
716      REAL heat0(klon,klev)   ! chauffage solaire ciel clair
717      REAL cool(klon,klev)    ! refroidissement infrarouge
718      REAL cool0(klon,klev)   ! refroidissement infrarouge ciel clair
719      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
720      real sollwdown(klon)    ! downward LW flux at surface
721cIM BEG
722      real sollwdownclr(klon)    ! downward CS LW flux at surface
723      real toplwdown(klon)       ! downward CS LW flux at TOA
724      real toplwdownclr(klon)    ! downward CS LW flux at TOA
725cIM END
726      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
727      REAL albpla(klon)
728      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface
729      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface
730c Le rayonnement n'est pas calcule tous les pas, il faut donc
731c                      sauvegarder les sorties du rayonnement
732      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
733      SAVE  sollwdownclr, toplwdown, toplwdownclr
734      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
735c
736      INTEGER itaprad
737      SAVE itaprad
738c
739      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
740      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
741c
742      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
743      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
744c
745      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
746c
747      REAL dist, rmu0(klon), fract(klon)
748      REAL zdtime, zlongi
749c
750      CHARACTER*2 str2
751      CHARACTER*2 iqn
752c
753      REAL qcheck
754      REAL z_avant(klon), z_apres(klon), z_factor(klon)
755      LOGICAL zx_ajustq
756c
757      REAL za, zb
758      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
759      real zqsat(klon,klev)
760      INTEGER i, k, iq, ig, j, nsrf, ll
761      REAL t_coup
762      PARAMETER (t_coup=234.0)
763c
764      REAL zphi(klon,klev)
765      REAL zx_tmp_x(iim), zx_tmp_yjjmp1
766      REAL zx_relief(iim,jjmp1)
767      REAL zx_aire(iim,jjmp1)
768cKE43
769c Variables locales pour la convection de K. Emanuel (sb):
770c
771      REAL upwd(klon,klev)      ! saturated updraft mass flux
772      REAL dnwd(klon,klev)      ! saturated downdraft mass flux
773      REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
774      REAL tvp(klon,klev)       ! virtual temp of lifted parcel
775      REAL cape(klon)           ! CAPE
776      SAVE cape
777      CHARACTER*40 capemaxcels  !max(CAPE)
778
779      REAL pbase(klon)          ! cloud base pressure
780      SAVE pbase
781      REAL bbase(klon)          ! cloud base buoyancy
782      SAVE bbase
783      REAL rflag(klon)          ! flag fonctionnement de convect
784      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
785c -- convect43:
786      INTEGER ntra              ! nb traceurs pour convect4.3
787      REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
788      REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
789      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
790      REAL dplcldt(klon), dplcldr(klon)
791c?     .     condm_con(klon,klev),conda_con(klon,klev),
792c?     .     mr_con(klon,klev),ep_con(klon,klev)
793c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
794c --
795c34EK
796c
797c Variables du changement
798c
799c con: convection
800c lsc: condensation a grande echelle (Large-Scale-Condensation)
801c ajs: ajustement sec
802c eva: evaporation de l'eau liquide nuageuse
803c vdf: couche limite (Vertical DiFfusion)
804      REAL d_t_con(klon,klev),d_q_con(klon,klev)
805      REAL d_u_con(klon,klev),d_v_con(klon,klev)
806      REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
807      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
808      REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
809      REAL rneb(klon,klev)
810c
811      REAL pmfu(klon,klev), pmfd(klon,klev)
812      REAL pen_u(klon,klev), pen_d(klon,klev)
813      REAL pde_u(klon,klev), pde_d(klon,klev)
814      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
815      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
816      REAL prfl(klon,klev+1), psfl(klon,klev+1)
817c
818      INTEGER ibas_con(klon), itop_con(klon)
819      REAL rain_con(klon), rain_lsc(klon)
820      REAL snow_con(klon), snow_lsc(klon)
821      REAL d_ts(klon,nbsrf)
822c
823      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
824      REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
825c
826      REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
827      REAL d_t_oro(klon,klev)
828      REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
829      REAL d_t_lif(klon,klev)
830      REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif
831
832      REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
833      real ratqsbas,ratqshaut
834      save ratqsbas,ratqshaut, ratqs
835      real zpt_conv(klon,klev)
836
837c Parametres lies au nouveau schema de nuages (SB, PDF)
838      real fact_cldcon
839      real facttemps
840      logical ok_newmicro
841      save ok_newmicro
842      save fact_cldcon,facttemps
843      real facteur
844
845      integer iflag_cldcon
846      save iflag_cldcon
847
848      logical ptconv(klon,klev)
849
850c
851c Variables liees a l'ecriture de la bande histoire physique
852c
853      INTEGER ecrit_mth
854      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
855c
856      INTEGER ecrit_day
857      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
858c
859      INTEGER ecrit_ins
860      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
861c
862      INTEGER ecrit_reg
863      SAVE ecrit_reg   ! frequence d'ecriture
864c
865      integer itau_w   ! pas de temps ecriture = itap + itau_phy
866c
867c
868c Variables locales pour effectuer les appels en serie
869c
870      REAL t_seri(klon,klev), q_seri(klon,klev)
871      REAL ql_seri(klon,klev),qs_seri(klon,klev)
872      REAL u_seri(klon,klev), v_seri(klon,klev)
873c
874      REAL tr_seri(klon,klev,nbtr)
875      REAL d_tr(klon,klev,nbtr)
876
877      REAL zx_rh(klon,klev)
878
879      INTEGER        length
880      PARAMETER    ( length = 100 )
881      REAL tabcntr0( length       )
882c
883      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
884      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
885      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
886      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
887      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
888c
889      INTEGER nid_day, nid_mth, nid_ins, nid_nmc
890      SAVE nid_day, nid_mth, nid_ins, nid_nmc
891c
892      INTEGER nhori, nvert
893      REAL zsto, zout, zsto1, zsto2
894      real zjulian
895      save zjulian
896
897      character*20 modname
898      character*80 abort_message
899      logical ok_sync
900      real date0
901      integer idayref
902
903C essai writephys
904      integer fid_day, fid_mth, fid_ins
905      parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
906      integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
907      parameter (prof2d_on = 1, prof3d_on = 2,
908     .           prof2d_av = 3, prof3d_av = 4)
909      character*30 nom_fichier
910      character*10 varname
911      character*40 vartitle
912      character*20 varunits
913C     Variables liees au bilan d'energie et d'enthalpi
914      REAL ztsol(klon)
915      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
916     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
917      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
918     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
919      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
920      REAL      d_h_vcol_phy
921      REAL      fs_bound, fq_bound
922      SAVE      d_h_vcol_phy
923      REAL      zero_v(klon)
924      CHARACTER*15 ztit
925      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
926      SAVE      ip_ebil
927      DATA      ip_ebil/0/
928      INTEGER   if_ebil ! level for energy conserv. dignostics
929      SAVE      if_ebil
930c+jld ec_conser
931      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
932      REAL ZRCPD
933c-jld ec_conser
934cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
935      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)   !temperature, humidite a 2m
936      REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
937      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
938      REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
939      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
940cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
941      REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
942      REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
943      SAVE sulfate_pi
944
945      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
946
947      REAL re(klon, klev)       ! Cloud droplet effective radius
948      REAL fl(klon, klev)  ! denominator of re
949
950      REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
951
952      ! Aerosol optical properties
953      REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
954      REAL cg_ae(klon,klev,2)
955
956      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
957      ! ok_ade=T -ADE=topswad-topsw
958
959      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
960      ! ok_aie=T ->
961      !        ok_ade=T -AIE=topswai-topswad
962      !        ok_ade=F -AIE=topswai-topsw
963
964      REAL aerindex(klon)       ! POLDER aerosol index
965     
966      ! Parameters
967      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
968      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
969cjq-end
970c
971c Declaration des constantes et des fonctions thermodynamiques
972c
973         SAVE fqfonte
974#include "YOMCST.h"
975#include "YOETHF.h"
976#include "FCTTRE.h"
977c======================================================================
978      modname = 'physiq'
979      IF (if_ebil.ge.1) THEN
980        DO i=1,klon
981          zero_v(i)=0.
982        END DO
983      END IF
984      ok_sync=.TRUE.
985      IF (nqmax .LT. 2) THEN
986         abort_message = 'eaux vapeur et liquide sont indispensables'
987         CALL abort_gcm (modname,abort_message,1)
988      ENDIF
989      IF (debut) THEN
990         CALL suphec ! initialiser constantes et parametres phys.
991c
992cIM 050204 BEG
993         DO i=1, klon
994          nday_rain(i)=0.
995         ENDDO
996cIM 050204 END
997c
998c======================================================================
999cIM BEG
1000        DO k=1, nlev
1001          DO l=1, nlevSTD
1002c
1003            bb=clevSTD(l)
1004c
1005            IF(l.GE.2) THEN
1006             aa=clevSTD(l)
1007             bb=aa(1:lnblnk1(aa))
1008            ENDIF
1009c
1010            IF(bb.EQ.clev(k)) THEN
1011c             print*,'k=',k,'l=',l,'clev=',clev(k)
1012              indENS(k)=l
1013c             print*,'k=',k,'l=',l,'clev=',clev(k),'indENS=',indENS(k)
1014            ENDIF
1015c
1016          ENDDO
1017        ENDDO
1018c
1019      ENDIF !debut
1020cIM END
1021      xjour = rjourvrai
1022c
1023c Si c'est le debut, il faut initialiser plusieurs choses
1024c          ********
1025c
1026       IF (debut) THEN
1027C
1028         fqfonte(:,:)=0.
1029         IF (if_ebil.ge.1) d_h_vcol_phy=0.
1030c
1031c appel a la lecture du run.def physique
1032c
1033         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1034     .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
1035     .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
1036     .                  ok_ade, ok_aie,
1037     .                  bl95_b0, bl95_b1)
1038cIM  .                  , RI0)
1039
1040c
1041c
1042c Initialiser les compteurs:
1043c
1044
1045         frugs = 0.
1046         itap    = 0
1047         itaprad = 0
1048         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1049     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
1050     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
1051     .       dlw,radsol,frugs,agesno,clesphy0,
1052     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1053     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1054     .       run_off_lic_0)
1055
1056c
1057         radpas = NINT( 86400./dtime/nbapp_rad)
1058c
1059C on remet le calendrier a zero
1060c
1061         IF (raz_date .eq. 1) THEN
1062           itau_phy = 0
1063         ENDIF
1064
1065c
1066         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
1067     ,                    ok_instan, ok_region )
1068c
1069         IF (ABS(dtime-pdtphys).GT.0.001) THEN
1070            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1071     .                        pdtphys
1072            abort_message='Pas physique n est pas correct '
1073            call abort_gcm(modname,abort_message,1)
1074         ENDIF
1075         IF (nlon .NE. klon) THEN
1076            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1077     .                      klon
1078            abort_message='nlon et klon ne sont pas coherents'
1079            call abort_gcm(modname,abort_message,1)
1080         ENDIF
1081         IF (nlev .NE. klev) THEN
1082            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1083     .                       klev
1084            abort_message='nlev et klev ne sont pas coherents'
1085            call abort_gcm(modname,abort_message,1)
1086         ENDIF
1087c
1088         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1089           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1090           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1091           abort_message='Nbre d appels au rayonnement insuffisant'
1092           call abort_gcm(modname,abort_message,1)
1093         ENDIF
1094         WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1095         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1096     .                   ok_cvl
1097c
1098cKE43
1099c Initialisation pour la convection de K.E. (sb):
1100         IF (iflag_con.GE.3) THEN
1101
1102         WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1103         WRITE(lunout,*)
1104     .      "On va utiliser le melange convectif des traceurs qui"
1105         WRITE(lunout,*)"est calcule dans convect4.3"
1106         WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1107
1108          DO i = 1, klon
1109           ema_cbmf(i) = 0.
1110           ema_pcb(i)  = 0.
1111           ema_pct(i)  = 0.
1112           ema_workcbmf(i) = 0.
1113          ENDDO
1114cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1115          DO i = 1, klon
1116           ibas_con(i) = 1
1117           itop_con(i) = klev+1
1118          ENDDO
1119cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1120
1121         ENDIF
1122
1123c34EK
1124         IF (ok_orodr) THEN
1125         DO i=1,klon
1126         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1127         ENDDO
1128         CALL SUGWD(klon,klev,paprs,pplay)
1129         DO i=1,klon
1130         zuthe(i)=0.
1131         zvthe(i)=0.
1132         if(zstd(i).gt.10.)then
1133           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1134           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1135         endif
1136         ENDDO
1137         ENDIF
1138c
1139c
1140         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1141         WRITE(lunout,*)'La frequence de lecture surface est de ',
1142     .                   lmt_pas
1143c
1144         ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
1145         IF (ok_mensuel) THEN
1146         WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1147     .                   ecrit_mth
1148         ENDIF
1149         ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
1150         IF (ok_journe) THEN
1151         WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1152     .                   ecrit_day
1153         ENDIF
1154ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
1155ccc         ecrit_ins = NINT(86400./dtime *0.25)  ! 4 fois par jour
1156         ecrit_ins = NINT(86400./dtime/48.)  ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1157         ecrit_ins = NINT(86400./dtime/12.)  ! toutes les deux heures
1158         IF (ok_instan) THEN
1159         WRITE(lunout,*)'La frequence de sortie instant. est de ',
1160     .                   ecrit_ins
1161         ENDIF
1162         ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
1163         IF (ok_region) THEN
1164         WRITE(lunout,*)'La frequence de sortie region est de ',
1165     .                   ecrit_reg
1166         ENDIF
1167
1168c
1169c Initialiser le couplage si necessaire
1170c
1171      npas = 0
1172      nexca = 0
1173      if (ocean == 'couple') then
1174        npas = itaufin/ iphysiq
1175        nexca = 86400 / dtime
1176        write(lunout,*)' ##### Ocean couple #####'
1177        write(lunout,*)' Valeurs des pas de temps'
1178        write(lunout,*)' npas = ', npas
1179        write(lunout,*)' nexca = ', nexca
1180      endif       
1181c
1182c
1183cIM
1184      capemaxcels = 't_max(X)'
1185      t2mincels = 't_min(X)'
1186      t2maxcels = 't_max(X)'
1187
1188c
1189c=============================================================
1190c   Initialisation des sorties
1191c=============================================================
1192
1193#ifdef CPP_IOIPSL
1194
1195#ifdef histhf
1196#include "ini_histhf.h"
1197#endif
1198
1199#ifdef histday
1200#include "ini_histday.h"
1201#endif
1202
1203#ifdef histmth
1204#include "ini_histmth.h"
1205#endif
1206
1207#ifdef histins
1208#include "ini_histins.h"
1209#endif
1210
1211#ifdef histISCCP
1212#include "ini_histISCCP.h"
1213#endif
1214
1215#ifdef histmthNMC
1216#include "ini_histmthNMC.h"
1217#endif
1218
1219#ifdef histREGDYN
1220#include "ini_histREGDYN.h"
1221#endif
1222
1223#ifdef histISCCP
1224#include "ini_histISCCP.h"
1225#endif
1226#endif
1227
1228cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1229      date0 = zjulian
1230C      date0 = day_ini
1231      WRITE(*,*) 'physiq date0 : ',date0
1232c
1233c
1234c
1235c Prescrire l'ozone dans l'atmosphere
1236c
1237c
1238cc         DO i = 1, klon
1239cc         DO k = 1, klev
1240cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1241cc         ENDDO
1242cc         ENDDO
1243c
1244#ifdef INCA
1245           iii = MOD(NINT(xjour),360)
1246           calday = FLOAT(iii) + gmtime
1247           WRITE(lunout,*) 'initial time ', xjour, calday
1248#ifdef INCAINFO
1249           WRITE(lunout,*) 'Appel CHEMINI ...'
1250#endif
1251           CALL chemini( rpi,
1252     $                   rg,
1253     $                   ra,
1254     $                   airephy,
1255     $                   rlat,
1256     $                   rlon,
1257     $                   presnivs,
1258     $                   calday,
1259     $                   tracnam,
1260     $                   natsnam,
1261c     $                   mxoutflds,
1262c     $                   outinst,
1263c     $                   outtimav,
1264     $                   klon,
1265     $                   nqmax,
1266     $                   pdtphys,
1267     $                   anne_ref,
1268     $                   day_ini)
1269#ifdef INCAINFO
1270           WRITE(lunout,*) 'OK.'
1271#endif
1272#endif
1273c
1274      ENDIF
1275c
1276c   ****************     Fin  de   IF ( debut  )   ***************
1277c
1278c
1279c Mettre a zero des variables de sortie (pour securite)
1280c
1281      DO i = 1, klon
1282         d_ps(i) = 0.0
1283      ENDDO
1284      DO k = 1, klev
1285      DO i = 1, klon
1286         d_t(i,k) = 0.0
1287         d_u(i,k) = 0.0
1288         d_v(i,k) = 0.0
1289      ENDDO
1290      ENDDO
1291      DO iq = 1, nqmax
1292      DO k = 1, klev
1293      DO i = 1, klon
1294         d_qx(i,k,iq) = 0.0
1295      ENDDO
1296      ENDDO
1297      ENDDO
1298c
1299c Ne pas affecter les valeurs entrees de u, v, h, et q
1300c
1301      DO k = 1, klev
1302      DO i = 1, klon
1303         t_seri(i,k)  = t(i,k)
1304         u_seri(i,k)  = u(i,k)
1305         v_seri(i,k)  = v(i,k)
1306         q_seri(i,k)  = qx(i,k,ivap)
1307         ql_seri(i,k) = qx(i,k,iliq)
1308         qs_seri(i,k) = 0.
1309      ENDDO
1310      ENDDO
1311      IF (nqmax.GE.3) THEN
1312      DO iq = 3, nqmax
1313      DO  k = 1, klev
1314      DO  i = 1, klon
1315         tr_seri(i,k,iq-2) = qx(i,k,iq)
1316      ENDDO
1317      ENDDO
1318      ENDDO
1319      ELSE
1320      DO k = 1, klev
1321      DO i = 1, klon
1322         tr_seri(i,k,1) = 0.0
1323      ENDDO
1324      ENDDO
1325      ENDIF
1326C
1327      DO i = 1, klon
1328        ztsol(i) = 0.
1329      ENDDO
1330      DO nsrf = 1, nbsrf
1331        DO i = 1, klon
1332          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1333        ENDDO
1334      ENDDO
1335C
1336      IF (if_ebil.ge.1) THEN
1337        ztit='after dynamic'
1338        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1339     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1340     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1341C     Comme les tendances de la physique sont ajoute dans la dynamique,
1342C     on devrait avoir que la variation d'entalpie par la dynamique
1343C     est egale a la variation de la physique au pas de temps precedent.
1344C     Donc la somme de ces 2 variations devrait etre nulle.
1345        call diagphy(airephy,ztit,ip_ebil
1346     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1347     e      , zero_v, zero_v, zero_v, ztsol
1348     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1349     s      , fs_bound, fq_bound )
1350      END IF
1351
1352c Diagnostiquer la tendance dynamique
1353c
1354      IF (ancien_ok) THEN
1355         DO k = 1, klev
1356         DO i = 1, klon
1357            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1358            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1359         ENDDO
1360         ENDDO
1361      ELSE
1362         DO k = 1, klev
1363         DO i = 1, klon
1364            d_t_dyn(i,k) = 0.0
1365            d_q_dyn(i,k) = 0.0
1366         ENDDO
1367         ENDDO
1368         ancien_ok = .TRUE.
1369      ENDIF
1370c
1371c Ajouter le geopotentiel du sol:
1372c
1373      DO k = 1, klev
1374      DO i = 1, klon
1375         zphi(i,k) = pphi(i,k) + pphis(i)
1376      ENDDO
1377      ENDDO
1378c
1379c Verifier les temperatures
1380c
1381      CALL hgardfou(t_seri,ftsol,'debutphy')
1382c
1383c Incrementer le compteur de la physique
1384c
1385      itap   = itap + 1
1386      julien = MOD(NINT(xjour),360)
1387      if (julien .eq. 0) julien = 360
1388c
1389c Mettre en action les conditions aux limites (albedo, sst, etc.).
1390c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1391c
1392      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1393         WRITE(lunout,*)' PHYS cond  julien ',julien
1394         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1395      ENDIF
1396c
1397c Re-evaporer l'eau liquide nuageuse
1398c
1399      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1400      DO i = 1, klon
1401         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1402c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1403         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1404         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1405         zb = MAX(0.0,ql_seri(i,k))
1406         za = - MAX(0.0,ql_seri(i,k))
1407     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1408         t_seri(i,k) = t_seri(i,k) + za
1409         q_seri(i,k) = q_seri(i,k) + zb
1410         ql_seri(i,k) = 0.0
1411         d_t_eva(i,k) = za
1412         d_q_eva(i,k) = zb
1413      ENDDO
1414      ENDDO
1415c
1416      IF (if_ebil.ge.2) THEN
1417        ztit='after reevap'
1418        CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
1419     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1420     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1421         call diagphy(airephy,ztit,ip_ebil
1422     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1423     e      , zero_v, zero_v, zero_v, ztsol
1424     e      , d_h_vcol, d_qt, d_ec
1425     s      , fs_bound, fq_bound )
1426C
1427      END IF
1428C
1429c
1430c Appeler la diffusion verticale (programme de couche limite)
1431c
1432      DO i = 1, klon
1433c       if (.not. ok_veget) then
1434c          frugs(i,is_ter) = SQRT(frugs(i,is_ter)**2+rugoro(i)**2)
1435c       endif
1436c         frugs(i,is_lic) = rugoro(i)
1437c         frugs(i,is_oce) = rugmer(i)
1438c         frugs(i,is_sic) = 0.001
1439         zxrugs(i) = 0.0
1440      ENDDO
1441      DO nsrf = 1, nbsrf
1442      DO i = 1, klon
1443c         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1444        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1445      ENDDO
1446      ENDDO
1447      DO nsrf = 1, nbsrf
1448      DO i = 1, klon
1449            zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1450      ENDDO
1451      ENDDO
1452c
1453C calculs necessaires au calcul de l'albedo dans l'interface
1454c
1455      CALL orbite(FLOAT(julien),zlongi,dist)
1456      IF (cycle_diurne) THEN
1457        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1458        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1459      ELSE
1460        rmu0 = -999.999
1461      ENDIF
1462cIM BEG
1463      DO i=1, klon
1464       sunlit(i)=1
1465       IF(rmu0(i).EQ.0.) sunlit(i)=0
1466       nbsunlit(1,i)=FLOAT(sunlit(i))
1467      ENDDO
1468cIM END
1469C     Calcul de l'abedo moyen par maille
1470      albsol(:)=0.
1471      albsollw(:)=0.
1472      DO nsrf = 1, nbsrf
1473      DO i = 1, klon
1474         albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1475         albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1476      ENDDO
1477      ENDDO
1478C
1479C     Repartition sous maille des flux LW et SW
1480C Modif OM+PASB+JLD
1481C Repartition du longwave par sous-surface linearisee
1482Cn
1483
1484       DO nsrf = 1, nbsrf
1485       DO i = 1, klon
1486c$$$        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
1487c$$$        fsollw(i,nsrf) = sollw(i)
1488         fsollw(i,nsrf) = sollw(i)
1489     $      + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf))
1490         fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1491       ENDDO
1492       ENDDO
1493
1494      fder = dlw
1495
1496
1497      CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
1498     e            t_seri,q_seri,u_seri,v_seri,
1499     e            julien, rmu0, co2_ppm,
1500     e            ok_veget, ocean, npas, nexca, ftsol,
1501     $            soil_model,cdmmax, cdhmax,
1502     $            ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
1503     $            paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1504     $            fluxlat,
1505cIM cf. JLD  e            rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1506     e            rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1507     e            rlon, rlat, cuphy, cvphy, frugs,
1508     e            debut, lafin, agesno,rugoro ,
1509     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1510     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1511     s            dsens, devap,
1512     s            ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1513     s            fqcalving, fqfonte,ffonte, run_off_lic_0)
1514c
1515CXXX PB
1516CXXX Incrementation des flux
1517CXXX
1518
1519      zxfluxt=0.
1520      zxfluxq=0.
1521      zxfluxu=0.
1522      zxfluxv=0.
1523      DO nsrf = 1, nbsrf
1524        DO k = 1, klev
1525          DO i = 1, klon
1526            zxfluxt(i,k) = zxfluxt(i,k) +
1527     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1528            zxfluxq(i,k) = zxfluxq(i,k) +
1529     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1530            zxfluxu(i,k) = zxfluxu(i,k) +
1531     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1532            zxfluxv(i,k) = zxfluxv(i,k) +
1533     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1534          END DO
1535        END DO
1536      END DO
1537      DO i = 1, klon
1538         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1539c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1540         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1541         fder(i) = dlw(i) + dsens(i) + devap(i)
1542      ENDDO
1543
1544
1545      DO k = 1, klev
1546      DO i = 1, klon
1547         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1548         q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1549         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1550         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1551      ENDDO
1552      ENDDO
1553c
1554      IF (if_ebil.ge.2) THEN
1555        ztit='after clmain'
1556        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1557     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1558     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1559         call diagphy(airephy,ztit,ip_ebil
1560     e      , zero_v, zero_v, zero_v, zero_v, sens
1561     e      , evap  , zero_v, zero_v, ztsol
1562     e      , d_h_vcol, d_qt, d_ec
1563     s      , fs_bound, fq_bound )
1564      END IF
1565C
1566c
1567c Incrementer la temperature du sol
1568c
1569      DO i = 1, klon
1570         zxtsol(i) = 0.0
1571         zxfluxlat(i) = 0.0
1572c
1573         zt2m(i) = 0.0
1574         zq2m(i) = 0.0
1575         zu10m(i) = 0.0
1576         zv10m(i) = 0.0
1577cIM cf JLD ??
1578         zxffonte(i) = 0.0
1579         zxfqcalving(i) = 0.0
1580         zxfqfonte(i) = 0.0
1581c
1582         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1583     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA)
1584     $       THEN
1585             WRITE(*,*) 'physiq : pb sous surface au point ', i,
1586     $           pctsrf(i, 1 : nbsrf)
1587         ENDIF
1588      ENDDO
1589      DO nsrf = 1, nbsrf
1590        DO i = 1, klon
1591c        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1592            ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1593cIM cf. JLD
1594            wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1595     $         + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1596            zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1597            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1598cccIM
1599            zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1600            zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1601            zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1602            zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1603cIM cf JLD ??
1604            zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1605            zxfqcalving(i) = zxfqcalving(i) +
1606     .                      fqcalving(i,nsrf)*pctsrf(i,nsrf)
1607            zxfqfonte(i) = zxfqfonte(i) +
1608     .                      fqfonte(i,nsrf)*pctsrf(i,nsrf)
1609c        ENDIF
1610        ENDDO
1611      ENDDO
1612
1613c
1614c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1615c
1616      DO nsrf = 1, nbsrf
1617        DO i = 1, klon
1618          IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1619cccIM
1620          IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1621          IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1622          IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1623          IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1624cIM cf JLD ??
1625          IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1626          IF (pctsrf(i,nsrf) .LT. epsfra)
1627     .    fqcalving(i,nsrf) = zxfqcalving(i)
1628          IF (pctsrf(i,nsrf) .LT. epsfra)
1629     .    fqfonte(i,nsrf) = zxfqfonte(i)
1630        ENDDO
1631      ENDDO
1632c
1633c
1634c Calculer la derive du flux infrarouge
1635c
1636cXXX      DO nsrf = 1, nbsrf
1637      DO i = 1, klon
1638cXXX        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1639            dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3
1640cXXX     .          *(ftsol(i,nsrf)-zxtsol(i))
1641cXXX     .          *pctsrf(i,nsrf)
1642cXXX        ENDIF
1643cXXX      ENDDO
1644      ENDDO
1645c
1646c Appeler la convection (au choix)
1647c
1648      DO k = 1, klev
1649      DO i = 1, klon
1650         conv_q(i,k) = d_q_dyn(i,k)
1651     .               + d_q_vdf(i,k)/dtime
1652         conv_t(i,k) = d_t_dyn(i,k)
1653     .               + d_t_vdf(i,k)/dtime
1654      ENDDO
1655      ENDDO
1656      IF (check) THEN
1657         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1658         WRITE(lunout,*) "avantcon=", za
1659      ENDIF
1660      zx_ajustq = .FALSE.
1661      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1662      IF (zx_ajustq) THEN
1663         DO i = 1, klon
1664            z_avant(i) = 0.0
1665         ENDDO
1666         DO k = 1, klev
1667         DO i = 1, klon
1668            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1669     .                        *(paprs(i,k)-paprs(i,k+1))/RG
1670         ENDDO
1671         ENDDO
1672      ENDIF
1673      IF (iflag_con.EQ.1) THEN
1674          stop'reactiver le call conlmd dans physiq.F'
1675c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1676c    .             d_t_con, d_q_con,
1677c    .             rain_con, snow_con, ibas_con, itop_con)
1678      ELSE IF (iflag_con.EQ.2) THEN
1679      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1680     e            conv_t, conv_q, zxfluxq(1,1), omega,
1681     s            d_t_con, d_q_con, rain_con, snow_con,
1682     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1683     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
1684      WHERE (rain_con < 0.) rain_con = 0.
1685      WHERE (snow_con < 0.) snow_con = 0.
1686      DO i = 1, klon
1687         ibas_con(i) = klev+1 - kcbot(i)
1688         itop_con(i) = klev+1 - kctop(i)
1689      ENDDO
1690      ELSE IF (iflag_con.GE.3) THEN
1691c nb of tracers for the KE convection:
1692          if (nqmax .GE. 4) then
1693              ntra = nbtr
1694          else
1695              ntra = 1
1696          endif
1697c
1698c sb, oct02:
1699c Schema de convection modularise et vectorise:
1700c (driver commun aux versions 3 et 4)
1701c
1702          IF (ok_cvl) THEN ! new driver for convectL
1703
1704          CALL concvl (iflag_con,
1705     .        dtime,paprs,pplay,t_seri,q_seri,
1706     .        u_seri,v_seri,tr_seri,nbtr,
1707     .        ema_work1,ema_work2,
1708     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1709     .        rain_con, snow_con, ibas_con, itop_con,
1710     .        upwd,dnwd,dnwd0,
1711     .        Ma,cape,tvp,iflagctrl,
1712     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
1713cIM cf. FH
1714              clwcon0=qcondc
1715
1716          ELSE ! ok_cvl
1717
1718          CALL conema3 (dtime,
1719     .        paprs,pplay,t_seri,q_seri,
1720     .        u_seri,v_seri,tr_seri,nbtr,
1721     .        ema_work1,ema_work2,
1722     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1723     .        rain_con, snow_con, ibas_con, itop_con,
1724     .        upwd,dnwd,dnwd0,bas,top,
1725     .        Ma,cape,tvp,rflag,
1726     .        pbase
1727     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1728     .        ,clwcon0)
1729
1730          ENDIF ! ok_cvl
1731
1732           IF (.NOT. ok_gust) THEN
1733           do i = 1, klon
1734            wd(i)=0.0
1735           enddo
1736           ENDIF
1737
1738c =================================================================== c
1739c Calcul des proprietes des nuages convectifs
1740c
1741      DO k = 1, klev
1742      DO i = 1, klon
1743         zx_t = t_seri(i,k)
1744         IF (thermcep) THEN
1745            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1746            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1747            zx_qs  = MIN(0.5,zx_qs)
1748            zcor   = 1./(1.-retv*zx_qs)
1749            zx_qs  = zx_qs*zcor
1750         ELSE
1751           IF (zx_t.LT.t_coup) THEN
1752              zx_qs = qsats(zx_t)/pplay(i,k)
1753           ELSE
1754              zx_qs = qsatl(zx_t)/pplay(i,k)
1755           ENDIF
1756         ENDIF
1757         zqsat(i,k)=zx_qs
1758      ENDDO
1759      ENDDO
1760
1761c   calcul des proprietes des nuages convectifs
1762             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1763             call clouds_gno
1764     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1765
1766c =================================================================== c
1767
1768          DO i = 1, klon
1769            ema_pcb(i)  = pbase(i)
1770          ENDDO
1771          DO i = 1, klon
1772            ema_pct(i)  = paprs(i,itop_con(i))
1773          ENDDO
1774          DO i = 1, klon
1775            ema_cbmf(i) = ema_workcbmf(i)
1776          ENDDO     
1777      ELSE
1778          WRITE(lunout,*) "iflag_con non-prevu", iflag_con
1779          CALL abort
1780      ENDIF
1781
1782c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1783c    .              d_u_con, d_v_con)
1784
1785      DO k = 1, klev
1786        DO i = 1, klon
1787         t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1788         q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1789         u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1790         v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1791        ENDDO
1792      ENDDO
1793c
1794      IF (if_ebil.ge.2) THEN
1795        ztit='after convect'
1796        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1797     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1798     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1799         call diagphy(airephy,ztit,ip_ebil
1800     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1801     e      , zero_v, rain_con, snow_con, ztsol
1802     e      , d_h_vcol, d_qt, d_ec
1803     s      , fs_bound, fq_bound )
1804      END IF
1805C
1806      IF (check) THEN
1807          za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1808          WRITE(lunout,*)"aprescon=", za
1809          zx_t = 0.0
1810          za = 0.0
1811          DO i = 1, klon
1812            za = za + airephy(i)/FLOAT(klon)
1813            zx_t = zx_t + (rain_con(i)+
1814     .                   snow_con(i))*airephy(i)/FLOAT(klon)
1815          ENDDO
1816          zx_t = zx_t/za*dtime
1817          WRITE(lunout,*)"Precip=", zx_t
1818      ENDIF
1819      IF (zx_ajustq) THEN
1820          DO i = 1, klon
1821            z_apres(i) = 0.0
1822          ENDDO
1823          DO k = 1, klev
1824            DO i = 1, klon
1825              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1826     .            *(paprs(i,k)-paprs(i,k+1))/RG
1827            ENDDO
1828          ENDDO
1829          DO i = 1, klon
1830            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1831     .          /z_apres(i)
1832          ENDDO
1833          DO k = 1, klev
1834            DO i = 1, klon
1835              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1836     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
1837                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
1838              ENDIF
1839            ENDDO
1840          ENDDO
1841      ENDIF
1842      zx_ajustq=.FALSE.
1843c
1844      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1845c
1846          IF (iflag_con .NE. 2 .AND. debut) THEN
1847              WRITE(lunout,*)'Pour l instant, seul conflx fonctionne ',
1848     $            'avec traceurs', iflag_con
1849              WRITE(lunout,*)' Mettre iflag_con',
1850     $            ' = 2 dans run.def et repasser'
1851c              CALL abort
1852              ENDIF
1853c
1854      ENDIF !--nqmax.GT.2
1855c
1856c Appeler l'ajustement sec
1857c
1858      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1859      DO k = 1, klev
1860      DO i = 1, klon
1861         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1862         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1863      ENDDO
1864      ENDDO
1865c
1866      IF (if_ebil.ge.2) THEN
1867        ztit='after dry_adjust'
1868        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1869     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1870     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1871      END IF
1872
1873
1874c-------------------------------------------------------------------------
1875c  Caclul des ratqs
1876c-------------------------------------------------------------------------
1877
1878c      print*,'calcul des ratqs'
1879c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1880c   ----------------
1881c   on ecrase le tableau ratqsc calcule par clouds_gno
1882      if (iflag_cldcon.eq.1) then
1883         do k=1,klev
1884         do i=1,klon
1885            if(ptconv(i,k)) then
1886              ratqsc(i,k)=ratqsbas
1887     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
1888            else
1889               ratqsc(i,k)=0.
1890            endif
1891         enddo
1892         enddo
1893      endif
1894
1895c   ratqs stables
1896c   -------------
1897      do k=1,klev
1898cIM RAJOUT boucle do=i
1899       do i=1, klon
1900cIM         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
1901cIM     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
1902         ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
1903     s   min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
1904cIM      print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
1905       enddo
1906      enddo
1907
1908
1909c  ratqs final
1910c  -----------
1911      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
1912c   les ratqs sont une conbinaison de ratqss et ratqsc
1913c   ratqs final
1914c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1915c   relaxation des ratqs
1916c         facttemps=exp(-pdtphys/1.e4)
1917         facteur=exp(-pdtphys*facttemps)
1918         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1919         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1920c         print*,'calcul des ratqs fini'
1921      else
1922c   on ne prend que le ratqs stable pour fisrtilp
1923         ratqs(:,:)=ratqss(:,:)
1924      endif
1925
1926
1927c
1928c Appeler le processus de condensation a grande echelle
1929c et le processus de precipitation
1930c-------------------------------------------------------------------------
1931      CALL fisrtilp(dtime,paprs,pplay,
1932     .           t_seri, q_seri,ptconv,ratqs,
1933     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1934     .           rain_lsc, snow_lsc,
1935     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
1936     .           frac_impa, frac_nucl,
1937     .           prfl, psfl, rhcl)
1938
1939      WHERE (rain_lsc < 0) rain_lsc = 0.
1940      WHERE (snow_lsc < 0) snow_lsc = 0.
1941      DO k = 1, klev
1942      DO i = 1, klon
1943         t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1944         q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1945         ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1946         cldfra(i,k) = rneb(i,k)
1947         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1948      ENDDO
1949      ENDDO
1950      IF (check) THEN
1951         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1952         WRITE(lunout,*)"apresilp=", za
1953         zx_t = 0.0
1954         za = 0.0
1955         DO i = 1, klon
1956            za = za + airephy(i)/FLOAT(klon)
1957            zx_t = zx_t + (rain_lsc(i)
1958     .                  + snow_lsc(i))*airephy(i)/FLOAT(klon)
1959        ENDDO
1960         zx_t = zx_t/za*dtime
1961         WRITE(lunout,*)"Precip=", zx_t
1962      ENDIF
1963c
1964      IF (if_ebil.ge.2) THEN
1965        ztit='after fisrt'
1966        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1967     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1968     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1969        call diagphy(airephy,ztit,ip_ebil
1970     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1971     e      , zero_v, rain_lsc, snow_lsc, ztsol
1972     e      , d_h_vcol, d_qt, d_ec
1973     s      , fs_bound, fq_bound )
1974      END IF
1975c
1976c-------------------------------------------------------------------
1977c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1978c-------------------------------------------------------------------
1979
1980c 1. NUAGES CONVECTIFS
1981c
1982      IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
1983
1984c Nuages diagnostiques pour Tiedtke
1985      CALL diagcld1(paprs,pplay,
1986     .             rain_con,snow_con,ibas_con,itop_con,
1987     .             diafra,dialiq)
1988      DO k = 1, klev
1989      DO i = 1, klon
1990      IF (diafra(i,k).GT.cldfra(i,k)) THEN
1991         cldliq(i,k) = dialiq(i,k)
1992         cldfra(i,k) = diafra(i,k)
1993      ENDIF
1994      ENDDO
1995      ENDDO
1996
1997      ELSE IF (iflag_cldcon.eq.3) THEN
1998c  On prend pour les nuages convectifs le max du calcul de la
1999c  convection et du calcul du pas de temps précédent diminué d'un facteur
2000c  facttemps
2001c      facttemps=pdtphys/1.e4
2002      facteur = pdtphys *facttemps
2003      do k=1,klev
2004         do i=1,klon
2005            rnebcon(i,k)=rnebcon(i,k)*facteur
2006            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
2007     s      then
2008                rnebcon(i,k)=rnebcon0(i,k)
2009                clwcon(i,k)=clwcon0(i,k)
2010            endif
2011         enddo
2012      enddo
2013
2014cIM calcul nuages par le simulateur ISCCP
2015      IF (ok_isccp) THEN
2016cIM calcul tau. emi nuages convectifs
2017      convfra(:,:)=rnebcon(:,:)
2018      convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
2019      CALL newmicro (paprs, pplay,ok_newmicro,
2020     .            t_seri, convliq, convfra, dtau_c, dem_c,
2021     .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
2022     .            flwp_c, fiwp_c, flwc_c, fiwc_c,
2023     e            ok_aie,
2024     e            sulfate, sulfate_pi,
2025     e            bl95_b0, bl95_b1,
2026     s            cldtaupi, re, fl)
2027c
2028cIM calcul tau. emi nuages startiformes
2029      CALL newmicro (paprs, pplay,ok_newmicro,
2030     .            t_seri, cldliq, cldfra, dtau_s, dem_s,
2031     .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
2032     .            flwp_s, fiwp_s, flwc_s, fiwc_s,
2033     e            ok_aie,
2034     e            sulfate, sulfate_pi,
2035     e            bl95_b0, bl95_b1,
2036     s            cldtaupi, re, fl)
2037c
2038      cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2039
2040cIM inversion des niveaux de pression ==> de haut en bas
2041      CALL haut2bas(klon, klev, pplay, pfull)
2042      CALL haut2bas(klon, klev, q_seri, qv)
2043      CALL haut2bas(klon, klev, cldtot, cc)
2044      CALL haut2bas(klon, klev, rnebcon, conv)
2045      CALL haut2bas(klon, klev, dtau_s, dtau_sH2B)
2046      CALL haut2bas(klon, klev, dtau_c, dtau_cH2B)
2047      CALL haut2bas(klon, klev, t_seri, at)
2048      CALL haut2bas(klon, klev, dem_s, dem_sH2B)
2049      CALL haut2bas(klon, klev, dem_c, dem_cH2B)
2050      CALL haut2bas(klon, klevp1, paprs, phalf)
2051
2052c     open(99,file='tautab.bin',access='sequential',
2053c    $     form='unformatted',status='old')
2054c     read(99) tautab
2055
2056cIM210503
2057      IF (debut) THEN
2058      open(99,file='tautab.formatted', FORM='FORMATTED')
2059      read(99,'(f30.20)') tautab
2060      close(99)
2061c
2062      open(99,file='invtau.formatted',form='FORMATTED')
2063      read(99,'(i10)') invtau
2064      close(99)
2065c
2066cIM: calcul coordonnees regions pour statistiques distribution
2067cIM: nuages en ftion du regime dynamique pour regions oceaniques
2068       IF (ok_regdyn) THEN !histREGDYN
2069       nsrf=3
2070       DO nreg=1, nbregdyn
2071       DO i=1, klon
2072
2073c       IF (debut) THEN
2074         IF(rlon(i).LT.0.) THEN
2075           rlonPOS(i)=rlon(i)+360.
2076         ELSE
2077           rlonPOS(i)=rlon(i) 
2078         ENDIF
2079c       ENDIF
2080
2081        pct_ocean(i,nreg)=0
2082
2083c test si c'est 1 point d'ocean
2084        IF(pctsrf(i,nsrf).EQ.1.) THEN
2085
2086         IF(nreg.EQ.1) THEN
2087
2088c TROP
2089          IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
2090           pct_ocean(i,nreg)=1
2091          ENDIF
2092
2093c PACIFIQUE NORD
2094          ELSEIF(nreg.EQ.2) THEN
2095           IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
2096            IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
2097             pct_ocean(i,nreg)=1
2098            ENDIF
2099           ENDIF
2100c CALIFORNIE ST-CU
2101         ELSEIF(nreg.EQ.3) THEN
2102          IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
2103           IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
2104            pct_ocean(i,nreg)=1
2105           ENDIF
2106          ENDIF
2107c HAWAI
2108        ELSEIF(nreg.EQ.4) THEN
2109         IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
2110          IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
2111           pct_ocean(i,nreg)=1
2112          ENDIF
2113         ENDIF
2114c WARM POOL
2115        ELSEIF(nreg.EQ.5) THEN
2116         IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
2117          IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
2118           pct_ocean(i,nreg)=1
2119          ENDIF
2120         ENDIF
2121        ENDIF !nbregdyn
2122c TROP
2123c        IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
2124c         pct_ocean(i)=.TRUE.
2125c         WRITE(*,*) 'pct_ocean =',i, rlon(i), rlat(i)
2126c          ENDIF !lon
2127c         ENDIF !lat
2128
2129        ENDIF !pctsrf
2130       ENDDO !klon
2131       ENDDO !nbregdyn
2132       ENDIF !ok_regdyn
2133 
2134cIM somme de toutes les nhistoW BEG
2135      DO nreg = 1, nbregdyn
2136       DO k = 1, kmaxm1
2137        DO l = 1, lmaxm1
2138         DO iw = 1, iwmax
2139          nhistoWt(k,l,iw,nreg)=0.
2140         ENDDO !iw
2141        ENDDO !l
2142       ENDDO !k
2143      ENDDO !nreg
2144cIM somme de toutes les nhistoW END
2145      ENDIF
2146cIM: initialisation de seed
2147        DO i=1, klon
2148          seed(i)=i+100
2149        ENDDO
2150     
2151cIM: pas de debug, debugcol
2152      debug=0
2153      debugcol=0
2154cIM260503
2155c o500 ==> distribution nuage ftion du regime dynamique a 500 hPa
2156        DO k=1, klevm1
2157        kp1=k+1
2158c       PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1)
2159        if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN
2160         DO i=1, klon
2161          o500(i)=omega(i,k)*RDAY/100.
2162c         if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1)
2163         ENDDO
2164         GOTO 1000
2165        endif
21661000  continue
2167      ENDDO
2168
2169      CALL ISCCP_CLOUD_TYPES(
2170     &     debug,
2171     &     debugcol,
2172     &     klon,
2173     &     sunlit,
2174     &     klev,
2175     &     ncol,
2176     &     seed,
2177     &     pfull,
2178     &     phalf,
2179     &     qv, cc, conv, dtau_sH2B, dtau_cH2B,
2180     &     top_height,
2181     &     overlap,
2182     &     tautab,
2183     &     invtau,
2184     &     ztsol,
2185     &     emsfc_lw,
2186     &     at, dem_sH2B, dem_cH2B,
2187     &     fq_isccp,
2188     &     totalcldarea,
2189     &     meanptop,
2190     &     meantaucld,
2191     &     boxtau,
2192     &     boxptop)
2193
2194
2195c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
2196      DO l=1, lmaxm1
2197       DO k=1, kmaxm1
2198        DO i=1, iim
2199         fq4d(i,1,k,l)=fq_isccp(1,k,l)
2200        ENDDO
2201        DO j=2, jjm
2202         DO i=1, iim
2203          ig=i+1+(j-2)*iim
2204          fq4d(i,j,k,l)=fq_isccp(ig,k,l)             
2205         ENDDO
2206        ENDDO
2207        DO i=1, iim
2208         fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
2209        ENDDO
2210       ENDDO
2211      ENDDO
2212c
2213      DO l=1, lmaxm1
2214       DO k=1, kmaxm1 
2215        DO j=1, jjmp1
2216         DO i=1, iim
2217           ni=(i-1)*lmaxm1+l
2218           nj=(j-1)*kmaxm1+k
2219           fq3d(ni,nj)=fq4d(i,j,k,l)
2220         ENDDO
2221        ENDDO
2222       ENDDO
2223      ENDDO
2224
2225c
2226c calculs statistiques distribution nuage ftion du regime dynamique
2227c
2228c Ce calcul doit etre fait a partir de valeurs mensuelles ??
2229      CALL histo_o500_pctau(nbregdyn,pct_ocean,o500,fq_isccp,
2230     &histoW,nhistoW)
2231c
2232c nhistoWt = somme de toutes les nhistoW
2233      DO nreg=1, nbregdyn
2234       DO k = 1, kmaxm1
2235        DO l = 1, lmaxm1
2236         DO iw = 1, iwmax
2237          nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
2238     &    nhistoW(k,l,iw,nreg)
2239         ENDDO
2240        ENDDO
2241       ENDDO
2242      ENDDO
2243c
2244      ENDIF !ok_isccp
2245
2246c   On prend la somme des fractions nuageuses et des contenus en eau
2247      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2248      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2249
2250      ENDIF
2251
2252c
2253c 2. NUAGES STARTIFORMES
2254c
2255      IF (ok_stratus) THEN
2256      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2257      DO k = 1, klev
2258      DO i = 1, klon
2259      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2260         cldliq(i,k) = dialiq(i,k)
2261         cldfra(i,k) = diafra(i,k)
2262      ENDIF
2263      ENDDO
2264      ENDDO
2265      ENDIF
2266c
2267c Precipitation totale
2268c
2269      DO i = 1, klon
2270         rain_fall(i) = rain_con(i) + rain_lsc(i)
2271         snow_fall(i) = snow_con(i) + snow_lsc(i)
2272      ENDDO
2273c
2274      IF (if_ebil.ge.2) THEN
2275        ztit="after diagcld"
2276        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2277     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2278     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2279      END IF
2280c
2281c Calculer l'humidite relative pour diagnostique
2282c
2283      DO k = 1, klev
2284      DO i = 1, klon
2285         zx_t = t_seri(i,k)
2286         IF (thermcep) THEN
2287            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2288            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2289            zx_qs  = MIN(0.5,zx_qs)
2290            zcor   = 1./(1.-retv*zx_qs)
2291            zx_qs  = zx_qs*zcor
2292         ELSE
2293           IF (zx_t.LT.t_coup) THEN
2294              zx_qs = qsats(zx_t)/pplay(i,k)
2295           ELSE
2296              zx_qs = qsatl(zx_t)/pplay(i,k)
2297           ENDIF
2298         ENDIF
2299         zx_rh(i,k) = q_seri(i,k)/zx_qs
2300         zqsat(i,k)=zx_qs
2301      ENDDO
2302      ENDDO
2303cjq - introduce the aerosol direct and first indirect radiative forcings
2304cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2305      IF (ok_ade.OR.ok_aie) THEN
2306         ! Get sulfate aerosol distribution
2307         CALL readsulfate(rjourvrai, debut, sulfate)
2308         CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2309
2310         ! Calculate aerosol optical properties (Olivier Boucher)
2311         CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2312     .        tau_ae, piz_ae, cg_ae, aerindex)
2313      ENDIF
2314
2315#ifdef INCA
2316           calday = FLOAT(julien) + gmtime
2317
2318#ifdef INCA_AER
2319      call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2320     &   prfl,psfl,pctsrf,airephy,xjour,rlat,rlon)
2321#endif
2322
2323#ifdef INCAINFO
2324           WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2325#endif
2326           CALL chemhook_begin (calday,
2327     $                          pctsrf(1,3),
2328     $                          rlat,
2329     $                          rlon,
2330     $                          airephy,
2331     $                          paprs,
2332     $                          pplay,
2333     $                          ycoefh,
2334     $                          pphi,
2335     $                          t_seri,
2336     $                          u,
2337     $                          v,
2338     $                          wo,
2339     $                          q_seri,
2340     $                          zxtsol,
2341     $                          zxsnow,
2342     $                          solsw,
2343     $                          albsol,
2344     $                          rain_fall,
2345     $                          snow_fall,
2346     $                          itop_con,
2347     $                          ibas_con,
2348     $                          cldfra,
2349     $                          iim,
2350     $                          jjm,
2351     $                          tr_seri)
2352#ifdef INCAINFO
2353           WRITE(lunout,*)'OK.'
2354#endif
2355#endif
2356c     
2357c Calculer les parametres optiques des nuages et quelques
2358c parametres pour diagnostiques:
2359c
2360      if (ok_newmicro) then
2361      CALL newmicro (paprs, pplay,ok_newmicro,
2362     .            t_seri, cldliq, cldfra, cldtau, cldemi,
2363     .            cldh, cldl, cldm, cldt, cldq,
2364     .            flwp, fiwp, flwc, fiwc,
2365     e            ok_aie,
2366     e            sulfate, sulfate_pi,
2367     e            bl95_b0, bl95_b1,
2368     s            cldtaupi, re, fl)
2369      else
2370      CALL nuage (paprs, pplay,
2371     .            t_seri, cldliq, cldfra, cldtau, cldemi,
2372     .            cldh, cldl, cldm, cldt, cldq,
2373     e            ok_aie,
2374     e            sulfate, sulfate_pi,
2375     e            bl95_b0, bl95_b1,
2376     s            cldtaupi, re, fl)
2377     
2378      endif
2379c
2380c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2381c
2382      IF (MOD(itaprad,radpas).EQ.0) THEN
2383      DO i = 1, klon
2384         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2385     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
2386     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
2387     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
2388         albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2389     .               + falblw(i,is_lic) * pctsrf(i,is_lic)
2390     .               + falblw(i,is_ter) * pctsrf(i,is_ter)
2391     .               + falblw(i,is_sic) * pctsrf(i,is_sic)
2392      ENDDO
2393      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2394     e            (dist, rmu0, fract,
2395     e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2396     e             wo,
2397     e             cldfra, cldemi, cldtau,
2398     s             heat,heat0,cool,cool0,radsol,albpla,
2399     s             topsw,toplw,solsw,sollw,
2400     s             sollwdown,
2401     s             topsw0,toplw0,solsw0,sollw0,
2402     s             lwdn0, lwdn, lwup0, lwup,
2403     s             swdn0, swdn, swup0, swup,
2404     e             ok_ade, ok_aie, ! new for aerosol radiative effects
2405     e             tau_ae, piz_ae, cg_ae, ! ="=
2406     s             topswad, solswad, ! ="=
2407     e             cldtaupi, ! ="=
2408     s             topswai, solswai) ! ="=
2409      itaprad = 0
2410      ENDIF
2411      itaprad = itaprad + 1
2412c
2413c Ajouter la tendance des rayonnements (tous les pas)
2414c
2415      DO k = 1, klev
2416      DO i = 1, klon
2417         t_seri(i,k) = t_seri(i,k)
2418     .               + (heat(i,k)-cool(i,k)) * dtime/86400.
2419      ENDDO
2420      ENDDO
2421c
2422      IF (if_ebil.ge.2) THEN
2423        ztit='after rad'
2424        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2425     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2426     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2427        call diagphy(airephy,ztit,ip_ebil
2428     e      , topsw, toplw, solsw, sollw, zero_v
2429     e      , zero_v, zero_v, zero_v, ztsol
2430     e      , d_h_vcol, d_qt, d_ec
2431     s      , fs_bound, fq_bound )
2432      END IF
2433c
2434c
2435c Calculer l'hydrologie de la surface
2436c
2437c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2438c     .            agesno, ftsol,fqsurf,fsnow, ruis)
2439c
2440      DO i = 1, klon
2441         zxqsurf(i) = 0.0
2442         zxsnow(i) = 0.0
2443      ENDDO
2444      DO nsrf = 1, nbsrf
2445      DO i = 1, klon
2446         zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2447         zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2448      ENDDO
2449      ENDDO
2450c
2451c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2452c
2453cXXX      DO nsrf = 1, nbsrf
2454cXXX      DO i = 1, klon
2455cXXX         IF (pctsrf(i,nsrf).LT.epsfra) THEN
2456cXXX            fqsurf(i,nsrf) = zxqsurf(i)
2457cXXX            fsnow(i,nsrf) = zxsnow(i)
2458cXXX         ENDIF
2459cXXX      ENDDO
2460cXXX      ENDDO
2461c
2462c Calculer le bilan du sol et la derive de temperature (couplage)
2463c
2464      DO i = 1, klon
2465c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2466c a la demande de JLD
2467         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2468      ENDDO
2469c
2470cmoddeblott(jan95)
2471c Appeler le programme de parametrisation de l'orographie
2472c a l'echelle sous-maille:
2473c
2474      IF (ok_orodr) THEN
2475c
2476c  selection des points pour lesquels le shema est actif:
2477        igwd=0
2478        DO i=1,klon
2479        itest(i)=0
2480c        IF ((zstd(i).gt.10.0)) THEN
2481        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2482          itest(i)=1
2483          igwd=igwd+1
2484          idx(igwd)=i
2485        ENDIF
2486        ENDDO
2487c        igwdim=MAX(1,igwd)
2488c
2489        CALL drag_noro(klon,klev,dtime,paprs,pplay,
2490     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
2491     e                   igwd,idx,itest,
2492     e                   t_seri, u_seri, v_seri,
2493     s                   zulow, zvlow, zustr, zvstr,
2494     s                   d_t_oro, d_u_oro, d_v_oro)
2495c
2496c  ajout des tendances
2497        DO k = 1, klev
2498        DO i = 1, klon
2499           t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2500           u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2501           v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2502        ENDDO
2503        ENDDO
2504c
2505      ENDIF ! fin de test sur ok_orodr
2506c
2507      IF (ok_orolf) THEN
2508c
2509c  selection des points pour lesquels le shema est actif:
2510        igwd=0
2511        DO i=1,klon
2512        itest(i)=0
2513        IF ((zpic(i)-zmea(i)).GT.100.) THEN
2514          itest(i)=1
2515          igwd=igwd+1
2516          idx(igwd)=i
2517        ENDIF
2518        ENDDO
2519c        igwdim=MAX(1,igwd)
2520c
2521        CALL lift_noro(klon,klev,dtime,paprs,pplay,
2522     e                   rlat,zmea,zstd,zpic,
2523     e                   itest,
2524     e                   t_seri, u_seri, v_seri,
2525     s                   zulow, zvlow, zustr, zvstr,
2526     s                   d_t_lif, d_u_lif, d_v_lif)
2527c
2528c  ajout des tendances
2529        DO k = 1, klev
2530        DO i = 1, klon
2531           t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2532           u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2533           v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2534        ENDDO
2535        ENDDO
2536c
2537      ENDIF ! fin de test sur ok_orolf
2538c
2539      IF (if_ebil.ge.2) THEN
2540        ztit='after orography'
2541        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2542     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2543     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2544      END IF
2545c
2546c
2547cAA
2548cAA Installation de l'interface online-offline pour traceurs
2549cAA
2550c====================================================================
2551c   Calcul  des tendances traceurs
2552c====================================================================
2553C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2554cKE43       des traceurs => il faut donc mettre des flags a .false.
2555      IF (iflag_con.GE.3) THEN
2556c           on ajoute les tendances calculees par KE43
2557cXXX OM on onhibe la convection sur les traceurs
2558        DO iq=1, nqmax-2 ! Sandrine a -3 ???
2559cXXX OM on inhibe la convection sur les traceur
2560cXXX        DO k = 1, nlev
2561cXXX        DO i = 1, klon
2562cXXX          tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2563cXXX        ENDDO
2564cXXX        ENDDO
2565        WRITE(iqn,'(i2.2)') iq
2566        CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2567        ENDDO
2568CMAF modif pour garder info du nombre de traceurs auxquels
2569C la physique s'applique
2570      ELSE
2571CMAF modif pour garder info du nombre de traceurs auxquels
2572C la physique s'applique
2573C
2574      call phytrac (rnpb,
2575     I                   itap, julien, gmtime,
2576     I                   debut,lafin,
2577     I                   nqmax-2,
2578     I                   nlon,nlev,dtime,
2579     I                   u,v,t,paprs,pplay,
2580     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2581     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2582     I                   frac_impa, frac_nucl,
2583     I                   rlon,presnivs,pphis,pphi,
2584     I                   albsol,
2585     I                   qx(1,1,1), rhcl,
2586     I                   cldfra, rneb, diafra, cldliq, itop_con,
2587     I                   ibas_con,
2588     I                   pmflxr,pmflxs,prfl,psfl,
2589#ifdef INCA_CH4
2590     I                   flxmass_w,
2591#endif
2592     O                   tr_seri)
2593      ENDIF
2594
2595      IF (offline) THEN
2596
2597         call phystokenc (
2598     I                   nlon,nlev,pdtphys,rlon,rlat,
2599     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2600     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
2601     I                   frac_impa, frac_nucl,
2602     I                   pphis,airephy,dtime,itap)
2603
2604
2605      ENDIF
2606
2607c
2608c Calculer le transport de l'eau et de l'energie (diagnostique)
2609c
2610      CALL transp (paprs,zxtsol,
2611     e                   t_seri, q_seri, u_seri, v_seri, zphi,
2612     s                   ve, vq, ue, uq)
2613c
2614c Accumuler les variables a stocker dans les fichiers histoire:
2615c
2616c
2617c
2618c+jld ec_conser
2619      DO k = 1, klev
2620      DO i = 1, klon
2621        ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i,k))
2622        d_t_ec(i,k)=0.5/ZRCPD
2623     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
2624        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2625        d_t_ec(i,k) = d_t_ec(i,k)/dtime
2626       END DO
2627      END DO
2628c-jld ec_conser
2629      IF (if_ebil.ge.1) THEN
2630        ztit='after physic'
2631        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
2632     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2633     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2634C     Comme les tendances de la physique sont ajoute dans la dynamique,
2635C     on devrait avoir que la variation d'entalpie par la dynamique
2636C     est egale a la variation de la physique au pas de temps precedent.
2637C     Donc la somme de ces 2 variations devrait etre nulle.
2638        call diagphy(airephy,ztit,ip_ebil
2639     e      , topsw, toplw, solsw, sollw, sens
2640     e      , evap, rain_fall, snow_fall, ztsol
2641     e      , d_h_vcol, d_qt, d_ec
2642     s      , fs_bound, fq_bound )
2643C
2644      d_h_vcol_phy=d_h_vcol
2645C
2646      END IF
2647C
2648c=======================================================================
2649c   SORTIES
2650c=======================================================================
2651
2652c   Interpollation sur quelques niveaux de pression
2653c   -----------------------------------------------
2654c
2655c on moyenne mensuellement les champs 3D et on les interpole sur les niveaux STD
2656c     if(itap.EQ.1.OR.itap.EQ.13.OR.itap.EQ.25.OR.itap.EQ.37) THEN
2657c     if(MOD(itap,12).EQ.1) THEN
2658cIM 120304 END
2659      DO k=1, nlevSTD
2660       call plevel(klon,klev,.true.,pplay,rlevSTD(k),
2661     .             t_seri,tlevSTD(:,k))
2662       call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2663     .             u_seri,ulevSTD(:,k))
2664       call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2665     .             v_seri,vlevSTD(:,k))
2666       call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2667     .             zphi,philevSTD(:,k))
2668       call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2669     .             qx(:,:,ivap),qlevSTD(:,k))
2670       call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2671     .             zx_rh,rhlevSTD(:,k))
2672      ENDDO !nlevSTD
2673c ENSEMBLES BEG
2674      DO k=1, nlevENS
2675cIM 170304
2676       tlev(:,k)=tlevSTD(:,indENS(k))
2677       ulev(:,k)=ulevSTD(:,indENS(k))
2678       vlev(:,k)=vlevSTD(:,indENS(k))
2679       philev(:,k)=philevSTD(:,indENS(k))
2680       qlev(:,k)=qlevSTD(:,indENS(k))
2681       rhlev(:,k)=rhlevSTD(:,indENS(k))
2682c
2683       call plevel(klon,klevp1,.true.,paprs,rlevENS(k),
2684     .             omega,wlev(:,k))
2685c
2686       ENDDO !k=1, nlevENS
2687cIM 100304 BEG
2688cIM interpolation a chaque pas de temps du SWup(clr) et SWdn(clr) a 200 hPa
2689      call plevel(klon,klevp1,.true.,paprs,20000.,
2690     $     swdn0,SWdn200clr)
2691      call plevel(klon,klevp1,.false.,paprs,20000.,
2692     $     swdn,SWdn200)
2693      call plevel(klon,klevp1,.false.,paprs,20000.,
2694     $     swup0,SWup200clr)
2695      call plevel(klon,klevp1,.false.,paprs,20000.,
2696     $     swup,SWup200)
2697c
2698      call plevel(klon,klevp1,.false.,paprs,20000.,
2699     $     lwdn0,LWdn200clr)
2700      call plevel(klon,klevp1,.false.,paprs,20000.,
2701     $     lwdn,LWdn200)
2702      call plevel(klon,klevp1,.false.,paprs,20000.,
2703     $     lwup0,LWup200clr)
2704      call plevel(klon,klevp1,.false.,paprs,20000.,
2705     $     lwup,LWup200)
2706c
2707cIM 100304 END
2708c     
2709c ENSEMBLES END
2710c
2711c slp sea level pressure
2712      slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
2713c
2714ccc prw = eau precipitable
2715      DO i = 1, klon
2716       prw(i) = 0.
2717       DO k = 1, klev
2718        prw(i) = prw(i) +
2719     .           q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2720       ENDDO
2721      ENDDO
2722c
2723cIM sorties bilans energie cinetique et potentielle MJO
2724      DO k = 1, klev
2725      DO i = 1, klon
2726        d_u_oli(i,k) = d_u_oro(i,k) + d_u_lif(i,k)
2727        d_v_oli(i,k) = d_v_oro(i,k) + d_v_lif(i,k)
2728      ENDDO
2729      ENDDO
2730c
2731      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
2732cIM      PRINT *,' PHYS cond  julien ',julien
2733c        CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
2734        DO i = 1, klon
2735         total_rain(i)=rain_fall(i)+snow_fall(i) 
2736         IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
2737        ENDDO
2738c
2739      ENDIF
2740c surface terre
2741      IF (debut) THEN
2742       DO i=1, klon
2743         IF(pctsrf_new(i,is_ter).GT.0.) THEN
2744            paire_ter(i)=airephy(i)*pctsrf_new(i,is_ter)
2745         ENDIF
2746       ENDDO
2747      ENDIF
2748cIM 050204 END
2749
2750c=============================================================
2751c
2752c Convertir les incrementations en tendances
2753c
2754      DO k = 1, klev
2755      DO i = 1, klon
2756         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2757         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2758         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2759         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2760         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2761      ENDDO
2762      ENDDO
2763c
2764      IF (nqmax.GE.3) THEN
2765      DO iq = 3, nqmax
2766      DO  k = 1, klev
2767      DO  i = 1, klon
2768         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2769      ENDDO
2770      ENDDO
2771      ENDDO
2772      ENDIF
2773c
2774c Sauvegarder les valeurs de t et q a la fin de la physique:
2775c
2776      DO k = 1, klev
2777      DO i = 1, klon
2778         t_ancien(i,k) = t_seri(i,k)
2779         q_ancien(i,k) = q_seri(i,k)
2780      ENDDO
2781      ENDDO
2782c
2783c 22.03.04 BEG
2784c=============================================================
2785c   Ecriture des sorties
2786c=============================================================
2787#ifdef CPP_IOIPSL
2788
2789#ifdef histhf
2790#include "write_histhf.h"
2791#endif
2792
2793#ifdef histday
2794#include "write_histday.h"
2795#endif
2796
2797#ifdef histmth
2798#include "write_histmth.h"
2799#endif
2800
2801#ifdef histins
2802#include "write_histins.h"
2803#endif
2804
2805#ifdef histREGDYN
2806#include "write_histREGDYN.h"
2807#endif
2808
2809#ifdef histISCCP
2810#include "write_histISCCP.h"
2811#endif
2812
2813#ifdef histmthNMC
2814#include "write_histmthNMC.h"
2815#endif
2816
2817#endif
2818
2819#ifdef INCA
2820#ifdef INCAINFO
2821           WRITE(lunout,*)'Appel CHEMHOOK_END ...'
2822#endif
2823           CALL chemhook_end (calday,
2824     $                        dtime,
2825     $                        pplay,
2826     $                        t_seri,
2827     $                        tr_seri,
2828     $                        nbtr,
2829     $                        paprs,
2830     $                        q_seri,
2831     $                        anne_ini,
2832     $                        day_ini,
2833     $                        xjour)
2834#ifdef INCAINFO
2835           WRITE(lunout,*)'OK.'
2836#endif
2837#endif
2838
2839c 22.03.04 END
2840c
2841c====================================================================
2842c Si c'est la fin, il faut conserver l'etat de redemarrage
2843c====================================================================
2844c
2845      IF (lafin) THEN
2846         itau_phy = itau_phy + itap
2847ccc         IF (ok_oasis) CALL quitcpl
2848         CALL phyredem ("restartphy.nc",dtime,radpas,
2849     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
2850     .      fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
2851     .      solsw, sollwdown,dlw,
2852     .      radsol,frugs,agesno,
2853     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2854     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
2855      ENDIF
2856     
2857
2858      RETURN
2859      END
2860      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2861      IMPLICIT none
2862c
2863c Calculer et imprimer l'eau totale. A utiliser pour verifier
2864c la conservation de l'eau
2865c
2866#include "YOMCST.h"
2867      INTEGER klon,klev
2868      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2869      REAL aire(klon)
2870      REAL qtotal, zx, qcheck
2871      INTEGER i, k
2872c
2873      zx = 0.0
2874      DO i = 1, klon
2875         zx = zx + aire(i)
2876      ENDDO
2877      qtotal = 0.0
2878      DO k = 1, klev
2879      DO i = 1, klon
2880         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2881     .                     *(paprs(i,k)-paprs(i,k+1))/RG
2882      ENDDO
2883      ENDDO
2884c
2885      qcheck = qtotal/zx
2886c
2887      RETURN
2888      END
2889      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2890      IMPLICIT none
2891c
2892c Tranformer une variable de la grille physique a
2893c la grille d'ecriture
2894c
2895      INTEGER nfield,nlon,iim,jjmp1, jjm
2896      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2897c
2898      INTEGER i, n, ig
2899c
2900      jjm = jjmp1 - 1
2901      DO n = 1, nfield
2902         DO i=1,iim
2903            ecrit(i,n) = fi(1,n)
2904            ecrit(i+jjm*iim,n) = fi(nlon,n)
2905         ENDDO
2906         DO ig = 1, nlon - 2
2907           ecrit(iim+ig,n) = fi(1+ig,n)
2908         ENDDO
2909      ENDDO
2910      RETURN
2911      END
Note: See TracBrowser for help on using the repository browser.