Ignore:
Timestamp:
May 23, 2011, 1:37:09 PM (13 years ago)
Author:
Ehouarn Millour
Message:

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Location:
LMDZ5/trunk/libf/dyn3d
Files:
2 added
12 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/comvert.h

    r1279 r1520  
    55!   INCLUDE 'comvert.h'
    66
    7       COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
    8      &               pa,preff,nivsigs(llm),nivsig(llm+1)
     7      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
     8     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
     9     &               aps(llm),bps(llm),scaleheight
    910
    10       REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
     11      common/comverti/disvert_type
     12
     13      real ap     ! hybrid pressure contribution at interlayers
     14      real bp     ! hybrid sigma contribution at interlayer
     15      real presnivs ! (reference) pressure at mid-layers
     16      real dpres
     17      real pa     ! reference pressure (Pa) at which hybrid coordinates
     18                  ! become purely pressure
     19      real preff  ! reference surface pressure (Pa)
     20      real nivsigs
     21      real nivsig
     22      real aps    ! hybrid pressure contribution at mid-layers
     23      real bps    ! hybrid sigma contribution at mid-layers
     24      real scaleheight ! atmospheric (reference) scale height (km)
     25
     26      integer disvert_type ! type of vertical discretization:
     27                           ! 1: Earth (default for planet_type==earth),
     28                           !     automatic generation
     29                           ! 2: Planets (default for planet_type!=earth),
     30                           !     using 'z2sig.def' (or 'esasig.def) file
    1131
    1232 !-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3d/disvert.F90

    r1492 r1520  
    11! $Id$
    22
    3 SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
     3SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
    44
    55  ! Auteur : P. Le Van
    … …  
    1717  ! ds(l) : distance entre les couches l et l-1 en coord.s
    1818
    19   REAL pa, preff
    20   REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
    21   REAL presnivs(llm)
     19  real,intent(in) :: pa, preff
     20  real,intent(out) :: ap(llmp1), bp(llmp1)
     21  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
     22  real,intent(out) :: presnivs(llm)
     23  real,intent(out) :: scaleheight
    2224
    2325  REAL sig(llm+1), dsig(llm)
    … …  
    2628  INTEGER l
    2729  REAL dsigmin
    28   REAL alpha, beta, deltaz, h
     30  REAL alpha, beta, deltaz,scaleheight
    2931  INTEGER iostat
    30   REAL pi, x
     32  REAL x
     33  character(len=*),parameter :: modname="disvert"
    3134
    3235  !-----------------------------------------------------------------------
    3336
    34   pi = 2 * ASIN(1.)
     37  ! default scaleheight is 8km for earth
     38  scaleheight=8.
    3539
    3640  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
    … …  
    3842  IF (iostat == 0) THEN
    3943     ! cas 1 on lit les options dans sigma.def:
    40      READ(99, *) h ! hauteur d'echelle 8.
     44     READ(99, *) scaleheight ! hauteur d'echelle 8.
    4145     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
    4246     READ(99, *) beta ! facteur d'acroissement en haut 1.3
    … …  
    4448     READ(99, *) k1 ! nombre de couches dans la transition haute
    4549     CLOSE(99)
    46      alpha=deltaz/(llm*h)
    47      write(lunout, *)'h, alpha, k0, k1, beta'
     50     alpha=deltaz/(llm*scaleheight)
     51     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
     52                               scaleheight, alpha, k0, k1, beta
    4853
    4954     alpha=deltaz/tanh(1./k0)*2.
    … …  
    5156     sig(1)=1.
    5257     do l=1, llm
    53         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
    54              *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    55         zk=-h*log(sig(l+1))
     58        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
     59             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
     60                  *beta**(l-(llm-k1))/log(beta))
     61        zk=-scaleheight*log(sig(l+1))
    5662
    5763        dzk1=alpha*tanh(l/k0)
    … …  
    7379           dsigmin=1.
    7480        else
    75            WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
     81           write(lunout,*) trim(modname), &
     82           ' ATTENTION discretisation z a ajuster'
    7683           dsigmin=1.
    7784        endif
    78         WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
     85        write(lunout,*) trim(modname), &
     86        ' Discretisation verticale DSIGMIN=',dsigmin
    7987     endif
    8088
    … …  
    119127  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    120128
    121   write(lunout, *)' BP '
     129  write(lunout, *)  trim(modname),': BP '
    122130  write(lunout, *) bp
    123   write(lunout, *)' AP '
     131  write(lunout, *)  trim(modname),': AP '
    124132  write(lunout, *) ap
    125133
    126134  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
    127135  write(lunout, *)'couches calcules pour une pression de surface =', preff
    128   write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
    129   write(lunout, *)'8km'
     136  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
     137  write(lunout, *) scaleheight,' km'
    130138  DO l = 1, llm
    131139     dpres(l) = bp(l) - bp(l+1)
    132140     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    133141     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
    134           log(preff/presnivs(l))*8. &
    135           , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     142          log(preff/presnivs(l))*scaleheight &
     143          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
    136144          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
    137145  ENDDO
    138146
    139   write(lunout, *) 'PRESNIVS '
     147  write(lunout, *) trim(modname),': PRESNIVS '
    140148  write(lunout, *) presnivs
    141149

  • LMDZ5/trunk/libf/dyn3d/etat0_netcdf.F90

    r1511 r1520  
    251251!*******************************************************************************
    252252  CALL pression(ip1jmp1, ap, bp, psol, p3d)
    253   CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
     253  if (disvert_type.eq.1) then
     254    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
     255  else ! we assume that we are in the disvert_type==2 case
     256    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
     257  endif
    254258  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
    255259!  WRITE(lunout,*) 'P3D :', p3d(10,20,:)

  • LMDZ5/trunk/libf/dyn3d/exner_hyb.F

    r1403 r1520  
    5151      REAL SSUM
    5252c
     53      logical,save :: firstcall=.true.
     54      character(len=*),parameter :: modname="exner_hyb"
     55     
     56      ! Sanity check
     57      if (firstcall) then
     58        ! check that vertical discretization is compatible
     59        ! with this routine
     60        if (disvert_type.ne.1) then
     61          call abort_gcm(modname,
     62     &     "this routine should only be called if disvert_type==1",42)
     63        endif
     64       
     65        ! sanity checks for Shallow Water case (1 vertical layer)
     66        if (llm.eq.1) then
     67          if (kappa.ne.1) then
     68            call abort_gcm(modname,
     69     &      "kappa!=1 , but running in Shallow Water mode!!",42)
     70          endif
     71          if (cpp.ne.r) then
     72            call abort_gcm(modname,
     73     &      "cpp!=r , but running in Shallow Water mode!!",42)
     74          endif
     75        endif ! of if (llm.eq.1)
     76
     77        firstcall=.false.
     78      endif ! of if (firstcall)
    5379
    5480      if (llm.eq.1) then
    55         ! Specific behaviour for Shallow Water (1 vertical layer) case
    56      
    57         ! Sanity checks
    58         if (kappa.ne.1) then
    59           call abort_gcm("exner_hyb",
    60      &    "kappa!=1 , but running in Shallow Water mode!!",42)
    61         endif
    62         if (cpp.ne.r) then
    63         call abort_gcm("exner_hyb",
    64      &    "cpp!=r , but running in Shallow Water mode!!",42)
    65         endif
    6681       
    6782        ! Compute pks(:),pk(:),pkf(:)
    … …  
    7792        ! our work is done, exit routine
    7893        return
     94
    7995      endif ! of if (llm.eq.1)
    8096
     97!!!! General case:
    8198     
    8299      unpl2k    = 1.+ 2.* kappa

  • LMDZ5/trunk/libf/dyn3d/guide_mod.F90

    r1403 r1520  
    644644! -----------------------------------------------------------------
    645645    CALL pression( ip1jmp1, ap, bp, psi, p )
    646     CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
    647 
     646    if (disvert_type==1) then
     647      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     648    else ! we assume that we are in the disvert_type==2 case
     649      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
     650    endif
    648651!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
    649652    unskap=1./kappa

  • LMDZ5/trunk/libf/dyn3d/iniacademic.F90

    r1505 r1520  
    7171
    7272  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73 
     74  character(len=*),parameter :: modname="iniacademic"
     75  character(len=80) :: abort_message
    7376
    7477  !-----------------------------------------------------------------------
    … …  
    210213
    211214        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    212         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     215        if (disvert_type.eq.1) then
     216          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     217        elseif (disvert_type.eq.2) then
     218          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
     219        else
     220          write(abort_message,*) "Wrong value for disvert_type: ", &
     221                              disvert_type
     222          call abort_gcm(modname,abort_message,0)
     223        endif
    213224        CALL massdair(p,masse)
    214225

  • LMDZ5/trunk/libf/dyn3d/iniconst.F

    r1502 r1520  
    55
    66      USE control_mod
     7#ifdef CPP_IOIPSL
     8      use IOIPSL
     9#else
     10! if not using IOIPSL, we still need to use (a local version of) getin
     11      use ioipsl_getincom
     12#endif
    713
    814      IMPLICIT NONE
    … …  
    2127#include "iniprint.h"
    2228
    23 
     29      character(len=*),parameter :: modname="iniconst"
     30      character(len=80) :: abort_message
    2431c
    2532c
    … …  
    4855      r       = cpp * kappa
    4956
    50       write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
     57      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
    5158c
    5259c-----------------------------------------------------------------------
    5360
    54        CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    55 c
    56 c
    57        RETURN
    58        END
     61! vertical discretization: default behavior depends on planet_type flag
     62      if (planet_type=="earth") then
     63        disvert_type=1
     64      else
     65        disvert_type=2
     66      endif
     67      ! but user can also specify using one or the other in run.def:
     68      call getin('disvert_type',disvert_type)
     69      write(lunout,*) trim(modname),': disvert_type=',disvert_type
     70     
     71      if (disvert_type==1) then
     72       ! standard case for Earth (automatic generation of levels)
     73       call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
     74     &              scaleheight)
     75      else if (disvert_type==2) then
     76        ! standard case for planets (levels generated using z2sig.def file)
     77        call disvert_noterre
     78      else
     79        write(abort_message,*) "Wrong value for disvert_type: ",
     80     &                        disvert_type
     81        call abort_gcm(modname,abort_message,0)
     82      endif
     83
     84      END

  • LMDZ5/trunk/libf/dyn3d/leapfrog.F

    r1505 r1520  
    172172
    173173      character*80 dynhist_file, dynhistave_file
    174       character(len=20) :: modname
     174      character(len=*),parameter :: modname="leapfrog"
    175175      character*80 abort_message
    176176
    … …  
    193193      itaufin   = nday*day_step
    194194      itaufinp1 = itaufin +1
    195       modname="leapfrog"
    196195     
    197196
    … …  
    211210      dq(:,:,:)=0.
    212211      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    213       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     212      if (disvert_type==1) then
     213        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     214      else ! we assume that we are in the disvert_type==2 case
     215        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     216      endif
    214217
    215218c-----------------------------------------------------------------------
    … …  
    359362
    360363         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    361          CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     364         if (disvert_type==1) then
     365           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     366         else ! we assume that we are in the disvert_type==2 case
     367           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     368         endif
    362369
    363370!           rdaym_ini  = itau * dtvr / daysec
    … …  
    463470
    464471        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    465         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     472        if (disvert_type==1) then
     473          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     474        else ! we assume that we are in the disvert_type==2 case
     475          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     476        endif
    466477
    467478

  • LMDZ5/trunk/libf/dyn3d/limy.F

    r524 r1520  
    1 !
    2 ! $Header$
    3 !
     1c
     2c $Id$
     3c
    44      SUBROUTINE limy(s0,sy,sm,pente_max)
    55c
    … …  
    4040      REAL qbyv(ip1jm,llm)
    4141
    42       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
     42      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
    4343      Logical extremum,first
    4444      save first
    … …  
    117117
    118118c     print*,dyqv(iip1+1)
    119 c     apn=abs(dyq(1)/dyqv(iip1+1))
     119c     appn=abs(dyq(1)/dyqv(iip1+1))
    120120c     print*,dyq(ip1jm+1)
    121121c     print*,dyqv(ip1jm-iip1+1)
    122 c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     122c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    123123c     do ij=2,iim
    124 c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    125 c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     124c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     125c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    126126c     enddo
    127 c     apn=min(pente_max/apn,1.)
    128 c     aps=min(pente_max/aps,1.)
     127c     appn=min(pente_max/appn,1.)
     128c     apps=min(pente_max/apps,1.)
    129129
    130130
    … …  
    132132
    133133c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    134 c    &   apn=0.
     134c    &   appn=0.
    135135c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    136136c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    137 c    &   aps=0.
     137c    &   apps=0.
    138138
    139139c   limitation des pentes aux poles
    140140c     do ij=1,iip1
    141 c        dyq(ij)=apn*dyq(ij)
    142 c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     141c        dyq(ij)=appn*dyq(ij)
     142c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    143143c     enddo
    144144

  • LMDZ5/trunk/libf/dyn3d/logic.h

    r1492 r1520  
    33!
    44!
    5 !
     5! NB: keep items of different kinds in seperate common blocs to avoid
     6!     "misaligned commons" issues
    67!-----------------------------------------------------------------------
    78! INCLUDE 'logic.h'
    89
    9       COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
     10      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
    1011     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
    1112     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    12      &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
     13     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
    1314
     15      COMMON/logici/ iflag_phys,iflag_trac
     16     
    1417      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
    1518     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
    1619     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    1720     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
     21      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
     22                     ! (only used if disvert_type==2)
    1823
    19       INTEGER iflag_phys
     24      integer iflag_phys,iflag_trac
    2025!-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3d/vlsplt.F

    r595 r1520  
    1 !
    2 ! $Header$
    3 !
    4 c
     1c
     2c $Id$
    53c
    64
    … …  
    478476      REAL qbyv(ip1jm,llm)
    479477
    480       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     478      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    481479c     REAL newq,oldmasse
    482480      Logical extremum,first,testcpu
    … …  
    602600C     PRINT*,dyq(1)
    603601C     PRINT*,dyqv(iip1+1)
    604 C     apn=abs(dyq(1)/dyqv(iip1+1))
     602C     appn=abs(dyq(1)/dyqv(iip1+1))
    605603C     PRINT*,dyq(ip1jm+1)
    606604C     PRINT*,dyqv(ip1jm-iip1+1)
    607 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     605C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    608606C     DO ij=2,iim
    609 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    610 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     607C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     608C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    611609C     ENDDO
    612 C     apn=min(pente_max/apn,1.)
    613 C     aps=min(pente_max/aps,1.)
     610C     appn=min(pente_max/appn,1.)
     611C     apps=min(pente_max/apps,1.)
    614612C
    615613C
    … …  
    617615C
    618616C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    619 C    &   apn=0.
     617C    &   appn=0.
    620618C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    621619C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    622 C    &   aps=0.
     620C    &   apps=0.
    623621C
    624622C   limitation des pentes aux poles
    625623C     DO ij=1,iip1
    626 C        dyq(ij)=apn*dyq(ij)
    627 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     624C        dyq(ij)=appn*dyq(ij)
     625C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    628626C     ENDDO
    629627C

  • LMDZ5/trunk/libf/dyn3d/vlspltqs.F

    r595 r1520  
    1 !
    2 ! $Header$
    3 !
     1c
     2c $Id$
     3c
    44       SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
    55     ,                                  p,pk,teta                 )
    … …  
    634634C     PRINT*,dyq(1)
    635635C     PRINT*,dyqv(iip1+1)
    636 C     apn=abs(dyq(1)/dyqv(iip1+1))
     636C     appn=abs(dyq(1)/dyqv(iip1+1))
    637637C     PRINT*,dyq(ip1jm+1)
    638638C     PRINT*,dyqv(ip1jm-iip1+1)
    639 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     639C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    640640C     DO ij=2,iim
    641 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    642 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     641C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     642C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    643643C     ENDDO
    644 C     apn=min(pente_max/apn,1.)
    645 C     aps=min(pente_max/aps,1.)
     644C     appn=min(pente_max/appn,1.)
     645C     apps=min(pente_max/apps,1.)
    646646C
    647647C
    … …  
    649649C
    650650C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    651 C    &   apn=0.
     651C    &   appn=0.
    652652C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    653653C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    654 C    &   aps=0.
     654C    &   apps=0.
    655655C
    656656C   limitation des pentes aux poles
    657657C     DO ij=1,iip1
    658 C        dyq(ij)=apn*dyq(ij)
    659 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     658C        dyq(ij)=appn*dyq(ij)
     659C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    660660C     ENDDO
    661661C
Note: See TracChangeset for help on using the changeset viewer.