/[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 36 by guez, Thu Dec 2 17:11:04 2010 UTC revision 42 by guez, Thu Mar 24 11:52:41 2011 UTC
# Line 27  contains Line 27  contains
27      ! "pa" should be defined before this procedure is called.      ! "pa" should be defined before this procedure is called.
28    
29      use nr_util, only: pi      use nr_util, only: pi
30        use new_unit_m, only: new_unit
31    
32      REAL s(llm+1)      REAL s(llm+1)
33      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at      ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
# Line 35  contains Line 36  contains
36      real ds(llm)      real ds(llm)
37      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"      ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
38    
39      INTEGER l      INTEGER l, unit
40      REAL alpha, x(llm)      REAL alpha, x(llm), trash
41    
42      character(len=7):: s_sampling = "LMD5"      character(len=7):: s_sampling = "LMD5"
43      ! (other allowed values are "param", "strato1" and "strato2")      ! (other allowed values are "param", "strato1", "strato2" and "read")
44    
45      real:: h = 7. ! scale height, in km      real:: h = 7. ! scale height, in km
46      ! (used only if "s_sampling" == "param" or "strato1")      ! (used only if "s_sampling" == "param" or "strato1")
# Line 63  contains Line 64  contains
64    
65      ! Compute "s":      ! Compute "s":
66    
     s(1) = 1.  
     s(llm+1) = 0.  
   
67      print *, "Enter namelist 'disvert_nml'."      print *, "Enter namelist 'disvert_nml'."
68      read(unit=*, nml=disvert_nml)      read(unit=*, nml=disvert_nml)
69      write(unit=*, nml=disvert_nml)      write(unit=*, nml=disvert_nml)
70    
71      select case (s_sampling)      select case (s_sampling)
72      case ("param")      case ("param")
73           s(1) = 1.
74           s(llm+1) = 0.
75         alpha = deltaz / tanh(1./k0) * 2.         alpha = deltaz / tanh(1./k0) * 2.
76         forall (l = 2: llm) s(l) &         forall (l = 2: llm) s(l) &
77              = cosh((l - 1) / k0) **(- alpha * k0 / h) &              = cosh((l - 1) / k0) **(- alpha * k0 / h) &
# Line 79  contains Line 79  contains
79              * beta **(l - 1 - (llm - k1)) / log(beta))              * beta **(l - 1 - (llm - k1)) / log(beta))
80      case ("LMD5")      case ("LMD5")
81         ! Ancienne discrétisation         ! Ancienne discrétisation
82           s(1) = 1.
83           s(llm+1) = 0.
84         forall (l = 1: llm) ds(l) &         forall (l = 1: llm) ds(l) &
85              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2              = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2
86         ds = ds / sum(ds)         ds = ds / sum(ds)
# Line 88  contains Line 90  contains
90         ENDDO         ENDDO
91      case ("strato1")      case ("strato1")
92         ! F. Lott 70 niveaux et plus         ! F. Lott 70 niveaux et plus
93           s(1) = 1.
94           s(llm+1) = 0.
95         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.) &
96                 + TANH(REAL(l - llm) / 10.) / 2.                 + TANH(REAL(l - llm) / 10.) / 2.
97    
# Line 100  contains Line 104  contains
104              / (1. - exp(- zz(llm + 1) / h))              / (1. - exp(- zz(llm + 1) / h))
105      case ("strato2")      case ("strato2")
106         ! Recommended by F. Lott for a domain including the stratosphere         ! Recommended by F. Lott for a domain including the stratosphere
107           s(1) = 1.
108           s(llm+1) = 0.
109         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)         forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
110    
111         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 114  contains
114         DO l = llm, 2, -1         DO l = llm, 2, -1
115            s(l) = s(l+1) + ds(l)            s(l) = s(l+1) + ds(l)
116         ENDDO         ENDDO
117        case("read")
118           call new_unit(unit)
119           open(unit, file="hybrid.csv", status="old", action="read", &
120                position="rewind")
121           read(unit, fmt=*) ! skip title line
122           do l = 1, llm + 1
123              read(unit, fmt=*) trash, s(l)
124           end do
125           close(unit)
126           ! Quick check:
127           if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then
128              print *, '"s" should be in descending order, from 1 to 0.'
129              stop 1
130           end if
131      case default      case default
132         print *, 'Wrong value for "s_sampling"'         print *, 'Wrong value for "s_sampling"'
133         stop 1         stop 1
Legend:
Removed from v.36  
changed lines
  Added in v.42
  ViewVC Help
Powered by ViewVC 1.1.21