Ignore:
Timestamp:
Aug 18, 2011, 12:09:27 PM (13 years ago)
Author:
emillour
Message:

Ehouarn: Mise a jour des dynamiques (seq et ) pour suivre
les changements et developpements de LMDZ5 terrestre
(mise a niveau avec LMDZ5 trunk, rev 1560). Ce qui ne devrais pas changer grand chose au fonctionnement de base du GCM).
Modification importante: correction du bug sur le cumul des flux de masse pour le transport des traceurs (mais dans les faits semble avoir peu d'impact).
(voir "commit_importants.log" pour les details des ajouts et modifications).

File:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90

    r269 r270  
    11!
    2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
    33!
    4       SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
    5      *                       tetagdiv,tetagrot,tetatemp             )
    6 c=======================================================================
    7 c   initialisation de la dissipation horizontale
    8 c=======================================================================
    9 c-----------------------------------------------------------------------
    10 c   declarations:
    11 c   -------------
    12 
    13       USE control_mod
    14 
    15       IMPLICIT NONE
    16 #include "dimensions.h"
    17 #include "paramet.h"
    18 #include "comdissipn.h"
    19 #include "comconst.h"
    20 #include "comvert.h"
    21 #include "logic.h"
    22 
    23       LOGICAL lstardis
    24       INTEGER nitergdiv,nitergrot,niterh
    25       REAL    tetagdiv,tetagrot,tetatemp
    26       REAL fact,zvert(llm),zz
    27       REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
    28       REAL ullm,vllm,umin,vmin,zhmin,zhmax
    29       REAL zllm,z1llm
    30 
    31       INTEGER l,ij,idum,ii
    32       REAL tetamin
    33       REAL Pup
    34 
    35       REAL ran1
    36 
    37 
    38 c-----------------------------------------------------------------------
    39 c
    40 c   calcul des valeurs propres des operateurs par methode iterrative:
    41 c   -----------------------------------------------------------------
    42 
    43       crot     = -1.
    44       cdivu    = -1.
    45       cdivh    = -1.
    46 
    47 c   calcul de la valeur propre de divgrad:
    48 c   --------------------------------------
    49       idum = 0
    50       DO l = 1, llm
    51        DO ij = 1, ip1jmp1
     4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
     5     tetagdiv,tetagrot,tetatemp             )
     6  !=======================================================================
     7  !   initialisation de la dissipation horizontale
     8  !=======================================================================
     9  !-----------------------------------------------------------------------
     10  !   declarations:
     11  !   -------------
     12
     13  USE control_mod, only : dissip_period,iperiod
     14
     15  IMPLICIT NONE
     16  include "dimensions.h"
     17  include "paramet.h"
     18  include "comdissipn.h"
     19  include "comconst.h"
     20  include "comvert.h"
     21  include "logic.h"
     22  include "iniprint.h"
     23
     24  LOGICAL,INTENT(in) :: lstardis
     25  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
     26  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
     27
     28! Local variables:
     29  REAL fact,zvert(llm),zz
     30  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
     31  REAL ullm,vllm,umin,vmin,zhmin,zhmax
     32  REAL zllm,z1llm
     33
     34  INTEGER l,ij,idum,ii
     35  REAL tetamin
     36  REAL Pup
     37  character (len=80) :: abort_message
     38
     39  REAL ran1
     40
     41
     42  !-----------------------------------------------------------------------
     43  !
     44  !   calcul des valeurs propres des operateurs par methode iterrative:
     45  !   -----------------------------------------------------------------
     46
     47  crot     = -1.
     48  cdivu    = -1.
     49  cdivh    = -1.
     50
     51  !   calcul de la valeur propre de divgrad:
     52  !   --------------------------------------
     53  idum = 0
     54  DO l = 1, llm
     55     DO ij = 1, ip1jmp1
    5256        deltap(ij,l) = 1.
    53        ENDDO
    54       ENDDO
    55 
    56       idum  = -1
    57       zh(1) = RAN1(idum)-.5
    58       idum  = 0
    59       DO ij = 2, ip1jmp1
    60         zh(ij) = RAN1(idum) -.5
    61       ENDDO
    62 
    63       CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
    64 
    65       CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    66 
    67       IF ( zhmin .GE. zhmax  )     THEN
    68          PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
    69          STOP'probleme generateur alleatoire dans inidissip'
    70       ENDIF
    71 
    72       zllm = ABS( zhmax )
    73       DO l = 1,50
    74          IF(lstardis) THEN
    75             CALL divgrad2(1,zh,deltap,niterh,zh)
    76          ELSE
    77             CALL divgrad (1,zh,niterh,zh)
    78          ENDIF
    79 
    80         CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    81 
    82          zllm  = ABS( zhmax )
    83          z1llm = 1./zllm
    84          DO ij = 1,ip1jmp1
    85             zh(ij) = zh(ij)* z1llm
    86          ENDDO
    87       ENDDO
    88 
    89       IF(lstardis) THEN
    90          cdivh = 1./ zllm
    91       ELSE
    92          cdivh = zllm ** ( -1./niterh )
    93       ENDIF
    94 
    95 c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
    96 c   -----------------------------------------------------------------
    97       print*,'calcul des valeurs propres'
    98 
    99       DO  20  ii = 1, 2
    100 c
    101          DO ij = 1, ip1jmp1
    102            zu(ij)  = RAN1(idum) -.5
    103          ENDDO
    104          CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
    105          DO ij = 1, ip1jm
    106             zv(ij) = RAN1(idum) -.5
    107          ENDDO
    108          CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
    109 
    110          CALL minmax(iip1*jjp1,zu,umin,ullm )
    111          CALL minmax(iip1*jjm, zv,vmin,vllm )
    112 
    113          ullm = ABS ( ullm )
    114          vllm = ABS ( vllm )
    115 
    116          DO  5  l = 1, 50
    117             IF(ii.EQ.1) THEN
    118 ccccc             CALL covcont( 1,zu,zv,zu,zv )
    119                IF(lstardis) THEN
    120                   CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
    121                ELSE
    122                   CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
    123                ENDIF
    124             ELSE
    125                IF(lstardis) THEN
    126                   CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
    127                ELSE
    128                   CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
    129                ENDIF
    130             ENDIF
    131 
    132             CALL minmax(iip1*jjp1,zu,umin,ullm )
    133             CALL minmax(iip1*jjm, zv,vmin,vllm )
    134 
    135             ullm = ABS  ( ullm )
    136             vllm = ABS  ( vllm )
    137 
    138             zllm  = MAX( ullm,vllm )
    139             z1llm = 1./ zllm
    140             DO ij = 1, ip1jmp1
    141               zu(ij) = zu(ij)* z1llm
    142             ENDDO
    143             DO ij = 1, ip1jm
    144                zv(ij) = zv(ij)* z1llm
    145             ENDDO
    146  5       CONTINUE
    147 
    148          IF ( ii.EQ.1 ) THEN
    149             IF(lstardis) THEN
    150                cdivu  = 1./zllm
    151             ELSE
    152                cdivu  = zllm **( -1./nitergdiv )
    153             ENDIF
    154          ELSE
    155             IF(lstardis) THEN
    156                crot   = 1./ zllm
    157             ELSE
    158                crot   = zllm **( -1./nitergrot )
    159             ENDIF
    160          ENDIF
    161 
    162  20   CONTINUE
    163 
    164 c   petit test pour les operateurs non star:
    165 c   ----------------------------------------
    166 
    167 c     IF(.NOT.lstardis) THEN
    168          fact    = rad*24./REAL(jjm)
    169          fact    = fact*fact
    170          PRINT*,'coef u ', fact/cdivu, 1./cdivu
    171          PRINT*,'coef r ', fact/crot , 1./crot
    172          PRINT*,'coef h ', fact/cdivh, 1./cdivh
    173 c     ENDIF
    174 
    175 c-----------------------------------------------------------------------
    176 c   variation verticale du coefficient de dissipation:
    177 c   --------------------------------------------------
    178 
    179 c First step: going from 1 to dissip_fac_mid (in gcm.def)
    180 c============
    181       DO l=1,llm
    182          zz      = 1. - preff/presnivs(l)
    183          zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
    184       ENDDO
    185 
    186          write(*,*) 'Dissipation : '
    187          write(*,*) 'Multiplication de la dissipation en altitude :'
    188          write(*,*) '  dissip_fac_mid =', dissip_fac_mid
    189 
    190 c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
    191 c==========================
    192 c Utilisation de la fonction tangente hyperbolique pour augmenter
    193 c arbitrairement la dissipation et donc la stabilite du modele a
    194 c partir d'une certaine altitude.
    195 
    196 c   Le facteur multiplicatif de basse atmosphere etant deja pris
    197 c   en compte, il faut diviser le facteur multiplicatif de haute
    198 c   atmosphere par celui-ci.
    199 
    200       if (ok_strato) then
    201 
    202        Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
    203        do l=1,llm
    204          zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
    205      &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
    206      &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
    207        enddo
    208 
    209          write(*,*) '  dissip_fac_up =', dissip_fac_up
    210          write(*,*) 'Transition mid /up:  Pupstart,delta =',
    211      &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
    212 
    213       endif
    214 
    215 
    216       PRINT*,'Constantes de temps de la diffusion horizontale'
    217 
    218       tetamin =  1.e+6
    219 
    220       DO l=1,llm
    221         tetaudiv(l)   = zvert(l)/tetagdiv
    222         tetaurot(l)   = zvert(l)/tetagrot
    223         tetah(l)      = zvert(l)/tetatemp
    224 
    225         IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
    226         IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
    227         IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
    228       ENDDO
    229 
    230       PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
    231       idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
    232       PRINT *,' INIDI tetamin idissip ',tetamin,idissip
    233       idissip = MAX(iperiod,idissip)
    234       dtdiss  = idissip * dtvr
    235       PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
    236 
    237       DO l = 1,llm
    238          PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
    239      *                   dtdiss*tetah(l)
    240       ENDDO
    241 
    242 c
    243       RETURN
    244       END
     57     ENDDO
     58  ENDDO
     59
     60  idum  = -1
     61  zh(1) = RAN1(idum)-.5
     62  idum  = 0
     63  DO ij = 2, ip1jmp1
     64     zh(ij) = RAN1(idum) -.5
     65  ENDDO
     66
     67  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
     68
     69  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     70
     71  IF ( zhmin .GE. zhmax  )     THEN
     72     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
     73     abort_message='probleme generateur alleatoire dans inidissip'
     74     call abort_gcm('inidissip',abort_message,1)
     75  ENDIF
     76
     77  zllm = ABS( zhmax )
     78  DO l = 1,50
     79     IF(lstardis) THEN
     80        CALL divgrad2(1,zh,deltap,niterh,zh)
     81     ELSE
     82        CALL divgrad (1,zh,niterh,zh)
     83     ENDIF
     84
     85     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     86
     87     zllm  = ABS( zhmax )
     88     z1llm = 1./zllm
     89     DO ij = 1,ip1jmp1
     90        zh(ij) = zh(ij)* z1llm
     91     ENDDO
     92  ENDDO
     93
     94  IF(lstardis) THEN
     95     cdivh = 1./ zllm
     96  ELSE
     97     cdivh = zllm ** ( -1./niterh )
     98  ENDIF
     99
     100  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
     101  !   -----------------------------------------------------------------
     102  write(lunout,*)'inidissip: calcul des valeurs propres'
     103
     104  DO    ii = 1, 2
     105     !
     106     DO ij = 1, ip1jmp1
     107        zu(ij)  = RAN1(idum) -.5
     108     ENDDO
     109     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
     110     DO ij = 1, ip1jm
     111        zv(ij) = RAN1(idum) -.5
     112     ENDDO
     113     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
     114
     115     CALL minmax(iip1*jjp1,zu,umin,ullm )
     116     CALL minmax(iip1*jjm, zv,vmin,vllm )
     117
     118     ullm = ABS ( ullm )
     119     vllm = ABS ( vllm )
     120
     121     DO    l = 1, 50
     122        IF(ii.EQ.1) THEN
     123           !cccc             CALL covcont( 1,zu,zv,zu,zv )
     124           IF(lstardis) THEN
     125              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
     126           ELSE
     127              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
     128           ENDIF
     129        ELSE
     130           IF(lstardis) THEN
     131              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
     132           ELSE
     133              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
     134           ENDIF
     135        ENDIF
     136
     137        CALL minmax(iip1*jjp1,zu,umin,ullm )
     138        CALL minmax(iip1*jjm, zv,vmin,vllm )
     139
     140        ullm = ABS  ( ullm )
     141        vllm = ABS  ( vllm )
     142
     143        zllm  = MAX( ullm,vllm )
     144        z1llm = 1./ zllm
     145        DO ij = 1, ip1jmp1
     146           zu(ij) = zu(ij)* z1llm
     147        ENDDO
     148        DO ij = 1, ip1jm
     149           zv(ij) = zv(ij)* z1llm
     150        ENDDO
     151     end DO
     152
     153     IF ( ii.EQ.1 ) THEN
     154        IF(lstardis) THEN
     155           cdivu  = 1./zllm
     156        ELSE
     157           cdivu  = zllm **( -1./nitergdiv )
     158        ENDIF
     159     ELSE
     160        IF(lstardis) THEN
     161           crot   = 1./ zllm
     162        ELSE
     163           crot   = zllm **( -1./nitergrot )
     164        ENDIF
     165     ENDIF
     166
     167  end DO
     168
     169  !   petit test pour les operateurs non star:
     170  !   ----------------------------------------
     171
     172  !     IF(.NOT.lstardis) THEN
     173  fact    = rad*24./REAL(jjm)
     174  fact    = fact*fact
     175  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
     176  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
     177  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
     178  !     ENDIF
     179
     180  !-----------------------------------------------------------------------
     181  !   variation verticale du coefficient de dissipation:
     182  !   --------------------------------------------------
     183
     184! First step: going from 1 to dissip_fac_mid (in gcm.def)
     185!============
     186  DO l=1,llm
     187     zz      = 1. - preff/presnivs(l)
     188     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     189  ENDDO
     190
     191  write(lunout,*) 'Dissipation : '
     192  write(lunout,*) 'Multiplication de la dissipation en altitude :'
     193  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
     194
     195! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
     196!==========================
     197! Utilisation de la fonction tangente hyperbolique pour augmenter
     198! arbitrairement la dissipation et donc la stabilite du modele a
     199! partir d'une certaine altitude.
     200
     201!   Le facteur multiplicatif de basse atmosphere etant deja pris
     202!   en compte, il faut diviser le facteur multiplicatif de haute
     203!   atmosphere par celui-ci.
     204
     205  if (ok_strato) then
     206
     207    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     208    do l=1,llm
     209      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
     210                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
     211               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
     212    enddo
     213
     214    write(*,*) '  dissip_fac_up =', dissip_fac_up
     215    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
     216                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
     217
     218  endif
     219
     220
     221  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
     222
     223  tetamin =  1.e+6
     224
     225  DO l=1,llm
     226     tetaudiv(l)   = zvert(l)/tetagdiv
     227     tetaurot(l)   = zvert(l)/tetagrot
     228     tetah(l)      = zvert(l)/tetatemp
     229
     230     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
     231     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
     232     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
     233  ENDDO
     234
     235  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
     236  IF (dissip_period == 0) THEN
     237     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
     238     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
     239     dissip_period = MAX(iperiod,dissip_period)
     240  END IF
     241
     242  dtdiss  = dissip_period * dtvr
     243  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
     244
     245  DO l = 1,llm
     246     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
     247          dtdiss*tetah(l)
     248  ENDDO
     249
     250END SUBROUTINE inidissip
Note: See TracChangeset for help on using the changeset viewer.