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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 178 - (hide annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 2 months ago) by guez
File size: 9717 byte(s) 
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

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