Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (31 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/bilan_dyn.F90

    r5245 r5246  
    22! $Id$
    33!
    4       SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
    5      s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
    6 
    7 c   AFAIRE
    8 c   Prevoir en champ nq+1 le diagnostique de l'energie
    9 c   en faisant Qzon=Cv T + L * ...
    10 c             vQ..A=Cp T + L * ...
     4SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum, &
     5        ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
     6
     7  !   AFAIRE
     8  !   Prevoir en champ nq+1 le diagnostique de l'energie
     9  !   en faisant Qzon=Cv T + L * ...
     10  !             vQ..A=Cp T + L * ...
    1111
    1212#ifdef CPP_IOIPSL
    13       USE IOIPSL
     13  USE IOIPSL
    1414#endif
    15       USE comconst_mod, ONLY: pi, cpp
    16       USE comvert_mod, ONLY: presnivs
    17       USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
    18 
    19       IMPLICIT NONE
    20 
    21       include "dimensions.h"
    22       include "paramet.h"
    23       include "comgeom2.h"
    24       include "iniprint.h"
    25 
    26 c====================================================================
    27 c
    28 c   Sous-programme consacre à des diagnostics dynamiques de base
    29 c
    30 c
    31 c   De facon generale, les moyennes des scalaires Q sont ponderees par
    32 c   la masse.
    33 c
    34 c   Les flux de masse sont eux simplement moyennes.
    35 c
    36 c====================================================================
    37 
    38 c   Arguments :
    39 c   ===========
    40 
    41       integer ntrac
    42       real dt_app,dt_cum
    43       real ps(iip1,jjp1)
    44       real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
    45       real flux_u(iip1,jjp1,llm)
    46       real flux_v(iip1,jjm,llm)
    47       real teta(iip1,jjp1,llm)
    48       real phi(iip1,jjp1,llm)
    49       real ucov(iip1,jjp1,llm)
    50       real vcov(iip1,jjm,llm)
    51       real trac(iip1,jjp1,llm,ntrac)
    52 
    53 c   Local :
    54 c   =======
    55 
    56       integer icum,ncum
    57       logical first
    58       real zz,zqy,zfactv(jjm,llm)
    59 
    60       integer nQ
    61       parameter (nQ=7)
    62 
    63 
    64 cym      character*6 nom(nQ)
    65 cym      character*6 unites(nQ)
    66       character*6,save :: nom(nQ)
    67       character*6,save :: unites(nQ)
    68 
    69       character*10 file
    70       integer ifile
    71       parameter (ifile=4)
    72 
    73       integer itemp,igeop,iecin,iang,iu,iovap,iun
    74       integer i_sortie
    75 
    76       save first,icum,ncum
    77       save itemp,igeop,iecin,iang,iu,iovap,iun
    78       save i_sortie
    79 
    80       real time
    81       integer itau
    82       save time,itau
    83       data time,itau/0.,0/
    84 
    85       data first/.true./
    86       data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
    87       data i_sortie/1/
    88 
    89       real ww
    90 
    91 c   variables dynamiques intermédiaires
    92       REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
    93       REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
    94       REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
    95       REAL vorpot(iip1,jjm,llm)
    96       REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
    97       REAL bern(iip1,jjp1,llm)
    98 
    99 c   champ contenant les scalaires advectés.
    100       real Q(iip1,jjp1,llm,nQ)
    101    
    102 c   champs cumulés
    103       real ps_cum(iip1,jjp1)
    104       real masse_cum(iip1,jjp1,llm)
    105       real flux_u_cum(iip1,jjp1,llm)
    106       real flux_v_cum(iip1,jjm,llm)
    107       real Q_cum(iip1,jjp1,llm,nQ)
    108       real flux_uQ_cum(iip1,jjp1,llm,nQ)
    109       real flux_vQ_cum(iip1,jjm,llm,nQ)
    110       real flux_wQ_cum(iip1,jjp1,llm,nQ)
    111       real dQ(iip1,jjp1,llm,nQ)
    112 
    113       save ps_cum,masse_cum,flux_u_cum,flux_v_cum
    114       save Q_cum,flux_uQ_cum,flux_vQ_cum
    115 
    116 c   champs de tansport en moyenne zonale
    117       integer ntr,itr
    118       parameter (ntr=5)
    119 
    120 cym      character*10 znom(ntr,nQ)
    121 cym      character*20 znoml(ntr,nQ)
    122 cym      character*10 zunites(ntr,nQ)
    123       character*10,save :: znom(ntr,nQ)
    124       character*20,save :: znoml(ntr,nQ)
    125       character*10,save :: zunites(ntr,nQ)
    126 
    127       integer iave,itot,immc,itrs,istn
    128       data iave,itot,immc,itrs,istn/1,2,3,4,5/
    129       character*3 ctrs(ntr)
    130       data ctrs/'  ','TOT','MMC','TRS','STN'/
    131 
    132       real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
    133       real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
    134       real zmasse(jjm,llm),zamasse(jjm)
    135 
    136       real zv(jjm,llm),psi(jjm,llm+1)
    137 
    138       integer i,j,l,iQ
    139 
    140 
    141 c   Initialisation du fichier contenant les moyennes zonales.
    142 c   ---------------------------------------------------------
    143 
    144       character*10 infile
    145 
    146       integer fileid
    147       integer thoriid, zvertiid
    148       save fileid
    149 
    150       integer ndex3d(jjm*llm)
    151 
    152 C   Variables locales
    153 C
    154       integer tau0
    155       real zjulian
    156       character*3 str
    157       character*10 ctrac
    158       integer ii,jj
    159       integer zan, dayref
    160 C
    161       real rlong(jjm),rlatg(jjm)
    162 
    163 
    164 
    165 c=====================================================================
    166 c   Initialisation
    167 c=====================================================================
    168 
    169       time=time+dt_app
    170       itau=itau+1
    171 cIM
    172       ndex3d=0
    173 
    174       if (first) then
    175 
    176 
    177         icum=0
    178 c       initialisation des fichiers
    179         first=.false.
    180 c   ncum est la frequence de stokage en pas de temps
    181         ncum=dt_cum/dt_app
    182         if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
    183            WRITE(lunout,*)
    184      .            'Pb : le pas de cumule doit etre multiple du pas'
    185            WRITE(lunout,*)'dt_app=',dt_app
    186            WRITE(lunout,*)'dt_cum=',dt_cum
    187            call abort_gcm('bilan_dyn','stopped',1)
     15  USE comconst_mod, ONLY: pi, cpp
     16  USE comvert_mod, ONLY: presnivs
     17  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
     18
     19  IMPLICIT NONE
     20
     21  include "dimensions.h"
     22  include "paramet.h"
     23  include "comgeom2.h"
     24  include "iniprint.h"
     25
     26  !====================================================================
     27  !
     28  !   Sous-programme consacre à des diagnostics dynamiques de base
     29  !
     30  !
     31  !   De facon generale, les moyennes des scalaires Q sont ponderees par
     32  !   la masse.
     33  !
     34  !   Les flux de masse sont eux simplement moyennes.
     35  !
     36  !====================================================================
     37
     38  !   Arguments :
     39  !   ===========
     40
     41  integer :: ntrac
     42  real :: dt_app,dt_cum
     43  real :: ps(iip1,jjp1)
     44  real :: masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
     45  real :: flux_u(iip1,jjp1,llm)
     46  real :: flux_v(iip1,jjm,llm)
     47  real :: teta(iip1,jjp1,llm)
     48  real :: phi(iip1,jjp1,llm)
     49  real :: ucov(iip1,jjp1,llm)
     50  real :: vcov(iip1,jjm,llm)
     51  real :: trac(iip1,jjp1,llm,ntrac)
     52
     53  !   Local :
     54  !   =======
     55
     56  integer :: icum,ncum
     57  logical :: first
     58  real :: zz,zqy,zfactv(jjm,llm)
     59
     60  integer :: nQ
     61  parameter (nQ=7)
     62
     63
     64  !ym      character*6 nom(nQ)
     65  !ym      character*6 unites(nQ)
     66  character*6,save :: nom(nQ)
     67  character*6,save :: unites(nQ)
     68
     69  character(len=10) :: file
     70  integer :: ifile
     71  parameter (ifile=4)
     72
     73  integer :: itemp,igeop,iecin,iang,iu,iovap,iun
     74  integer :: i_sortie
     75
     76  save first,icum,ncum
     77  save itemp,igeop,iecin,iang,iu,iovap,iun
     78  save i_sortie
     79
     80  real :: time
     81  integer :: itau
     82  save time,itau
     83  data time,itau/0.,0/
     84
     85  data first/.true./
     86  data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
     87  data i_sortie/1/
     88
     89  real :: ww
     90
     91  !   variables dynamiques intermédiaires
     92  REAL :: vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
     93  REAL :: ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
     94  REAL :: massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
     95  REAL :: vorpot(iip1,jjm,llm)
     96  REAL :: w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
     97  REAL :: bern(iip1,jjp1,llm)
     98
     99  !   champ contenant les scalaires advectés.
     100  real :: Q(iip1,jjp1,llm,nQ)
     101
     102  !   champs cumulés
     103  real :: ps_cum(iip1,jjp1)
     104  real :: masse_cum(iip1,jjp1,llm)
     105  real :: flux_u_cum(iip1,jjp1,llm)
     106  real :: flux_v_cum(iip1,jjm,llm)
     107  real :: Q_cum(iip1,jjp1,llm,nQ)
     108  real :: flux_uQ_cum(iip1,jjp1,llm,nQ)
     109  real :: flux_vQ_cum(iip1,jjm,llm,nQ)
     110  real :: flux_wQ_cum(iip1,jjp1,llm,nQ)
     111  real :: dQ(iip1,jjp1,llm,nQ)
     112
     113  save ps_cum,masse_cum,flux_u_cum,flux_v_cum
     114  save Q_cum,flux_uQ_cum,flux_vQ_cum
     115
     116  !   champs de tansport en moyenne zonale
     117  integer :: ntr,itr
     118  parameter (ntr=5)
     119
     120  !ym      character*10 znom(ntr,nQ)
     121  !ym      character*20 znoml(ntr,nQ)
     122  !ym      character*10 zunites(ntr,nQ)
     123  character*10,save :: znom(ntr,nQ)
     124  character*20,save :: znoml(ntr,nQ)
     125  character*10,save :: zunites(ntr,nQ)
     126
     127  integer :: iave,itot,immc,itrs,istn
     128  data iave,itot,immc,itrs,istn/1,2,3,4,5/
     129  character(len=3) :: ctrs(ntr)
     130  data ctrs/'  ','TOT','MMC','TRS','STN'/
     131
     132  real :: zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
     133  real :: zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
     134  real :: zmasse(jjm,llm),zamasse(jjm)
     135
     136  real :: zv(jjm,llm),psi(jjm,llm+1)
     137
     138  integer :: i,j,l,iQ
     139
     140
     141  !   Initialisation du fichier contenant les moyennes zonales.
     142  !   ---------------------------------------------------------
     143
     144  character(len=10) :: infile
     145
     146  integer :: fileid
     147  integer :: thoriid, zvertiid
     148  save fileid
     149
     150  integer :: ndex3d(jjm*llm)
     151
     152  !   Variables locales
     153  !
     154  integer :: tau0
     155  real :: zjulian
     156  character(len=3) :: str
     157  character(len=10) :: ctrac
     158  integer :: ii,jj
     159  integer :: zan, dayref
     160  !
     161  real :: rlong(jjm),rlatg(jjm)
     162
     163
     164
     165  !=====================================================================
     166  !   Initialisation
     167  !=====================================================================
     168
     169  time=time+dt_app
     170  itau=itau+1
     171  !IM
     172  ndex3d=0
     173
     174  if (first) then
     175
     176
     177    icum=0
     178    ! initialisation des fichiers
     179    first=.false.
     180  !   ncum est la frequence de stokage en pas de temps
     181    ncum=dt_cum/dt_app
     182    if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
     183       WRITE(lunout,*) &
     184             'Pb : le pas de cumule doit etre multiple du pas'
     185       WRITE(lunout,*)'dt_app=',dt_app
     186       WRITE(lunout,*)'dt_cum=',dt_cum
     187       call abort_gcm('bilan_dyn','stopped',1)
     188    endif
     189
     190    if (i_sortie.eq.1) then
     191     file='dynzon'
     192     call inigrads(ifile,1 &
     193           ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi &
     194           ,llm,presnivs,1. &
     195           ,dt_cum,file,'dyn_zon ')
     196    endif
     197
     198    nom(itemp)='T'
     199    nom(igeop)='gz'
     200    nom(iecin)='K'
     201    nom(iang)='ang'
     202    nom(iu)='u'
     203    nom(iovap)='ovap'
     204    nom(iun)='un'
     205
     206    unites(itemp)='K'
     207    unites(igeop)='m2/s2'
     208    unites(iecin)='m2/s2'
     209    unites(iang)='ang'
     210    unites(iu)='m/s'
     211    unites(iovap)='kg/kg'
     212    unites(iun)='un'
     213
     214
     215  !   Initialisation du fichier contenant les moyennes zonales.
     216  !   ---------------------------------------------------------
     217
     218  infile='dynzon'
     219
     220  zan = annee_ref
     221  dayref = day_ref
     222  CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
     223  tau0 = itau_dyn
     224
     225  rlong=0.
     226  rlatg=rlatv*180./pi
     227
     228  call histbeg(infile, 1, rlong, jjm, rlatg, &
     229        1, 1, 1, jjm, &
     230        tau0, zjulian, dt_cum, thoriid, fileid)
     231
     232  !
     233  !  Appel a histvert pour la grille verticale
     234  !
     235  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
     236        llm, presnivs, zvertiid)
     237  !
     238  !  Appels a histdef pour la definition des variables a sauvegarder
     239  do iQ=1,nQ
     240     do itr=1,ntr
     241        if(itr.eq.1) then
     242           znom(itr,iQ)=nom(iQ)
     243           znoml(itr,iQ)=nom(iQ)
     244           zunites(itr,iQ)=unites(iQ)
     245        else
     246           znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
     247           znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
     248           zunites(itr,iQ)='m/s * '//unites(iQ)
    188249        endif
    189 
    190         if (i_sortie.eq.1) then
    191          file='dynzon'
    192          call inigrads(ifile,1
    193      s  ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
    194      s  ,llm,presnivs,1.
    195      s  ,dt_cum,file,'dyn_zon ')
    196         endif
    197 
    198         nom(itemp)='T'
    199         nom(igeop)='gz'
    200         nom(iecin)='K'
    201         nom(iang)='ang'
    202         nom(iu)='u'
    203         nom(iovap)='ovap'
    204         nom(iun)='un'
    205 
    206         unites(itemp)='K'
    207         unites(igeop)='m2/s2'
    208         unites(iecin)='m2/s2'
    209         unites(iang)='ang'
    210         unites(iu)='m/s'
    211         unites(iovap)='kg/kg'
    212         unites(iun)='un'
    213 
    214 
    215 c   Initialisation du fichier contenant les moyennes zonales.
    216 c   ---------------------------------------------------------
    217 
    218       infile='dynzon'
    219 
    220       zan = annee_ref
    221       dayref = day_ref
    222       CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
    223       tau0 = itau_dyn
    224      
    225       rlong=0.
    226       rlatg=rlatv*180./pi
    227        
    228       call histbeg(infile, 1, rlong, jjm, rlatg,
    229      .             1, 1, 1, jjm,
    230      .             tau0, zjulian, dt_cum, thoriid, fileid)
    231 
    232 C
    233 C  Appel a histvert pour la grille verticale
    234 C
    235       call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
    236      .              llm, presnivs, zvertiid)
    237 C
    238 C  Appels a histdef pour la definition des variables a sauvegarder
    239       do iQ=1,nQ
    240          do itr=1,ntr
    241             if(itr.eq.1) then
    242                znom(itr,iQ)=nom(iQ)
    243                znoml(itr,iQ)=nom(iQ)
    244                zunites(itr,iQ)=unites(iQ)
    245             else
    246                znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
    247                znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
    248                zunites(itr,iQ)='m/s * '//unites(iQ)
    249             endif
    250          enddo
    251       enddo
    252 
    253 c   Declarations des champs avec dimension verticale
    254 c      print*,'1HISTDEF'
    255       do iQ=1,nQ
    256          do itr=1,ntr
    257       IF (prt_level > 5)
    258      . WRITE(lunout,*)'var ',itr,iQ
    259      .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
    260             call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
    261      .        zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
    262      .        32,'ave(X)',dt_cum,dt_cum)
    263          enddo
    264 c   Declarations pour les fonctions de courant
    265 c      print*,'2HISTDEF'
    266           call histdef(fileid,'psi'//nom(iQ)
    267      .      ,'stream fn. '//znoml(itot,iQ),
    268      .      zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
    269      .      32,'ave(X)',dt_cum,dt_cum)
    270       enddo
    271 
    272 
    273 c   Declarations pour les champs de transport d'air
    274 c      print*,'3HISTDEF'
    275       call histdef(fileid, 'masse', 'masse',
    276      .             'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
    277      .             32, 'ave(X)', dt_cum, dt_cum)
    278       call histdef(fileid, 'v', 'v',
    279      .             'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
    280      .             32, 'ave(X)', dt_cum, dt_cum)
    281 c   Declarations pour les fonctions de courant
    282 c      print*,'4HISTDEF'
    283           call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
    284      .      1,jjm,thoriid,llm,1,llm,zvertiid,
    285      .      32,'ave(X)',dt_cum,dt_cum)
    286 
    287 
    288 c   Declaration des champs 1D de transport en latitude
    289 c      print*,'5HISTDEF'
    290       do iQ=1,nQ
    291          do itr=2,ntr
    292             call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
    293      .        zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
    294      .        32,'ave(X)',dt_cum,dt_cum)
    295          enddo
    296       enddo
    297 
    298 
    299 c      print*,'8HISTDEF'
    300                CALL histend(fileid)
    301 
    302 
    303       endif
    304 
    305 
    306 c=====================================================================
    307 c   Calcul des champs dynamiques
    308 c   ----------------------------
    309 
    310 c   énergie cinétique
    311       ucont(:,:,:)=0
    312       CALL covcont(llm,ucov,vcov,ucont,vcont)
    313       CALL enercin(vcov,ucov,vcont,ucont,ecin)
    314 
    315 c   moment cinétique
    316       do l=1,llm
    317          ang(:,:,l)=ucov(:,:,l)+constang(:,:)
    318          unat(:,:,l)=ucont(:,:,l)*cu(:,:)
    319       enddo
    320 
    321       Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
    322       Q(:,:,:,igeop)=phi(:,:,:)
    323       Q(:,:,:,iecin)=ecin(:,:,:)
    324       Q(:,:,:,iang)=ang(:,:,:)
    325       Q(:,:,:,iu)=unat(:,:,:)
    326       Q(:,:,:,iovap)=trac(:,:,:,1)
    327       Q(:,:,:,iun)=1.
    328 
    329 
    330 c=====================================================================
    331 c   Cumul
    332 c=====================================================================
    333 c
    334       if(icum.EQ.0) then
    335          ps_cum=0.
    336          masse_cum=0.
    337          flux_u_cum=0.
    338          flux_v_cum=0.
    339          Q_cum=0.
    340          flux_vQ_cum=0.
    341          flux_uQ_cum=0.
    342       endif
    343 
    344       IF (prt_level > 5)
    345      . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
    346       icum=icum+1
    347 
    348 c   accumulation des flux de masse horizontaux
    349       ps_cum=ps_cum+ps
    350       masse_cum=masse_cum+masse
    351       flux_u_cum=flux_u_cum+flux_u
    352       flux_v_cum=flux_v_cum+flux_v
    353       do iQ=1,nQ
    354       Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
    355       enddo
    356 
    357 c=====================================================================
    358 c  FLUX ET TENDANCES
    359 c=====================================================================
    360 
    361 c   Flux longitudinal
    362 c   -----------------
    363       do iQ=1,nQ
    364          do l=1,llm
    365             do j=1,jjp1
    366                do i=1,iim
    367                   flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
    368      s            +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
    369                enddo
    370                flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
    371             enddo
    372          enddo
    373       enddo
    374 
    375 c    flux méridien
    376 c    -------------
    377       do iQ=1,nQ
    378          do l=1,llm
    379             do j=1,jjm
    380                do i=1,iip1
    381                   flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
    382      s            +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
    383                enddo
    384             enddo
    385          enddo
    386       enddo
    387 
    388 
    389 c    tendances
    390 c    ---------
    391 
    392 c   convergence horizontale
    393       call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
    394 
    395 c   calcul de la vitesse verticale
    396       call convmas(flux_u_cum,flux_v_cum,convm)
    397       CALL vitvert(convm,w)
    398 
    399       do iQ=1,nQ
    400          do l=1,llm-1
    401             do j=1,jjp1
    402                do i=1,iip1
    403                   ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
    404                   dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
    405                   dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
    406                enddo
    407             enddo
    408          enddo
    409       enddo
    410       IF (prt_level > 5)
    411      . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
    412 c=====================================================================
    413 c   PAS DE TEMPS D'ECRITURE
    414 c=====================================================================
    415       if (icum.eq.ncum) then
    416 c=====================================================================
    417 
    418       IF (prt_level > 5)
    419      . WRITE(lunout,*)'Pas d ecriture'
    420 
    421 c   Normalisation
    422       do iQ=1,nQ
    423          Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
    424       enddo
    425       zz=1./REAL(ncum)
    426       ps_cum=ps_cum*zz
    427       masse_cum=masse_cum*zz
    428       flux_u_cum=flux_u_cum*zz
    429       flux_v_cum=flux_v_cum*zz
    430       flux_uQ_cum=flux_uQ_cum*zz
    431       flux_vQ_cum=flux_vQ_cum*zz
    432       dQ=dQ*zz
    433 
    434 
    435 c   A retravailler eventuellement
    436 c   division de dQ par la masse pour revenir aux bonnes grandeurs
    437       do iQ=1,nQ
    438          dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
    439       enddo
    440  
    441 c=====================================================================
    442 c   Transport méridien
    443 c=====================================================================
    444 
    445 c   cumul zonal des masses des mailles
    446 c   ----------------------------------
    447       zv=0.
    448       zmasse=0.
    449       call massbar(masse_cum,massebx,masseby)
    450       do l=1,llm
    451          do j=1,jjm
    452             do i=1,iim
    453                zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
    454                zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
    455             enddo
    456             zfactv(j,l)=cv(1,j)/zmasse(j,l)
    457          enddo
    458       enddo
    459 
    460 c     print*,'3OK'
    461 c   --------------------------------------------------------------
    462 c   calcul de la moyenne zonale du transport :
    463 c   ------------------------------------------
    464 c
    465 c                                     --
    466 c TOT : la circulation totale       [ vq ]
    467 c
    468 c                                      -     -
    469 c MMC : mean meridional circulation [ v ] [ q ]
    470 c
    471 c                                     ----      --       - -
    472 c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
    473 c
    474 c                                     - * - *       - -       -     -
    475 c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
    476 c
    477 c                                              - -
    478 c    on utilise aussi l'intermediaire TMP :  [ v q ]
    479 c
    480 c    la variable zfactv transforme un transport meridien cumule
    481 c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
    482 c
    483 c   --------------------------------------------------------------
    484 
    485 
    486 c   ----------------------------------------
    487 c   Transport dans le plan latitude-altitude
    488 c   ----------------------------------------
    489 
    490       zvQ=0.
    491       psiQ=0.
    492       do iQ=1,nQ
    493          zvQtmp=0.
    494          do l=1,llm
    495             do j=1,jjm
    496 c              print*,'j,l,iQ=',j,l,iQ
    497 c   Calcul des moyennes zonales du transort total et de zvQtmp
    498                do i=1,iim
    499                   zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
    500      s                            +flux_vQ_cum(i,j,l,iQ)
    501                   zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
    502      s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
    503                   zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
    504      s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
    505                   zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
    506                enddo
    507 c              print*,'aOK'
    508 c   Decomposition
    509                zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
    510                zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
    511                zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
    512                zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
    513                zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
    514                zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
    515             enddo
    516          enddo
    517 c   fonction de courant meridienne pour la quantite Q
    518          do l=llm,1,-1
    519             do j=1,jjm
    520                psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
    521             enddo
    522          enddo
    523       enddo
    524 
    525 c   fonction de courant pour la circulation meridienne moyenne
    526       psi=0.
    527       do l=llm,1,-1
    528          do j=1,jjm
    529             psi(j,l)=psi(j,l+1)+zv(j,l)
    530             zv(j,l)=zv(j,l)*zfactv(j,l)
    531          enddo
    532       enddo
    533 
    534 c     print*,'4OK'
    535 c   sorties proprement dites
    536       if (i_sortie.eq.1) then
    537       do iQ=1,nQ
    538          do itr=1,ntr
    539             call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ)
    540      s      ,jjm*llm,ndex3d)
    541          enddo
    542          call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ)
    543      s      ,jjm*llm,ndex3d)
    544       enddo
    545 
    546       call histwrite(fileid,'masse',itau,zmasse
    547      s   ,jjm*llm,ndex3d)
    548       call histwrite(fileid,'v',itau,zv
    549      s   ,jjm*llm,ndex3d)
    550       psi=psi*1.e-9
    551       call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
    552 
    553       endif
    554 
    555 
    556 c   -----------------
    557 c   Moyenne verticale
    558 c   -----------------
    559 
    560       zamasse=0.
    561       do l=1,llm
    562          zamasse(:)=zamasse(:)+zmasse(:,l)
    563       enddo
    564       zavQ=0.
    565       do iQ=1,nQ
    566          do itr=2,ntr
    567             do l=1,llm
    568                zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
    569             enddo
    570             zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
    571             call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ)
    572      s      ,jjm*llm,ndex3d)
    573          enddo
    574       enddo
    575 
    576 c     on doit pouvoir tracer systematiquement la fonction de courant.
    577 
    578 c=====================================================================
    579 c/////////////////////////////////////////////////////////////////////
    580       icum=0                  !///////////////////////////////////////
    581       endif ! icum.eq.ncum    !///////////////////////////////////////
    582 c/////////////////////////////////////////////////////////////////////
    583 c=====================================================================
    584 
    585       return
    586       end
     250     enddo
     251  enddo
     252
     253  !   Declarations des champs avec dimension verticale
     254   ! print*,'1HISTDEF'
     255  do iQ=1,nQ
     256     do itr=1,ntr
     257  IF (prt_level > 5) &
     258        WRITE(lunout,*)'var ',itr,iQ &
     259        ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
     260        call histdef(fileid,znom(itr,iQ),znoml(itr,iQ), &
     261              zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid, &
     262              32,'ave(X)',dt_cum,dt_cum)
     263     enddo
     264  !   Declarations pour les fonctions de courant
     265   ! print*,'2HISTDEF'
     266      call histdef(fileid,'psi'//nom(iQ) &
     267            ,'stream fn. '//znoml(itot,iQ), &
     268            zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid, &
     269            32,'ave(X)',dt_cum,dt_cum)
     270  enddo
     271
     272
     273  !   Declarations pour les champs de transport d'air
     274   ! print*,'3HISTDEF'
     275  call histdef(fileid, 'masse', 'masse', &
     276        'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
     277        32, 'ave(X)', dt_cum, dt_cum)
     278  call histdef(fileid, 'v', 'v', &
     279        'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
     280        32, 'ave(X)', dt_cum, dt_cum)
     281  !   Declarations pour les fonctions de courant
     282   ! print*,'4HISTDEF'
     283      call histdef(fileid,'psi','stream fn. MMC ','mega t/s', &
     284            1,jjm,thoriid,llm,1,llm,zvertiid, &
     285            32,'ave(X)',dt_cum,dt_cum)
     286
     287
     288  !   Declaration des champs 1D de transport en latitude
     289   ! print*,'5HISTDEF'
     290  do iQ=1,nQ
     291     do itr=2,ntr
     292        call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ), &
     293              zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99, &
     294              32,'ave(X)',dt_cum,dt_cum)
     295     enddo
     296  enddo
     297
     298
     299   ! print*,'8HISTDEF'
     300           CALL histend(fileid)
     301
     302
     303  endif
     304
     305
     306  !=====================================================================
     307  !   Calcul des champs dynamiques
     308  !   ----------------------------
     309
     310  !   énergie cinétique
     311  ucont(:,:,:)=0
     312  CALL covcont(llm,ucov,vcov,ucont,vcont)
     313  CALL enercin(vcov,ucov,vcont,ucont,ecin)
     314
     315  !   moment cinétique
     316  do l=1,llm
     317     ang(:,:,l)=ucov(:,:,l)+constang(:,:)
     318     unat(:,:,l)=ucont(:,:,l)*cu(:,:)
     319  enddo
     320
     321  Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
     322  Q(:,:,:,igeop)=phi(:,:,:)
     323  Q(:,:,:,iecin)=ecin(:,:,:)
     324  Q(:,:,:,iang)=ang(:,:,:)
     325  Q(:,:,:,iu)=unat(:,:,:)
     326  Q(:,:,:,iovap)=trac(:,:,:,1)
     327  Q(:,:,:,iun)=1.
     328
     329
     330  !=====================================================================
     331  !   Cumul
     332  !=====================================================================
     333  !
     334  if(icum.EQ.0) then
     335     ps_cum=0.
     336     masse_cum=0.
     337     flux_u_cum=0.
     338     flux_v_cum=0.
     339     Q_cum=0.
     340     flux_vQ_cum=0.
     341     flux_uQ_cum=0.
     342  endif
     343
     344  IF (prt_level > 5) &
     345        WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
     346  icum=icum+1
     347
     348  !   accumulation des flux de masse horizontaux
     349  ps_cum=ps_cum+ps
     350  masse_cum=masse_cum+masse
     351  flux_u_cum=flux_u_cum+flux_u
     352  flux_v_cum=flux_v_cum+flux_v
     353  do iQ=1,nQ
     354  Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
     355  enddo
     356
     357  !=====================================================================
     358  !  FLUX ET TENDANCES
     359  !=====================================================================
     360
     361  !   Flux longitudinal
     362  !   -----------------
     363  do iQ=1,nQ
     364     do l=1,llm
     365        do j=1,jjp1
     366           do i=1,iim
     367              flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ) &
     368                    +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
     369           enddo
     370           flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
     371        enddo
     372     enddo
     373  enddo
     374
     375  !    flux méridien
     376  !    -------------
     377  do iQ=1,nQ
     378     do l=1,llm
     379        do j=1,jjm
     380           do i=1,iip1
     381              flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ) &
     382                    +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
     383           enddo
     384        enddo
     385     enddo
     386  enddo
     387
     388
     389  !    tendances
     390  !    ---------
     391
     392  !   convergence horizontale
     393  call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
     394
     395  !   calcul de la vitesse verticale
     396  call convmas(flux_u_cum,flux_v_cum,convm)
     397  CALL vitvert(convm,w)
     398
     399  do iQ=1,nQ
     400     do l=1,llm-1
     401        do j=1,jjp1
     402           do i=1,iip1
     403              ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
     404              dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
     405              dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
     406           enddo
     407        enddo
     408     enddo
     409  enddo
     410  IF (prt_level > 5) &
     411        WRITE(lunout,*)'Apres les calculs fait a chaque pas'
     412  !=====================================================================
     413  !   PAS DE TEMPS D'ECRITURE
     414  !=====================================================================
     415  if (icum.eq.ncum) then
     416  !=====================================================================
     417
     418  IF (prt_level > 5) &
     419        WRITE(lunout,*)'Pas d ecriture'
     420
     421  !   Normalisation
     422  do iQ=1,nQ
     423     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
     424  enddo
     425  zz=1./REAL(ncum)
     426  ps_cum=ps_cum*zz
     427  masse_cum=masse_cum*zz
     428  flux_u_cum=flux_u_cum*zz
     429  flux_v_cum=flux_v_cum*zz
     430  flux_uQ_cum=flux_uQ_cum*zz
     431  flux_vQ_cum=flux_vQ_cum*zz
     432  dQ=dQ*zz
     433
     434
     435  !   A retravailler eventuellement
     436  !   division de dQ par la masse pour revenir aux bonnes grandeurs
     437  do iQ=1,nQ
     438     dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
     439  enddo
     440
     441  !=====================================================================
     442  !   Transport méridien
     443  !=====================================================================
     444
     445  !   cumul zonal des masses des mailles
     446  !   ----------------------------------
     447  zv=0.
     448  zmasse=0.
     449  call massbar(masse_cum,massebx,masseby)
     450  do l=1,llm
     451     do j=1,jjm
     452        do i=1,iim
     453           zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
     454           zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
     455        enddo
     456        zfactv(j,l)=cv(1,j)/zmasse(j,l)
     457     enddo
     458  enddo
     459
     460  ! print*,'3OK'
     461  !   --------------------------------------------------------------
     462  !   calcul de la moyenne zonale du transport :
     463  !   ------------------------------------------
     464  !
     465  !                                 --
     466  ! TOT : la circulation totale       [ vq ]
     467  !
     468  !                                  -     -
     469  ! MMC : mean meridional circulation [ v ] [ q ]
     470  !
     471  !                                 ----      --       - -
     472  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
     473  !
     474  !                                 - * - *       - -       -     -
     475  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
     476  !
     477  !                                          - -
     478  !    on utilise aussi l'intermediaire TMP :  [ v q ]
     479  !
     480  !    la variable zfactv transforme un transport meridien cumule
     481  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
     482  !
     483  !   --------------------------------------------------------------
     484
     485
     486  !   ----------------------------------------
     487  !   Transport dans le plan latitude-altitude
     488  !   ----------------------------------------
     489
     490  zvQ=0.
     491  psiQ=0.
     492  do iQ=1,nQ
     493     zvQtmp=0.
     494     do l=1,llm
     495        do j=1,jjm
     496           ! print*,'j,l,iQ=',j,l,iQ
     497  !   Calcul des moyennes zonales du transort total et de zvQtmp
     498           do i=1,iim
     499              zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ) &
     500                    +flux_vQ_cum(i,j,l,iQ)
     501              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
     502                    Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
     503              zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy &
     504                    /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
     505              zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
     506           enddo
     507           ! print*,'aOK'
     508  !   Decomposition
     509           zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
     510           zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
     511           zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
     512           zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
     513           zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
     514           zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
     515        enddo
     516     enddo
     517  !   fonction de courant meridienne pour la quantite Q
     518     do l=llm,1,-1
     519        do j=1,jjm
     520           psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
     521        enddo
     522     enddo
     523  enddo
     524
     525  !   fonction de courant pour la circulation meridienne moyenne
     526  psi=0.
     527  do l=llm,1,-1
     528     do j=1,jjm
     529        psi(j,l)=psi(j,l+1)+zv(j,l)
     530        zv(j,l)=zv(j,l)*zfactv(j,l)
     531     enddo
     532  enddo
     533
     534  ! print*,'4OK'
     535  !   sorties proprement dites
     536  if (i_sortie.eq.1) then
     537  do iQ=1,nQ
     538     do itr=1,ntr
     539        call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ) &
     540              ,jjm*llm,ndex3d)
     541     enddo
     542     call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ) &
     543           ,jjm*llm,ndex3d)
     544  enddo
     545
     546  call histwrite(fileid,'masse',itau,zmasse &
     547        ,jjm*llm,ndex3d)
     548  call histwrite(fileid,'v',itau,zv &
     549        ,jjm*llm,ndex3d)
     550  psi=psi*1.e-9
     551  call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
     552
     553  endif
     554
     555
     556  !   -----------------
     557  !   Moyenne verticale
     558  !   -----------------
     559
     560  zamasse=0.
     561  do l=1,llm
     562     zamasse(:)=zamasse(:)+zmasse(:,l)
     563  enddo
     564  zavQ=0.
     565  do iQ=1,nQ
     566     do itr=2,ntr
     567        do l=1,llm
     568           zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
     569        enddo
     570        zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
     571        call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ) &
     572              ,jjm*llm,ndex3d)
     573     enddo
     574  enddo
     575
     576  ! on doit pouvoir tracer systematiquement la fonction de courant.
     577
     578  !=====================================================================
     579  !/////////////////////////////////////////////////////////////////////
     580  icum=0                  !///////////////////////////////////////
     581  endif ! icum.eq.ncum    !///////////////////////////////////////
     582  !/////////////////////////////////////////////////////////////////////
     583  !=====================================================================
     584
     585  return
     586end subroutine bilan_dyn
Note: See TracChangeset for help on using the changeset viewer.