/[lmdze]/trunk/libf/phylmd/Orography/gwprofil.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/Orography/gwprofil.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (show annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
File size: 5209 byte(s) 
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".
1 module gwprofil_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE gwprofil(nlon, nlev, kgwd, kdx, ktest, kkcrith, kcrit, paphm1, &
8 prho, pstab, pvph, pri, ptau, pdmod, psig, pvar)
9
10 ! Method. The stress profile for gravity waves is computed as
11 ! follows: it is constant (no gwd) at the levels between the
12 ! ground and the top of the blocked layer (kkenvh). It decreases
13 ! linearly with height from the top of the blocked layer to
14 ! 3*varor (kknu), to simulate lee waves or nonlinear gravity wave
15 ! breaking. Above it is constant, except when the wave encounters
16 ! a critical level (kcrit) or when it breaks.
17
18 ! Reference.
19 ! See ECMWF research department documentation of the "I.F.S."
20
21 ! Modifications.
22 ! Passage of the new gwdrag TO I.F.S. (F. LOTT, 22/11/93)
23
24 USE dimphy, ONLY : klev, klon
25 USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra
26
27 ! 0.1 ARGUMENTS
28
29 INTEGER nlon, nlev
30 INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)
31
32 REAL paphm1(nlon, nlev+1), pstab(nlon, nlev+1), prho(nlon, nlev+1), &
33 pvph(nlon, nlev+1), pri(nlon, nlev+1), ptau(nlon, nlev+1)
34
35 REAL pdmod(nlon)
36 REAL, INTENT (IN) :: psig(nlon)
37 REAL, INTENT (IN) :: pvar(nlon)
38
39 ! 0.2 LOCAL ARRAYS
40
41 INTEGER ilevh, ji, kgwd, jl, jk
42 REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
43 REAL zdelp, zdelpt
44 REAL zdz2(klon, klev), znorm(klon), zoro(klon)
45 REAL ztau(klon, klev+1)
46
47 !-----------------------------------------------------------------------
48
49 ! 1. INITIALIZATION
50
51 ! COMPUTATIONAL CONSTANTS.
52
53 ilevh = klev/3
54
55 DO jl = 1, klon
56 IF (ktest(jl)==1) THEN
57 zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
58 ztau(jl, klev+1) = ptau(jl, klev+1)
59 END IF
60 end DO
61
62 DO jk = klev, 2, -1
63 ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
64 ! BLOCKING LAYER.
65 DO jl = 1, klon
66 IF (ktest(jl)==1) THEN
67 IF (jk>kkcrith(jl)) THEN
68 ptau(jl, jk) = ztau(jl, klev+1)
69 ELSE
70 ptau(jl, jk) = grahilo*ztau(jl, klev+1)
71 END IF
72 END IF
73 end DO
74
75 ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
76 DO jl = 1, klon
77 IF (ktest(jl)==1) THEN
78 IF (jk<kkcrith(jl)) THEN
79 znorm(jl) = gkdrag * prho(jl, jk) * sqrt(pstab(jl, jk)) &
80 * pvph(jl, jk)* zoro(jl)
81 zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
82 END IF
83 END IF
84 end DO
85
86 ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
87 ! AND STRESS: BREAKING EVALUATION AND CRITICAL
88 ! LEVEL
89 DO jl = 1, klon
90 IF (ktest(jl)==1) THEN
91 IF (jk<kkcrith(jl)) THEN
92 IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
93 ptau(jl, jk) = 0.0
94 ELSE
95 zsqr = sqrt(pri(jl, jk))
96 zalfa = sqrt(pstab(jl, jk)*zdz2(jl, jk))/pvph(jl, jk)
97 zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
98 IF (zriw<grcrit) THEN
99 zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
100 zb = 1./grcrit + 2./zsqr
101 zalpha = 0.5*(-zb+sqrt(zdel))
102 zdz2n = (pvph(jl, jk)*zalpha)**2/pstab(jl, jk)
103 ptau(jl, jk) = znorm(jl)*zdz2n
104 ELSE
105 ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
106 END IF
107 ptau(jl, jk) = min(ptau(jl, jk), ptau(jl, jk+1))
108 END IF
109 END IF
110 END IF
111 end DO
112 end DO
113
114 ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
115
116 DO jl = 1, klon
117 IF (ktest(jl)==1) THEN
118 ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
119 ztau(jl, nstra) = ptau(jl, nstra)
120 END IF
121 end DO
122
123 DO jk = 1, klev
124 DO jl = 1, klon
125 IF (ktest(jl)==1) THEN
126 IF (jk>kkcrith(jl)) THEN
127 zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
128 zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
129 ptau(jl, jk) = ztau(jl, klev+1) &
130 + (ztau(jl, kkcrith(jl)) - ztau(jl, klev+1))*zdelp/zdelpt
131 END IF
132 END IF
133 end DO
134
135 ! REORGANISATION IN THE STRATOSPHERE
136 DO jl = 1, klon
137 IF (ktest(jl)==1) THEN
138 IF (jk < nstra) THEN
139 zdelp = paphm1(jl, nstra)
140 zdelpt = paphm1(jl, jk)
141 ptau(jl, jk) = ztau(jl, nstra) * zdelpt / zdelp
142 END IF
143 END IF
144 end DO
145
146 ! REORGANISATION IN THE TROPOSPHERE
147 DO jl = 1, klon
148 IF (ktest(jl)==1) THEN
149 IF (jk<kkcrith(jl) .AND. jk > nstra) THEN
150 zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
151 zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
152 ptau(jl, jk) = ztau(jl, kkcrith(jl)) &
153 + (ztau(jl, nstra) - ztau(jl, kkcrith(jl)))*zdelp/zdelpt
154 END IF
155 END IF
156 end DO
157 end DO
158
159 END SUBROUTINE gwprofil
160
161 end module gwprofil_m
  ViewVC Help
Powered by ViewVC 1.1.21