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

Last change on this file since 856 was 849, checked in by emillour, 13 years ago

Work on common dynamics and interfacing with different physics:

  • Put calls to PVtheta in dynamics between CPP_EARTH flags (because it calls tetalevel, which is supposed to be in the physics; only OK for Earth...).
  • Adapted makelmdz script so that one can compile main programs in dyn* or phy* (makelmdz_fcm already capable of doing that).
  • Moved start2archive-VT.F and start2archive-VT.F to phyvenus (as start2archive.F and newstart.F); leave adapting them to Titan for later.
  • Small correction to phyvenus/testphys1d.F (use module control_mod instead of control.h and call disvert_noterre).
  • removed "use histcom" in phyvenus/physiq.F and phytitan/physiq.F ; it is not needed since there is already a "use ioipsl" (and it moreover confused makelmdz_fcm...)
  • Had to add declaration of variable "zlsm1" in phytitan/physiq.F because it is used in "write_hist.h"; but note that it is used while not initialized (but what should it be initialized to?).

EM

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