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

Last change on this file since 1198 was 1160, checked in by slebonnois, 11 years ago

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