/[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 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC
# Line 6  module comvert Line 6  module comvert
6    
7    private llm    private llm
8    
9    real ap(llm+1) ! in Pa    real ap(llm+1), pa ! in Pa
10    real bp(llm+1), dpres(llm)    real bp(llm+1)
11    real presnivs(llm), pa ! in Pa    real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa
12    real, parameter:: preff = 101325. ! in Pa    real, parameter:: preff = 101325. ! in Pa
13    real nivsigs(llm), nivsig(llm+1)    real nivsigs(llm), nivsig(llm+1)
14    
# Line 22  contains Line 22  contains
22      ! Auteur : P. Le Van      ! Auteur : P. Le Van
23    
24      ! This procedure sets the vertical grid.      ! This procedure sets the vertical grid.
25      ! It defines the host variables "ap", "bp", "dpres", "presnivs",      ! It defines the host variables "ap", "bp", "presnivs", "nivsigs"
26      ! "nivsigs" and "nivsig".      ! and "nivsig".
27      ! "pa" should be defined before this procedure is called.      ! "pa" should be defined before this procedure is called.
28    
29      use comconst, only: pi      use nr_util, only: pi
30        use jumble, only: new_unit
31        use unit_nml_m, only: unit_nml
32    
33      REAL s(llm+1)      REAL s(llm+1)
34      ! (atmospheric hybrid sigma-pressure coordinate at the interface      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
35      ! between layers "l" and "l-1")      ! the interface between layers "l" and "l-1"
36    
37      real ds(llm)      real ds(llm)
38      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
39    
40      INTEGER l      INTEGER l, unit
41      REAL alpha, x(llm)      REAL alpha, x(llm)
42    
43      character(len=7):: s_sampling = "LMD5"      character(len=7):: s_sampling = "LMD5"
44      ! (other allowed values are "param", "strato1" and "strato2")      ! (other allowed values are "param", "strato1", "strato2" and "read")
45    
46      real:: h = 7. ! scale height, in km      real:: h = 7. ! scale height, in km
47      ! (used only if "s_sampling" == "param" or "strato1")      ! (used only if "s_sampling" == "param" or "strato1")
# Line 63  contains Line 65  contains
65    
66      ! Compute "s":      ! Compute "s":
67    
     s(1) = 1.  
     s(llm+1) = 0.  
   
68      print *, "Enter namelist 'disvert_nml'."      print *, "Enter namelist 'disvert_nml'."
69      read(unit=*, nml=disvert_nml)      read(unit=*, nml=disvert_nml)
70      write(unit=*, nml=disvert_nml)      write(unit_nml, nml=disvert_nml)
71    
72      select case (s_sampling)      select case (s_sampling)
73      case ("param")      case ("param")
74           s(1) = 1.
75           s(llm+1) = 0.
76         alpha = deltaz / tanh(1./k0) * 2.         alpha = deltaz / tanh(1./k0) * 2.
77         forall (l = 2: llm) s(l) &         forall (l = 2: llm) s(l) &
78              = cosh((l - 1) / k0) **(- alpha * k0 / h) &              = cosh((l - 1) / k0) **(- alpha * k0 / h) &
# Line 79  contains Line 80  contains
80              * beta **(l - 1 - (llm - k1)) / log(beta))              * beta **(l - 1 - (llm - k1)) / log(beta))
81      case ("LMD5")      case ("LMD5")
82         ! Ancienne discrétisation         ! Ancienne discrétisation
83           s(1) = 1.
84           s(llm+1) = 0.
85         forall (l = 1: llm) ds(l) &         forall (l = 1: llm) ds(l) &
86              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2
87         ds = ds / sum(ds)         ds = ds / sum(ds)
# Line 88  contains Line 91  contains
91         ENDDO         ENDDO
92      case ("strato1")      case ("strato1")
93         ! F. Lott 70 niveaux et plus         ! F. Lott 70 niveaux et plus
94           s(1) = 1.
95           s(llm+1) = 0.
96         forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &         forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &
97                 + TANH(REAL(l - llm) / 10.) / 2.                 + TANH(REAL(l - llm) / 10.) / 2.
98    
# Line 100  contains Line 105  contains
105              / (1. - exp(- zz(llm + 1) / h))              / (1. - exp(- zz(llm + 1) / h))
106      case ("strato2")      case ("strato2")
107         ! Recommended by F. Lott for a domain including the stratosphere         ! Recommended by F. Lott for a domain including the stratosphere
108           s(1) = 1.
109           s(llm+1) = 0.
110         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
111    
112         ds = (1. + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.         ds = (1. + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
# Line 108  contains Line 115  contains
115         DO l = llm, 2, -1         DO l = llm, 2, -1
116            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
117         ENDDO         ENDDO
118        case("read")
119           call new_unit(unit)
120           open(unit, file="hybrid.csv", status="old", action="read", &
121                position="rewind")
122           read(unit, fmt=*) ! skip title line
123           do l = 1, llm + 1
124              read(unit, fmt=*) s(l)
125           end do
126           close(unit)
127           ! Quick check:
128           if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then
129              print *, '"s" should be in descending order, from 1 to 0.'
130              stop 1
131           end if
132      case default      case default
133         print *, 'Wrong value for "s_sampling"'         print *, 'Wrong value for "s_sampling"'
134         stop 1         stop 1
# Line 118  contains Line 139  contains
139      bp(llm + 1) = 0.      bp(llm + 1) = 0.
140      ap = pa * (s - bp)      ap = pa * (s - bp)
141    
142      forall (l = 1: llm)      forall (l = 1: llm) presnivs(l) = 0.5 &
143         dpres(l) = bp(l) - bp(l+1)           * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
        presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)  
     end forall  
144    
145    END SUBROUTINE disvert    END SUBROUTINE disvert
146    
Legend:
Removed from v.3  
changed lines
  Added in v.57
  ViewVC Help
Powered by ViewVC 1.1.21