/[lmdze]/trunk/dyn3d/disvert.f
ViewVC logotype

Diff of /trunk/dyn3d/disvert.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/comvert.f90 revision 48 by guez, Tue Jul 19 12:54:20 2011 UTC trunk/dyn3d/disvert.f revision 90 by guez, Wed Mar 12 21:16:36 2014 UTC
# Line 1  Line 1 
1  module comvert  module disvert_m
2    
3    use dimens_m, only: llm    use dimens_m, only: llm
4    
# Line 6  module comvert Line 6  module comvert
6    
7    private llm    private llm
8    
9    real ap(llm+1), pa ! in Pa    real, save:: ap(llm+1), pa ! in Pa
10    real bp(llm+1)    real, save:: bp(llm+1)
11    real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa  
12    real, parameter:: preff = 101325. ! in Pa    real, save:: presnivs(llm)
13    real nivsigs(llm), nivsig(llm+1)    ! pressions approximatives des milieux de couches, en Pa
14    
15    save    real, parameter:: preff = 101325. ! in Pa
16    
17  contains  contains
18    
19    SUBROUTINE disvert    SUBROUTINE disvert
20    
21      ! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05      ! From dyn3d/disvert.F, version 1.1.1.1, 2004/05/19 12:53:05
22      ! Auteur : P. Le Van      ! Author: P. Le Van
23    
24      ! This procedure sets the vertical grid.      ! This procedure sets the vertical grid. It defines the host
25      ! It defines the host variables "ap", "bp", "presnivs", "nivsigs"      ! variables "ap", "bp", "presnivs". "pa" should be defined before
26      ! and "nivsig".      ! this procedure is called.
     ! "pa" should be defined before this procedure is called.  
27    
     use nr_util, only: pi  
28      use jumble, only: new_unit      use jumble, only: new_unit
29        use nr_util, only: pi, assert
30        use unit_nml_m, only: unit_nml
31    
32        ! Local:
33    
34      REAL s(llm+1)      REAL s(llm+1)
35      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
36      ! the interface between layers "l" and "l-1"      ! the interface between layers "l" and "l-1"
37    
38      real ds(llm)      real ds(llm)
39      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"      ! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s"
40    
41      INTEGER l, unit      INTEGER l, unit
42      REAL alpha, x(llm), trash      REAL alpha, x(llm)
43    
44      character(len=7):: s_sampling = "LMD5"      character(len=7):: vert_sampling = "tropo"
45      ! (other allowed values are "param", "strato1", "strato2" and "read")      ! other allowed values are "param", "strato" and "read"
46    
47      real:: h = 7. ! scale height, in km      real:: h = 7. ! scale height, in km
48      ! (used only if "s_sampling" == "param" or "strato1")      ! used only if vert_sampling == "param"
49    
50      ! These variables are used only in the case 's_sampling == "param"':      ! These variables are used only in the case vert_sampling == "param":
51      real:: deltaz = 0.04 ! épaisseur de la première couche      real:: deltaz = 0.04 ! \'epaisseur de la premi\`ere couche
52      real:: beta = 1.3 ! facteur d'accroissement en haut      real:: beta = 1.3 ! facteur d'accroissement en haut
53      real:: k0 = 20. ! nombre de couches dans la transition surface      real:: k0 = 20. ! nombre de couches dans la transition surface
54      real:: k1 = 1.2 ! nombre de couches dans la transition haute      real:: k1 = 1.2 ! nombre de couches dans la transition haute
55    
56      REAL ZZ(llm + 1), DZ(llm) ! in km      namelist /disvert_nml/h, deltaz, beta, k0, k1, vert_sampling
   
     namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling  
57    
58      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
59    
60      print *, "Call sequence information: disvert"      print *, "Call sequence information: disvert"
61    
     forall (l = 1: llm) nivsigs(l) = REAL(l)  
     forall (l = 1: llm + 1) nivsig(l) = REAL(l)  
   
     ! Compute "s":  
   
62      print *, "Enter namelist 'disvert_nml'."      print *, "Enter namelist 'disvert_nml'."
63      read(unit=*, nml=disvert_nml)      read(unit=*, nml=disvert_nml)
64      write(unit=*, nml=disvert_nml)      write(unit_nml, nml=disvert_nml)
65    
66      select case (s_sampling)      select case (vert_sampling)
67      case ("param")      case ("param")
68         s(1) = 1.         s(1) = 1.
69         s(llm+1) = 0.         s(llm+1) = 0.
# Line 77  contains Line 72  contains
72              = cosh((l - 1) / k0) **(- alpha * k0 / h) &              = cosh((l - 1) / k0) **(- alpha * k0 / h) &
73              * exp(- alpha / h * tanh((llm - k1) / k0) &              * exp(- alpha / h * tanh((llm - k1) / k0) &
74              * beta **(l - 1 - (llm - k1)) / log(beta))              * beta **(l - 1 - (llm - k1)) / log(beta))
75      case ("LMD5")         call compute_ab
76         ! Ancienne discrétisation      case ("tropo")
77         s(1) = 1.         s(1) = 1.
78         s(llm+1) = 0.         s(llm+1) = 0.
79         forall (l = 1: llm) ds(l) &         forall (l = 1: llm) ds(l) &
80              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2              = 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
81         ds = ds / sum(ds)         ds = ds / sum(ds)
82    
83         DO l = llm, 2, -1         DO l = llm, 2, -1
84            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
85         ENDDO         ENDDO
     case ("strato1")  
        ! F. Lott 70 niveaux et plus  
        s(1) = 1.  
        s(llm+1) = 0.  
        forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &  
                + TANH(REAL(l - llm) / 10.) / 2.  
86    
87         zz(1) = 0.         call compute_ab
88         DO l = 2, llm + 1      case ("strato")
           zz(l) = zz(l - 1) + dz(l - 1)  
        end DO  
   
        s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &  
             / (1. - exp(- zz(llm + 1) / h))  
     case ("strato2")  
89         ! Recommended by F. Lott for a domain including the stratosphere         ! Recommended by F. Lott for a domain including the stratosphere
90         s(1) = 1.         s(1) = 1.
91         s(llm+1) = 0.         s(llm+1) = 0.
92         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
93    
94         ds = (1. + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.         ds = (0.3 + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
95         ds = ds / sum(ds)         ds = ds / sum(ds)
96    
97         DO l = llm, 2, -1         DO l = llm, 2, -1
98            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
99         ENDDO         ENDDO
100    
101           call compute_ab
102      case("read")      case("read")
103           ! Read "ap" and "bp". First line is skipped (title line). "ap"
104           ! should be in Pa. First couple of values should correspond to
105           ! the surface, that is : "bp" should be in descending order.
106         call new_unit(unit)         call new_unit(unit)
107         open(unit, file="hybrid.csv", status="old", action="read", &         open(unit, file="hybrid.csv", status="old", action="read", &
108              position="rewind")              position="rewind")
109         read(unit, fmt=*) ! skip title line         read(unit, fmt=*) ! skip title line
110         do l = 1, llm + 1         do l = 1, llm + 1
111            read(unit, fmt=*) trash, s(l)            read(unit, fmt=*) ap(l), bp(l)
112         end do         end do
113         close(unit)         close(unit)
114         ! Quick check:         ! Quick check:
115         if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then         call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
116            print *, '"s" should be in descending order, from 1 to 0.'              bp(llm + 1) == 0., "disvert: bad ap or bp values")
           stop 1  
        end if  
117      case default      case default
118         print *, 'Wrong value for "s_sampling"'         print *, 'Wrong value for "vert_sampling"'
119         stop 1         stop 1
120      END select      END select
121    
     ! Calcul de "ap" et "bp" :  
     bp(:llm) = EXP(1. - 1. / s(:llm)**2)  
     bp(llm + 1) = 0.  
     ap = pa * (s - bp)  
   
122      forall (l = 1: llm) presnivs(l) = 0.5 &      forall (l = 1: llm) presnivs(l) = 0.5 &
123           * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)           * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
124    
125      contains
126    
127        subroutine compute_ab
128    
129          ! Calcul de "ap" et "bp" :
130          bp(:llm) = EXP(1. - 1. / s(:llm)**2)
131          bp(llm + 1) = 0.
132          ap = pa * (s - bp)
133    
134        end subroutine compute_ab
135    
136    END SUBROUTINE disvert    END SUBROUTINE disvert
137    
138  end module comvert  end module disvert_m
Legend:
Removed from v.48  
changed lines
  Added in v.90
  ViewVC Help
Powered by ViewVC 1.1.21