--- trunk/libf/phylmd/aaam_bud.f90 2012/01/10 19:02:02 56 +++ trunk/phylmd/aaam_bud.f90 2019/06/13 14:40:06 328 @@ -4,24 +4,21 @@ contains - subroutine aaam_bud(rea, rg, ome, plat, plon, phis, dragu, liftu, phyu, & - dragv, liftv, phyv, p, u, v, aam, torsfc) + subroutine aaam_bud(rg, ome, phis, dragu, liftu, phyu, dragv, liftv, phyv, & + p, u, v, aam, torsfc) ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute - ! different terms of the axial AAAM Budget and mountain torque. + ! different terms of the axial AAAM budget and mountain torque. ! Only valid for regular rectangular grids. Should be called after ! "lift_noro". - USE dimens_m, ONLY : iim, jjm + USE dimensions, ONLY : iim, jjm use nr_util, only: assert_eq, assert, pi + use phyetat0_m, only: rlat, rlon + USE suphec_m, ONLY: ra - real, intent(in):: rea ! Earth radius real, intent(in):: rg ! gravity constant real, intent(in):: ome ! Earth rotation rate - - REAL, intent(in):: plat(:), plon(:) - ! (nlon) latitude and longitude in degrees - real, intent(in):: phis(:) ! (nlon) Geopotential at the ground REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal) REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal) @@ -41,7 +38,6 @@ INTEGER nlev ! number of vertical levels INTEGER i, j, k, l - REAL hadley, hadday REAL dlat, dlon ! latitude and longitude increments (radians) REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial) @@ -49,7 +45,6 @@ REAL tmou(3) ! resolved mountain torque (3 components) REAL tsso(3) ! parameterised moutain drag torque (3 components) REAL tbls(3) ! parameterised boundary layer torque (3 components) - integer iax REAL ZS(801, 401) ! topographic height REAL PS(801, 401) ! surface pressure @@ -60,9 +55,9 @@ !------------------------------------------------------------------- - call assert(size(plat) == (/size(plon), size(phis), size(dragu), & - size(liftu), size(phyu), size(dragv), size(liftv), size(phyv), & - size(p, 1), size(u, 1), size(v, 1)/), "aaam_bud nlon") + call assert(size(phis) == (/size(dragu), size(liftu), size(phyu), & + size(dragv), size(liftv), size(phyv), size(p, 1), size(u, 1), & + size(v, 1)/), "aaam_bud nlon") nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev") if (iim + 1 > 801 .or. jjm + 1 > 401) then @@ -70,8 +65,6 @@ stop 1 endif - hadley = 1e18 - hadday = 1e18 * 24. * 3600. dlat = pi / jjm dlon = 2 * pi / real(iim) @@ -92,7 +85,7 @@ vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg enddo - zlat(1) = plat(1) * pi / 180. + zlat(1) = rlat(1) * pi / 180. do i = 1, iim + 1 zs(i, 1) = phis(1) / rg @@ -115,8 +108,8 @@ ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1) blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1) blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1) - zlon(iim + 1) = - plon(l + 1) * pi / 180. - zlat(j) = plat(l + 1) * pi / 180. + zlon(iim + 1) = - rlon(l + 1) * pi / 180. + zlat(j) = rlat(l + 1) * pi / 180. ub(iim + 1, j) = 0. vb(iim + 1, j) = 0. @@ -135,7 +128,7 @@ ssov(i, j) = dragv(l) + liftv(l) blsu(i, j) = phyu(l) - dragu(l) - liftu(l) blsv(i, j) = phyv(l) - dragv(l) - liftv(l) - zlon(i) = plon(l) * pi / 180. + zlon(i) = rlon(l) * pi / 180. ub(i, j) = 0. vb(i, j) = 0. @@ -155,7 +148,7 @@ ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg enddo - zlat(jjm + 1) = plat(l) * pi / 180. + zlat(jjm + 1) = rlat(l) * pi / 180. do i = 1, iim + 1 zs(i, jjm + 1) = phis(l) / rg @@ -172,32 +165,32 @@ DO j = 1, jjm DO i = 1, iim - raam(1) = raam(1) - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) & + raam(1) = raam(1) - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) & * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) & * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) & - + rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) & + + ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) & * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1)) - oaam(1) = oaam(1) - ome * rea**4 * dlon * dlat / rg * 0.5 & + oaam(1) = oaam(1) - ome * ra**4 * dlon * dlat / rg * 0.5 & * (cos(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) & * ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 & * sin(zlat(j + 1)) * ps(i , j + 1)) - raam(2) = raam(2) - rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) & + raam(2) = raam(2) - ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) & * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) & * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) & - - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) & + - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) & * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1)) - oaam(2) = oaam(2) - ome * rea**4 * dlon * dlat / rg * 0.5 & + oaam(2) = oaam(2) - ome * ra**4 * dlon * dlat / rg * 0.5 & * (sin(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) & * ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 & * sin(zlat(j + 1)) * ps(i , j + 1)) - raam(3) = raam(3) + rea**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 & + raam(3) = raam(3) + ra**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 & * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1)) - oaam(3) = oaam(3) + ome * rea**4 * dlon * dlat / rg * 0.5 & + oaam(3) = oaam(3) + ome * ra**4 * dlon * dlat / rg * 0.5 & * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 & * ps(i, j + 1)) ENDDO @@ -207,10 +200,10 @@ DO j = 1, jjm DO i = 1, iim - tmou(1) = tmou(1) - rea**2 * dlon * 0.5 * sin(zlon(i)) & + tmou(1) = tmou(1) - ra**2 * dlon * 0.5 * sin(zlon(i)) & * (zs(i, j) - zs(i, j + 1)) & * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j)) - tmou(2) = tmou(2) + rea**2 * dlon * 0.5 * cos(zlon(i)) & + tmou(2) = tmou(2) + ra**2 * dlon * 0.5 * cos(zlon(i)) & * (zs(i, j) - zs(i, j + 1)) & * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j)) ENDDO @@ -218,13 +211,13 @@ DO j = 2, jjm DO i = 1, iim - tmou(1) = tmou(1) + rea**2 * dlat * 0.5 * sin(zlat(j)) & + tmou(1) = tmou(1) + ra**2 * dlat * 0.5 * sin(zlat(j)) & * (zs(i + 1, j) - zs(i, j)) & * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j)) - tmou(2) = tmou(2) + rea**2 * dlat * 0.5 * sin(zlat(j)) & + tmou(2) = tmou(2) + ra**2 * dlat * 0.5 * sin(zlat(j)) & * (zs(i + 1, j) - zs(i, j)) & * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j)) - tmou(3) = tmou(3) - rea**2 * dlat * 0.5* cos(zlat(j)) & + tmou(3) = tmou(3) - ra**2 * dlat * 0.5* cos(zlat(j)) & * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j)) ENDDO ENDDO @@ -233,30 +226,30 @@ DO j = 2, jjm DO i = 1, iim - tsso(1) = tsso(1) - rea**3 * cos(zlat(j)) * dlon * dlat* & + tsso(1) = tsso(1) - ra**3 * cos(zlat(j)) * dlon * dlat* & ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) & - + rea**3 * cos(zlat(j)) * dlon * dlat* & + + ra**3 * cos(zlat(j)) * dlon * dlat* & ssov(i, j) * sin(zlon(i)) - tsso(2) = tsso(2) - rea**3 * cos(zlat(j)) * dlon * dlat* & + tsso(2) = tsso(2) - ra**3 * cos(zlat(j)) * dlon * dlat* & ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) & - - rea**3 * cos(zlat(j)) * dlon * dlat* & + - ra**3 * cos(zlat(j)) * dlon * dlat* & ssov(i, j) * cos(zlon(i)) - tsso(3) = tsso(3) + rea**3 * cos(zlat(j)) * dlon * dlat* & + tsso(3) = tsso(3) + ra**3 * cos(zlat(j)) * dlon * dlat* & ssou(i, j) * cos(zlat(j)) - tbls(1) = tbls(1) - rea**3 * cos(zlat(j)) * dlon * dlat* & + tbls(1) = tbls(1) - ra**3 * cos(zlat(j)) * dlon * dlat* & blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) & - + rea**3 * cos(zlat(j)) * dlon * dlat* & + + ra**3 * cos(zlat(j)) * dlon * dlat* & blsv(i, j) * sin(zlon(i)) - tbls(2) = tbls(2) - rea**3 * cos(zlat(j)) * dlon * dlat* & + tbls(2) = tbls(2) - ra**3 * cos(zlat(j)) * dlon * dlat* & blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) & - - rea**3 * cos(zlat(j)) * dlon * dlat* & + - ra**3 * cos(zlat(j)) * dlon * dlat* & blsv(i, j) * cos(zlon(i)) - tbls(3) = tbls(3) + rea**3 * cos(zlat(j)) * dlon * dlat* & + tbls(3) = tbls(3) + ra**3 * cos(zlat(j)) * dlon * dlat* & blsu(i, j) * cos(zlat(j)) ENDDO ENDDO