Line data Source code
1 : module getapex
2 : !
3 : ! Calculate quantities needed to transform scalar fields between geographic
4 : ! and geomagnetic coordinate systems.
5 : !
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use cam_logfile, only: iulog
8 : use cam_abortutils, only: endrun
9 : use edyn_geogrid, only: nlon,nlonp1,jspole,jnpole
10 : use edyn_maggrid, only: nmlonp1,nmlat,ylatm,ylonm
11 : use infnan, only: nan, assignment(=)
12 :
13 : implicit none
14 : save
15 :
16 : private
17 :
18 : public :: get_apex ! Allocate and initialize apex data
19 : public :: magfield, bx, by, bz
20 : public :: bmod2, bmod
21 : public :: xb, yb, zb
22 : public :: be3arr, dddarr, dvec
23 : public :: alatm, alonm
24 : public :: gdlatdeg, gdlondeg
25 : public :: rjac
26 :
27 : real(r8),dimension(:,:), allocatable :: & ! geo lat,lon coords on mag grid
28 : gdlatdeg, & ! geographic latitude of each magnetic grid point (deg)
29 : gdlondeg ! geographic longitude of each magnetic grid point (deg)
30 : !
31 : ! Variables on geographic grid needed by other modules must
32 : ! be allocated dynamically to be grid-independent (sub alloc_apex):
33 : !
34 : real(r8),allocatable :: & ! (nlonp1,jspole:jnpole,3,2)
35 : dvec(:,:,:,:) ! vectors from apxmall
36 :
37 : real(r8),allocatable :: & ! (nlonp1,jspole:jnpole)
38 : dddarr(:,:), & ! from apxmall
39 : be3arr(:,:) ! from apxmall
40 :
41 : real(r8),allocatable :: & ! (nlonp1,jspole:jnpole)
42 : alatm(:,:), & ! geomagnetic latitude at each geographic grid point (radians)
43 : alonm(:,:), & ! geomagnetic longitude at each geographic grid point (radians)
44 : xb(:,:), & ! northward component of magnetic field
45 : yb(:,:), & ! eastward component of magnetic field
46 : zb(:,:), & ! downward component of magnetic field (gauss)
47 : bmod(:,:) ! magnitude of magnetic field (gauss)
48 : !
49 : ! rjac: scaled derivatives of geomagnetic coords wrt geographic coordinates.
50 : ! rjac(1,1) = cos(thetas)/cos(theta)*d(lamdas)/d(lamda)
51 : ! rjac(1,2) = cos(thetas)*d(lamdas)/d(theta)
52 : ! rjac(2,1) = 1./cos(theta)*d(thetas)/d(lamda)
53 : ! rjac(2,2) = d(thetas)/d(theta)
54 : ! where (lamda,theta) are geographic coordinates
55 : ! (lamdas,thetas) are geomagnetic coordinates
56 : !
57 : real(r8),allocatable :: &
58 : rjac(:,:,:,:) ! (nlon+1,jspole:jnpole,2,2)
59 : !
60 : ! Parameters defined by sub magfield (allocated in alloc_magfield):
61 : !
62 : real(r8),allocatable,dimension(:,:) :: & ! (0:nlon+1,jspole-1:jnpole+1)
63 : bx,by,bz,bmod2
64 :
65 : contains
66 : !-----------------------------------------------------------------------
67 768 : subroutine get_apex( )
68 : !
69 : ! This is called once per run from main.
70 : !
71 : use edyn_params, only: re_dyn, h0, hs, dtr, rtd, cm2km
72 : use apex, only: apex_mall, apex_q2g
73 : use edyn_geogrid, only: glat_edyn_geo=>glat, glon_edyn_geo=>glon
74 :
75 : !
76 : ! Local:
77 : integer :: i, j, ier
78 : real(r8) :: rekm, h0km, alt, hr, ror03, glat, glon
79 : real(r8) :: qdlon, qdlat, gdlon, gdlat
80 : integer, parameter :: nalt=2
81 :
82 : !
83 : ! Non-scalar arguments returned by apxmall:
84 : real(r8) :: &
85 : b(3), bhat(3), &
86 : d1(3), d2(3),d3(3), &
87 : e1(3), e2(3),e3(3), &
88 : f1(2), f2(2)
89 : real(r8) :: bmag, alon, xlatm, vmp, w, d, be3, si, sim, xlatqd, f
90 :
91 : !
92 : ! Allocate arrays that are needed by other modules:
93 768 : call alloc_apex
94 768 : call alloc_magfield
95 :
96 768 : dddarr = nan
97 768 : dvec = nan
98 768 : rekm = re_dyn*cm2km ! earth radius (km)
99 768 : h0km = h0*cm2km
100 768 : alt = hs*cm2km ! modified apex reference altitude (km)
101 768 : hr = alt
102 768 : ror03= ((rekm + alt)/(rekm + h0km))**3
103 : !
104 : ! Loop over 2d geographic grid:
105 : !
106 74496 : do j = jspole, jnpole
107 73728 : glat = glat_edyn_geo(j)
108 10765056 : do i = 1, nlonp1
109 10690560 : if (i == nlonp1) then
110 73728 : glon = glon_edyn_geo(1)
111 : else
112 10616832 : glon = glon_edyn_geo(i)
113 : end if
114 :
115 : call apex_mall ( &
116 : glat,glon,alt,hr, & !Inputs
117 : b,bhat,bmag,si, & !Mag Fld
118 : alon, & !Apx Lon
119 : xlatm,vmp,w,d,be3,sim,d1,d2,d3,e1,e2,e3, & !Mod Apx
120 10690560 : xlatqd,f,f1,f2, ier) !Qsi-Dpl
121 :
122 10690560 : if (ier /= 0) then
123 0 : call endrun('get_apex: apxmall error')
124 : end if
125 :
126 10690560 : alatm(i,j) = xlatm*dtr
127 10690560 : alonm(i,j) = alon *dtr
128 10690560 : xb (i,j) = b(2)*1.e-5_r8 ! nT -> gauss
129 10690560 : yb (i,j) = b(1)*1.e-5_r8 ! nT -> gauss
130 10690560 : zb (i,j) = -b(3)*1.e-5_r8 ! nT -> gauss
131 10690560 : bmod (i,j) = bmag*1.e-5_r8 ! nT -> gauss
132 :
133 10690560 : rjac (i,j,1,1) = f2(2)
134 10690560 : rjac (i,j,1,2) = -f2(1)
135 10690560 : rjac (i,j,2,1) = -f1(2)
136 10690560 : rjac (i,j,2,2) = f1(1)
137 : !
138 : ! Set up parameters for magnetic to geographic interpolation.
139 : !
140 10690560 : dvec(i,j,1,1) = d1(1)
141 10690560 : dvec(i,j,2,1) = d1(2)
142 10690560 : dvec(i,j,3,1) = d1(3)
143 10690560 : dvec(i,j,1,2) = d2(1)
144 10690560 : dvec(i,j,2,2) = d2(2)
145 10690560 : dvec(i,j,3,2) = d2(3)
146 10690560 : dddarr(i,j) = d
147 : !
148 : ! Scale be3 from 130 km to a reference height of 90 km.
149 21454848 : be3arr(i,j) = be3 * ror03
150 : end do ! i=1,nlonp1
151 : end do ! j=jspole,jnpole
152 : !
153 : ! Set up parameters for geographic to magnetic interpolation
154 62976 : do i = 1, nmlonp1
155 62208 : qdlon = ylonm(i)*rtd
156 6097152 : do j = 1, nmlat
157 6034176 : qdlat = ylatm(j)*rtd
158 : !
159 : ! Convert from Quasi-Dipole to geographic coordinates.
160 : ! gdlat,gdlon are returned by apxq2g.
161 : !
162 6034176 : call apex_q2g(qdlat, qdlon, alt, gdlat, gdlon, ier)
163 6034176 : if (ier /= 0) then
164 0 : write(iulog,"(i3,i3,i3)") '>>> Error from apex_q2g: ier=',ier, &
165 0 : ' i=',i,' j=',j
166 0 : call endrun('get_apex: apex_q2g ier')
167 : end if
168 6034176 : gdlat = gdlat * dtr
169 6034176 : gdlon = gdlon * dtr
170 : !
171 : ! gdlatdeg,gdlondeg will be coordY,coordX of the mag grid for ESMF
172 : ! regridding (see edyn_esmf.F)
173 : !
174 6034176 : gdlatdeg(i,j) = gdlat*rtd
175 12130560 : gdlondeg(i,j) = gdlon*rtd
176 : enddo ! j=1,nmlat
177 : enddo ! i=1,nmlonp1
178 768 : end subroutine get_apex
179 : !-----------------------------------------------------------------------
180 768 : subroutine magfield
181 : !
182 : ! Calculate magnetic field parameters (bx,by,bz)
183 : ! (see also TIEGCM magfield.F)
184 : ! This is called once per run and when crossing year boundary from edyn_init, after get_apex.
185 : ! All arrays are on the global domain, all processors execute.
186 : !
187 : ! Local:
188 : integer :: i,j
189 : !
190 : ! QUESTION: in TIEGCM, dipmin is resolution dependent - how do we
191 : ! handle this for different resolutions in WACCM?
192 : !
193 : ! real(r8),parameter :: dipmin=0.17 ! set for 5.0-deg TIEGCM (also known as sin10)
194 : real(r8),parameter :: dipmin=0.24_r8 ! set for 2.5-deg TIEGCM (also known as sin10)
195 : real(r8) :: cos10
196 :
197 768 : cos10 = sqrt(1._r8-dipmin**2)
198 74496 : do j=jspole,jnpole ! 1,nlat
199 10691328 : do i=1,nlon
200 10616832 : bx(i,j) = yb(i,j)/bmod(i,j)
201 10616832 : by(i,j) = xb(i,j)/bmod(i,j)
202 10616832 : bz(i,j) = -zb(i,j)/bmod(i,j)
203 10616832 : bmod2(i,j) = bmod(i,j)
204 : !
205 : ! Set minimum dip to 10 degrees
206 10690560 : if (abs(bz(i,j))-dipmin < 0._r8) then
207 771840 : bx(i,j) = bx(i,j)*(cos10/sqrt(1._r8-bz(i,j)**2))
208 771840 : by(i,j) = by(i,j)*(cos10/sqrt(1._r8-bz(i,j)**2))
209 771840 : bz(i,j) = sign(dipmin,bz(i,j))
210 : endif
211 : enddo ! i=1,nlon
212 : enddo ! j=jspole,jnpole
213 :
214 : !
215 : ! Values at jspole-1:
216 768 : j=jspole-1 ! j=0
217 111360 : do i=1,nlon
218 110592 : bx(i,j) = -bx(1+mod(i-1+nlon/2,nlon),jspole)
219 110592 : by(i,j) = -by(1+mod(i-1+nlon/2,nlon),jspole)
220 110592 : bz(i,j) = bz(1+mod(i-1+nlon/2,nlon),jspole)
221 111360 : bmod2(i,j) = bmod2(1+mod(i-1+nlon/2,nlon),jspole)
222 : enddo
223 : !
224 : ! Values at jnpole+1:
225 768 : j=jnpole+1 ! j=nlat+1
226 111360 : do i=1,nlon
227 110592 : bx(i,j) = -bx(1+mod(i-1+nlon/2,nlon),jnpole)
228 110592 : by(i,j) = -by(1+mod(i-1+nlon/2,nlon),jnpole)
229 110592 : bz(i,j) = bz(1+mod(i-1+nlon/2,nlon),jnpole)
230 111360 : bmod2(i,j) = bmod2(1+mod(i-1+nlon/2,nlon),jnpole)
231 : enddo
232 : !
233 : ! Periodic points:
234 : ! FIX: not sure about this, but
235 : ! I am following tiegcm, but with a single point on each end instead of 2
236 : !
237 76032 : do j=jspole-1,jnpole+1
238 75264 : bx(nlonp1,j) = bx(1,j)
239 75264 : by(nlonp1,j) = by(1,j)
240 75264 : bz(nlonp1,j) = bz(1,j)
241 75264 : bmod2(nlonp1,j) = bmod2(1,j)
242 :
243 75264 : bx(0,j) = bx(nlon,j)
244 75264 : by(0,j) = by(nlon,j)
245 75264 : bz(0,j) = bz(nlon,j)
246 76032 : bmod2(0,j) = bmod2(nlon,j)
247 : enddo
248 :
249 768 : end subroutine magfield
250 : !-----------------------------------------------------------------------
251 768 : subroutine alloc_magfield
252 :
253 : !------------------------------------------------------------------------------------------
254 : ! Do allocations, checking if previously allocated in case of year boundary crossing
255 : !------------------------------------------------------------------------------------------
256 3072 : if (.not.allocated(bx)) allocate(bx(0:nlonp1,jspole-1:jnpole+1))
257 3072 : if (.not.allocated(by)) allocate(by(0:nlonp1,jspole-1:jnpole+1))
258 3072 : if (.not.allocated(bz)) allocate(bz(0:nlonp1,jspole-1:jnpole+1))
259 3072 : if (.not.allocated(bmod2)) allocate(bmod2(0:nlonp1,jspole-1:jnpole+1))
260 :
261 768 : end subroutine alloc_magfield
262 : !-----------------------------------------------------------------------
263 :
264 768 : subroutine alloc_apex
265 :
266 : !------------------------------------------------------------------------------------------
267 : ! Do allocations, checking if previously allocated in case of year boundary crossing
268 : !------------------------------------------------------------------------------------------
269 :
270 3072 : if (.not.allocated(xb)) allocate(xb (nlonp1,jspole:jnpole))
271 3072 : if (.not.allocated(yb)) allocate(yb (nlonp1,jspole:jnpole))
272 3072 : if (.not.allocated(zb)) allocate(zb (nlonp1,jspole:jnpole))
273 3072 : if (.not.allocated(bmod)) allocate(bmod (nlonp1,jspole:jnpole))
274 3072 : if (.not.allocated(alatm)) allocate(alatm(nlonp1,jspole:jnpole))
275 3072 : if (.not.allocated(alonm))allocate(alonm(nlonp1,jspole:jnpole))
276 :
277 4608 : if (.not.allocated(dvec)) allocate(dvec (nlonp1,jspole:jnpole,3,2))
278 3072 : if (.not.allocated(dddarr)) allocate(dddarr(nlonp1,jspole:jnpole))
279 3072 : if (.not.allocated(be3arr)) allocate(be3arr(nlonp1,jspole:jnpole))
280 :
281 4608 : if (.not.allocated(rjac)) allocate(rjac(nlon+1,jspole:jnpole,2,2))
282 :
283 3072 : if (.not.allocated(gdlatdeg)) allocate(gdlatdeg(nmlonp1,nmlat))
284 3072 : if (.not.allocated(gdlondeg)) allocate(gdlondeg(nmlonp1,nmlat))
285 :
286 768 : end subroutine alloc_apex
287 : !-----------------------------------------------------------------------
288 : end module getapex
|
