source: trunk/libf/phytitan/physiq.F @ 53

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