source: trunk/LMDZ.VENUS/libf/phyvenus/physiq.F @ 402

Last change on this file since 402 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

  • Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
  • Quelques ajouts dans la doc.
File size: 44.0 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
3!
4c
5      SUBROUTINE physiq (nlon,nlev,nqmax,
6     .            debut,lafin,rjourvrai,gmtime,pdtphys,
7     .            paprs,pplay,ppk,pphi,pphis,presnivs,
8     .            u,v,t,qx,
9     .            omega,
10     .            d_u, d_v, d_t, d_qx, d_ps)
11
12c======================================================================
13c
14c Modifications pour la physique de Venus
15c     S. Lebonnois (LMD/CNRS) Septembre 2005
16c
17c ---------------------------------------------------------------------
18c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
19c
20c Objet: Moniteur general de la physique du modele
21cAA      Modifications quant aux traceurs :
22cAA                  -  uniformisation des parametrisations ds phytrac
23cAA                  -  stockage des moyennes des champs necessaires
24cAA                     en mode traceur off-line
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
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 RDAY 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 ppk  ---input-R-fonction d'Exner au milieu de couche
40c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
41c pphis---input-R-geopotentiel du sol
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-mass mixing ratio traceurs (kg/kg)
47c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
48c omega---input-R-vitesse verticale en Pa/s
49c
50c d_u-----output-R-tendance physique de "u" (m/s/s)
51c d_v-----output-R-tendance physique de "v" (m/s/s)
52c d_t-----output-R-tendance physique de "t" (K/s)
53c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
54c d_ps----output-R-tendance physique de la pression au sol
55c======================================================================
56      USE ioipsl
57      USE histcom
58      USE infotrac
59      USE control_mod
60      use dimphy
61      USE comgeomphy
62      IMPLICIT none
63c======================================================================
64c   CLEFS CPP POUR LES IO
65c   =====================
66c#define histhf
67#define histday
68#define histmth
69#define histins
70c======================================================================
71#include "dimensions.h"
72      integer jjmp1
73      parameter (jjmp1=jjm+1-1/jjm)
74#include "dimsoil.h"
75#include "clesphys.h"
76#include "temps.h"
77#include "iniprint.h"
78#include "timerad.h"
79#include "logic.h"
80c======================================================================
81      LOGICAL ok_journe ! sortir le fichier journalier
82      save ok_journe
83c      PARAMETER (ok_journe=.true.)
84c
85      LOGICAL ok_mensuel ! sortir le fichier mensuel
86      save ok_mensuel
87c      PARAMETER (ok_mensuel=.true.)
88c
89      LOGICAL ok_instan ! sortir le fichier instantane
90      save ok_instan
91c      PARAMETER (ok_instan=.true.)
92c
93c======================================================================
94c
95c Variables argument:
96c
97      INTEGER nlon
98      INTEGER nlev
99      INTEGER nqmax
100      REAL rjourvrai
101      REAL gmtime
102      REAL pdtphys
103      LOGICAL debut, lafin
104      REAL paprs(klon,klev+1)
105      REAL pplay(klon,klev)
106      REAL pphi(klon,klev)
107      REAL pphis(klon)
108      REAL presnivs(klev)
109
110! ADAPTATION GCM POUR CP(T)
111      REAL ppk(klon,klev)
112
113      REAL u(klon,klev)
114      REAL v(klon,klev)
115      REAL t(klon,klev)
116      REAL qx(klon,klev,nqmax)
117
118      REAL,save,allocatable :: t_ancien(:,:)
119      REAL,save,allocatable :: u_ancien(:,:)
120      LOGICAL ancien_ok
121      SAVE ancien_ok
122
123      REAL d_u_dyn(klon,klev)
124      REAL d_t_dyn(klon,klev)
125
126      REAL omega(klon,klev)
127
128      REAL d_u(klon,klev)
129      REAL d_v(klon,klev)
130      REAL d_t(klon,klev)
131      REAL d_qx(klon,klev,nqmax)
132      REAL d_ps(klon)
133
134      REAL,save,allocatable :: swnet(:,:)
135      REAL,save,allocatable :: lwnet(:,:)
136c
137      logical ok_hf
138      real ecrit_hf
139      integer nid_hf
140      save ok_hf, ecrit_hf, nid_hf
141
142#ifdef histhf
143      data ok_hf,ecrit_hf/.true.,0.25/
144#else
145      data ok_hf/.false./
146#endif
147
148c Variables propres a la physique
149c
150      REAL dtime
151      SAVE dtime                  ! pas temporel de la physique
152c
153      INTEGER radpas
154      SAVE radpas                 ! frequence d'appel rayonnement
155c
156      REAL,save,allocatable :: radsol(:) ! bilan radiatif au sol calcule par code radiatif
157      REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche)
158      INTEGER,save :: itap        ! compteur pour la physique
159      REAL,save,allocatable :: ftsol(:)    ! temperature du sol
160      REAL,save,allocatable :: ftsoil(:,:) ! temperature dans le sol
161      REAL,save,allocatable :: falbe(:)    ! albedo
162      REAL delp(klon,klev)        ! epaisseur d'une couche
163     
164CMODDEB FLOTT
165c
166c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
167c
168      REAL,save,allocatable :: zmea(:)   ! orographie moyenne
169      REAL,save,allocatable :: zstd(:)   ! deviation standard de l'OESM
170      REAL,save,allocatable :: zsig(:)   ! pente de l'OESM
171      REAL,save,allocatable :: zgam(:)   ! anisotropie de l'OESM
172      REAL,save,allocatable :: zthe(:)   ! orientation de l'OESM
173      REAL,save,allocatable :: zpic(:)   ! Maximum de l'OESM
174      REAL,save,allocatable :: zval(:)   ! Minimum de l'OESM
175      REAL,save,allocatable :: rugoro(:) ! longueur de rugosite de l'OESM
176
177      INTEGER igwd,idx(klon),itest(klon)
178c
179c  Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
180
181      REAL zulow(klon),zvlow(klon)
182      REAL zustrdr(klon), zvstrdr(klon)
183      REAL zustrli(klon), zvstrli(klon)
184      REAL zustrhi(klon), zvstrhi(klon)
185
186c Pour calcul GW drag oro et nonoro: CALCUL de N2:
187      real zdzlev(klon,klev)
188      real ztlev(klon,klev),zpklev(klon,klev)
189      real ztetalay(klon,klev),ztetalev(klon,klev)
190      real zdtetalev(klon,klev)
191      real zn2(klon,klev) ! BV^2 at plev
192
193c Pour les bilans de moment angulaire,
194      integer bilansmc
195c Pour le transport de ballons
196      integer ballons
197c j'ai aussi besoin
198c du stress de couche limite a la surface:
199
200      REAL zustrcl(klon),zvstrcl(klon)
201
202c et du stress total c de la physique:
203
204      REAL zustrph(klon),zvstrph(klon)
205      REAL,save,allocatable :: zuthe(:),zvthe(:)
206
207c Variables locales:
208c
209      REAL cdragh(klon) ! drag coefficient pour T and Q
210      REAL cdragm(klon) ! drag coefficient pour vent
211c
212cAA  Pour  TRACEURS
213cAA
214      REAL,save,allocatable :: source(:,:)
215      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
216      REAL yu1(klon)            ! vents dans la premiere couche U
217      REAL yv1(klon)            ! vents dans la premiere couche V
218
219      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
220      REAL,save,allocatable :: dlw(:)  ! derivee infra rouge
221      REAL,save,allocatable :: fder(:) ! Derive de flux (sensible et latente)
222      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
223      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
224      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
225      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
226c
227
228c======================================================================
229c
230c Declaration des procedures appelees
231c
232      EXTERNAL ajsec     ! ajustement sec
233      EXTERNAL clmain    ! couche limite
234      EXTERNAL hgardfou  ! verifier les temperatures
235c     EXTERNAL orbite    ! calculer l'orbite
236      EXTERNAL phyetat0  ! lire l'etat initial de la physique
237      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
238      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
239      EXTERNAL suphec    ! initialiser certaines constantes
240      EXTERNAL transp    ! transport total de l'eau et de l'energie
241      INTEGER  lnblnk1
242      EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
243                         !caracter
244      EXTERNAL abort_gcm
245      EXTERNAL printflag
246      EXTERNAL zenang
247      EXTERNAL diagetpq
248      EXTERNAL conf_phys
249      EXTERNAL diagphy
250      EXTERNAL mucorr
251      EXTERNAL phytrac
252c
253c Variables locales
254c
255CXXX PB
256      REAL fluxt(klon,klev)   ! flux turbulent de chaleur
257      REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
258      REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
259c
260      REAL flux_dyn(klon,klev)  ! flux de chaleur produit par la dynamique
261      REAL flux_ajs(klon,klev)  ! flux de chaleur ajustement sec
262      REAL flux_ec(klon,klev)   ! flux de chaleur Ec
263c
264c Le rayonnement n'est pas calcule tous les pas, il faut donc
265c                      sauvegarder les sorties du rayonnement
266      REAL,save,allocatable :: heat(:,:)    ! chauffage solaire
267      REAL,save,allocatable :: cool(:,:)    ! refroidissement infrarouge
268      REAL,save,allocatable :: dtrad(:,:)   ! K s-1
269      REAL,save,allocatable :: topsw(:), toplw(:)
270      REAL,save,allocatable :: solsw(:), sollw(:)
271      REAL,save,allocatable :: sollwdown(:) ! downward LW flux at surface
272      REAL tmpout(klon,klev)  ! K s-1
273
274      INTEGER itaprad
275      SAVE itaprad
276c
277      REAL dist, rmu0(klon), fract(klon)
278      REAL zdtime, zlongi
279c
280      INTEGER i, k, iq, ig, j, ll
281c
282      REAL zphi(klon,klev)
283c
284c Variables du changement
285c
286c ajs: ajustement sec
287c vdf: couche limite (Vertical DiFfusion)
288      REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
289      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
290c
291      REAL d_ts(klon)
292c
293      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
294      REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
295c
296CMOD LOTT: Tendances Orography Sous-maille
297      REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
298      REAL d_t_oro(klon,klev)
299      REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
300      REAL d_t_lif(klon,klev)
301C          Tendances Ondes de G non oro (runs strato).
302      REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
303      REAL d_t_hin(klon,klev)
304
305c
306c Variables liees a l'ecriture de la bande histoire physique
307c
308      INTEGER ecrit_mth
309      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
310c
311      INTEGER ecrit_day
312      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
313c
314      INTEGER ecrit_ins
315      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
316c
317      integer itau_w   ! pas de temps ecriture = itap + itau_phy
318
319c Variables locales pour effectuer les appels en serie
320c
321      REAL t_seri(klon,klev)
322      REAL u_seri(klon,klev), v_seri(klon,klev)
323c
324      REAL tr_seri(klon,klev,nqmax)
325      REAL d_tr(klon,klev,nqmax)
326
327      INTEGER        length
328      PARAMETER    ( length = 100 )
329      REAL tabcntr0( length       )
330c
331c pour ioipsl
332      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
333      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
334      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
335      REAL zx_tmp_2d(iim,jjmp1),zx_tmp_3d(iim,jjmp1,klev)
336      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
337
338      INTEGER nid_day, nid_mth, nid_ins
339      SAVE nid_day, nid_mth, nid_ins
340c
341      INTEGER nhori, nvert, idayref
342      REAL zsto, zout, zsto1, zsto2, zero
343      parameter (zero=0.0e0)
344      real zjulian
345      save zjulian
346
347      CHARACTER*2  str2
348      character*20 modname
349      character*80 abort_message
350      logical ok_sync
351
352      character*30 nom_fichier
353      character*10 varname
354      character*40 vartitle
355      character*20 varunits
356C     Variables liees au bilan d'energie et d'enthalpi
357      REAL ztsol(klon)
358      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
359     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
360      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
361     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
362      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
363      REAL      d_h_vcol_phy
364      REAL      fs_bound, fq_bound
365      SAVE      d_h_vcol_phy
366      REAL      zero_v(klon),zero_v2(klon,klev)
367      CHARACTER*15 ztit
368      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
369      SAVE      ip_ebil
370      DATA      ip_ebil/2/
371      INTEGER   if_ebil ! level for energy conserv. dignostics
372      SAVE      if_ebil
373c+jld ec_conser
374      REAL d_t_ec(klon,klev)    ! tendance du a la conversion Ec -> E thermique
375c-jld ec_conser
376
377c TEST VENUS...
378      REAL mang(klon,klev)    ! moment cinetique
379      REAL mangtot            ! moment cinetique total
380
381c Declaration des constantes et des fonctions thermodynamiques
382c
383#include "YOMCST.h"
384
385c======================================================================
386c INITIALISATIONS
387c================
388
389      modname = 'physiq'
390      ok_sync=.TRUE.
391
392      bilansmc = 0
393      ballons  = 0
394
395      IF (if_ebil.ge.1) THEN
396        DO i=1,klon
397          zero_v(i)=0.
398        END DO
399        DO i=1,klon
400         DO j=1,klev
401          zero_v2(i,j)=0.
402         END DO
403        END DO
404      END IF
405     
406c PREMIER APPEL SEULEMENT
407c========================
408      IF (debut) THEN
409         allocate(t_ancien(klon,klev),u_ancien(klon,klev))
410         allocate(swnet(klon,klevp1),lwnet(klon,klevp1))
411         allocate(radsol(klon),ftsol(klon),falbe(klon))
412         allocate(rlev(klon,klevp1),ftsoil(klon,nsoilmx))
413         allocate(zmea(klon),zstd(klon),zsig(klon),zgam(klon))
414         allocate(zthe(klon),zpic(klon),zval(klon),rugoro(klon))
415         allocate(zuthe(klon),zvthe(klon),dlw(klon),fder(klon))
416         allocate(heat(klon,klev),cool(klon,klev))
417         allocate(dtrad(klon,klev),topsw(klon),toplw(klon))
418         allocate(solsw(klon),sollw(klon),sollwdown(klon))
419         allocate(source(klon,nqmax))
420
421         CALL suphec ! initialiser constantes et parametres phys.
422
423         IF (if_ebil.ge.1) d_h_vcol_phy=0.
424c
425c appel a la lecture du physiq.def
426c
427         call conf_phys(ok_journe, ok_mensuel,
428     .                  ok_instan,
429     .                  if_ebil)
430
431c
432c Initialiser les compteurs:
433c
434         itap    = 0
435         itaprad = 0
436c         
437c Lecture startphy.nc :
438c
439c REMETTRE TOUS LES PARAMETRES POUR OROGW...
440         CALL phyetat0 ("startphy.nc",dtime,
441     .       rlatd,rlond,ftsol,ftsoil,
442     .       falbe, solsw, sollw,
443     .       dlw,radsol,
444     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
445     .       tabcntr0,
446     .       t_ancien, ancien_ok)
447
448c dtime est lu dans startphy
449c pdtphys est calcule a partir des nouvelles conditions:
450c Reinitialisation du pas de temps physique quand changement
451         IF (ABS(dtime-pdtphys).GT.0.001) THEN
452            WRITE(lunout,*) 'Pas physique a change',dtime,
453     .                        pdtphys
454c           abort_message='Pas physique n est pas correct '
455c           call abort_gcm(modname,abort_message,1)
456c----------------
457c pour initialiser convenablement le time_counter, il faut tenir compte
458c du changement de dtime en changeant itau_phy (point de depart)
459            itau_phy = NINT(itau_phy*dtime/pdtphys)
460c----------------
461            dtime=pdtphys
462         ENDIF
463
464         radpas = NINT( RDAY/pdtphys/nbapp_rad)
465
466         CALL printflag( tabcntr0,radpas,ok_journe,ok_instan )
467c
468c---------
469c FLOTT
470       IF (ok_orodr) THEN
471         DO i=1,klon
472         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
473         ENDDO
474         CALL SUGWD(klon,klev,paprs,pplay)
475         DO i=1,klon
476         zuthe(i)=0.
477         zvthe(i)=0.
478         if(zstd(i).gt.10.)then
479           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
480           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
481         endif
482         ENDDO
483       ENDIF
484
485      if (bilansmc.eq.1) then
486C  OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
487C  DU BILAN DE MOMENT ANGULAIRE.
488      open(27,file='aaam_bud.out',form='formatted')
489      open(28,file='fields_2d.out',form='formatted')
490      write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)'
491      write(*,*)'Ouverture de fields_2d.out (FL Vous parle)'
492      endif !bilansmc
493
494c--------------SLEBONNOIS
495C  OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
496C  DES BALLONS
497      if (ballons.eq.1) then
498      open(30,file='ballons-lat.out',form='formatted')
499      open(31,file='ballons-lon.out',form='formatted')
500      open(32,file='ballons-u.out',form='formatted')
501      open(33,file='ballons-v.out',form='formatted')
502      open(34,file='ballons-alt.out',form='formatted')
503      write(*,*)'Ouverture des ballons*.out'
504      endif !ballons
505c-------------
506
507c---------
508C TRACEURS
509C source dans couche limite
510         source = 0.0 ! pas de source, pour l'instant
511C
512c---------
513c
514c Verifications:
515c
516         IF (nlon .NE. klon) THEN
517            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
518     .                      klon
519            abort_message='nlon et klon ne sont pas coherents'
520            call abort_gcm(modname,abort_message,1)
521         ENDIF
522         IF (nlev .NE. klev) THEN
523            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
524     .                       klev
525            abort_message='nlev et klev ne sont pas coherents'
526            call abort_gcm(modname,abort_message,1)
527         ENDIF
528c
529         IF (dtime*FLOAT(radpas).GT.(RDAY*0.25).AND.cycle_diurne)
530     $    THEN
531           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
532           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
533           abort_message='Nbre d appels au rayonnement insuffisant'
534           call abort_gcm(modname,abort_message,1)
535         ENDIF
536c
537         WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
538     .                   iflag_ajs
539c
540         ecrit_mth = NINT(RDAY/dtime*ecritphy)  ! tous les ecritphy jours
541         IF (ok_mensuel) THEN
542         WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
543     .                   ecrit_mth
544         ENDIF
545
546         ecrit_day = NINT(RDAY/dtime *1.0)  ! tous les jours
547         IF (ok_journe) THEN
548         WRITE(lunout,*)'La frequence de sortie journaliere est de ',
549     .                   ecrit_day
550         ENDIF
551
552         ecrit_ins = NINT(RDAY/dtime*ecritphy)  ! Fraction de jour reglable
553         IF (ok_instan) THEN
554         WRITE(lunout,*)'La frequence de sortie instant. est de ',
555     .                   ecrit_ins
556         ENDIF
557
558c Initialisation des sorties
559c===========================
560
561#ifdef CPP_IOIPSL
562
563#ifdef histhf
564#include "ini_histhf.h"
565#endif
566
567#ifdef histday
568#include "ini_histday.h"
569#endif
570
571#ifdef histmth
572#include "ini_histmth.h"
573#endif
574
575#ifdef histins
576#include "ini_histins.h"
577#endif
578
579#endif
580
581c
582c Initialiser les valeurs de u pour calculs tendances
583c (pour T, c'est fait dans phyetat0)
584c
585      DO k = 1, klev
586      DO i = 1, klon
587         u_ancien(i,k) = u(i,k)
588      ENDDO
589      ENDDO
590
591         
592      ENDIF ! debut
593c====================================================================
594c======================================================================
595
596c Mettre a zero des variables de sortie (pour securite)
597c
598      DO i = 1, klon
599         d_ps(i) = 0.0
600      ENDDO
601      DO k = 1, klev
602      DO i = 1, klon
603         d_t(i,k) = 0.0
604         d_u(i,k) = 0.0
605         d_v(i,k) = 0.0
606      ENDDO
607      ENDDO
608      DO iq = 1, nqmax
609      DO k = 1, klev
610      DO i = 1, klon
611         d_qx(i,k,iq) = 0.0
612      ENDDO
613      ENDDO
614      ENDDO
615c
616c Ne pas affecter les valeurs entrees de u, v, h, et q
617c
618      DO k = 1, klev
619      DO i = 1, klon
620         t_seri(i,k)  = t(i,k)
621         u_seri(i,k)  = u(i,k)
622         v_seri(i,k)  = v(i,k)
623      ENDDO
624      ENDDO
625      DO iq = 1, nqmax
626      DO  k = 1, klev
627      DO  i = 1, klon
628         tr_seri(i,k,iq) = qx(i,k,iq)
629      ENDDO
630      ENDDO
631      ENDDO
632C
633      DO i = 1, klon
634          ztsol(i) = ftsol(i)
635      ENDDO
636C
637      IF (if_ebil.ge.1) THEN
638        ztit='after dynamic'
639        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
640     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
641     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
642C     Comme les tendances de la physique sont ajoute dans la dynamique,
643C     on devrait avoir que la variation d'entalpie par la dynamique
644C     est egale a la variation de la physique au pas de temps precedent.
645C     Donc la somme de ces 2 variations devrait etre nulle.
646        call diagphy(airephy,ztit,ip_ebil
647     e      , zero_v, zero_v, zero_v, zero_v, zero_v
648     e      , zero_v, zero_v, zero_v, ztsol
649     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
650     s      , fs_bound, fq_bound )
651      END IF
652
653c====================================================================
654c Diagnostiquer la tendance dynamique
655c
656      IF (ancien_ok) THEN
657         DO k = 1, klev
658         DO i = 1, klon
659            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
660            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
661         ENDDO
662         ENDDO
663
664! ADAPTATION GCM POUR CP(T)
665         do i=1,klon
666          flux_dyn(i,1) = 0.0
667          do j=2,klev
668            flux_dyn(i,j) = flux_dyn(i,j-1)
669     . +cpdet(t_seri(i,j-1))/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j))
670          enddo
671         enddo
672         
673      ELSE
674         DO k = 1, klev
675         DO i = 1, klon
676            d_u_dyn(i,k) = 0.0
677            d_t_dyn(i,k) = 0.0
678         ENDDO
679         ENDDO
680         ancien_ok = .TRUE.
681      ENDIF
682c====================================================================
683c
684c Ajouter le geopotentiel du sol:
685c
686      DO k = 1, klev
687      DO i = 1, klon
688         zphi(i,k) = pphi(i,k) + pphis(i)
689      ENDDO
690      ENDDO
691c====================================================================
692c
693c Verifier les temperatures
694c
695      CALL hgardfou(t_seri,ftsol,'debutphy')
696c====================================================================
697c
698c Incrementer le compteur de la physique
699c
700      itap   = itap + 1
701
702c====================================================================
703c Orbite et eclairement
704c====================================================================
705
706c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
707c donc pas de variations de Ls, ni de dist.
708c La seule chose qui compte, c'est la rotation de la planete devant
709c le Soleil...
710     
711      zlongi = 0.0
712      dist   = 0.72  ! en UA
713
714c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
715c a sa valeur, et prendre en compte leur evolution,
716c il faudra refaire un orbite.F...
717c     CALL orbite(zlongi,dist)
718
719      IF (cycle_diurne) THEN
720        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
721        CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
722      ELSE
723        call mucorr(klon,zlongi,rlatd,rmu0,fract)
724      ENDIF
725     
726c====================================================================
727c   Calcul  des tendances traceurs
728c====================================================================
729
730      if (iflag_trac.eq.1) then
731      call phytrac (
732     I                   itap, gmtime,
733     I                   debut,lafin,
734     I                   nqmax,
735     I                   nlon,nlev,dtime,
736     I                   u,v,t,paprs,pplay,
737     I                   rlatd,
738     I                   rlond,presnivs,pphis,pphi,
739     I                   falbe,
740     O                   tr_seri)
741      endif
742
743c
744c====================================================================
745c Appeler la diffusion verticale (programme de couche limite)
746c====================================================================
747
748c-------------------------------
749c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
750c l'equilibre radiatif du sol
751      if (1.eq.0) then
752              if (debut) then
753                print*,"ATTENTION, CLMAIN SHUNTEE..."
754              endif
755
756      DO i = 1, klon
757         sens(i) = 0.0e0 ! flux de chaleur sensible au sol
758         fder(i) = 0.0e0
759         dlw(i)  = 0.0e0
760      ENDDO
761
762c Incrementer la temperature du sol
763c
764      DO i = 1, klon
765         d_ts(i)  = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
766         ftsol(i) = ftsol(i) + d_ts(i)
767         do j=1,nsoilmx
768           ftsoil(i,j)=ftsol(i)
769         enddo
770      ENDDO
771
772c-------------------------------
773      else
774c-------------------------------
775
776      fder = dlw
777
778c     print*,"radsol avant clmain=",radsol(klon/2)
779c     print*,"solsw avant clmain=",solsw(klon/2)
780c     print*,"sollw avant clmain=",sollw(klon/2)
781
782! ADAPTATION GCM POUR CP(T)
783      CALL clmain(dtime,itap,
784     e            t_seri,u_seri,v_seri,
785     e            rmu0,
786     e            ftsol,
787     $            ftsoil,
788     $            paprs,pplay,ppk,radsol,falbe,
789     e            solsw, sollw, sollwdown, fder,
790     e            rlond, rlatd, cuphy, cvphy,   
791     e            debut, lafin,
792     s            d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
793     s            fluxt,fluxu,fluxv,cdragh,cdragm,
794     s            dsens,
795     s            ycoefh,yu1,yv1)
796
797c     print*,"radsol apres clmain=",radsol(klon/2)
798c     print*,"solsw apres clmain=",solsw(klon/2)
799c     print*,"sollw apres clmain=",sollw(klon/2)
800
801CXXX Incrementation des flux
802      DO i = 1, klon
803         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
804         fder(i) = dlw(i) + dsens(i)
805      ENDDO
806CXXX
807
808      DO k = 1, klev
809      DO i = 1, klon
810         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
811         d_t_vdf(i,k)= d_t_vdf(i,k)/dtime          ! K/s
812         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
813         d_u_vdf(i,k)= d_u_vdf(i,k)/dtime          ! (m/s)/s
814         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
815         d_v_vdf(i,k)= d_v_vdf(i,k)/dtime          ! (m/s)/s
816      ENDDO
817      ENDDO
818c
819c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
820c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
821c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
822c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
823c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime
824
825C TRACEURS
826
827      if (iflag_trac.eq.1) then
828         DO k = 1, klev
829         DO i = 1, klon
830            delp(i,k) = paprs(i,k)-paprs(i,k+1)
831         ENDDO
832         ENDDO
833         DO iq=1, nqmax
834             CALL cltrac(dtime,ycoefh,t_seri,
835     s               tr_seri(1,1,iq),source,
836     e               paprs, pplay,delp,
837     s               d_tr_vdf(1,1,iq))
838             tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
839             d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime          ! /s
840         ENDDO
841      endif
842
843      IF (if_ebil.ge.2) THEN
844        ztit='after clmain'
845        CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
846     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
847     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
848         call diagphy(airephy,ztit,ip_ebil
849     e      , zero_v, zero_v, zero_v, zero_v, sens
850     e      , zero_v, zero_v, zero_v, ztsol
851     e      , d_h_vcol, d_qt, d_ec
852     s      , fs_bound, fq_bound )
853      END IF
854C
855c
856c Incrementer la temperature du sol
857c
858c     print*,'Tsol avant clmain:',ftsol(1)
859      DO i = 1, klon
860         ftsol(i) = ftsol(i) + d_ts(i)
861      ENDDO
862c     print*,'Tsol apres clmain:',ftsol(1)
863
864c Calculer la derive du flux infrarouge
865c
866      DO i = 1, klon
867            dlw(i) = - 4.0*RSIGMA*ftsol(i)**3
868      ENDDO
869
870c-------------------------------
871      endif  ! fin du VENUS TEST
872
873c
874c Appeler l'ajustement sec
875c
876c===================================================================
877c Convection seche
878c===================================================================
879c
880      d_t_ajs(:,:)=0.
881      d_u_ajs(:,:)=0.
882      d_v_ajs(:,:)=0.
883      d_tr_ajs(:,:,:)=0.
884c
885      IF(prt_level>9)WRITE(lunout,*)
886     .    'AVANT LA CONVECTION SECHE , iflag_ajs='
887     s   ,iflag_ajs
888
889      if(iflag_ajs.eq.0) then
890c  Rien
891c  ====
892         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
893
894      else if(iflag_ajs.eq.1) then
895
896c  Ajustement sec
897c  ==============
898         IF(prt_level>9)WRITE(lunout,*)'ajsec'
899
900! ADAPTATION GCM POUR CP(T)
901         CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
902     .              tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
903
904! ADAPTATION GCM POUR CP(T)
905         do i=1,klon
906          flux_ajs(i,1) = 0.0
907          do j=2,klev
908            flux_ajs(i,j) = flux_ajs(i,j-1)
909     .        + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
910     .                                 *(paprs(i,j-1)-paprs(i,j))
911          enddo
912         enddo
913         
914         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
915         d_t_ajs(:,:)= d_t_ajs(:,:)/dtime          ! K/s
916         u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
917         d_u_ajs(:,:)= d_u_ajs(:,:)/dtime          ! (m/s)/s
918         v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
919         d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
920      if (iflag_trac.eq.1) then
921           tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
922           d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
923      endif
924
925c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
926c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
927c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
928c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
929c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime
930
931      endif
932c
933      IF (if_ebil.ge.2) THEN
934        ztit='after dry_adjust'
935        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
936     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
937     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
938        call diagphy(airephy,ztit,ip_ebil
939     e      , zero_v, zero_v, zero_v, zero_v, sens
940     e      , zero_v, zero_v, zero_v, ztsol
941     e      , d_h_vcol, d_qt, d_ec
942     s      , fs_bound, fq_bound )
943      END IF
944
945c====================================================================
946c RAYONNEMENT
947c====================================================================
948
949      IF (MOD(itaprad,radpas).EQ.0) THEN
950c             print*,'RAYONNEMENT ',
951c    $             ' (itaprad=',itaprad,'/radpas=',radpas,')'
952
953      dtimerad = dtime*FLOAT(radpas)  ! pas de temps du rayonnement (s)
954c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
955           
956c     print*,"radsol avant radlwsw=",radsol(klon/2)
957c     print*,"solsw avant radlwsw=",solsw(klon/2)
958c     print*,"sollw avant radlwsw=",sollw(klon/2)
959c     print*,"avant radlwsw"
960      CALL radlwsw
961     e            (dist, rmu0, fract,
962     e             paprs, pplay,ftsol, t_seri,
963     s             heat,cool,radsol,
964     s             topsw,toplw,solsw,sollw,
965     s             sollwdown,
966     s             lwnet, swnet)
967c     print*,"apres radlwsw"
968
969c     print*,"radsol apres radlwsw=",radsol(klon/2)
970c     print*,"solsw apres radlwsw=",solsw(klon/2)
971c     print*,"sollw apres radlwsw=",sollw(klon/2)
972      itaprad = 0
973      DO k = 1, klev
974      DO i = 1, klon
975         dtrad(i,k) = heat(i,k)-cool(i,k)
976         dtrad(i,k) = dtrad(i,k)/RDAY  !K/s
977      ENDDO
978      ENDDO
979c        print*,"dtrad1=",dtrad(1,:)*dtime
980c        print*,"dtrad2=",dtrad(klon/2,:)*dtime
981c        print*,"dtrad3=",dtrad(klon,:)*dtime
982     
983      ENDIF
984      itaprad = itaprad + 1
985c====================================================================
986c
987c Ajouter la tendance des rayonnements (tous les pas)
988c
989      DO k = 1, klev
990      DO i = 1, klon
991         t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
992      ENDDO
993      ENDDO
994 
995      IF (if_ebil.ge.2) THEN
996        ztit='after rad'
997        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
998     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
999     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1000        call diagphy(airephy,ztit,ip_ebil
1001     e      , topsw, toplw, solsw, sollw, zero_v
1002     e      , zero_v, zero_v, zero_v, ztsol
1003     e      , d_h_vcol, d_qt, d_ec
1004     s      , fs_bound, fq_bound )
1005      END IF
1006c
1007
1008c====================================================================
1009c   Calcul  des gravity waves  FLOTT
1010c====================================================================
1011c
1012      if (ok_orodr.or.ok_gw_nonoro) then
1013c  CALCUL DE N2
1014       do i=1,klon
1015        do k=2,klev
1016          ztlev(i,k)  = (t_seri(i,k)+t_seri(i,k-1))/2.
1017          zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1018        enddo
1019       enddo
1020       call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1021       call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1022       do i=1,klon
1023        do k=2,klev
1024          zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1025          zdzlev(i,k)    = (zphi(i,k)-zphi(i,k-1))/RG
1026          zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k))
1027          zn2(i,k) = max(zn2(i,k),1.e-12)  ! securite
1028        enddo
1029       enddo
1030
1031      endif
1032     
1033c ----------------------------ORODRAG
1034      IF (ok_orodr) THEN
1035c
1036c  selection des points pour lesquels le shema est actif:
1037        igwd=0
1038        DO i=1,klon
1039        itest(i)=0
1040c        IF ((zstd(i).gt.10.0)) THEN
1041        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1042          itest(i)=1
1043          igwd=igwd+1
1044          idx(igwd)=i
1045        ENDIF
1046        ENDDO
1047c        igwdim=MAX(1,igwd)
1048c
1049c A ADAPTER POUR VENUS!!!
1050        CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
1051     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
1052     e                   igwd,idx,itest,
1053     e                   t_seri, u_seri, v_seri,
1054     s                   zulow, zvlow, zustrdr, zvstrdr,
1055     s                   d_t_oro, d_u_oro, d_v_oro)
1056
1057c  ajout des tendances
1058           t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1059           d_t_oro(:,:)= d_t_oro(:,:)/dtime          ! K/s
1060           u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1061           d_u_oro(:,:)= d_u_oro(:,:)/dtime          ! (m/s)/s
1062           v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1063           d_v_oro(:,:)= d_v_oro(:,:)/dtime          ! (m/s)/s
1064c   
1065      ELSE
1066         d_t_oro = 0.
1067         d_u_oro = 0.
1068         d_v_oro = 0.
1069         zustrdr = 0.
1070         zvstrdr = 0.
1071c
1072      ENDIF ! fin de test sur ok_orodr
1073c
1074c ----------------------------OROLIFT
1075      IF (ok_orolf) THEN
1076c
1077c  selection des points pour lesquels le shema est actif:
1078        igwd=0
1079        DO i=1,klon
1080        itest(i)=0
1081        IF ((zpic(i)-zmea(i)).GT.100.) THEN
1082          itest(i)=1
1083          igwd=igwd+1
1084          idx(igwd)=i
1085        ENDIF
1086        ENDDO
1087c        igwdim=MAX(1,igwd)
1088c
1089c A ADAPTER POUR VENUS!!!
1090c            CALL lift_noro(klon,klev,dtime,paprs,pplay,
1091c     e                   rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1092c     e                   igwd,idx,itest,
1093c     e                   t_seri, u_seri, v_seri,
1094c     s                   zulow, zvlow, zustrli, zvstrli,
1095c     s                   d_t_lif, d_u_lif, d_v_lif               )
1096
1097c
1098c  ajout des tendances
1099           t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1100           d_t_lif(:,:)= d_t_lif(:,:)/dtime          ! K/s
1101           u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1102           d_u_lif(:,:)= d_u_lif(:,:)/dtime          ! (m/s)/s
1103           v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1104           d_v_lif(:,:)= d_v_lif(:,:)/dtime          ! (m/s)/s
1105c
1106      ELSE
1107         d_t_lif = 0.
1108         d_u_lif = 0.
1109         d_v_lif = 0.
1110         zustrli = 0.
1111         zvstrli = 0.
1112c
1113      ENDIF ! fin de test sur ok_orolf
1114
1115c ---------------------------- NON-ORO GRAVITY WAVES
1116       IF(ok_gw_nonoro) then
1117
1118      call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1119     e               t_seri, u_seri, v_seri,
1120     o               zustrhi,zvstrhi,
1121     o               d_t_hin, d_u_hin, d_v_hin)
1122
1123c  ajout des tendances
1124
1125         t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1126         d_t_hin(:,:)= d_t_hin(:,:)/dtime          ! K/s
1127         u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1128         d_u_hin(:,:)= d_u_hin(:,:)/dtime          ! (m/s)/s
1129         v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1130         d_v_hin(:,:)= d_v_hin(:,:)/dtime          ! (m/s)/s
1131
1132      ELSE
1133         d_t_hin = 0.
1134         d_u_hin = 0.
1135         d_v_hin = 0.
1136         zustrhi = 0.
1137         zvstrhi = 0.
1138
1139      ENDIF ! fin de test sur ok_gw_nonoro
1140
1141c====================================================================
1142c Transport de ballons
1143c====================================================================
1144      if (ballons.eq.1) then
1145         CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
1146c    C               t,pplay,u,v,pphi)   ! alt above surface (smoothed for GCM)
1147     C               t,pplay,u,v,zphi)   ! alt above planet average radius
1148      endif !ballons
1149
1150c====================================================================
1151c Bilan de mmt angulaire
1152c====================================================================
1153      if (bilansmc.eq.1) then
1154CMODDEB FLOTT
1155C  CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1156C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1157
1158      DO i = 1, klon
1159        zustrph(i)=0.
1160        zvstrph(i)=0.
1161        zustrcl(i)=0.
1162        zvstrcl(i)=0.
1163      ENDDO
1164      DO k = 1, klev
1165      DO i = 1, klon
1166       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1167     c            (paprs(i,k)-paprs(i,k+1))/rg
1168       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1169     c            (paprs(i,k)-paprs(i,k+1))/rg
1170       zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1171     c            (paprs(i,k)-paprs(i,k+1))/rg
1172       zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1173     c            (paprs(i,k)-paprs(i,k+1))/rg
1174      ENDDO
1175      ENDDO
1176
1177      CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1178     C               ra,rg,romega,
1179     C               rlatd,rlond,pphis,
1180     C               zustrdr,zustrli,zustrcl,
1181     C               zvstrdr,zvstrli,zvstrcl,
1182     C               paprs,u,v)
1183                     
1184CCMODFIN FLOTT
1185      endif !bilansmc
1186
1187c====================================================================
1188c====================================================================
1189c Calculer le transport de l'eau et de l'energie (diagnostique)
1190c
1191c  A REVOIR POUR VENUS...
1192c
1193c     CALL transp (paprs,ftsol,
1194c    e                   t_seri, q_seri, u_seri, v_seri, zphi,
1195c    s                   ve, vq, ue, uq)
1196c
1197c====================================================================
1198c+jld ec_conser
1199      DO k = 1, klev
1200      DO i = 1, klon
1201        d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1202     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
1203        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1204        d_t_ec(i,k) = d_t_ec(i,k)/dtime
1205       END DO
1206      END DO
1207         do i=1,klon
1208          flux_ec(i,1) = 0.0
1209          do j=2,klev
1210            flux_ec(i,j) = flux_ec(i,j-1)
1211     . +cpdet(t_seri(i,j-1))/RG*d_t_ec(i,j-1)*(paprs(i,j-1)-paprs(i,j))
1212          enddo
1213         enddo
1214         
1215c-jld ec_conser
1216c====================================================================
1217      IF (if_ebil.ge.1) THEN
1218        ztit='after physic'
1219        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1220     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1221     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1222C     Comme les tendances de la physique sont ajoute dans la dynamique,
1223C     on devrait avoir que la variation d'entalpie par la dynamique
1224C     est egale a la variation de la physique au pas de temps precedent.
1225C     Donc la somme de ces 2 variations devrait etre nulle.
1226        call diagphy(airephy,ztit,ip_ebil
1227     e      , topsw, toplw, solsw, sollw, sens
1228     e      , zero_v, zero_v, zero_v, ztsol
1229     e      , d_h_vcol, d_qt, d_ec
1230     s      , fs_bound, fq_bound )
1231C
1232      d_h_vcol_phy=d_h_vcol
1233C
1234      END IF
1235C
1236c=======================================================================
1237c   SORTIES
1238c=======================================================================
1239
1240c Convertir les incrementations en tendances
1241c
1242      DO k = 1, klev
1243      DO i = 1, klon
1244         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1245         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1246         d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1247      ENDDO
1248      ENDDO
1249c
1250      DO iq = 1, nqmax
1251      DO  k = 1, klev
1252      DO  i = 1, klon
1253         d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1254      ENDDO
1255      ENDDO
1256      ENDDO
1257     
1258c------------------------
1259c Calcul moment cinetique
1260c------------------------
1261c TEST VENUS...
1262c     mangtot = 0.0
1263c     DO k = 1, klev
1264c     DO i = 1, klon
1265c       mang(i,k) = RA*cos(rlatd(i)*RPI/180.)
1266c    .     *(u_seri(i,k)+RA*cos(rlatd(i)*RPI/180.)*ROMEGA)
1267c    .     *airephy(i)*(paprs(i,k)-paprs(i,k+1))/RG
1268c       mangtot=mangtot+mang(i,k)
1269c     ENDDO
1270c     ENDDO
1271c     print*,"Moment cinetique total = ",mangtot
1272c
1273c------------------------
1274c
1275c Sauvegarder les valeurs de t et u a la fin de la physique:
1276c
1277      DO k = 1, klev
1278      DO i = 1, klon
1279         u_ancien(i,k) = u_seri(i,k)
1280         t_ancien(i,k) = t_seri(i,k)
1281      ENDDO
1282      ENDDO
1283c
1284c=============================================================
1285c   Ecriture des sorties
1286c=============================================================
1287     
1288#ifdef CPP_IOIPSL
1289
1290#ifdef histhf
1291#include "write_histhf.h"
1292#endif
1293
1294#ifdef histday
1295#include "write_histday.h"
1296#endif
1297
1298#ifdef histmth
1299#include "write_histmth.h"
1300#endif
1301
1302#ifdef histins
1303#include "write_histins.h"
1304#endif
1305
1306#endif
1307
1308c ---------- Sorties testphys1d -------------
1309     
1310      if (klon.eq.1) then
1311        call writeg1d(klon,klev,t_seri,"Temp","Temperature")
1312        call writeg1d(klon,1,ftsol,"tsurf","Surface Temp")
1313      DO k = 1, klev
1314      DO i = 1, klon
1315c        tmpout(i,k) = real(heat(i,k))/RDAY
1316         tmpout(i,k) = heat(i,k)/RDAY
1317      ENDDO
1318      ENDDO
1319        call writeg1d(klon,klev,tmpout,
1320     .                 "heat","Solar heating")
1321      DO k = 1, klev
1322      DO i = 1, klon
1323c        tmpout(i,k) = real(dtrad(i,k))
1324         tmpout(i,k) = dtrad(i,k)
1325      ENDDO
1326      ENDDO
1327        call writeg1d(klon,klev,tmpout,
1328     .                 "dtrad","IR cooling")
1329        call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux")
1330        call writeg1d(klon,klev,swnet,"swnet","Net SW flux")
1331        call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux")
1332        call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux")
1333        call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux")
1334c       call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface")
1335c       call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface")
1336        call writeg1d(klon,1,radsol,"surfnet","Net flux at surface")
1337        call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain")
1338        call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec")
1339        call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec")
1340      endif
1341     
1342c====================================================================
1343c Si c'est la fin, il faut conserver l'etat de redemarrage
1344c====================================================================
1345c
1346      IF (lafin) THEN
1347         itau_phy = itau_phy + itap
1348c REMETTRE TOUS LES PARAMETRES POUR OROGW...
1349         CALL phyredem ("restartphy.nc",dtime,radpas,
1350     .      rlatd, rlond, ftsol, ftsoil,
1351     .      falbe,
1352     .      solsw, sollw,dlw,
1353     .      radsol,
1354     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
1355     .      t_ancien)
1356     
1357c--------------FLOTT
1358CMODEB LOTT
1359C  FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1360C  DU BILAN DE MOMENT ANGULAIRE.
1361      if (bilansmc.eq.1) then
1362        write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)'
1363        close(27)                                     
1364        close(28)                                     
1365      endif !bilansmc
1366CMODFIN
1367c-------------
1368c--------------SLEBONNOIS
1369C  FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1370C  DES BALLONS
1371      if (ballons.eq.1) then
1372        write(*,*)'Fermeture des ballons*.out'
1373        close(30)                                     
1374        close(31)                                     
1375        close(32)                                     
1376        close(33)                                     
1377        close(34)                                     
1378      endif !ballons
1379c-------------
1380      ENDIF
1381     
1382
1383      RETURN
1384      END
1385
1386
1387
1388***********************************************************************
1389***********************************************************************
1390***********************************************************************
1391***********************************************************************
1392***********************************************************************
1393***********************************************************************
1394***********************************************************************
1395***********************************************************************
1396
1397      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
1398      IMPLICIT none
1399c
1400c Tranformer une variable de la grille physique a
1401c la grille d'ecriture
1402c
1403      INTEGER nfield,nlon,iim,jjmp1, jjm
1404      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
1405c
1406      INTEGER i, n, ig
1407c
1408      jjm = jjmp1 - 1
1409      DO n = 1, nfield
1410         DO i=1,iim
1411            ecrit(i,n) = fi(1,n)
1412            ecrit(i+jjm*iim,n) = fi(nlon,n)
1413         ENDDO
1414         DO ig = 1, nlon - 2
1415           ecrit(iim+ig,n) = fi(1+ig,n)
1416         ENDDO
1417      ENDDO
1418      RETURN
1419      END
Note: See TracBrowser for help on using the repository browser.