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

Last change on this file since 973 was 973, checked in by slebonnois, 12 years ago

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