GCC Code Coverage Report

Directory: ./
File: phys/aaam_bud.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 90 92 97.8%
Branches: 34 38 89.5%
Line Branch Exec Source
1
2 ! $Id: aaam_bud.F90 2350 2015-08-25 11:40:19Z emillour $
3
4 480 SUBROUTINE aaam_bud(iam, nlon, nlev, rjour, rsec, rea, rg, ome, plat, plon, &
5 480 phis, dragu, liftu, phyu, dragv, liftv, phyv, p, u, v, aam, torsfc)
6
7 USE dimphy
8 USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo
9 IMPLICIT NONE
10 ! ======================================================================
11 ! Auteur(s): F.Lott (LMD/CNRS) date: 20031020
12 ! Object: Compute different terms of the axial AAAM Budget.
13 ! No outputs, every AAM quantities are written on the IAM
14 ! File.
15
16 ! Modif : I.Musat (LMD/CNRS) date : 20041020
17 ! Outputs : axial components of wind AAM "aam" and total surface torque
18 ! "torsfc",
19 ! but no write in the iam file.
20
21 ! WARNING: Only valid for regular rectangular grids.
22 ! REMARK: CALL DANS PHYSIQ AFTER lift_noro:
23 ! CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
24 ! C ra,rg,romega,
25 ! C rlat,rlon,pphis,
26 ! C zustrdr,zustrli,zustrph,
27 ! C zvstrdr,zvstrli,zvstrph,
28 ! C paprs,u,v)
29
30 ! ======================================================================
31 ! Explicit Arguments:
32 ! ==================
33 ! iam-----input-I-File number where AAMs and torques are written
34 ! It is a formatted file that has been opened
35 ! in physiq.F
36 ! nlon----input-I-Total number of horizontal points that get into physics
37 ! nlev----input-I-Number of vertical levels
38 ! rjour -R-Jour compte depuis le debut de la simu (run.def)
39 ! rsec -R-Seconde de la journee
40 ! rea -R-Earth radius
41 ! rg -R-gravity constant
42 ! ome -R-Earth rotation rate
43 ! plat ---input-R-Latitude en degres
44 ! plon ---input-R-Longitude en degres
45 ! phis ---input-R-Geopotential at the ground
46 ! dragu---input-R-orodrag stress (zonal)
47 ! liftu---input-R-orolift stress (zonal)
48 ! phyu----input-R-Stress total de la physique (zonal)
49 ! dragv---input-R-orodrag stress (Meridional)
50 ! liftv---input-R-orolift stress (Meridional)
51 ! phyv----input-R-Stress total de la physique (Meridional)
52 ! p-------input-R-Pressure (Pa) at model half levels
53 ! u-------input-R-Horizontal wind (m/s)
54 ! v-------input-R-Meridional wind (m/s)
55 ! aam-----output-R-Axial Wind AAM (=raam(3))
56 ! torsfc--output-R-Total surface torque (=tmou(3)+tsso(3)+tbls(3))
57
58 ! Implicit Arguments:
59 ! ===================
60
61 ! nbp_lon--common-I: Number of longitude intervals
62 ! (nbp_lat-1)--common-I: Number of latitude intervals
63 ! klon-common-I: Number of points seen by the physics
64 ! nbp_lon*(nbp_lat-2)+2 for instance
65 ! klev-common-I: Number of vertical layers
66 ! ======================================================================
67 ! Local Variables:
68 ! ================
69 ! dlat-----R: Latitude increment (Radians)
70 ! dlon-----R: Longitude increment (Radians)
71 ! raam ---R: Wind AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
72 ! oaam ---R: Mass AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
73 ! tmou-----R: Resolved Mountain torque (3 components)
74 ! tsso-----R: Parameterised Moutain drag torque (3 components)
75 ! tbls-----R: Parameterised Boundary layer torque (3 components)
76
77 ! LOCAL ARRAY:
78 ! ===========
79 ! zs ---R: Topographic height
80 ! ps ---R: Surface Pressure
81 ! ub ---R: Barotropic wind zonal
82 ! vb ---R: Barotropic wind meridional
83 ! zlat ---R: Latitude in radians
84 ! zlon ---R: Longitude in radians
85 ! ======================================================================
86
87 ! ARGUMENTS
88
89 INTEGER iam, nlon, nlev
90 REAL, INTENT (IN) :: rjour, rsec, rea, rg, ome
91 REAL plat(nlon), plon(nlon), phis(nlon)
92 REAL dragu(nlon), liftu(nlon), phyu(nlon)
93 REAL dragv(nlon), liftv(nlon), phyv(nlon)
94 REAL p(nlon, nlev+1), u(nlon, nlev), v(nlon, nlev)
95
96 ! Variables locales:
97
98 INTEGER i, j, k, l
99 REAL xpi, hadley, hadday
100 REAL dlat, dlon
101 REAL raam(3), oaam(3), tmou(3), tsso(3), tbls(3)
102 INTEGER iax
103 ! IM ajout aam, torsfc
104 ! aam = composante axiale du Wind AAM raam
105 ! torsfc = composante axiale de (tmou+tsso+tbls)
106 REAL aam, torsfc
107
108 REAL zs(801, 401), ps(801, 401)
109 REAL ub(801, 401), vb(801, 401)
110 REAL ssou(801, 401), ssov(801, 401)
111 REAL blsu(801, 401), blsv(801, 401)
112 REAL zlon(801), zlat(401)
113
114 CHARACTER (LEN=20) :: modname = 'aaam_bud'
115 CHARACTER (LEN=80) :: abort_message
116
117
118
119 ! PUT AAM QUANTITIES AT ZERO:
120
121
2/4
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 480 times.
 480 IF (nbp_lon+1>801 .OR. nbp_lat>401) THEN
122 ✗ abort_message = 'Pb de dimension dans aaam_bud'
123 ✗ CALL abort_physic(modname, abort_message, 1)
124 END IF
125
126 xpi = acos(-1.)
127 hadley = 1.E18
128 hadday = 1.E18*24.*3600.
129
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
 480 IF(klon_glo.EQ.1) THEN
130 dlat = xpi
131 ELSE
132 480 dlat = xpi/real(nbp_lat-1)
133 ENDIF
134 480 dlon = 2.*xpi/real(nbp_lon)
135
136
2/2
✓ Branch 0 taken 1440 times.
✓ Branch 1 taken 480 times.
 1920 DO iax = 1, 3
137 1440 oaam(iax) = 0.
138 1440 raam(iax) = 0.
139 1440 tmou(iax) = 0.
140 1440 tsso(iax) = 0.
141 1920 tbls(iax) = 0.
142 END DO
143
144 ! MOUNTAIN HEIGHT, PRESSURE AND BAROTROPIC WIND:
145
146 ! North pole values (j=1):
147
148 l = 1
149
150 480 ub(1, 1) = 0.
151 480 vb(1, 1) = 0.
152
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO k = 1, nlev
153 18720 ub(1, 1) = ub(1, 1) + u(l, k)*(p(l,k)-p(l,k+1))/rg
154 19200 vb(1, 1) = vb(1, 1) + v(l, k)*(p(l,k)-p(l,k+1))/rg
155 END DO
156
157 480 zlat(1) = plat(l)*xpi/180.
158
159
2/2
✓ Branch 0 taken 15840 times.
✓ Branch 1 taken 480 times.
 16320 DO i = 1, nbp_lon + 1
160
161 15840 zs(i, 1) = phis(l)/rg
162 15840 ps(i, 1) = p(l, 1)
163 15840 ub(i, 1) = ub(1, 1)
164 15840 vb(i, 1) = vb(1, 1)
165 15840 ssou(i, 1) = dragu(l) + liftu(l)
166 15840 ssov(i, 1) = dragv(l) + liftv(l)
167 15840 blsu(i, 1) = phyu(l) - dragu(l) - liftu(l)
168 16320 blsv(i, 1) = phyv(l) - dragv(l) - liftv(l)
169
170 END DO
171
172
173
2/2
✓ Branch 0 taken 14880 times.
✓ Branch 1 taken 480 times.
 15360 DO j = 2, nbp_lat-1
174
175 ! Values at Greenwich (Periodicity)
176
177 14880 zs(nbp_lon+1, j) = phis(l+1)/rg
178 14880 ps(nbp_lon+1, j) = p(l+1, 1)
179 14880 ssou(nbp_lon+1, j) = dragu(l+1) + liftu(l+1)
180 14880 ssov(nbp_lon+1, j) = dragv(l+1) + liftv(l+1)
181 14880 blsu(nbp_lon+1, j) = phyu(l+1) - dragu(l+1) - liftu(l+1)
182 14880 blsv(nbp_lon+1, j) = phyv(l+1) - dragv(l+1) - liftv(l+1)
183 zlon(nbp_lon+1) = -plon(l+1)*xpi/180.
184 14880 zlat(j) = plat(l+1)*xpi/180.
185
186 14880 ub(nbp_lon+1, j) = 0.
187 14880 vb(nbp_lon+1, j) = 0.
188
2/2
✓ Branch 0 taken 580320 times.
✓ Branch 1 taken 14880 times.
 595200 DO k = 1, nlev
189 580320 ub(nbp_lon+1, j) = ub(nbp_lon+1, j) + u(l+1, k)*(p(l+1,k)-p(l+1,k+1))/rg
190 595200 vb(nbp_lon+1, j) = vb(nbp_lon+1, j) + v(l+1, k)*(p(l+1,k)-p(l+1,k+1))/rg
191 END DO
192
193
194
2/2
✓ Branch 0 taken 476160 times.
✓ Branch 1 taken 14880 times.
 491520 DO i = 1, nbp_lon
195
196 476160 l = l + 1
197 476160 zs(i, j) = phis(l)/rg
198 476160 ps(i, j) = p(l, 1)
199 476160 ssou(i, j) = dragu(l) + liftu(l)
200 476160 ssov(i, j) = dragv(l) + liftv(l)
201 476160 blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
202 476160 blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
203 zlon(i) = plon(l)*xpi/180.
204
205 476160 ub(i, j) = 0.
206 476160 vb(i, j) = 0.
207
2/2
✓ Branch 0 taken 18570240 times.
✓ Branch 1 taken 476160 times.
 19061280 DO k = 1, nlev
208 18570240 ub(i, j) = ub(i, j) + u(l, k)*(p(l,k)-p(l,k+1))/rg
209 19046400 vb(i, j) = vb(i, j) + v(l, k)*(p(l,k)-p(l,k+1))/rg
210 END DO
211
212 END DO
213
214 END DO
215
216
217 ! South Pole
218
219
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
 480 IF (nbp_lat-1>1) THEN
220 480 l = l + 1
221 480 ub(1, nbp_lat) = 0.
222 480 vb(1, nbp_lat) = 0.
223
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO k = 1, nlev
224 18720 ub(1, nbp_lat) = ub(1, nbp_lat) + u(l, k)*(p(l,k)-p(l,k+1))/rg
225 19200 vb(1, nbp_lat) = vb(1, nbp_lat) + v(l, k)*(p(l,k)-p(l,k+1))/rg
226 END DO
227 480 zlat(nbp_lat) = plat(l)*xpi/180.
228
229
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 15840 times.
 16320 DO i = 1, nbp_lon + 1
230 15840 zs(i, nbp_lat) = phis(l)/rg
231 15840 ps(i, nbp_lat) = p(l, 1)
232 15840 ssou(i, nbp_lat) = dragu(l) + liftu(l)
233 15840 ssov(i, nbp_lat) = dragv(l) + liftv(l)
234 15840 blsu(i, nbp_lat) = phyu(l) - dragu(l) - liftu(l)
235 15840 blsv(i, nbp_lat) = phyv(l) - dragv(l) - liftv(l)
236 15840 ub(i, nbp_lat) = ub(1, nbp_lat)
237 16320 vb(i, nbp_lat) = vb(1, nbp_lat)
238 END DO
239 END IF
240
241
242 ! MOMENT ANGULAIRE
243
244
2/2
✓ Branch 0 taken 15360 times.
✓ Branch 1 taken 480 times.
 15840 DO j = 1, nbp_lat-1
245
2/2
✓ Branch 0 taken 491520 times.
✓ Branch 1 taken 15360 times.
 507360 DO i = 1, nbp_lon
246
247 raam(1) = raam(1) - rea**3*dlon*dlat*0.5*(cos(zlon(i))*sin(zlat(j))*cos &
248 (zlat(j))*ub(i,j)+cos(zlon(i))*sin(zlat(j+1))*cos(zlat(j+ &
249 1))*ub(i,j+1)) + rea**3*dlon*dlat*0.5*(sin(zlon(i))*cos(zlat(j))*vb(i &
250 491520 ,j)+sin(zlon(i))*cos(zlat(j+1))*vb(i,j+1))
251
252 oaam(1) = oaam(1) - ome*rea**4*dlon*dlat/rg*0.5*(cos(zlon(i))*cos(zlat( &
253 j))**2*sin(zlat(j))*ps(i,j)+cos(zlon(i))*cos(zlat(j+ &
254 491520 1))**2*sin(zlat(j+1))*ps(i,j+1))
255
256 raam(2) = raam(2) - rea**3*dlon*dlat*0.5*(sin(zlon(i))*sin(zlat(j))*cos &
257 (zlat(j))*ub(i,j)+sin(zlon(i))*sin(zlat(j+1))*cos(zlat(j+ &
258 1))*ub(i,j+1)) - rea**3*dlon*dlat*0.5*(cos(zlon(i))*cos(zlat(j))*vb(i &
259 ,j)+cos(zlon(i))*cos(zlat(j+1))*vb(i,j+1))
260
261 oaam(2) = oaam(2) - ome*rea**4*dlon*dlat/rg*0.5*(sin(zlon(i))*cos(zlat( &
262 j))**2*sin(zlat(j))*ps(i,j)+sin(zlon(i))*cos(zlat(j+ &
263 1))**2*sin(zlat(j+1))*ps(i,j+1))
264
265 raam(3) = raam(3) + rea**3*dlon*dlat*0.5*(cos(zlat(j))**2*ub(i,j)+cos( &
266 491520 zlat(j+1))**2*ub(i,j+1))
267
268 oaam(3) = oaam(3) + ome*rea**4*dlon*dlat/rg*0.5*(cos(zlat(j))**3*ps(i,j &
269 15360 )+cos(zlat(j+1))**3*ps(i,j+1))
270
271 END DO
272 END DO
273
274
275 ! COUPLE DES MONTAGNES:
276
277
278 DO j = 1, nbp_lat-1
279 480 DO i = 1, nbp_lon
280 tmou(1) = tmou(1) - rea**2*dlon*0.5*sin(zlon(i))*(zs(i,j)-zs(i,j+1))*( &
281 cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
282 tmou(2) = tmou(2) + rea**2*dlon*0.5*cos(zlon(i))*(zs(i,j)-zs(i,j+1))*( &
283 cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
284 END DO
285 END DO
286
287
2/2
✓ Branch 0 taken 14880 times.
✓ Branch 1 taken 480 times.
 15360 DO j = 2, nbp_lat-1
288
2/2
✓ Branch 0 taken 476160 times.
✓ Branch 1 taken 14880 times.
 491520 DO i = 1, nbp_lon
289 tmou(1) = tmou(1) + rea**2*dlat*0.5*sin(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
290 476160 cos(zlon(i+1))*ps(i+1,j)+cos(zlon(i))*ps(i,j))
291 tmou(2) = tmou(2) + rea**2*dlat*0.5*sin(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
292 sin(zlon(i+1))*ps(i+1,j)+sin(zlon(i))*ps(i,j))
293 tmou(3) = tmou(3) - rea**2*dlat*0.5*cos(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
294 491040 ps(i+1,j)+ps(i,j))
295 END DO
296 END DO
297
298
299 ! COUPLES DES DIFFERENTES FRICTION AU SOL:
300
301 l = 1
302
2/2
✓ Branch 0 taken 14880 times.
✓ Branch 1 taken 480 times.
 15360 DO j = 2, nbp_lat-1
303
2/2
✓ Branch 0 taken 476160 times.
✓ Branch 1 taken 14880 times.
 491520 DO i = 1, nbp_lon
304 l = l + 1
305 tsso(1) = tsso(1) - rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*sin(zlat(j &
306 ))*cos(zlon(i)) + rea**3*cos(zlat(j))*dlon*dlat*ssov(i, j)*sin(zlon(i &
307 476160 ))
308
309 tsso(2) = tsso(2) - rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*sin(zlat(j &
310 ))*sin(zlon(i)) - rea**3*cos(zlat(j))*dlon*dlat*ssov(i, j)*cos(zlon(i &
311 ))
312
313 tsso(3) = tsso(3) + rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*cos(zlat(j &
314 476160 ))
315
316 tbls(1) = tbls(1) - rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*sin(zlat(j &
317 ))*cos(zlon(i)) + rea**3*cos(zlat(j))*dlon*dlat*blsv(i, j)*sin(zlon(i &
318 476160 ))
319
320 tbls(2) = tbls(2) - rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*sin(zlat(j &
321 ))*sin(zlon(i)) - rea**3*cos(zlat(j))*dlon*dlat*blsv(i, j)*cos(zlon(i &
322 ))
323
324 tbls(3) = tbls(3) + rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*cos(zlat(j &
325 491040 ))
326
327 END DO
328 END DO
329
330
331 ! write(*,*) 'AAM',rsec,
332 ! write(*,*) 'AAM',rjour+rsec/86400.,
333 ! c raam(3)/hadday,oaam(3)/hadday,
334 ! c tmou(3)/hadley,tsso(3)/hadley,tbls(3)/hadley
335
336 ! write(iam,100)rjour+rsec/86400.,
337 ! c raam(1)/hadday,oaam(1)/hadday,
338 ! c tmou(1)/hadley,tsso(1)/hadley,tbls(1)/hadley,
339 ! c raam(2)/hadday,oaam(2)/hadday,
340 ! c tmou(2)/hadley,tsso(2)/hadley,tbls(2)/hadley,
341 ! c raam(3)/hadday,oaam(3)/hadday,
342 ! c tmou(3)/hadley,tsso(3)/hadley,tbls(3)/hadley
343 100 FORMAT (F12.5, 15(1X,F12.5))
344
345 ! write(iam+1,*)((zs(i,j),i=1,nbp_lon),j=1,nbp_lat)
346 ! write(iam+1,*)((ps(i,j),i=1,nbp_lon),j=1,nbp_lat)
347 ! write(iam+1,*)((ub(i,j),i=1,nbp_lon),j=1,nbp_lat)
348 ! write(iam+1,*)((vb(i,j),i=1,nbp_lon),j=1,nbp_lat)
349 ! write(iam+1,*)((ssou(i,j),i=1,nbp_lon),j=1,nbp_lat)
350 ! write(iam+1,*)((ssov(i,j),i=1,nbp_lon),j=1,nbp_lat)
351 ! write(iam+1,*)((blsu(i,j),i=1,nbp_lon),j=1,nbp_lat)
352 ! write(iam+1,*)((blsv(i,j),i=1,nbp_lon),j=1,nbp_lat)
353
354 480 aam = raam(3)
355 480 torsfc = tmou(3) + tsso(3) + tbls(3)
356
357 480 RETURN
358 END SUBROUTINE aaam_bud
359