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

Last change on this file since 1301 was 1301, checked in by slebonnois, 10 years ago

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