Line data Source code
1 : module mo_apex
2 :
3 : !-------------------------------------------------------------------------------
4 : ! Purpose:
5 : !
6 : ! Calculate apex coordinates and magnetic field magnitudes
7 : ! at global geographic grid for year of current model run.
8 : !
9 : ! Method:
10 : !
11 : ! The magnetic field parameters output by this module are time and height
12 : ! independent. They are chunked for waccm physics, i.e., allocated as
13 : ! (pcols,begchunk:endchunk)
14 : ! Interface sub apexmag is called once per run from sub inti.
15 : ! Sub apexmag may be called for years 1900 through 2005.
16 : ! This module is dependent on routines in apex_subs.F (modified IGRF model).
17 : ! Apex_subs has several authors, but has been modified and maintained
18 : ! in recent years by Roy Barnes (bozo@ucar.edu).
19 : ! Subs apxmka and apxmall are called with the current lat x lon grid
20 : ! resolution.
21 : !
22 : ! Author: Ben Foster, foster@ucar.edu (Nov, 2003)
23 : !-------------------------------------------------------------------------------
24 :
25 : use shr_kind_mod, only: r8 => shr_kind_r8
26 : use ppgrid, only: pcols, begchunk, endchunk ! physics grid
27 : use cam_abortutils, only: endrun
28 : use cam_logfile, only: iulog
29 : use spmd_utils, only: masterproc
30 : use apex, only: apex_mka, apex_mall, apex_dypol, apex_set_igrf
31 : use apex, only: apex_beg_yr, apex_end_yr
32 : implicit none
33 :
34 : private
35 : public :: mo_apex_readnl
36 : public :: mo_apex_init
37 : public :: mo_apex_init1
38 : public :: alatm, alonm, bnorth, beast, bdown, bmag
39 : public :: d1vec, d2vec, colatp, elonp
40 : public :: maglon0 ! geographic longitude at the equator where geomagnetic longitude is zero (radians)
41 :
42 : ! year to initialize apex
43 : real(r8), public, protected :: geomag_year = -1._r8
44 : logical, public, protected :: geomag_year_updated = .true.
45 :
46 : integer :: fixed_geomag_year = -1
47 :
48 : !-------------------------------------------------------------------------------
49 : ! Magnetic field output arrays, chunked for physics:
50 : ! (these are allocated (pcols,begchunk:endchunk) by sub allocate_arrays)
51 : !-------------------------------------------------------------------------------
52 : real(r8), protected, allocatable, dimension(:,:) :: & ! (pcols,begchunk:endchunk)
53 : alatm, & ! apex mag latitude at each geographic grid point (radians)
54 : alonm, & ! apex mag longitude at each geographic grid point (radians)
55 : bnorth, & ! northward component of magnetic field
56 : beast, & ! eastward component of magnetic field
57 : bdown, & ! downward component of magnetic field
58 : bmag ! magnitude of magnetic field
59 : real(r8), protected, allocatable, dimension(:,:,:) :: & ! (3,pcols,begchunk:endchunk)
60 : d1vec, & ! base vectors more-or-less magnetic eastward direction
61 : d2vec ! base vectors more-or-less magnetic downward/equatorward direction
62 : real(r8), protected :: &
63 : colatp, & ! geocentric colatitude of geomagnetic dipole north pole (deg)
64 : elonp ! East longitude of geomagnetic dipole north pole (deg)
65 :
66 : real(r8), protected :: maglon0
67 :
68 : character(len=256) :: igrf_geomag_coefs_file = 'igrf_geomag_coefs_file'
69 :
70 : contains
71 :
72 : !======================================================================
73 : !======================================================================
74 1536 : subroutine mo_apex_readnl(nlfile)
75 :
76 : use namelist_utils, only : find_group_name
77 : use units, only : getunit, freeunit
78 : use spmd_utils, only : mpicom, masterprocid, mpi_integer, mpi_character
79 :
80 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
81 :
82 : ! Local variables
83 : integer :: unitn, ierr
84 : character(len=*), parameter :: subname = 'mo_apex_readnl'
85 :
86 : namelist /geomag_nl/ fixed_geomag_year, igrf_geomag_coefs_file
87 :
88 : ! Read namelist
89 1536 : if (masterproc) then
90 2 : unitn = getunit()
91 2 : open( unitn, file=trim(nlfile), status='old' )
92 2 : call find_group_name(unitn, 'geomag_nl', status=ierr)
93 2 : if (ierr == 0) then
94 0 : read(unitn, geomag_nl, iostat=ierr)
95 0 : if (ierr /= 0) then
96 0 : call endrun(subname // ':: ERROR reading namelist')
97 : end if
98 : end if
99 2 : close(unitn)
100 2 : call freeunit(unitn)
101 : end if
102 :
103 : ! Broadcast namelist variables
104 1536 : call mpi_bcast(fixed_geomag_year, 1, mpi_integer, masterprocid, mpicom, ierr)
105 1536 : call mpi_bcast(igrf_geomag_coefs_file, len(igrf_geomag_coefs_file), mpi_character, masterprocid, mpicom, ierr)
106 :
107 1536 : end subroutine mo_apex_readnl
108 :
109 : !======================================================================
110 : !======================================================================
111 0 : subroutine mo_apex_init1()
112 : use time_manager, only: get_curr_date
113 : use phys_grid, only: get_grid_dims
114 :
115 : integer :: i, j, ist ! indices
116 :
117 : integer :: nglats
118 : integer :: nglons
119 : integer, parameter :: ngalts = 2 ! number of altitudes
120 :
121 0 : real(r8), allocatable :: gridlats(:)
122 0 : real(r8), allocatable :: gridlons(:)
123 : real(r8) :: gridalts(ngalts) ! altitudes passed to apxmka
124 :
125 : integer :: ngcols, hdim1_d, hdim2_d
126 : integer :: yr, mon, day, sec
127 :
128 : ! read the IGRF coefs from file
129 0 : call apex_set_igrf( igrf_geomag_coefs_file )
130 :
131 0 : if (fixed_geomag_year>0) then
132 0 : yr = fixed_geomag_year
133 : else
134 0 : call get_curr_date(yr, mon, day, sec)
135 : end if
136 :
137 0 : if ( yr < apex_beg_yr ) yr = apex_beg_yr
138 0 : if ( yr > apex_end_yr-1 ) yr = apex_end_yr-1
139 :
140 0 : if (.not.(yr > geomag_year)) then
141 0 : geomag_year_updated = .false.
142 : return
143 : else
144 0 : geomag_year_updated = .true.
145 : endif
146 :
147 0 : geomag_year = dble(yr)+0.5_r8
148 :
149 : !-------------------------------------------------------------------------------
150 : ! Center min, max altitudes about 130 km
151 : !-------------------------------------------------------------------------------
152 0 : gridalts(:ngalts) = (/ 90._r8, 170._r8 /)
153 :
154 : !-------------------------------------------------------------------------------
155 : ! Initialize APEX with a regular lat/lon grid ...
156 : ! (Note apex_mka expects longitudes in -180 -> +180)
157 : !-------------------------------------------------------------------------------
158 0 : call get_grid_dims(hdim1_d,hdim2_d)
159 0 : ngcols = hdim1_d*hdim2_d
160 0 : if ( ngcols < 1000 ) then ! 10-degrees
161 0 : nglats = 19
162 0 : nglons = 37
163 0 : elseif ( ngcols < 10000 ) then ! 5-degrees
164 0 : nglats = 37
165 0 : nglons = 73
166 0 : elseif ( ngcols < 20000 ) then ! 2-degree
167 0 : nglats = 91
168 0 : nglons = 181
169 0 : elseif ( ngcols < 100000 ) then ! 1-degree
170 0 : nglats = 181
171 0 : nglons = 361
172 : else ! half-degee
173 0 : nglats = 361
174 0 : nglons = 721
175 : endif
176 :
177 0 : allocate ( gridlats(nglats), gridlons(nglons) )
178 0 : do i = 1,nglons
179 0 : gridlons(i) = -180._r8 + dble(i-1)*360._r8/(nglons-1)
180 : enddo
181 0 : do j = 1,nglats
182 0 : gridlats(j) = -90._r8 + dble(j-1)*180._r8/(nglats-1)
183 : enddo
184 :
185 : call apex_mka( geomag_year, gridlats, gridlons, gridalts, &
186 0 : nglats, nglons, ngalts, ist )
187 :
188 0 : if( ist /= 0 ) then
189 0 : write(iulog,"(/,'>>> mo_apex_init: Error from apxmka: ist=',i5)") ist
190 0 : call endrun("mo_apex_init: Error from apxmka")
191 : end if
192 :
193 0 : deallocate( gridlats, gridlons )
194 :
195 0 : if (masterproc) then
196 0 : if (fixed_geomag_year<1) then
197 0 : write(iulog, "('mo_apex_init: model yr,mon,day,sec ',4i6)") yr, mon, day, sec
198 : endif
199 0 : write(iulog, "('mo_apex_init: nglons,nglats ', 2i6)") nglons, nglats
200 : endif
201 :
202 0 : end subroutine mo_apex_init1
203 :
204 : !======================================================================
205 : !======================================================================
206 0 : subroutine mo_apex_init(phys_state)
207 : !-------------------------------------------------------------------------------
208 : ! Driver for apex code to calculate apex magnetic coordinates at
209 : ! current geographic spatial resolution for given year. This calls
210 : ! routines in apex_subs.F.
211 : !
212 : ! This is called once per run from sub inti.
213 : !-------------------------------------------------------------------------------
214 :
215 0 : use physconst,only : pi
216 : use physics_types, only: physics_state
217 : use epp_ionization,only: epp_ionization_setmag
218 :
219 : ! Input/output arguments
220 : type(physics_state), intent(in), dimension(begchunk:endchunk) :: phys_state
221 :
222 : !-------------------------------------------------------------------------------
223 : ! Local variables
224 : !-------------------------------------------------------------------------------
225 : real(r8), parameter :: re = 6.378165e8_r8 ! earth radius (cm)
226 : real(r8), parameter :: h0 = 9.0e6_r8 ! base height (90 km)
227 : real(r8), parameter :: hs = 1.3e7_r8
228 : real(r8), parameter :: eps = 1.e-6_r8 ! epsilon
229 : real(r8), parameter :: cm2km = 1.e-5_r8
230 :
231 : integer :: c, i, ist ! indices
232 : integer :: ncol
233 :
234 : real(r8) :: alt, hr, alon, alat, & ! apxmall args
235 : vmp, w, d, be3, sim, xlatqd, f, si, collat, collon
236 :
237 : !-------------------------------------------------------------------------------
238 : ! Non-scalar arguments returned by APXMALL:
239 : !-------------------------------------------------------------------------------
240 : real(r8) :: bhat(3)
241 : real(r8) :: d3(3)
242 : real(r8) :: e1(3), e2(3), e3(3)
243 : real(r8) :: f1(2), f2(2)
244 :
245 : real(r8) :: bg(3), d1g(3), d2g(3), bmg
246 :
247 : real(r8) :: rdum
248 :
249 0 : real(r8) :: maglat(pcols,begchunk:endchunk)
250 :
251 : real(r8), parameter :: rtd = 180._r8/pi ! radians to degrees
252 : real(r8), parameter :: dtr = pi/180._r8 ! degrees to radians
253 :
254 0 : call mo_apex_init1()
255 0 : if ((.not.geomag_year_updated) .and. (allocated(alatm))) return
256 :
257 : !-------------------------------------------------------------------------------
258 : ! Allocate output arrays
259 : !-------------------------------------------------------------------------------
260 0 : call allocate_arrays()
261 :
262 0 : alt = hs*cm2km ! altitude for apxmall (km)
263 0 : hr = alt ! reference altitude (km)
264 :
265 : !------------------------------------------------------------------------------
266 : ! Apex coords alon, alat are returned for each geographic grid point:
267 : ! first form global arrays
268 : !------------------------------------------------------------------------------
269 0 : do c = begchunk, endchunk
270 0 : ncol = phys_state(c)%ncol
271 0 : do i = 1,ncol
272 0 : collat = phys_state(c)%lat(i)*rtd ! latitude of current column (deg)
273 0 : collon = phys_state(c)%lon(i)*rtd ! latitude of current column (deg)
274 0 : if ( collon < -180._r8 ) collon = collon+360._r8
275 0 : if ( collon > 180._r8 ) collon = collon-360._r8
276 : call apex_mall( &
277 : collat, collon, alt, hr, & ! Inputs
278 0 : bg, bhat, bmag(i,c), si, & ! Mag Fld
279 : alon, alat, & ! Apex lon,lat output
280 : vmp, w, d, be3, sim, d1vec(:,i,c), d2vec(:,i,c), d3, e1, e2, e3, & ! Mod Apex
281 0 : xlatqd, f, f1, f2, ist ) ! Qsi-Dpl
282 0 : if( ist /= 0 ) then
283 0 : write(iulog,"(/,'>>> mo_apex_init: Error from apxmall: ist=',i4)") ist
284 0 : call endrun('mo_apex_init: Error from apxmall')
285 : end if
286 0 : beast (i,c) = bg(1)
287 0 : bnorth(i,c) = bg(2)
288 0 : bdown (i,c) = -bg(3)
289 0 : alonm (i,c) = alon*dtr ! mag lons (radians)
290 0 : alatm (i,c) = alat*dtr ! mag lats (radians)
291 0 : maglat(i,c) = alat ! mag lats (degrees)
292 : enddo
293 : enddo
294 :
295 : ! find geograghic latitude ( maglon0 ) where the geomagnetic latitude is zero at the equator
296 : ! by first extracting the geographic coordinates at zero degrees longitude ...
297 0 : collat = 0._r8
298 0 : collon = 0._r8
299 : call apex_mall( &
300 : collat, collon, alt, hr, & ! Inputs
301 : bg, bhat, bmg, si, & ! Mag Fld
302 : alon, alat, & ! Apex lon,lat output
303 : vmp, w, d, be3, sim, d1g, d2g, d3, e1, e2, e3, & ! Mod Apex
304 0 : xlatqd, f, f1, f2, ist ) ! Qsi-Dpl
305 :
306 0 : if( ist /= 0 ) then
307 0 : write(iulog,"(/,'>>> mo_apex_init: Error from apxmall: ist=',i4)") ist
308 0 : call endrun('mo_apex_init: Error from apxmall')
309 : end if
310 :
311 0 : maglon0 = -alon*dtr ! (radians) geograghic latitude where the geomagnetic latitude is zero
312 : ! where longitude ranges from -180E to 180E
313 :
314 0 : call apex_dypol( colatp, elonp, rdum ) ! get geomagnetic dipole north pole
315 :
316 0 : if (masterproc) then
317 0 : write(iulog, "('mo_apex_init: colatp,elonp ', 2f12.6)") colatp, elonp
318 0 : write(iulog, "('mo_apex_init: Calculated apex magnetic coordinates for year AD ',f8.2)") geomag_year
319 : endif
320 :
321 0 : call epp_ionization_setmag(maglat)
322 :
323 0 : end subroutine mo_apex_init
324 :
325 0 : subroutine allocate_arrays
326 : !------------------------------------------------------------------------------
327 : ! Allocate module output arrays for chunked physics grid.
328 : !------------------------------------------------------------------------------
329 :
330 : !------------------------------------------------------------------------------
331 : ! local variables
332 : !------------------------------------------------------------------------------
333 : integer :: istat ! status of allocate statements
334 :
335 0 : if (.not.allocated(alatm)) then
336 0 : allocate(alatm(pcols,begchunk:endchunk),stat=istat)
337 0 : if (istat /= 0) then
338 0 : write(iulog,"('>>> allocate_arrays: allocate of alatm failed: istat=',i5)") istat
339 0 : call endrun
340 : end if
341 : end if
342 :
343 0 : if (.not.allocated(alonm)) then
344 0 : allocate(alonm(pcols,begchunk:endchunk),stat=istat)
345 0 : if (istat /= 0) then
346 0 : write(iulog,"('>>> allocate_arrays: allocate of alonm failed: istat=',i5)") istat
347 0 : call endrun
348 : end if
349 : end if
350 :
351 0 : if (.not.allocated(bnorth)) then
352 0 : allocate(bnorth(pcols,begchunk:endchunk),stat=istat)
353 0 : if (istat /= 0) then
354 0 : write(iulog,"('>>> allocate_arrays: allocate of bnorth failed: istat=',i5)") istat
355 0 : call endrun
356 : end if
357 : end if
358 :
359 0 : if (.not.allocated(beast)) then
360 0 : allocate(beast(pcols,begchunk:endchunk),stat=istat)
361 0 : if (istat /= 0) then
362 0 : write(iulog,"('>>> allocate_arrays: allocate of beast failed: istat=',i5)") istat
363 0 : call endrun
364 : end if
365 : end if
366 :
367 0 : if (.not.allocated(bdown)) then
368 0 : allocate(bdown(pcols,begchunk:endchunk),stat=istat)
369 0 : if (istat /= 0) then
370 0 : write(iulog,"('>>> allocate_arrays: allocate of bdown failed: istat=',i5)") istat
371 0 : call endrun
372 : end if
373 : end if
374 :
375 0 : if (.not.allocated(bmag)) then
376 0 : allocate(bmag(pcols,begchunk:endchunk),stat=istat)
377 0 : if (istat /= 0) then
378 0 : write(iulog,"('>>> allocate_arrays: allocate of bmag failed: istat=',i5)") istat
379 0 : call endrun
380 : end if
381 : end if
382 0 : if (.not.allocated(d1vec)) then
383 0 : allocate(d1vec(3,pcols,begchunk:endchunk),stat=istat)
384 0 : if (istat /= 0) then
385 0 : write(iulog,"('>>> allocate_arrays: allocate of d1vec failed: istat=',i5)") istat
386 0 : call endrun
387 : endif
388 : endif
389 :
390 0 : if (.not.allocated(d2vec)) then
391 0 : allocate(d2vec(3,pcols,begchunk:endchunk),stat=istat)
392 0 : if (istat /= 0) then
393 0 : write(iulog,"('>>> allocate_arrays: allocate of d2vec failed: istat=',i5)") istat
394 0 : call endrun
395 : endif
396 : endif
397 :
398 0 : end subroutine allocate_arrays
399 :
400 : end module mo_apex
|
