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

Annotation of /trunk/phylmd/Orography/gwprofil.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (hide annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/Orography/gwprofil.f90
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 guez 54 module gwprofil_m
2 guez 23
3 guez 54 IMPLICIT NONE
4 guez 23
5 guez 54 contains
6 guez 23
7 guez 54 SUBROUTINE gwprofil(nlon, nlev, kgwd, kdx, ktest, kkcrith, kcrit, paphm1, &
8     prho, pstab, pvph, pri, ptau, pdmod, psig, pvar)
9 guez 23
10 guez 54 ! 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 guez 23
18 guez 54 ! Reference.
19     ! See ECMWF research department documentation of the "I.F.S."
20 guez 23
21 guez 54 ! Modifications.
22     ! Passage of the new gwdrag TO I.F.S. (F. LOTT, 22/11/93)
23 guez 23
24 guez 54 USE dimphy, ONLY : klev, klon
25     USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra
26 guez 23
27 guez 54 ! 0.1 ARGUMENTS
28 guez 23
29 guez 54 INTEGER nlon, nlev
30     INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)
31 guez 23
32 guez 54 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 guez 23
35 guez 54 REAL pdmod(nlon)
36     REAL, INTENT (IN) :: psig(nlon)
37     REAL, INTENT (IN) :: pvar(nlon)
38 guez 23
39 guez 54 ! 0.2 LOCAL ARRAYS
40 guez 23
41 guez 54 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 guez 23
47 guez 54 !-----------------------------------------------------------------------
48 guez 23
49 guez 54 ! 1. INITIALIZATION
50 guez 23
51 guez 54 ! COMPUTATIONAL CONSTANTS.
52 guez 23
53 guez 54 ilevh = klev/3
54 guez 23
55 guez 54 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 guez 23
62 guez 54 DO jk = klev, 2, -1
63     ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
64     ! BLOCKING LAYER.
65     DO jl = 1, klon
66 guez 23 IF (ktest(jl)==1) THEN
67 guez 54 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 guez 23 END IF
73 guez 54 end DO
74 guez 23
75 guez 54 ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
76     DO jl = 1, klon
77 guez 23 IF (ktest(jl)==1) THEN
78 guez 54 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 guez 23 END IF
84 guez 54 end DO
85 guez 23
86 guez 54 ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
87     ! AND STRESS: BREAKING EVALUATION AND CRITICAL
88     ! LEVEL
89     DO jl = 1, klon
90 guez 23 IF (ktest(jl)==1) THEN
91 guez 54 IF (jk<kkcrith(jl)) THEN
92     IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
93     ptau(jl, jk) = 0.0
94 guez 23 ELSE
95 guez 54 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 guez 23 END IF
109 guez 54 END IF
110 guez 23 END IF
111 guez 54 end DO
112     end DO
113 guez 23
114 guez 54 ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
115 guez 23
116 guez 54 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 guez 23
123 guez 54 DO jk = 1, klev
124     DO jl = 1, klon
125 guez 23 IF (ktest(jl)==1) THEN
126 guez 54 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 guez 23 END IF
133 guez 54 end DO
134 guez 23
135 guez 54 ! REORGANISATION IN THE STRATOSPHERE
136     DO jl = 1, klon
137 guez 23 IF (ktest(jl)==1) THEN
138 guez 54 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 guez 23 END IF
144 guez 54 end DO
145 guez 23
146 guez 54 ! REORGANISATION IN THE TROPOSPHERE
147     DO jl = 1, klon
148 guez 23 IF (ktest(jl)==1) THEN
149 guez 54 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 guez 23 END IF
156 guez 54 end DO
157     end DO
158 guez 23
159 guez 54 END SUBROUTINE gwprofil
160 guez 23
161 guez 54 end module gwprofil_m
  ViewVC Help
Powered by ViewVC 1.1.21