Line data Source code
1 : module edyn_maggrid
2 : use shr_kind_mod, only : r8 => shr_kind_r8 ! 8-byte reals
3 : use cam_logfile, only: iulog
4 : use edyn_params, only: finit
5 :
6 : implicit none
7 :
8 : !
9 : ! Global geomagnetic grid:
10 : !
11 : integer, protected :: &
12 : nmlat, & ! number of mag latitudes
13 : nmlath, & ! index of magnetic equator
14 : nmlon, & ! number of mag longitudes
15 : nmlonp1 ! number of longitudes plus periodic point
16 :
17 : !
18 : ! geomagnetic grid resolution parameters:
19 : !
20 : integer, protected :: res_nlev
21 : integer, protected :: res_ngrid
22 :
23 : !
24 : ! Mag grid coordinates:
25 : !
26 : real(r8), allocatable, protected :: &
27 : ylatm(:), & ! magnetic latitudes (radians)
28 : ylonm(:), & ! magnetic longitudes (radians)
29 : gmlat(:), & ! magnetic latitudes (degrees)
30 : gmlon(:) ! magnetic longitudes (degrees)
31 : real(r8), protected :: dlonm,dlatm
32 : !
33 : ! Level coordinates will be same as geographic levels:
34 : !
35 : integer, protected :: nmlev ! number of levels (same as nlev in geographic)
36 :
37 : real(r8), allocatable, protected :: &
38 : rcos0s(:), & ! cos(theta0)/cos(thetas)
39 : dt0dts(:), & ! d(theta0)/d(thetas)
40 : dt1dts(:) ! dt0dts/abs(sinim) (non-zero at equator)
41 :
42 :
43 : real(r8), protected :: table(91,2) = finit
44 :
45 : logical, private :: debug = .false. ! set true for prints to stdout at each call
46 :
47 : contains
48 :
49 : !-----------------------------------------------------------------------
50 768 : subroutine alloc_maggrid( mag_nlon, mag_nlat, mag_nlev, mag_ngrid )
51 :
52 : integer, intent(in) :: mag_nlon, mag_nlat, mag_nlev, mag_ngrid
53 :
54 768 : res_nlev = mag_nlev
55 768 : res_ngrid = mag_ngrid
56 :
57 768 : nmlat = mag_nlat ! number of mag latitudes
58 768 : nmlath = (nmlat+1)/2 ! index of magnetic equator
59 768 : nmlon = mag_nlon ! number of mag longitudes
60 768 : nmlonp1 = nmlon+1 ! number of longitudes plus periodic point
61 :
62 2304 : allocate(ylatm(nmlat))
63 2304 : allocate(ylonm(nmlonp1))
64 1536 : allocate(gmlat(nmlat))
65 1536 : allocate(gmlon(nmlonp1))
66 1536 : allocate(rcos0s(nmlat))
67 1536 : allocate(dt0dts(nmlat))
68 1536 : allocate(dt1dts(nmlat))
69 :
70 768 : end subroutine alloc_maggrid
71 :
72 : !-----------------------------------------------------------------------
73 768 : subroutine set_maggrid()
74 : use edyn_params, only: pi, pi_dyn, rtd, r0
75 : use edyn_mpi, only: nlev => nlev_geo
76 : !
77 : ! Local:
78 : integer :: i, j, n
79 : real(r8) :: tanths2, dtheta, real8
80 1536 : real(r8) :: tanth0(nmlat)
81 1536 : real(r8) :: tanths(nmlat)
82 1536 : real(r8) :: theta0(nmlat)
83 1536 : real(r8) :: hamh0(nmlat)
84 :
85 : real(r8), parameter :: e = 1.e-6_r8
86 : real(r8), parameter :: r1 = 1.06e7_r8
87 : real(r8), parameter :: alfa = 1.668_r8
88 :
89 : real(r8) :: table2(91, 3:5)
90 :
91 768 : real8 = real(nmlat-1, r8)
92 768 : dlatm = pi_dyn / real8
93 768 : real8 = real(nmlon, r8)
94 768 : dlonm = 2._r8 * pi_dyn / real8
95 : !
96 : ! ylatm is equally spaced in theta0, but holds the corresponding
97 : ! value of thetas.
98 : !
99 75264 : do j = 1, nmlat
100 74496 : real8 = real(j-1, r8)
101 75264 : theta0(j) = -pi_dyn/2._r8+real8*dlatm ! note use of pi_dyn
102 : end do ! j=1,nmlat
103 73728 : do j=2,nmlat-1
104 72960 : tanth0(j) = abs(tan(theta0(j)))
105 : hamh0(j) = r1*tanth0(j)+r0*tanth0(j)**(2._r8+2._r8*alfa)/ &
106 72960 : (1._r8+tanth0(j)**2)**alfa
107 72960 : tanths(j) = sqrt(hamh0(j)/r0)
108 72960 : ylatm(j) = sign(atan(tanths(j)),theta0(j))
109 72960 : rcos0s(j) = sqrt((1._r8+tanths(j)**2)/(1._r8+tanth0(j)**2))
110 : !
111 : ! Timegcm has an alternate calculation for dt1dts and dt0dts if dynamo
112 : ! is not called.
113 : !
114 72960 : tanths2 = tanths(j)**2
115 0 : dt1dts(j) = &
116 : (r0*sqrt(1._r8+4._r8*tanths2)*(1._r8+tanths2))/ &
117 : (r1*(1._r8+tanth0(j)**2)+2._r8*r0*tanth0(j)**(2._r8*alfa+1._r8)* &
118 72960 : (1._r8+alfa+tanth0(j)**2)/(1._r8+tanth0(j)**2)**alfa)
119 73728 : dt0dts(j) = dt1dts(j)*2._r8*tanths(j)/sqrt(1._r8+4._r8*tanths2)
120 : end do ! j=2,nmlat-1
121 : !
122 : ! Magnetic poles:
123 : !
124 768 : ylatm(1) = theta0(1)
125 768 : ylatm(nmlat) = theta0(nmlat)
126 768 : rcos0s(1) = 1._r8
127 768 : rcos0s(nmlat) = 1._r8
128 768 : dt0dts(1) = 1._r8
129 768 : dt0dts(nmlat) = 1._r8
130 : !
131 : ! Magnetic longitudes:
132 : !
133 62976 : do i=1,nmlonp1
134 62208 : real8 = real(i-1, r8)
135 62976 : ylonm(i) = -pi+real8*dlonm
136 : ! ylonm(i) = real8*dlonm
137 : end do ! i=1,nmlonp1
138 : !
139 : ! Define mag grid in degrees, and mag levels:
140 : !
141 75264 : gmlat(:) = ylatm(:)*rtd
142 62976 : gmlon(:) = ylonm(:)*rtd
143 : !
144 : ! Magnetic levels are same as midpoint geographic levels:
145 : !
146 768 : nmlev = nlev
147 :
148 : !
149 : ! Calculate table:
150 : !
151 768 : table(1,1) = 0._r8
152 768 : table(1,2) = 0._r8
153 768 : dtheta = pi / 180._r8
154 69888 : do i = 2, 91
155 69888 : table(i,1) = table(i-1,1)+dtheta
156 : end do
157 69120 : do i=2,90
158 68352 : table2(i,4) = tan(table(i,1))
159 69120 : table(i,2) = table(i,1)
160 : end do ! i=2,90
161 768 : table(91,2) = table(91,1)
162 6144 : do n=1,7
163 484608 : do i=2,90
164 478464 : table2(i,3) = table(i,2)
165 478464 : table(i,2) = tan(table2(i,3))
166 : table2(i,5) = sqrt(r1/r0*table(i,2)+table(i,2)**(2._r8*(1._r8+alfa))/ &
167 478464 : (1._r8+table(i,2)**2)**alfa)
168 : table(i,2) = table2(i,3)-(table2(i,5)-table2(i,4))*2._r8* &
169 : table2(i,5)/(r1/r0*(1._r8+table(i,2)**2)+2._r8*table(i,2)** &
170 : (2._r8*alfa+1._r8)*(1._r8+alfa+table(i,2)**2)/ &
171 483840 : (1._r8+table(i,2)**2)**alfa)
172 : end do ! i=2,90
173 : end do ! n=1,7
174 :
175 768 : if (debug) then
176 0 : write(iulog,"('set_maggrid: table= ',/,(6e12.4))") table
177 0 : write(iulog,"('set_maggrid: table2=',/,(6e12.4))") table2
178 : end if
179 :
180 768 : end subroutine set_maggrid
181 : !-----------------------------------------------------------------------
182 : end module edyn_maggrid
|
