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

Last change on this file since 104 was 101, checked in by slebonnois, 14 years ago

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