source: trunk/libf/phyvenus/physiq.F @ 4

Last change on this file since 4 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

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