source: trunk/libf/phyvenus/physiq.F @ 124

Last change on this file since 124 was 119, checked in by slebonnois, 14 years ago

Sebastien Lebonnois: apres validation des versions Venus et Titan,
correction d'un certain nombre de bugs.
File size: 43.7 KB
RevLine 
[3]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,
[97]7     .            paprs,pplay,ppk,pphi,pphis,presnivs,
[3]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======================================================================
[101]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======================================================================
[3]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
[101]118      REAL,save,allocatable :: t_ancien(:,:)
119      REAL,save,allocatable :: u_ancien(:,:)
[3]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
[101]134      REAL,save,allocatable :: swnet(:,:)
135      REAL,save,allocatable :: lwnet(:,:)
[3]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
[101]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
[3]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
[101]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
[3]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)
[101]205      REAL,save,allocatable :: zuthe(:),zvthe(:)
[3]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
[105]214      REAL,save,allocatable :: source(:,:)
[3]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
[101]220      REAL,save,allocatable :: dlw(:)  ! derivee infra rouge
221      REAL,save,allocatable :: fder(:) ! Derive de flux (sensible et latente)
[3]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
[101]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
[3]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
[97]342      REAL zsto, zout, zsto1, zsto2, zero
[3]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
[119]392      bilansmc = 0
[3]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
[101]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))
[105]419         allocate(source(klon,nqmax))
[3]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,
[97]443     .       dlw,radsol,
[3]444     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
445     .       tabcntr0,
446     .       t_ancien, ancien_ok)
447
448         radpas = NINT( RDAY/pdtphys/nbapp_rad)
449
450         CALL printflag( tabcntr0,radpas,ok_journe,ok_instan )
451c
452c---------
453c FLOTT
454       IF (ok_orodr) THEN
455         DO i=1,klon
456         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
457         ENDDO
458         CALL SUGWD(klon,klev,paprs,pplay)
459         DO i=1,klon
460         zuthe(i)=0.
461         zvthe(i)=0.
462         if(zstd(i).gt.10.)then
463           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
464           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
465         endif
466         ENDDO
467       ENDIF
468
469      if (bilansmc.eq.1) then
470C  OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
471C  DU BILAN DE MOMENT ANGULAIRE.
472      open(27,file='aaam_bud.out',form='formatted')
473      open(28,file='fields_2d.out',form='formatted')
474      write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)'
475      write(*,*)'Ouverture de fields_2d.out (FL Vous parle)'
476      endif !bilansmc
477
478c--------------SLEBONNOIS
479C  OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
480C  DES BALLONS
481      if (ballons.eq.1) then
482      open(30,file='ballons-lat.out',form='formatted')
483      open(31,file='ballons-lon.out',form='formatted')
484      open(32,file='ballons-u.out',form='formatted')
485      open(33,file='ballons-v.out',form='formatted')
486      open(34,file='ballons-alt.out',form='formatted')
487      write(*,*)'Ouverture des ballons*.out'
488      endif !ballons
489c-------------
490
491c---------
492C TRACEURS
493C source dans couche limite
494         source = 0.0 ! pas de source, pour l'instant
495C
496c---------
497c
498c Verifications:
499c
500         IF (ABS(dtime-pdtphys).GT.0.001) THEN
501            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
502     .                        pdtphys
503c           abort_message='Pas physique n est pas correct '
504c           call abort_gcm(modname,abort_message,1)
505            dtime=pdtphys
506         ENDIF
507         IF (nlon .NE. klon) THEN
508            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
509     .                      klon
510            abort_message='nlon et klon ne sont pas coherents'
511            call abort_gcm(modname,abort_message,1)
512         ENDIF
513         IF (nlev .NE. klev) THEN
514            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
515     .                       klev
516            abort_message='nlev et klev ne sont pas coherents'
517            call abort_gcm(modname,abort_message,1)
518         ENDIF
519c
520         IF (dtime*FLOAT(radpas).GT.(RDAY*0.25).AND.cycle_diurne)
521     $    THEN
522           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
523           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
524           abort_message='Nbre d appels au rayonnement insuffisant'
525           call abort_gcm(modname,abort_message,1)
526         ENDIF
527c
528         WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
529     .                   iflag_ajs
530c
[97]531         ecrit_mth = NINT(RDAY/dtime*ecritphy)  ! tous les ecritphy jours
[3]532         IF (ok_mensuel) THEN
533         WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
534     .                   ecrit_mth
535         ENDIF
[97]536
[3]537         ecrit_day = NINT(RDAY/dtime *1.0)  ! tous les jours
538         IF (ok_journe) THEN
539         WRITE(lunout,*)'La frequence de sortie journaliere est de ',
540     .                   ecrit_day
541         ENDIF
[97]542
[101]543         ecrit_ins = NINT(RDAY/dtime*ecritphy)  ! Fraction de jour reglable
[3]544         IF (ok_instan) THEN
545         WRITE(lunout,*)'La frequence de sortie instant. est de ',
546     .                   ecrit_ins
547         ENDIF
548
549c Initialisation des sorties
550c===========================
551
552#ifdef CPP_IOIPSL
553
554#ifdef histhf
555#include "ini_histhf.h"
556#endif
557
558#ifdef histday
559#include "ini_histday.h"
560#endif
561
562#ifdef histmth
563#include "ini_histmth.h"
564#endif
565
566#ifdef histins
567#include "ini_histins.h"
568#endif
569
570#endif
571
572c
573c Initialiser les valeurs de u pour calculs tendances
574c (pour T, c'est fait dans phyetat0)
575c
576      DO k = 1, klev
577      DO i = 1, klon
578         u_ancien(i,k) = u(i,k)
579      ENDDO
580      ENDDO
581
582         
583      ENDIF ! debut
584c====================================================================
585c======================================================================
586
587c Mettre a zero des variables de sortie (pour securite)
588c
589      DO i = 1, klon
590         d_ps(i) = 0.0
591      ENDDO
592      DO k = 1, klev
593      DO i = 1, klon
594         d_t(i,k) = 0.0
595         d_u(i,k) = 0.0
596         d_v(i,k) = 0.0
597      ENDDO
598      ENDDO
599      DO iq = 1, nqmax
600      DO k = 1, klev
601      DO i = 1, klon
602         d_qx(i,k,iq) = 0.0
603      ENDDO
604      ENDDO
605      ENDDO
606c
607c Ne pas affecter les valeurs entrees de u, v, h, et q
608c
609      DO k = 1, klev
610      DO i = 1, klon
611         t_seri(i,k)  = t(i,k)
612         u_seri(i,k)  = u(i,k)
613         v_seri(i,k)  = v(i,k)
614      ENDDO
615      ENDDO
616      DO iq = 1, nqmax
617      DO  k = 1, klev
618      DO  i = 1, klon
619         tr_seri(i,k,iq) = qx(i,k,iq)
620      ENDDO
621      ENDDO
622      ENDDO
623C
624      DO i = 1, klon
625          ztsol(i) = ftsol(i)
626      ENDDO
627C
628      IF (if_ebil.ge.1) THEN
629        ztit='after dynamic'
630        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
631     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
632     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
633C     Comme les tendances de la physique sont ajoute dans la dynamique,
634C     on devrait avoir que la variation d'entalpie par la dynamique
635C     est egale a la variation de la physique au pas de temps precedent.
636C     Donc la somme de ces 2 variations devrait etre nulle.
637        call diagphy(airephy,ztit,ip_ebil
638     e      , zero_v, zero_v, zero_v, zero_v, zero_v
639     e      , zero_v, zero_v, zero_v, ztsol
640     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
641     s      , fs_bound, fq_bound )
642      END IF
643
644c====================================================================
645c Diagnostiquer la tendance dynamique
646c
647      IF (ancien_ok) THEN
648         DO k = 1, klev
649         DO i = 1, klon
650            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
651            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
652         ENDDO
653         ENDDO
654
655! ADAPTATION GCM POUR CP(T)
656         do i=1,klon
657          flux_dyn(i,1) = 0.0
658          do j=2,klev
659            flux_dyn(i,j) = flux_dyn(i,j-1)
660     . +cpdet(t_seri(i,j-1))/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j))
661          enddo
662         enddo
663         
664      ELSE
665         DO k = 1, klev
666         DO i = 1, klon
667            d_u_dyn(i,k) = 0.0
668            d_t_dyn(i,k) = 0.0
669         ENDDO
670         ENDDO
671         ancien_ok = .TRUE.
672      ENDIF
673c====================================================================
674c
675c Ajouter le geopotentiel du sol:
676c
677      DO k = 1, klev
678      DO i = 1, klon
679         zphi(i,k) = pphi(i,k) + pphis(i)
680      ENDDO
681      ENDDO
682c====================================================================
683c
684c Verifier les temperatures
685c
686      CALL hgardfou(t_seri,ftsol,'debutphy')
687c====================================================================
688c
689c Incrementer le compteur de la physique
690c
691      itap   = itap + 1
692
693c====================================================================
694c Orbite et eclairement
695c====================================================================
696
697c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
698c donc pas de variations de Ls, ni de dist.
699c La seule chose qui compte, c'est la rotation de la planete devant
700c le Soleil...
701     
702      zlongi = 0.0
703      dist   = 0.72  ! en UA
704
705c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
706c a sa valeur, et prendre en compte leur evolution,
707c il faudra refaire un orbite.F...
708c     CALL orbite(zlongi,dist)
709
710      IF (cycle_diurne) THEN
711        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
712        CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
713      ELSE
714        call mucorr(klon,zlongi,rlatd,rmu0,fract)
715      ENDIF
716     
717c====================================================================
718c   Calcul  des tendances traceurs
719c====================================================================
720
721      if (iflag_trac.eq.1) then
722      call phytrac (
723     I                   itap, gmtime,
724     I                   debut,lafin,
725     I                   nqmax,
726     I                   nlon,nlev,dtime,
727     I                   u,v,t,paprs,pplay,
728     I                   rlatd,
729     I                   rlond,presnivs,pphis,pphi,
730     I                   falbe,
731     O                   tr_seri)
732      endif
733
734c
735c====================================================================
736c Appeler la diffusion verticale (programme de couche limite)
737c====================================================================
738
739c-------------------------------
740c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
741c l'equilibre radiatif du sol
742      if (1.eq.0) then
743              if (debut) then
744                print*,"ATTENTION, CLMAIN SHUNTEE..."
745              endif
746
747      DO i = 1, klon
748         sens(i) = 0.0e0 ! flux de chaleur sensible au sol
749         fder(i) = 0.0e0
750         dlw(i)  = 0.0e0
751      ENDDO
752
753c Incrementer la temperature du sol
754c
755      DO i = 1, klon
756         d_ts(i)  = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
757         ftsol(i) = ftsol(i) + d_ts(i)
758         do j=1,nsoilmx
759           ftsoil(i,j)=ftsol(i)
760         enddo
761      ENDDO
762
763c-------------------------------
764      else
765c-------------------------------
766
767      fder = dlw
768
769c     print*,"radsol avant clmain=",radsol(klon/2)
770c     print*,"solsw avant clmain=",solsw(klon/2)
771c     print*,"sollw avant clmain=",sollw(klon/2)
772
773! ADAPTATION GCM POUR CP(T)
774      CALL clmain(dtime,itap,
775     e            t_seri,u_seri,v_seri,
776     e            rmu0,
777     e            ftsol,
778     $            ftsoil,
779     $            paprs,pplay,ppk,radsol,falbe,
780     e            solsw, sollw, sollwdown, fder,
781     e            rlond, rlatd, cuphy, cvphy,   
782     e            debut, lafin,
783     s            d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
784     s            fluxt,fluxu,fluxv,cdragh,cdragm,
785     s            dsens,
786     s            ycoefh,yu1,yv1)
787
788c     print*,"radsol apres clmain=",radsol(klon/2)
789c     print*,"solsw apres clmain=",solsw(klon/2)
790c     print*,"sollw apres clmain=",sollw(klon/2)
791
792CXXX Incrementation des flux
793      DO i = 1, klon
794         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
795         fder(i) = dlw(i) + dsens(i)
796      ENDDO
797CXXX
798
799      DO k = 1, klev
800      DO i = 1, klon
801         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
802         d_t_vdf(i,k)= d_t_vdf(i,k)/dtime          ! K/s
803         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
804         d_u_vdf(i,k)= d_u_vdf(i,k)/dtime          ! (m/s)/s
805         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
806         d_v_vdf(i,k)= d_v_vdf(i,k)/dtime          ! (m/s)/s
807      ENDDO
808      ENDDO
809c
810c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
811c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
812c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
813c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
814c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime
815
816C TRACEURS
817
818      if (iflag_trac.eq.1) then
[101]819         DO k = 1, klev
820         DO i = 1, klon
821            delp(i,k) = paprs(i,k)-paprs(i,k+1)
822         ENDDO
823         ENDDO
[3]824         DO iq=1, nqmax
825             CALL cltrac(dtime,ycoefh,t_seri,
826     s               tr_seri(1,1,iq),source,
827     e               paprs, pplay,delp,
828     s               d_tr_vdf(1,1,iq))
829             tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
830             d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime          ! /s
831         ENDDO
832      endif
833
834      IF (if_ebil.ge.2) THEN
835        ztit='after clmain'
836        CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
837     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
838     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
839         call diagphy(airephy,ztit,ip_ebil
840     e      , zero_v, zero_v, zero_v, zero_v, sens
841     e      , zero_v, zero_v, zero_v, ztsol
842     e      , d_h_vcol, d_qt, d_ec
843     s      , fs_bound, fq_bound )
844      END IF
845C
846c
847c Incrementer la temperature du sol
848c
849c     print*,'Tsol avant clmain:',ftsol(1)
850      DO i = 1, klon
851         ftsol(i) = ftsol(i) + d_ts(i)
852      ENDDO
853c     print*,'Tsol apres clmain:',ftsol(1)
854
855c Calculer la derive du flux infrarouge
856c
857      DO i = 1, klon
858            dlw(i) = - 4.0*RSIGMA*ftsol(i)**3
859      ENDDO
860
861c-------------------------------
862      endif  ! fin du VENUS TEST
863
864c
865c Appeler l'ajustement sec
866c
867c===================================================================
868c Convection seche
869c===================================================================
870c
871      d_t_ajs(:,:)=0.
872      d_u_ajs(:,:)=0.
873      d_v_ajs(:,:)=0.
874      d_tr_ajs(:,:,:)=0.
875c
876      IF(prt_level>9)WRITE(lunout,*)
877     .    'AVANT LA CONVECTION SECHE , iflag_ajs='
878     s   ,iflag_ajs
879
880      if(iflag_ajs.eq.0) then
881c  Rien
882c  ====
883         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
884
885      else if(iflag_ajs.eq.1) then
886
887c  Ajustement sec
888c  ==============
889         IF(prt_level>9)WRITE(lunout,*)'ajsec'
890
891! ADAPTATION GCM POUR CP(T)
892         CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
893     .              tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
894
895! ADAPTATION GCM POUR CP(T)
896         do i=1,klon
897          flux_ajs(i,1) = 0.0
898          do j=2,klev
899            flux_ajs(i,j) = flux_ajs(i,j-1)
900     .        + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
901     .                                 *(paprs(i,j-1)-paprs(i,j))
902          enddo
903         enddo
904         
905         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
906         d_t_ajs(:,:)= d_t_ajs(:,:)/dtime          ! K/s
907         u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
908         d_u_ajs(:,:)= d_u_ajs(:,:)/dtime          ! (m/s)/s
909         v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
910         d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
911      if (iflag_trac.eq.1) then
912           tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
913           d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
914      endif
915
916c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
917c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
918c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
919c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
920c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime
921
922      endif
923c
924      IF (if_ebil.ge.2) THEN
925        ztit='after dry_adjust'
926        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
927     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
928     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
929        call diagphy(airephy,ztit,ip_ebil
930     e      , zero_v, zero_v, zero_v, zero_v, sens
931     e      , zero_v, zero_v, zero_v, ztsol
932     e      , d_h_vcol, d_qt, d_ec
933     s      , fs_bound, fq_bound )
934      END IF
935
936c====================================================================
937c RAYONNEMENT
938c====================================================================
939
940      IF (MOD(itaprad,radpas).EQ.0) THEN
941c             print*,'RAYONNEMENT ',
942c    $             ' (itaprad=',itaprad,'/radpas=',radpas,')'
943
944      dtimerad = dtime*FLOAT(radpas)  ! pas de temps du rayonnement (s)
945c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
946           
947c     print*,"radsol avant radlwsw=",radsol(klon/2)
948c     print*,"solsw avant radlwsw=",solsw(klon/2)
949c     print*,"sollw avant radlwsw=",sollw(klon/2)
950c     print*,"avant radlwsw"
951      CALL radlwsw
952     e            (dist, rmu0, fract,
953     e             paprs, pplay,ftsol, t_seri,
954     s             heat,cool,radsol,
955     s             topsw,toplw,solsw,sollw,
956     s             sollwdown,
957     s             lwnet, swnet)
958c     print*,"apres radlwsw"
959
960c     print*,"radsol apres radlwsw=",radsol(klon/2)
961c     print*,"solsw apres radlwsw=",solsw(klon/2)
962c     print*,"sollw apres radlwsw=",sollw(klon/2)
963      itaprad = 0
964      DO k = 1, klev
965      DO i = 1, klon
966         dtrad(i,k) = heat(i,k)-cool(i,k)
967         dtrad(i,k) = dtrad(i,k)/RDAY  !K/s
968      ENDDO
969      ENDDO
970c        print*,"dtrad1=",dtrad(1,:)*dtime
971c        print*,"dtrad2=",dtrad(klon/2,:)*dtime
972c        print*,"dtrad3=",dtrad(klon,:)*dtime
973     
974      ENDIF
975      itaprad = itaprad + 1
976c====================================================================
977c
978c Ajouter la tendance des rayonnements (tous les pas)
979c
980      DO k = 1, klev
981      DO i = 1, klon
982         t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
983      ENDDO
984      ENDDO
985 
986      IF (if_ebil.ge.2) THEN
987        ztit='after rad'
988        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
989     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
990     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
991        call diagphy(airephy,ztit,ip_ebil
992     e      , topsw, toplw, solsw, sollw, zero_v
993     e      , zero_v, zero_v, zero_v, ztsol
994     e      , d_h_vcol, d_qt, d_ec
995     s      , fs_bound, fq_bound )
996      END IF
997c
998
999c====================================================================
1000c   Calcul  des gravity waves  FLOTT
1001c====================================================================
1002c
1003      if (ok_orodr.or.ok_gw_nonoro) then
1004c  CALCUL DE N2
1005       do i=1,klon
1006        do k=2,klev
1007          ztlev(i,k)  = (t_seri(i,k)+t_seri(i,k-1))/2.
1008          zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1009        enddo
1010       enddo
1011       call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1012       call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1013       do i=1,klon
1014        do k=2,klev
1015          zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1016          zdzlev(i,k)    = (zphi(i,k)-zphi(i,k-1))/RG
1017          zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k))
1018          zn2(i,k) = max(zn2(i,k),1.e-12)  ! securite
1019        enddo
1020       enddo
1021
1022      endif
1023     
1024c ----------------------------ORODRAG
1025      IF (ok_orodr) THEN
1026c
1027c  selection des points pour lesquels le shema est actif:
1028        igwd=0
1029        DO i=1,klon
1030        itest(i)=0
1031c        IF ((zstd(i).gt.10.0)) THEN
1032        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1033          itest(i)=1
1034          igwd=igwd+1
1035          idx(igwd)=i
1036        ENDIF
1037        ENDDO
1038c        igwdim=MAX(1,igwd)
1039c
1040c A ADAPTER POUR VENUS!!!
1041        CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
1042     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
1043     e                   igwd,idx,itest,
1044     e                   t_seri, u_seri, v_seri,
1045     s                   zulow, zvlow, zustrdr, zvstrdr,
1046     s                   d_t_oro, d_u_oro, d_v_oro)
1047
1048c  ajout des tendances
1049           t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1050           d_t_oro(:,:)= d_t_oro(:,:)/dtime          ! K/s
1051           u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1052           d_u_oro(:,:)= d_u_oro(:,:)/dtime          ! (m/s)/s
1053           v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1054           d_v_oro(:,:)= d_v_oro(:,:)/dtime          ! (m/s)/s
1055c   
1056      ELSE
1057         d_t_oro = 0.
1058         d_u_oro = 0.
1059         d_v_oro = 0.
1060         zustrdr = 0.
1061         zvstrdr = 0.
1062c
1063      ENDIF ! fin de test sur ok_orodr
1064c
1065c ----------------------------OROLIFT
1066      IF (ok_orolf) THEN
1067c
1068c  selection des points pour lesquels le shema est actif:
1069        igwd=0
1070        DO i=1,klon
1071        itest(i)=0
1072        IF ((zpic(i)-zmea(i)).GT.100.) THEN
1073          itest(i)=1
1074          igwd=igwd+1
1075          idx(igwd)=i
1076        ENDIF
1077        ENDDO
1078c        igwdim=MAX(1,igwd)
1079c
1080c A ADAPTER POUR VENUS!!!
1081c            CALL lift_noro(klon,klev,dtime,paprs,pplay,
1082c     e                   rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1083c     e                   igwd,idx,itest,
1084c     e                   t_seri, u_seri, v_seri,
1085c     s                   zulow, zvlow, zustrli, zvstrli,
1086c     s                   d_t_lif, d_u_lif, d_v_lif               )
1087
1088c
1089c  ajout des tendances
1090           t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1091           d_t_lif(:,:)= d_t_lif(:,:)/dtime          ! K/s
1092           u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1093           d_u_lif(:,:)= d_u_lif(:,:)/dtime          ! (m/s)/s
1094           v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1095           d_v_lif(:,:)= d_v_lif(:,:)/dtime          ! (m/s)/s
1096c
1097      ELSE
1098         d_t_lif = 0.
1099         d_u_lif = 0.
1100         d_v_lif = 0.
1101         zustrli = 0.
1102         zvstrli = 0.
1103c
1104      ENDIF ! fin de test sur ok_orolf
1105
1106c ---------------------------- NON-ORO GRAVITY WAVES
1107       IF(ok_gw_nonoro) then
1108
1109      call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1110     e               t_seri, u_seri, v_seri,
1111     o               zustrhi,zvstrhi,
1112     o               d_t_hin, d_u_hin, d_v_hin)
1113
1114c  ajout des tendances
1115
1116         t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1117         d_t_hin(:,:)= d_t_hin(:,:)/dtime          ! K/s
1118         u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1119         d_u_hin(:,:)= d_u_hin(:,:)/dtime          ! (m/s)/s
1120         v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1121         d_v_hin(:,:)= d_v_hin(:,:)/dtime          ! (m/s)/s
1122
1123      ELSE
1124         d_t_hin = 0.
1125         d_u_hin = 0.
1126         d_v_hin = 0.
1127         zustrhi = 0.
1128         zvstrhi = 0.
1129
1130      ENDIF ! fin de test sur ok_gw_nonoro
1131
1132c====================================================================
1133c Transport de ballons
1134c====================================================================
1135      if (ballons.eq.1) then
1136         CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
1137c    C               t,pplay,u,v,pphi)   ! alt above surface (smoothed for GCM)
1138     C               t,pplay,u,v,zphi)   ! alt above planet average radius
1139      endif !ballons
1140
1141c====================================================================
1142c Bilan de mmt angulaire
1143c====================================================================
1144      if (bilansmc.eq.1) then
1145CMODDEB FLOTT
1146C  CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1147C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1148
1149      DO i = 1, klon
1150        zustrph(i)=0.
1151        zvstrph(i)=0.
1152        zustrcl(i)=0.
1153        zvstrcl(i)=0.
1154      ENDDO
1155      DO k = 1, klev
1156      DO i = 1, klon
1157       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1158     c            (paprs(i,k)-paprs(i,k+1))/rg
1159       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1160     c            (paprs(i,k)-paprs(i,k+1))/rg
1161       zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1162     c            (paprs(i,k)-paprs(i,k+1))/rg
1163       zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1164     c            (paprs(i,k)-paprs(i,k+1))/rg
1165      ENDDO
1166      ENDDO
1167
1168      CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1169     C               ra,rg,romega,
1170     C               rlatd,rlond,pphis,
1171     C               zustrdr,zustrli,zustrcl,
1172     C               zvstrdr,zvstrli,zvstrcl,
1173     C               paprs,u,v)
1174                     
1175CCMODFIN FLOTT
1176      endif !bilansmc
1177
1178c====================================================================
1179c====================================================================
1180c Calculer le transport de l'eau et de l'energie (diagnostique)
1181c
1182c  A REVOIR POUR VENUS...
1183c
1184c     CALL transp (paprs,ftsol,
1185c    e                   t_seri, q_seri, u_seri, v_seri, zphi,
1186c    s                   ve, vq, ue, uq)
1187c
1188c====================================================================
1189c+jld ec_conser
1190      DO k = 1, klev
1191      DO i = 1, klon
1192        d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1193     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
1194        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1195        d_t_ec(i,k) = d_t_ec(i,k)/dtime
1196       END DO
1197      END DO
1198         do i=1,klon
1199          flux_ec(i,1) = 0.0
1200          do j=2,klev
1201            flux_ec(i,j) = flux_ec(i,j-1)
1202     . +cpdet(t_seri(i,j-1))/RG*d_t_ec(i,j-1)*(paprs(i,j-1)-paprs(i,j))
1203          enddo
1204         enddo
1205         
1206c-jld ec_conser
1207c====================================================================
1208      IF (if_ebil.ge.1) THEN
1209        ztit='after physic'
1210        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1211     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1212     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1213C     Comme les tendances de la physique sont ajoute dans la dynamique,
1214C     on devrait avoir que la variation d'entalpie par la dynamique
1215C     est egale a la variation de la physique au pas de temps precedent.
1216C     Donc la somme de ces 2 variations devrait etre nulle.
1217        call diagphy(airephy,ztit,ip_ebil
1218     e      , topsw, toplw, solsw, sollw, sens
1219     e      , zero_v, zero_v, zero_v, ztsol
1220     e      , d_h_vcol, d_qt, d_ec
1221     s      , fs_bound, fq_bound )
1222C
1223      d_h_vcol_phy=d_h_vcol
1224C
1225      END IF
1226C
1227c=======================================================================
1228c   SORTIES
1229c=======================================================================
1230
1231c Convertir les incrementations en tendances
1232c
1233      DO k = 1, klev
1234      DO i = 1, klon
1235         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1236         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1237         d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1238      ENDDO
1239      ENDDO
1240c
1241      DO iq = 1, nqmax
1242      DO  k = 1, klev
1243      DO  i = 1, klon
1244         d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1245      ENDDO
1246      ENDDO
1247      ENDDO
1248     
1249c------------------------
1250c Calcul moment cinetique
1251c------------------------
1252c TEST VENUS...
1253c     mangtot = 0.0
1254c     DO k = 1, klev
1255c     DO i = 1, klon
1256c       mang(i,k) = RA*cos(rlatd(i)*RPI/180.)
1257c    .     *(u_seri(i,k)+RA*cos(rlatd(i)*RPI/180.)*ROMEGA)
1258c    .     *airephy(i)*(paprs(i,k)-paprs(i,k+1))/RG
1259c       mangtot=mangtot+mang(i,k)
1260c     ENDDO
1261c     ENDDO
1262c     print*,"Moment cinetique total = ",mangtot
1263c
1264c------------------------
1265c
1266c Sauvegarder les valeurs de t et u a la fin de la physique:
1267c
1268      DO k = 1, klev
1269      DO i = 1, klon
1270         u_ancien(i,k) = u_seri(i,k)
1271         t_ancien(i,k) = t_seri(i,k)
1272      ENDDO
1273      ENDDO
1274c
1275c=============================================================
1276c   Ecriture des sorties
1277c=============================================================
1278     
1279#ifdef CPP_IOIPSL
1280
1281#ifdef histhf
1282#include "write_histhf.h"
1283#endif
1284
1285#ifdef histday
1286#include "write_histday.h"
1287#endif
1288
1289#ifdef histmth
1290#include "write_histmth.h"
1291#endif
1292
1293#ifdef histins
1294#include "write_histins.h"
1295#endif
1296
1297#endif
1298
1299c ---------- Sorties testphys1d -------------
1300     
1301      if (klon.eq.1) then
1302        call writeg1d(klon,klev,t_seri,"Temp","Temperature")
1303        call writeg1d(klon,1,ftsol,"tsurf","Surface Temp")
1304      DO k = 1, klev
1305      DO i = 1, klon
1306c        tmpout(i,k) = real(heat(i,k))/RDAY
1307         tmpout(i,k) = heat(i,k)/RDAY
1308      ENDDO
1309      ENDDO
1310        call writeg1d(klon,klev,tmpout,
1311     .                 "heat","Solar heating")
1312      DO k = 1, klev
1313      DO i = 1, klon
1314c        tmpout(i,k) = real(dtrad(i,k))
1315         tmpout(i,k) = dtrad(i,k)
1316      ENDDO
1317      ENDDO
1318        call writeg1d(klon,klev,tmpout,
1319     .                 "dtrad","IR cooling")
1320        call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux")
1321        call writeg1d(klon,klev,swnet,"swnet","Net SW flux")
1322        call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux")
1323        call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux")
1324        call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux")
1325c       call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface")
1326c       call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface")
1327        call writeg1d(klon,1,radsol,"surfnet","Net flux at surface")
1328        call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain")
1329        call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec")
1330        call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec")
1331      endif
1332     
1333c====================================================================
1334c Si c'est la fin, il faut conserver l'etat de redemarrage
1335c====================================================================
1336c
1337      IF (lafin) THEN
1338         itau_phy = itau_phy + itap
1339c REMETTRE TOUS LES PARAMETRES POUR OROGW...
1340         CALL phyredem ("restartphy.nc",dtime,radpas,
1341     .      rlatd, rlond, ftsol, ftsoil,
1342     .      falbe,
1343     .      solsw, sollw,dlw,
1344     .      radsol,
1345     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
1346     .      t_ancien)
1347     
1348c--------------FLOTT
1349CMODEB LOTT
1350C  FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1351C  DU BILAN DE MOMENT ANGULAIRE.
1352      if (bilansmc.eq.1) then
1353        write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)'
1354        close(27)                                     
1355        close(28)                                     
1356      endif !bilansmc
1357CMODFIN
1358c-------------
1359c--------------SLEBONNOIS
1360C  FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1361C  DES BALLONS
1362      if (ballons.eq.1) then
1363        write(*,*)'Fermeture des ballons*.out'
1364        close(30)                                     
1365        close(31)                                     
1366        close(32)                                     
1367        close(33)                                     
1368        close(34)                                     
1369      endif !ballons
1370c-------------
1371      ENDIF
1372     
1373
1374      RETURN
1375      END
1376
1377
1378
1379***********************************************************************
1380***********************************************************************
1381***********************************************************************
1382***********************************************************************
1383***********************************************************************
1384***********************************************************************
1385***********************************************************************
1386***********************************************************************
1387
1388      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
1389      IMPLICIT none
1390c
1391c Tranformer une variable de la grille physique a
1392c la grille d'ecriture
1393c
1394      INTEGER nfield,nlon,iim,jjmp1, jjm
1395      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
1396c
1397      INTEGER i, n, ig
1398c
1399      jjm = jjmp1 - 1
1400      DO n = 1, nfield
1401         DO i=1,iim
1402            ecrit(i,n) = fi(1,n)
1403            ecrit(i+jjm*iim,n) = fi(nlon,n)
1404         ENDDO
1405         DO ig = 1, nlon - 2
1406           ecrit(iim+ig,n) = fi(1+ig,n)
1407         ENDDO
1408      ENDDO
1409      RETURN
1410      END
Note: See TracBrowser for help on using the repository browser.