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

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

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