source: LMDZ5/trunk/libf/dyn3dpar/disvert.F90 @ 1934

Last change on this file since 1934 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.3 KB
RevLine 
[1279]1! $Id: disvert.F90 1907 2013-11-26 13:10:46Z jghattas $
[630]2
[1520]3SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
[630]4
[1472]5  ! Auteur : P. Le Van
[630]6
[1620]7  use new_unit_m, only: new_unit
8  use ioipsl, only: getin
[1645]9  use assert_m, only: assert
[1620]10
[1472]11  IMPLICIT NONE
[630]12
[1472]13  include "dimensions.h"
14  include "paramet.h"
15  include "iniprint.h"
16  include "logic.h"
[630]17
[1472]18  ! s = sigma ** kappa : coordonnee verticale
19  ! dsig(l) : epaisseur de la couche l ds la coord. s
20  ! sig(l) : sigma a l'interface des couches l et l-1
21  ! ds(l) : distance entre les couches l et l-1 en coord.s
[630]22
[1520]23  real,intent(in) :: pa, preff
[1645]24  real,intent(out) :: ap(llmp1) ! in Pa
25  real, intent(out):: bp(llmp1)
[1520]26  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
27  real,intent(out) :: presnivs(llm)
28  real,intent(out) :: scaleheight
[630]29
[1472]30  REAL sig(llm+1), dsig(llm)
31  real zk, zkm1, dzk1, dzk2, k0, k1
[630]32
[1620]33  INTEGER l, unit
[1472]34  REAL dsigmin
[1524]35  REAL alpha, beta, deltaz
[1520]36  REAL x
37  character(len=*),parameter :: modname="disvert"
[630]38
[1620]39  character(len=6):: vert_sampling
40  ! (allowed values are "param", "tropo", "strato" and "read")
41
[1472]42  !-----------------------------------------------------------------------
43
[1620]44  print *, "Call sequence information: disvert"
45
[1520]46  ! default scaleheight is 8km for earth
47  scaleheight=8.
[1472]48
[1623]49  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
[1620]50  call getin('vert_sampling', vert_sampling)
51  print *, 'vert_sampling = ' // vert_sampling
[1472]52
[1620]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')
[1520]57     READ(99, *) scaleheight ! hauteur d'echelle 8.
[1472]58     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
59     READ(99, *) beta ! facteur d'acroissement en haut 1.3
60     READ(99, *) k0 ! nombre de couches dans la transition surf
61     READ(99, *) k1 ! nombre de couches dans la transition haute
62     CLOSE(99)
[1520]63     alpha=deltaz/(llm*scaleheight)
64     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
65                               scaleheight, alpha, k0, k1, beta
[1472]66
67     alpha=deltaz/tanh(1./k0)*2.
68     zkm1=0.
69     sig(1)=1.
70     do l=1, llm
[1520]71        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
72             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
73                  *beta**(l-(llm-k1))/log(beta))
74        zk=-scaleheight*log(sig(l+1))
[630]75
76        dzk1=alpha*tanh(l/k0)
77        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
[1472]78        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
[630]79        zkm1=zk
[1472]80     enddo
[630]81
[1472]82     sig(llm+1)=0.
[1645]83
84     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
85     bp(llmp1) = 0.
86
87     ap = pa * (sig - bp)
[1620]88  case("tropo")
[1472]89     DO l = 1, llm
[1620]90        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
91        dsig(l) = 1.0 + 7.0 * SIN(x)**2
92     ENDDO
93     dsig = dsig / sum(dsig)
94     sig(llm+1) = 0.
95     DO l = llm, 1, -1
96        sig(l) = sig(l+1) + dsig(l)
97     ENDDO
[1645]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))
[1620]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.
[1472]113     endif
[1620]114     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
[630]115
[1472]116     DO l = 1, llm
[1480]117        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
[1620]118        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
119             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
[1472]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
[1645]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))
[1620]133  case("read")
[1645]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.
[1620]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
[1645]142        read(unit, fmt=*) ap(l), bp(l)
[1620]143     end do
144     close(unit)
[1645]145     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
146          bp(llm + 1) == 0., "disvert: bad ap or bp values")
[1620]147  case default
148     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
149  END select
[1000]150
[1472]151  DO l=1, llm
152     nivsigs(l) = REAL(l)
153  ENDDO
[630]154
[1472]155  DO l=1, llmp1
156     nivsig(l)= REAL(l)
157  ENDDO
[630]158
[1520]159  write(lunout, *)  trim(modname),': BP '
[1472]160  write(lunout, *) bp
[1520]161  write(lunout, *)  trim(modname),': AP '
[1472]162  write(lunout, *) ap
[630]163
[1472]164  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
165  write(lunout, *)'couches calcules pour une pression de surface =', preff
[1520]166  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
167  write(lunout, *) scaleheight,' km'
[1472]168  DO l = 1, llm
169     dpres(l) = bp(l) - bp(l+1)
170     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
171     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
[1520]172          log(preff/presnivs(l))*scaleheight &
173          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
[1472]174          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
175  ENDDO
[1279]176
[1520]177  write(lunout, *) trim(modname),': PRESNIVS '
[1472]178  write(lunout, *) presnivs
[1279]179
[1472]180END SUBROUTINE disvert
Note: See TracBrowser for help on using the repository browser.