source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/phylmd/physiq.F @ 3421

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

Pour la conservation du flux d'eau OM, JLD
LF

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