Line data Source code
1 : module adotv_mod
2 : use shr_kind_mod,only: r8 => shr_kind_r8 ! 8-byte reals
3 :
4 : implicit none
5 :
6 : contains
7 :
8 7296 : subroutine calc_adotv(z, un, vn, wn, adotv1, adotv2, adota1, adota2, &
9 7296 : a1dta2, be3, sini, lev0, lev1, lon0, lon1, lat0, lat1)
10 : !
11 : ! Calculate adotv1,2, adota1,2, a1dta2 and be3.
12 : ! All fields should be on O+ grid
13 : !
14 : use edyn_params, only: r0,h0
15 : use edyn_geogrid, only: jspole, jnpole
16 : use getapex, only: &
17 : zb, & ! downward component of magnetic field
18 : bmod, & ! magnitude of magnetic field (gauss)
19 : dvec, & ! (nlonp1,nlat,3,2)
20 : dddarr, & ! (nlonp1,nlat)
21 : be3arr, & ! (nlonp1,nlat)
22 : alatm ! (nlonp1,0:nlatp1)
23 : !
24 : ! Args:
25 : integer,intent(in) :: lev0, lev1, lon0, lon1, lat0, lat1
26 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1), intent(in) :: &
27 : z, & ! geopotential height (cm)
28 : un, & ! neutral zonal velocity (cm/s)
29 : vn ! neutral meridional velocity (cm/s)
30 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1), intent(in) :: &
31 : wn ! vertical velocity (cm/s)
32 :
33 : real(r8), dimension(lon0:lon1,lat0:lat1,lev0:lev1), intent(out) :: &
34 : adotv1, adotv2
35 : real(r8), dimension(lon0:lon1,lat0:lat1), intent(out) :: &
36 : adota1, adota2, a1dta2, be3, sini
37 : !
38 : ! Local:
39 : integer :: k, i, j
40 : real(r8) :: r0or, rat, sinalat
41 7296 : real(r8) :: clm2(lon0:lon1,lat0:lat1)
42 : !
43 37946496 : adotv1 = 0.0_r8
44 37946496 : adotv2 = 0.0_r8
45 291840 : adota1 = 0.0_r8
46 291840 : adota2 = 0.0_r8
47 291840 : a1dta2 = 0.0_r8
48 291840 : be3 = 0.0_r8
49 291840 : sini = 0.0_r8
50 :
51 29184 : do j = lat0, lat1
52 291840 : do i = lon0, lon1
53 262656 : sinalat = sin(alatm(i,j)) ! sin(lam)
54 262656 : clm2(i,j) = 1._r8 - (sinalat * sinalat) ! cos^2(lam)
55 262656 : be3(i,j) = 1.e-9_r8*be3arr(i,j) ! be3 is in T (be3arr in nT)
56 262656 : sini(i,j) = zb(i,j)/bmod(i,j) ! sin(I_m)
57 :
58 34145280 : do k=lev0,lev1-1
59 : !
60 : ! d_1 = (R_0/R)^1.5
61 33882624 : r0or = r0/(r0 + 0.5_r8* (z(k,i,j) + z(k+1,i,j)) - h0)
62 33882624 : rat = 1.e-2_r8*r0or**1.5_r8 ! 1/100 conversion in cm
63 : !
64 : ! A_1 dot V = fac( d_1(1) u + d_1(2) v + d_1(3) w
65 : adotv1(i,j,k) = rat*( &
66 0 : dvec(i,j,1,1) * un(k,i,j) + &
67 33882624 : dvec(i,j,2,1) * vn(k,i,j) + &
68 67765248 : dvec(i,j,3,1) * wn(k,i,j))
69 :
70 : !
71 : ! Note: clm2 is being used here to represent the squared cosine
72 : ! of the quasi-dipole latitude, not of the M(90) latitude,
73 : ! since the wind values are aligned vertically,
74 : ! not along the field line.
75 : !
76 : rat = rat * sqrt((4._r8 - (3._r8 * clm2(i,j))) / &
77 33882624 : (4._r8 - (3._r8 * r0or * clm2(i,j))))
78 : !
79 : ! A_2 dot V = fac( d_2(1) u + d_2(2) v + d_2(3) w
80 : adotv2(i,j,k) = rat * ( &
81 0 : dvec(i,j,1,2) * un(k,i,j) + &
82 33882624 : dvec(i,j,2,2) * vn(k,i,j) + &
83 68027904 : dvec(i,j,3,2) * wn(k,i,j))
84 : end do ! k=lev0,lev1-1
85 :
86 : !
87 : ! Calculation of adota(n) = d(n)**2/D
88 : ! a1dta2 = (d(1) dot d(2)) /D
89 : !
90 262656 : adota1(i,j) = (dvec(i,j,1,1)**2 + dvec(i,j,2,1)**2 + &
91 525312 : dvec(i,j,3,1)**2) / dddarr(i,j)
92 262656 : adota2(i,j) = (dvec(i,j,1,2)**2 + dvec(i,j,2,2)**2 + &
93 525312 : dvec(i,j,3,2)**2) / dddarr(i,j)
94 262656 : a1dta2(i,j) = (dvec(i,j,1,1) * dvec(i,j,1,2) + &
95 262656 : dvec(i,j,2,1) * dvec(i,j,2,2) + &
96 547200 : dvec(i,j,3,1) * dvec(i,j,3,2)) / dddarr(i,j)
97 : end do ! i=lon0,lon1
98 :
99 : end do ! j=lat0,lat1
100 :
101 7296 : end subroutine calc_adotv
102 :
103 :
104 : end module adotv_mod
|
