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