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

Last change on this file was 556, checked in by lmdzadmin, 20 years ago

Initialisations diverses necessaires au portage Prism AC
LF

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