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

Diff of /trunk/dyn3d/disvert.f90

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

trunk/libf/dyn3d/comvert.f90 revision 42 by guez, Thu Mar 24 11:52:41 2011 UTC trunk/dyn3d/disvert.f revision 313 by guez, Mon Dec 10 15:54:30 2018 UTC
# Line 1  Line 1 
1  module comvert  module disvert_m
2    
3    use dimens_m, only: llm    use dimensions, only: llm
4    
5    implicit none    implicit none
6    
7    private llm    private llm, hybrid, funcd, y, ya, compute_ab
8    
9    real ap(llm+1), pa ! in Pa    real, save:: ap(llm+1) ! 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 s(llm+1)
13    real nivsigs(llm), nivsig(llm+1)    ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
14      ! half-level, between layers "l" and "l-1"
15    
16    save    real, save:: presnivs(llm)
17      ! approximate full level pressure for a reference surface pressure, in Pa
18    
19      real, parameter:: preff = 101325. ! in Pa
20      real y, ya ! for the hybrid function
21    
22  contains  contains
23    
24    SUBROUTINE disvert    SUBROUTINE disvert
25    
26      ! 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
27      ! Auteur : P. Le Van      ! Author: P. Le Van
28    
29      ! This procedure sets the vertical grid.      ! This procedure sets the vertical grid. It defines the host
30      ! It defines the host variables "ap", "bp", "presnivs", "nivsigs"      ! variables "ap", "bp", "presnivs". "pa" should be defined before
31      ! and "nivsig".      ! this procedure is called.
32      ! "pa" should be defined before this procedure is called.  
33        ! Libraries:
34      use nr_util, only: pi      use jumble, only: read_column, new_unit
35      use new_unit_m, only: new_unit      use nr_util, only: pi, assert
36        
37      REAL s(llm+1)      use dynetat0_chosen_m, only: pa
38      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at      use unit_nml_m, only: unit_nml
39      ! the interface between layers "l" and "l-1"  
40        ! Local:
41    
42      real ds(llm)      real ds(llm)
43      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"      ! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s"
44    
45      INTEGER l, unit      INTEGER l, unit
46      REAL alpha, x(llm), trash      REAL x(llm)
47        real:: dsigmin = 1.
48      character(len=7):: s_sampling = "LMD5"      real zz(llm) ! in km
49      ! (other allowed values are "param", "strato1", "strato2" and "read")  
50        character(len=20):: vert_sampling = "tropo"
51      real:: h = 7. ! scale height, in km      ! Allowed values: "tropo", "strato_custom", "strato",
52      ! (used only if "s_sampling" == "param" or "strato1")      ! "read_hybrid", "read_pressure".
53    
54      ! These variables are used only in the case 's_sampling == "param"':      ! These variables are used only in the case vert_sampling ==
55      real:: deltaz = 0.04 ! épaisseur de la première couche      ! "strato_custom", and all are in km:
56      real:: beta = 1.3 ! facteur d'accroissement en haut      real:: vert_scale = 7. ! scale height
57      real:: k0 = 20. ! nombre de couches dans la transition surface      real:: vert_dzmin = 0.017 ! width of first layer
58      real:: k1 = 1.2 ! nombre de couches dans la transition haute      real:: vert_dzlow = 1. ! dz in the low atmosphere
59        real:: vert_z0low = 8.7 ! height at which resolution reaches dzlow
60      REAL ZZ(llm + 1), DZ(llm) ! in km      real:: vert_dzmid = 2. ! dz in the mid atmosphere
61        real:: vert_z0mid = 70. ! height at which resolution reaches dzmid
62      namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling      real:: vert_h_mid = 20. ! width of the transition
63        real:: vert_dzhig = 11. ! dz in the high atmosphere
64        real:: vert_z0hig = 75. ! height at which resolution reaches dz
65        real:: vert_h_hig = 20. ! width of the transition
66    
67        real, allocatable:: p(:) ! (2:llm or llm + 1) pressure (in hPa)
68    
69        namelist /disvert_nml/vert_sampling, vert_scale, vert_dzmin, vert_dzlow, &
70             vert_z0low, vert_dzmid, vert_z0mid, vert_h_mid, vert_dzhig, &
71             vert_z0hig, vert_h_hig, dsigmin
72    
73      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
74    
75      print *, "Call sequence information: disvert"      print *, "Call sequence information: disvert"
76    
77      forall (l = 1: llm) nivsigs(l) = REAL(l)      write(unit=*, nml=disvert_nml)
     forall (l = 1: llm + 1) nivsig(l) = REAL(l)  
   
     ! Compute "s":  
   
78      print *, "Enter namelist 'disvert_nml'."      print *, "Enter namelist 'disvert_nml'."
79      read(unit=*, nml=disvert_nml)      read(unit=*, nml=disvert_nml)
80      write(unit=*, nml=disvert_nml)      write(unit_nml, nml=disvert_nml)
81    
82        s(1) = 1.
83        s(llm+1) = 0.
84    
85        select case (vert_sampling)
86    
87        case ("tropo")
88           ! with llm = 19 and dsigmin = 1 for CMIP 3
89    
     select case (s_sampling)  
     case ("param")  
        s(1) = 1.  
        s(llm+1) = 0.  
        alpha = deltaz / tanh(1./k0) * 2.  
        forall (l = 2: llm) s(l) &  
             = cosh((l - 1) / k0) **(- alpha * k0 / h) &  
             * exp(- alpha / h * tanh((llm - k1) / k0) &  
             * beta **(l - 1 - (llm - k1)) / log(beta))  
     case ("LMD5")  
        ! Ancienne discrétisation  
        s(1) = 1.  
        s(llm+1) = 0.  
90         forall (l = 1: llm) ds(l) &         forall (l = 1: llm) ds(l) &
91              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2              = dsigmin + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
92         ds = ds / sum(ds)         ds = ds / sum(ds)
93    
94         DO l = llm, 2, -1         DO l = llm, 2, -1
95            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
96         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.  
97    
98         zz(1) = 0.         call compute_ab
99         DO l = 2, llm + 1  
100            zz(l) = zz(l - 1) + dz(l - 1)      case ("strato")
101         end DO         ! with llm = 39 and dsigmin = 0.3 for CMIP5
102    
        s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &  
             / (1. - exp(- zz(llm + 1) / h))  
     case ("strato2")  
        ! Recommended by F. Lott for a domain including the stratosphere  
        s(1) = 1.  
        s(llm+1) = 0.  
103         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
104    
105         ds = (1. + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.         ds = (dsigmin + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
106         ds = ds / sum(ds)         ds = ds / sum(ds)
107    
108         DO l = llm, 2, -1         DO l = llm, 2, -1
109            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
110         ENDDO         ENDDO
111      case("read")  
112           call compute_ab
113    
114        case ("strato_custom")
115           ! with llm = 79 for CMIP 6
116    
117           zz(1) = 0.
118           DO l = 1, llm - 1
119              zz(l + 1) = zz(l) + vert_dzmin + vert_dzlow &
120                   * TANH(zz(l) / vert_z0low) + (vert_dzmid - vert_dzlow) &
121                   * (TANH((zz(l) - vert_z0mid) / vert_h_mid) &
122                   - TANH(- vert_z0mid / vert_h_mid)) &
123                   + (vert_dzhig - vert_dzmid - vert_dzlow) &
124                   * (TANH((zz(l) - vert_z0hig) / vert_h_hig) &
125                   - TANH(- vert_z0hig / vert_h_hig))
126           ENDDO
127    
128           allocate(p(2: llm))
129           p = preff * EXP(- zz(2:) / vert_scale)
130           ya = pa / preff
131           s(2: llm) = hybrid(p)
132    
133           call compute_ab
134    
135        case("read_hybrid")
136           ! Read "ap" and "bp". First line is skipped (title line). "ap"
137           ! should be in Pa. First couple of values should correspond to
138           ! the surface, that is : "bp" should be in descending order.
139         call new_unit(unit)         call new_unit(unit)
140         open(unit, file="hybrid.csv", status="old", action="read", &         open(unit, file="hybrid.csv", status="old", action="read", &
141              position="rewind")              position="rewind")
142         read(unit, fmt=*) ! skip title line         read(unit, fmt=*) ! skip title line
143         do l = 1, llm + 1         do l = 1, llm + 1
144            read(unit, fmt=*) trash, s(l)            read(unit, fmt=*) ap(l), bp(l)
145         end do         end do
146         close(unit)         close(unit)
147         ! Quick check:         ! Quick check:
148         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., &
149            print *, '"s" should be in descending order, from 1 to 0.'              bp(llm + 1) == 0., "disvert: bad ap or bp values")
150            stop 1         s(2: llm) = ap(2: llm) / pa + bp(2: llm)
151         end if  
152        case("read_pressure")
153           ! Read pressure values, in Pa, in descending order, from preff
154           ! to 0. First line is skipped (title line).
155           call read_column("pressure.txt", p, first=2)
156           call assert(size(p) == llm + 1, "disvert: bad number of pressure values")
157           ! Quick check:
158           call assert(p(1) == preff, p(llm + 1) == 0., &
159                "disvert: bad pressure values")
160           ya = pa / preff
161           s(2: llm) = hybrid(p(2: llm))
162           call compute_ab
163    
164      case default      case default
165         print *, 'Wrong value for "s_sampling"'         print *, 'Wrong value for "vert_sampling"'
166         stop 1         stop 1
     END select  
167    
168      ! Calcul de "ap" et "bp" :      END select
     bp(:llm) = EXP(1. - 1. / s(:llm)**2)  
     bp(llm + 1) = 0.  
     ap = pa * (s - bp)  
169    
170      forall (l = 1: llm) presnivs(l) = 0.5 &      forall (l = 1: llm) presnivs(l) = 0.5 &
171           * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)           * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
172    
173    END SUBROUTINE disvert    END SUBROUTINE disvert
174    
175  end module comvert    !**********************************************************
176    
177      subroutine compute_ab
178    
179        use dynetat0_chosen_m, only: pa
180    
181        ! Calcul de "ap" et "bp".
182    
183        where (s >= 1. / sqrt(1. - log(tiny(0.))))
184           bp = exp(1. - 1. / s**2)
185        elsewhere
186           bp = 0.
187        end where
188    
189        ap = pa * (s - bp)
190    
191      end subroutine compute_ab
192    
193      !**********************************************************
194    
195      function hybrid(p)
196    
197        ! This procedure computes the hybrid sigma-pressure coordinate
198        ! from pressure values, assuming some reference surface
199        ! pressure. The procedure assumes, and does not check, that
200        ! pressure values are given in descending order.
201    
202        use numer_rec_95, only: rtsafe
203    
204        real, intent(in):: p(:) ! pressure (in hPa)
205        real hybrid(size(p)) ! hybrid sigma-pressure coordinate
206    
207        ! Local:
208        integer l
209    
210        !-------------------------------------------------------
211    
212        y = p(1) / preff
213        hybrid(1) = rtsafe(funcd, x1 = 0., x2 = 1., xacc = 1e-4)
214    
215        do l = 2, size(p)
216           y = p(l) / preff
217           ! Assuming descending order in pressure:
218           hybrid(l) = rtsafe(funcd, x1 = 0., x2 = hybrid(l - 1), &
219                xacc = hybrid(l - 1) * 1e-4)
220        end do
221    
222      end function hybrid
223    
224      !**********************************************************
225    
226      SUBROUTINE funcd(s, fval, fderiv)
227    
228        REAL, INTENT(IN):: s
229        REAL, INTENT(OUT):: fval, fderiv
230    
231        ! Local:
232        real b
233    
234        !------------------------------------
235    
236        if (s**3 > 1. / huge(0.)) then
237           b = exp(1. - 1. / s**2)
238           fval = ya * s + (1. - ya) * b - y
239           fderiv = ya + 2 * (1. - ya) * b / s**3
240        else
241           fval = ya * s - y
242           fderiv = ya
243        end if
244    
245      END SUBROUTINE funcd
246    
247    end module disvert_m
Legend:
Removed from v.42  
changed lines
  Added in v.313
  ViewVC Help
Powered by ViewVC 1.1.21