/[lmdze]/trunk/phylmd/aaam_bud.f90
ViewVC logotype

Annotation of /trunk/phylmd/aaam_bud.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (hide annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
Original Path: trunk/Sources/phylmd/aaam_bud.f
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 guez 56 module aaam_bud_m
2 guez 3
3 guez 56 implicit none
4 guez 3
5 guez 56 contains
6 guez 3
7 guez 227 subroutine aaam_bud(rg, ome, phis, dragu, liftu, phyu, dragv, liftv, phyv, &
8     p, u, v, aam, torsfc)
9 guez 3
10 guez 56 ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
11 guez 61 ! different terms of the axial AAAM budget and mountain torque.
12 guez 56 ! Only valid for regular rectangular grids. Should be called after
13     ! "lift_noro".
14 guez 3
15 guez 56 USE dimens_m, ONLY : iim, jjm
16     use nr_util, only: assert_eq, assert, pi
17 guez 227 use phyetat0_m, only: rlat, rlon
18 guez 171 USE suphec_m, ONLY: ra
19 guez 3
20 guez 56 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 guez 3
30 guez 56 REAL, intent(in):: p(:, :)
31     ! (nlon, nlev + 1) pressure (Pa) at model half levels
32 guez 3
33 guez 56 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 guez 3
37 guez 56 ! Local Variables:
38 guez 3
39 guez 56 INTEGER nlev ! number of vertical levels
40     INTEGER i, j, k, l
41     REAL dlat, dlon ! latitude and longitude increments (radians)
42 guez 3
43 guez 56 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 guez 3
49 guez 56 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 guez 3
56 guez 56 !-------------------------------------------------------------------
57 guez 3
58 guez 227 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 guez 56 nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")
62 guez 3
63 guez 56 if (iim + 1 > 801 .or. jjm + 1 > 401) then
64     print *, ' Problème de dimension dans aaam_bud'
65     stop 1
66     endif
67 guez 3
68 guez 56 dlat = pi / jjm
69     dlon = 2 * pi / real(iim)
70 guez 3
71 guez 56 oaam = 0.
72     raam = 0.
73     tmou = 0.
74     tsso = 0.
75     tbls = 0.
76 guez 3
77 guez 56 ! Mountain height, pressure and barotropic wind:
78 guez 3
79 guez 56 ! North pole values (j = 1):
80 guez 3
81 guez 56 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 guez 3
88 guez 227 zlat(1) = rlat(1) * pi / 180.
89 guez 3
90 guez 56 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 guez 3
101 guez 56 l = 1
102     do j = 2, jjm
103     ! Values at Greenwich (Periodicity)
104 guez 3
105 guez 56 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 guez 227 zlon(iim + 1) = - rlon(l + 1) * pi / 180.
112     zlat(j) = rlat(l + 1) * pi / 180.
113 guez 3
114 guez 56 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 guez 3
123 guez 56 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 guez 227 zlon(i) = rlon(l) * pi / 180.
132 guez 3
133 guez 56 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 guez 3
142 guez 56 ! South Pole
143 guez 3
144 guez 56 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 guez 227 zlat(jjm + 1) = rlat(l) * pi / 180.
152 guez 3
153 guez 56 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 guez 3
164 guez 56 ! Moment angulaire
165 guez 3
166 guez 56 DO j = 1, jjm
167     DO i = 1, iim
168 guez 171 raam(1) = raam(1) - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) &
169 guez 56 * 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 guez 171 + ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
172 guez 56 * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
173 guez 3
174 guez 171 oaam(1) = oaam(1) - ome * ra**4 * dlon * dlat / rg * 0.5 &
175 guez 56 * (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 guez 3
179 guez 171 raam(2) = raam(2) - ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) &
180 guez 56 * 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 guez 171 - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
183 guez 56 * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
184 guez 3
185 guez 171 oaam(2) = oaam(2) - ome * ra**4 * dlon * dlat / rg * 0.5 &
186 guez 56 * (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 guez 3
190 guez 171 raam(3) = raam(3) + ra**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 &
191 guez 56 * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1))
192 guez 3
193 guez 171 oaam(3) = oaam(3) + ome * ra**4 * dlon * dlat / rg * 0.5 &
194 guez 56 * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 &
195     * ps(i, j + 1))
196     ENDDO
197     ENDDO
198 guez 3
199 guez 56 ! Couple des montagnes :
200 guez 3
201 guez 56 DO j = 1, jjm
202     DO i = 1, iim
203 guez 171 tmou(1) = tmou(1) - ra**2 * dlon * 0.5 * sin(zlon(i)) &
204 guez 56 * (zs(i, j) - zs(i, j + 1)) &
205     * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
206 guez 171 tmou(2) = tmou(2) + ra**2 * dlon * 0.5 * cos(zlon(i)) &
207 guez 56 * (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 guez 3
212 guez 56 DO j = 2, jjm
213     DO i = 1, iim
214 guez 171 tmou(1) = tmou(1) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
215 guez 56 * (zs(i + 1, j) - zs(i, j)) &
216     * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
217 guez 171 tmou(2) = tmou(2) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
218 guez 56 * (zs(i + 1, j) - zs(i, j)) &
219     * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
220 guez 171 tmou(3) = tmou(3) - ra**2 * dlat * 0.5* cos(zlat(j)) &
221 guez 56 * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
222     ENDDO
223     ENDDO
224 guez 3
225 guez 56 ! Couples des differentes friction au sol :
226    
227     DO j = 2, jjm
228     DO i = 1, iim
229 guez 171 tsso(1) = tsso(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
230 guez 56 ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
231 guez 171 + ra**3 * cos(zlat(j)) * dlon * dlat* &
232 guez 56 ssov(i, j) * sin(zlon(i))
233    
234 guez 171 tsso(2) = tsso(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
235 guez 56 ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
236 guez 171 - ra**3 * cos(zlat(j)) * dlon * dlat* &
237 guez 56 ssov(i, j) * cos(zlon(i))
238    
239 guez 171 tsso(3) = tsso(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
240 guez 56 ssou(i, j) * cos(zlat(j))
241    
242 guez 171 tbls(1) = tbls(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
243 guez 56 blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
244 guez 171 + ra**3 * cos(zlat(j)) * dlon * dlat* &
245 guez 56 blsv(i, j) * sin(zlon(i))
246    
247 guez 171 tbls(2) = tbls(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
248 guez 56 blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
249 guez 171 - ra**3 * cos(zlat(j)) * dlon * dlat* &
250 guez 56 blsv(i, j) * cos(zlon(i))
251    
252 guez 171 tbls(3) = tbls(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
253 guez 56 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