Changeset 1520


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
Files:
5 added
26 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

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

    r1279 r1520  
    11!
    2 ! $Header$
     2! $Id$
    33!
    44!-----------------------------------------------------------------------
    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
    1112
    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
     31
     32 !-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3dpar/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/dyn3dpar/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/dyn3dpar/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/dyn3dpar/exner_hyb_p.F

    r1403 r1520  
    5353      EXTERNAL SSUM
    5454      INTEGER ije,ijb,jje,jjb
    55 c
    56 c$OMP BARRIER
    57 
    58       if (llm.eq.1) then
    59         ! Specific behaviour for Shallow Water (1 vertical layer) case
    60      
    61         ! Sanity checks
    62         if (kappa.ne.1) then
    63           call abort_gcm("exner_hyb",
    64      &    "kappa!=1 , but running in Shallow Water mode!!",42)
    65         endif
    66         if (cpp.ne.r) then
    67         call abort_gcm("exner_hyb",
    68      &    "cpp!=r , but running in Shallow Water mode!!",42)
     55      logical,save :: firstcall=.true.
     56!$OMP THREADPRIVATE(firstcall)
     57      character(len=*),parameter :: modname="exner_hyb_p"
     58c
     59
     60      ! Sanity check
     61      if (firstcall) then
     62        ! check that vertical discretization is compatible
     63        ! with this routine
     64        if (disvert_type.ne.1) then
     65          call abort_gcm(modname,
     66     &     "this routine should only be called if disvert_type==1",42)
    6967        endif
    7068       
     69        ! sanity checks for Shallow Water case (1 vertical layer)
     70        if (llm.eq.1) then
     71          if (kappa.ne.1) then
     72            call abort_gcm(modname,
     73     &      "kappa!=1 , but running in Shallow Water mode!!",42)
     74          endif
     75          if (cpp.ne.r) then
     76            call abort_gcm(modname,
     77     &      "cpp!=r , but running in Shallow Water mode!!",42)
     78          endif
     79        endif ! of if (llm.eq.1)
     80
     81        firstcall=.false.
     82      endif ! of if (firstcall)
     83
     84c$OMP BARRIER
     85
     86! Specific behaviour for Shallow Water (1 vertical layer) case
     87      if (llm.eq.1) then
     88     
    7189        ! Compute pks(:),pk(:),pkf(:)
    7290        ijb=ij_begin
    … …  
    116134      endif ! of if (llm.eq.1)
    117135
     136!!!! General case:
    118137
    119138      unpl2k    = 1.+ 2.* kappa

  • LMDZ5/trunk/libf/dyn3dpar/gcm.F

    r1454 r1520  
    538538c       write(78,*) 'q',q
    539539
    540 c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
     540c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
    541541      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
    542542     .              time_0)

  • LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90

    r1512 r1520  
    455455!       Calcul niveaux pression milieu de couches
    456456        CALL pression_p( ip1jmp1, ap, bp, ps, p )
    457         CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
    458         unskap=1./kappa
     457        if (disvert_type==1) then
     458          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     459        else
     460          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
     461        endif
     462        unskap=1./kappa
    459463        DO l = 1, llm
    460464            DO j=jjb_u,jje_u
    … …  
    751755    ELSE
    752756        CALL pression_p( ip1jmp1, ap, bp, psi, p )
    753         CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     757        if (disvert_type==1) then
     758          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     759        else ! we assume that we are in the disvert_type==2 case
     760          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
     761        endif
    754762        unskap=1./kappa
    755763        DO l = 1, llm

  • LMDZ5/trunk/libf/dyn3dpar/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/dyn3dpar/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/dyn3dpar/leapfrog_p.F

    r1505 r1520  
    164164
    165165      character*80 dynhist_file, dynhistave_file
    166       character*20 modname
     166      character(len=*),parameter :: modname="leapfrog"
    167167      character*80 abort_message
    168168
    … …  
    208208      itaufin   = nday*day_step
    209209      itaufinp1 = itaufin +1
    210       modname="leapfrog_p"
    211210
    212211      itau = 0
    … …  
    236235      dq(:,:,:)=0.
    237236      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    238       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     237      if (disvert_type==1) then
     238        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     239      else ! we assume that we are in the disvert_type==2 case
     240        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     241      endif
    239242c$OMP END MASTER
    240243c-----------------------------------------------------------------------
    … …  
    691694
    692695c$OMP BARRIER
    693          CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     696         if (disvert_type==1) then
     697           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     698         else ! we assume that we are in the disvert_type==2 case
     699           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     700         endif
    694701c$OMP BARRIER
    695702           jD_cur = jD_ref + day_ini - day_ref
    … …  
    10301037        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
    10311038c$OMP BARRIER
    1032         CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1039        if (disvert_type==1) then
     1040          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1041        else ! we assume that we are in the disvert_type==2 case
     1042          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     1043        endif
    10331044c$OMP BARRIER
    10341045

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

    r764 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
    … …  
    5252      REAL      SSUM
    5353      integer ismax,ismin
    54       EXTERNAL  SSUM, ismin,ismax
     54      EXTERNAL  SSUM, convflu,ismin,ismax
     55      EXTERNAL filtreg
    5556
    5657      data first/.true./
    … …  
    116117
    117118c     print*,dyqv(iip1+1)
    118 c     apn=abs(dyq(1)/dyqv(iip1+1))
     119c     appn=abs(dyq(1)/dyqv(iip1+1))
    119120c     print*,dyq(ip1jm+1)
    120121c     print*,dyqv(ip1jm-iip1+1)
    121 c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     122c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    122123c     do ij=2,iim
    123 c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    124 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)
    125126c     enddo
    126 c     apn=min(pente_max/apn,1.)
    127 c     aps=min(pente_max/aps,1.)
     127c     appn=min(pente_max/appn,1.)
     128c     apps=min(pente_max/apps,1.)
    128129
    129130
    … …  
    131132
    132133c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    133 c    &   apn=0.
     134c    &   appn=0.
    134135c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    135136c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    136 c    &   aps=0.
     137c    &   apps=0.
    137138
    138139c   limitation des pentes aux poles
    139140c     do ij=1,iip1
    140 c        dyq(ij)=apn*dyq(ij)
    141 c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     141c        dyq(ij)=appn*dyq(ij)
     142c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    142143c     enddo
    143144

  • LMDZ5/trunk/libf/dyn3dpar/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
    20 !$OMP THREADPRIVATE(/logic/)
     24      integer iflag_phys,iflag_trac
     25!$OMP THREADPRIVATE(/logicl/)
     26!$OMP THREADPRIVATE(/logici/)
    2127!-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3dpar/vlsplt_p.F

    r764 r1520  
     1c
     2c $Id$
     3c
     4
    15      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
    26c
    … …  
    561565      REAL qbyv(ip1jm,llm)
    562566
    563       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     567      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    564568c     REAL newq,oldmasse
    565569      Logical extremum,first,testcpu
    … …  
    732736C     PRINT*,dyq(1)
    733737C     PRINT*,dyqv(iip1+1)
    734 C     apn=abs(dyq(1)/dyqv(iip1+1))
     738C     appn=abs(dyq(1)/dyqv(iip1+1))
    735739C     PRINT*,dyq(ip1jm+1)
    736740C     PRINT*,dyqv(ip1jm-iip1+1)
    737 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     741C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    738742C     DO ij=2,iim
    739 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    740 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     743C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     744C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    741745C     ENDDO
    742 C     apn=min(pente_max/apn,1.)
    743 C     aps=min(pente_max/aps,1.)
     746C     appn=min(pente_max/appn,1.)
     747C     apps=min(pente_max/apps,1.)
    744748C
    745749C
    … …  
    747751C
    748752C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    749 C    &   apn=0.
     753C    &   appn=0.
    750754C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    751755C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    752 C    &   aps=0.
     756C    &   apps=0.
    753757C
    754758C   limitation des pentes aux poles
    755759C     DO ij=1,iip1
    756 C        dyq(ij)=apn*dyq(ij)
    757 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     760C        dyq(ij)=appn*dyq(ij)
     761C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    758762C     ENDDO
    759763C

  • LMDZ5/trunk/libf/dyn3dpar/vlspltqs_p.F

    r764 r1520  
    1 !
    2 ! $Header$
    3 !
     1c
     2c $Id$
     3c
    44       SUBROUTINE vlspltqs_p ( q,pente_max,masse,w,pbaru,pbarv,pdt,
    55     ,                                  p,pk,teta                 )
    … …  
    772772C     PRINT*,dyq(1)
    773773C     PRINT*,dyqv(iip1+1)
    774 C     apn=abs(dyq(1)/dyqv(iip1+1))
     774C     appn=abs(dyq(1)/dyqv(iip1+1))
    775775C     PRINT*,dyq(ip1jm+1)
    776776C     PRINT*,dyqv(ip1jm-iip1+1)
    777 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     777C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    778778C     DO ij=2,iim
    779 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    780 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     779C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     780C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    781781C     ENDDO
    782 C     apn=min(pente_max/apn,1.)
    783 C     aps=min(pente_max/aps,1.)
     782C     appn=min(pente_max/appn,1.)
     783C     apps=min(pente_max/apps,1.)
    784784C
    785785C
    … …  
    787787C
    788788C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    789 C    &   apn=0.
     789C    &   appn=0.
    790790C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    791791C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    792 C    &   aps=0.
     792C    &   apps=0.
    793793C
    794794C   limitation des pentes aux poles
    795795C     DO ij=1,iip1
    796 C        dyq(ij)=apn*dyq(ij)
    797 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     796C        dyq(ij)=appn*dyq(ij)
     797C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    798798C     ENDDO
    799799C
Note: See TracChangeset for help on using the changeset viewer.