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

Diff of /trunk/dyn3d/disvert.f

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

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