Ignore:
Timestamp:
Sep 7, 2012, 2:49:58 PM (12 years ago)
Author:
emillour
Message:

Common dynamics: updates to keep up with LMDZ5 Earth (rev 1649)
See file "DOC/chantiers/commit_importants.log" for details.
EM

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r128 r776  
    1 ! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
     1! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
    22
    33SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
    44
    55  ! Auteur : P. Le Van
     6
     7  use new_unit_m, only: new_unit
     8  use ioipsl, only: getin
     9  use assert_m, only: assert
    610
    711  IMPLICIT NONE
    … …  
    1822
    1923  real,intent(in) :: pa, preff
    20   real,intent(out) :: ap(llmp1), bp(llmp1)
     24  real,intent(out) :: ap(llmp1) ! in Pa
     25  real, intent(out):: bp(llmp1)
    2126  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
    2227  real,intent(out) :: presnivs(llm)
    … …  
    2631  real zk, zkm1, dzk1, dzk2, k0, k1
    2732
    28   INTEGER l
     33  INTEGER l, unit
    2934  REAL dsigmin
    3035  REAL alpha, beta, deltaz
    31   INTEGER iostat
    3236  REAL x
    3337  character(len=*),parameter :: modname="disvert"
    3438
     39  character(len=6):: vert_sampling
     40  ! (allowed values are "param", "tropo", "strato" and "read")
     41
    3542  !-----------------------------------------------------------------------
     43
     44  print *, "Call sequence information: disvert"
    3645
    3746  ! default scaleheight is 8km for earth
    3847  scaleheight=8.
    3948
    40   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     49  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
     50  call getin('vert_sampling', vert_sampling)
     51  print *, 'vert_sampling = ' // vert_sampling
    4152
    42   IF (iostat == 0) THEN
    43      ! cas 1 on lit les options dans sigma.def:
     53  select case (vert_sampling)
     54  case ("param")
     55     ! On lit les options dans sigma.def:
     56     OPEN(99, file='sigma.def', status='old', form='formatted')
    4457     READ(99, *) scaleheight ! hauteur d'echelle 8.
    4558     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
    … …  
    6982     sig(llm+1)=0.
    7083
    71      DO l = 1, llm
    72         dsig(l) = sig(l)-sig(l+1)
    73      end DO
    74   ELSE
    75      if (ok_strato) then
    76         if (llm==39) then
    77            dsigmin=0.3
    78         else if (llm==50) then
    79            dsigmin=1.
    80         else
    81            write(lunout,*) trim(modname), &
    82            ' ATTENTION discretisation z a ajuster'
    83            dsigmin=1.
    84         endif
    85         write(lunout,*) trim(modname), &
    86         ' Discretisation verticale DSIGMIN=',dsigmin
    87      endif
     84     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
     85     bp(llmp1) = 0.
    8886
     87     ap = pa * (sig - bp)
     88  case("tropo")
    8989     DO l = 1, llm
    9090        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    91 
    92         IF (ok_strato) THEN
    93            dsig(l) =(dsigmin + 7. * SIN(x)**2) &
    94                 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
    95         ELSE
    96            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    97         ENDIF
     91        dsig(l) = 1.0 + 7.0 * SIN(x)**2
    9892     ENDDO
    9993     dsig = dsig / sum(dsig)
    … …  
    10296        sig(l) = sig(l+1) + dsig(l)
    10397     ENDDO
    104   ENDIF
     98
     99     bp(1)=1.
     100     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     101     bp(llmp1) = 0.
     102
     103     ap(1)=0.
     104     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
     105  case("strato")
     106     if (llm==39) then
     107        dsigmin=0.3
     108     else if (llm==50) then
     109        dsigmin=1.
     110     else
     111        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
     112        dsigmin=1.
     113     endif
     114     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
     115
     116     DO l = 1, llm
     117        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     118        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     119             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     120     ENDDO
     121     dsig = dsig / sum(dsig)
     122     sig(llm+1) = 0.
     123     DO l = llm, 1, -1
     124        sig(l) = sig(l+1) + dsig(l)
     125     ENDDO
     126
     127     bp(1)=1.
     128     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     129     bp(llmp1) = 0.
     130
     131     ap(1)=0.
     132     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
     133  case("read")
     134     ! Read "ap" and "bp". First line is skipped (title line). "ap"
     135     ! should be in Pa. First couple of values should correspond to
     136     ! the surface, that is : "bp" should be in descending order.
     137     call new_unit(unit)
     138     open(unit, file="hybrid.txt", status="old", action="read", &
     139          position="rewind")
     140     read(unit, fmt=*) ! skip title line
     141     do l = 1, llm + 1
     142        read(unit, fmt=*) ap(l), bp(l)
     143     end do
     144     close(unit)
     145     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
     146          bp(llm + 1) == 0., "disvert: bad ap or bp values")
     147  case default
     148     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
     149  END select
    105150
    106151  DO l=1, llm
    … …  
    111156     nivsig(l)= REAL(l)
    112157  ENDDO
    113 
    114   ! .... Calculs de ap(l) et de bp(l) ....
    115   ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    116 
    117   bp(llmp1) = 0.
    118 
    119   DO l = 1, llm
    120      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    121      ap(l) = pa * ( sig(l) - bp(l) )
    122   ENDDO
    123 
    124   bp(1)=1.
    125   ap(1)=0.
    126 
    127   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    128158
    129159  write(lunout, *)  trim(modname),': BP '
Note: See TracChangeset for help on using the changeset viewer.