/[lmdze]/trunk/Sources/phylmd/aaam_bud.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/aaam_bud.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (show annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
File size: 9629 byte(s)
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).

1 module aaam_bud_m
2
3 implicit none
4
5 contains
6
7 subroutine aaam_bud(rg, ome, phis, dragu, liftu, phyu, dragv, liftv, phyv, &
8 p, u, v, aam, torsfc)
9
10 ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
11 ! different terms of the axial AAAM budget and mountain torque.
12 ! Only valid for regular rectangular grids. Should be called after
13 ! "lift_noro".
14
15 USE dimens_m, ONLY : iim, jjm
16 use nr_util, only: assert_eq, assert, pi
17 use phyetat0_m, only: rlat, rlon
18 USE suphec_m, ONLY: ra
19
20 real, intent(in):: rg ! gravity constant
21 real, intent(in):: ome ! Earth rotation rate
22 real, intent(in):: phis(:) ! (nlon) Geopotential at the ground
23 REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal)
24 REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal)
25 REAL, intent(in):: phyu(:) ! (nlon) Stress total de la physique (zonal)
26 REAL, intent(in):: dragv(:) ! (nlon) orodrag stress (Meridional)
27 REAL, intent(in):: liftv(:) ! (nlon) orolift stress (Meridional)
28 REAL, intent(in):: phyv(:) ! (nlon) Stress total de la physique (Meridional)
29
30 REAL, intent(in):: p(:, :)
31 ! (nlon, nlev + 1) pressure (Pa) at model half levels
32
33 real, intent(in):: u(:, :), v(:, :) ! (nlon, nlev) horizontal wind (m/s)
34 REAL, intent(out):: aam ! axial component of wind AAM
35 REAL, intent(out):: torsfc ! axial component of total surface torque
36
37 ! Local Variables:
38
39 INTEGER nlev ! number of vertical levels
40 INTEGER i, j, k, l
41 REAL dlat, dlon ! latitude and longitude increments (radians)
42
43 REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial)
44 REAL oaam(3) ! mass AAM (components 1 & 2: equatorial; component 3: axial)
45 REAL tmou(3) ! resolved mountain torque (3 components)
46 REAL tsso(3) ! parameterised moutain drag torque (3 components)
47 REAL tbls(3) ! parameterised boundary layer torque (3 components)
48
49 REAL ZS(801, 401) ! topographic height
50 REAL PS(801, 401) ! surface pressure
51 REAL UB(801, 401), VB(801, 401) ! barotropic wind, zonal and meridional
52 REAL SSOU(801, 401), SSOV(801, 401)
53 REAL BLSU(801, 401), BLSV(801, 401)
54 REAL ZLON(801), ZLAT(401) ! longitude and latitude in radians
55
56 !-------------------------------------------------------------------
57
58 call assert(size(phis) == (/size(dragu), size(liftu), size(phyu), &
59 size(dragv), size(liftv), size(phyv), size(p, 1), size(u, 1), &
60 size(v, 1)/), "aaam_bud nlon")
61 nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")
62
63 if (iim + 1 > 801 .or. jjm + 1 > 401) then
64 print *, ' Problème de dimension dans aaam_bud'
65 stop 1
66 endif
67
68 dlat = pi / jjm
69 dlon = 2 * pi / real(iim)
70
71 oaam = 0.
72 raam = 0.
73 tmou = 0.
74 tsso = 0.
75 tbls = 0.
76
77 ! Mountain height, pressure and barotropic wind:
78
79 ! North pole values (j = 1):
80
81 ub(1, 1) = 0.
82 vb(1, 1) = 0.
83 do k = 1, nlev
84 ub(1, 1) = ub(1, 1) + u(1, k) * (p(1, k) - p(1, k + 1)) / rg
85 vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg
86 enddo
87
88 zlat(1) = rlat(1) * pi / 180.
89
90 do i = 1, iim + 1
91 zs(i, 1) = phis(1) / rg
92 ps(i, 1) = p(1, 1)
93 ub(i, 1) = ub(1, 1)
94 vb(i, 1) = vb(1, 1)
95 ssou(i, 1) = dragu(1) + liftu(1)
96 ssov(i, 1) = dragv(1) + liftv(1)
97 blsu(i, 1) = phyu(1) - dragu(1) - liftu(1)
98 blsv(i, 1) = phyv(1) - dragv(1) - liftv(1)
99 enddo
100
101 l = 1
102 do j = 2, jjm
103 ! Values at Greenwich (Periodicity)
104
105 zs(iim + 1, j) = phis(l + 1) / rg
106 ps(iim + 1, j) = p(l + 1, 1)
107 ssou(iim + 1, j) = dragu(l + 1) + liftu(l + 1)
108 ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1)
109 blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1)
110 blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1)
111 zlon(iim + 1) = - rlon(l + 1) * pi / 180.
112 zlat(j) = rlat(l + 1) * pi / 180.
113
114 ub(iim + 1, j) = 0.
115 vb(iim + 1, j) = 0.
116 do k = 1, nlev
117 ub(iim + 1, j) = ub(iim + 1, j) &
118 + u(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
119 vb(iim + 1, j) = vb(iim + 1, j) &
120 + v(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
121 enddo
122
123 do i = 1, iim
124 l = l + 1
125 zs(i, j) = phis(l) / rg
126 ps(i, j) = p(l, 1)
127 ssou(i, j) = dragu(l) + liftu(l)
128 ssov(i, j) = dragv(l) + liftv(l)
129 blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
130 blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
131 zlon(i) = rlon(l) * pi / 180.
132
133 ub(i, j) = 0.
134 vb(i, j) = 0.
135 do k = 1, nlev
136 ub(i, j) = ub(i, j) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
137 vb(i, j) = vb(i, j) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
138 enddo
139 enddo
140 enddo
141
142 ! South Pole
143
144 l = l + 1
145 ub(1, jjm + 1) = 0.
146 vb(1, jjm + 1) = 0.
147 do k = 1, nlev
148 ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
149 vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
150 enddo
151 zlat(jjm + 1) = rlat(l) * pi / 180.
152
153 do i = 1, iim + 1
154 zs(i, jjm + 1) = phis(l) / rg
155 ps(i, jjm + 1) = p(l, 1)
156 ssou(i, jjm + 1) = dragu(l) + liftu(l)
157 ssov(i, jjm + 1) = dragv(l) + liftv(l)
158 blsu(i, jjm + 1) = phyu(l) - dragu(l) - liftu(l)
159 blsv(i, jjm + 1) = phyv(l) - dragv(l) - liftv(l)
160 ub(i, jjm + 1) = ub(1, jjm + 1)
161 vb(i, jjm + 1) = vb(1, jjm + 1)
162 enddo
163
164 ! Moment angulaire
165
166 DO j = 1, jjm
167 DO i = 1, iim
168 raam(1) = raam(1) - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) &
169 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) &
170 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
171 + ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
172 * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
173
174 oaam(1) = oaam(1) - ome * ra**4 * dlon * dlat / rg * 0.5 &
175 * (cos(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
176 * ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 &
177 * sin(zlat(j + 1)) * ps(i , j + 1))
178
179 raam(2) = raam(2) - ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) &
180 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) &
181 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
182 - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
183 * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
184
185 oaam(2) = oaam(2) - ome * ra**4 * dlon * dlat / rg * 0.5 &
186 * (sin(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
187 * ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 &
188 * sin(zlat(j + 1)) * ps(i , j + 1))
189
190 raam(3) = raam(3) + ra**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 &
191 * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1))
192
193 oaam(3) = oaam(3) + ome * ra**4 * dlon * dlat / rg * 0.5 &
194 * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 &
195 * ps(i, j + 1))
196 ENDDO
197 ENDDO
198
199 ! Couple des montagnes :
200
201 DO j = 1, jjm
202 DO i = 1, iim
203 tmou(1) = tmou(1) - ra**2 * dlon * 0.5 * sin(zlon(i)) &
204 * (zs(i, j) - zs(i, j + 1)) &
205 * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
206 tmou(2) = tmou(2) + ra**2 * dlon * 0.5 * cos(zlon(i)) &
207 * (zs(i, j) - zs(i, j + 1)) &
208 * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
209 ENDDO
210 ENDDO
211
212 DO j = 2, jjm
213 DO i = 1, iim
214 tmou(1) = tmou(1) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
215 * (zs(i + 1, j) - zs(i, j)) &
216 * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
217 tmou(2) = tmou(2) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
218 * (zs(i + 1, j) - zs(i, j)) &
219 * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
220 tmou(3) = tmou(3) - ra**2 * dlat * 0.5* cos(zlat(j)) &
221 * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
222 ENDDO
223 ENDDO
224
225 ! Couples des differentes friction au sol :
226
227 DO j = 2, jjm
228 DO i = 1, iim
229 tsso(1) = tsso(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
230 ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
231 + ra**3 * cos(zlat(j)) * dlon * dlat* &
232 ssov(i, j) * sin(zlon(i))
233
234 tsso(2) = tsso(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
235 ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
236 - ra**3 * cos(zlat(j)) * dlon * dlat* &
237 ssov(i, j) * cos(zlon(i))
238
239 tsso(3) = tsso(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
240 ssou(i, j) * cos(zlat(j))
241
242 tbls(1) = tbls(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
243 blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
244 + ra**3 * cos(zlat(j)) * dlon * dlat* &
245 blsv(i, j) * sin(zlon(i))
246
247 tbls(2) = tbls(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
248 blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
249 - ra**3 * cos(zlat(j)) * dlon * dlat* &
250 blsv(i, j) * cos(zlon(i))
251
252 tbls(3) = tbls(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
253 blsu(i, j) * cos(zlat(j))
254 ENDDO
255 ENDDO
256
257 aam = raam(3)
258 torsfc = tmou(3) + tsso(3) + tbls(3)
259
260 END subroutine aaam_bud
261
262 end module aaam_bud_m

  ViewVC Help
Powered by ViewVC 1.1.21