Line data Source code
1 : module amie_module
2 : !
3 : ! Module used to read data from the AMIE outputs (POT,mean energy,
4 : ! and energy flux).
5 : !
6 :
7 : use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl
8 : use cam_logfile, only: iulog
9 : use spmd_utils, only: masterproc
10 : use edyn_maggrid, only: nmlat, nmlonp1
11 : use edyn_maggrid, only: ylonm ! magnetic latitudes (nmlat) (radians)
12 : use edyn_maggrid, only: ylatm ! magnetic longtitudes (nmlonp1) (radians)
13 : use cam_pio_utils, only: cam_pio_openfile, cam_pio_closefile
14 : use pio, only: pio_inq_dimid, pio_inquire_dimension
15 : use pio, only: pio_inquire, pio_inq_varid
16 : use pio, only: file_desc_t, pio_noerr, pio_nowrite, pio_get_var
17 : use utils_mod, only: check_ncerr, check_alloc, boxcar_ave
18 : use edyn_mpi, only: ntask, mytid
19 : use edyn_params, only: pi, dtr, rtd
20 : use input_data_utils, only: time_coordinate
21 :
22 : implicit none
23 :
24 : private
25 : public :: init_amie
26 : public :: getamie
27 :
28 : ! Define parameters for AMIE input data file:
29 : integer, parameter :: &
30 : ithmx = 55, & ! maximum number of latitudes of AMIE data
31 : jmxm = 2*ithmx-1, & ! maximum number of global latitudes
32 : lonmx = 36 ! maximum number of longitudes of AMIE data
33 : integer :: lonp1,latp1
34 : !
35 : ! Define fields for AMIE input data file:
36 : ! electric potential in Volt
37 : ! mean energy in KeV
38 : ! energy flux in W/m^2
39 : ! cusplat_nh_input(sh) and cuspmlt_nh_input(sh) are
40 : ! AMIE cusp latitude and MLT in NH and SH
41 : ! hpi_nh(sh) are AMIE hemi-integrated power
42 : ! pcp_nh(sh) are AMIE polar-cap potential drop
43 : ! Time interpolated AMIE outputs with suffix _amie
44 : !
45 : real(r8), allocatable, dimension(:,:,:), save :: & ! (lonp1,latp1,ntimes)
46 : pot_nh_input, pot_sh_input, &
47 : ekv_nh_input, ekv_sh_input, &
48 : efx_nh_input, efx_sh_input
49 : real(r8), allocatable, dimension(:,:), save :: & ! (lonp1,latp1)
50 : pot_nh_amie, pot_sh_amie, ekv_nh_amie, ekv_sh_amie, &
51 : efx_nh_amie, efx_sh_amie
52 : integer, allocatable, dimension(:), save :: & ! (ntimes)
53 : year, month, day, jday
54 : real(r8), allocatable, dimension(:), save :: & ! (ntimes)
55 : cusplat_nh_input, cuspmlt_nh_input, hpi_nh_input, &
56 : pcp_nh_input, amie_nh_ut, &
57 : cusplat_sh_input, cuspmlt_sh_input, hpi_sh_input, &
58 : pcp_sh_input, amie_sh_ut
59 : real(r8) :: &
60 : cusplat_nh_amie, cuspmlt_nh_amie, cusplat_sh_amie, &
61 : cuspmlt_sh_amie, hpi_sh_amie, hpi_nh_amie, pcp_sh_amie, &
62 : pcp_nh_amie
63 : !
64 : type(file_desc_t) :: ncid_nh
65 : type(file_desc_t) :: ncid_sh
66 :
67 : character(len=cl), allocatable :: amienh_files(:)
68 : character(len=cl), allocatable :: amiesh_files(:)
69 : integer :: num_files, file_ndx
70 :
71 : type(time_coordinate) :: time_coord_nh
72 : type(time_coordinate) :: time_coord_sh
73 :
74 : contains
75 :
76 : !-----------------------------------------------------------------------
77 0 : subroutine init_amie(amienh_list,amiesh_list)
78 :
79 : character(len=*),intent(in) :: amienh_list(:)
80 : character(len=*),intent(in) :: amiesh_list(:)
81 :
82 : integer :: n, nfiles
83 :
84 0 : nfiles = min( size(amienh_list), size(amiesh_list) )
85 0 : num_files = 0
86 :
87 0 : count_files: do n = 1,nfiles
88 0 : if (len_trim(amienh_list(n))<1 .or. len_trim(amiesh_list(n))<1 .or. &
89 0 : trim(amienh_list(n))=='NONE' .or. trim(amiesh_list(n))=='NONE') then
90 : exit count_files
91 : else
92 0 : num_files = num_files + 1
93 : end if
94 : end do count_files
95 :
96 0 : allocate(amienh_files(num_files), amiesh_files(num_files))
97 0 : amienh_files(:num_files) = amienh_list(:num_files)
98 0 : amiesh_files(:num_files) = amiesh_list(:num_files)
99 0 : file_ndx = 1
100 0 : call open_files()
101 :
102 0 : end subroutine init_amie
103 :
104 : !-----------------------------------------------------------------------
105 0 : subroutine rdamie_nh(amienh)
106 : !
107 : ! Read AMIE data for the northern hemisphere from amienh
108 : !
109 :
110 : ! Dummy argument
111 : character(len=*), intent(in) :: amienh
112 : ! Local variables:
113 : integer :: istat, ntimes, ndims, nvars, ngatts
114 : integer :: idunlim, ier
115 : integer :: id_lon, id_lat, id_time
116 : integer :: idv_year, idv_mon, idv_day, idv_jday
117 : integer :: idv_ut, idv_cusplat, idv_cuspmlt
118 : integer :: idv_hpi, idv_pcp
119 : character(len=*), parameter :: subname = 'rdamie_nh'
120 : !
121 0 : if (masterproc) then
122 0 : write(iulog, "(/,72('-'))")
123 0 : write(iulog, "(a,': read AMIE data for northern hemisphere:')") subname
124 : end if
125 : !
126 : ! Open netcdf file:
127 0 : call cam_pio_openfile(ncid_nh, amienh, pio_nowrite)
128 : !
129 : ! Get AMIE grid dimension:
130 0 : istat = pio_inq_dimid(ncid_nh, 'lon', id_lon)
131 0 : istat = pio_inquire_dimension(ncid_nh, id_lon, len=lonp1)
132 0 : call check_ncerr(istat, subname, 'AMIE longitude dimension')
133 :
134 0 : istat = pio_inq_dimid(ncid_nh, 'lat', id_lat)
135 0 : istat = pio_inquire_dimension(ncid_nh, id_lat, len=latp1)
136 0 : call check_ncerr(istat, subname, 'AMIE latitude dimension')
137 :
138 0 : call time_coord_nh%initialize( amienh, set_weights=.false. )
139 :
140 : !
141 : ! Get time dimension:
142 0 : istat = pio_inquire(ncid_nh, unlimiteddimid=id_time)
143 0 : istat = pio_inquire_dimension(ncid_nh, id_time, len=ntimes)
144 0 : call check_ncerr(istat, subname, 'AMIE time dimension')
145 : !
146 : ! Search for requested AMIE output fields
147 0 : istat = pio_inquire(ncid_nh, ndims, nvars, ngatts, idunlim)
148 : !
149 : ! Get 1-D AMIE fields (ntimes)
150 0 : if (.not. allocated(year)) then
151 0 : allocate(year(ntimes), stat=ier)
152 0 : call check_alloc(ier, subname, 'year', ntimes=ntimes)
153 : end if
154 0 : istat = pio_inq_varid(ncid_nh, 'year', idv_year)
155 0 : call check_ncerr(istat, subname, 'AMIE year id')
156 0 : istat = pio_get_var(ncid_nh, idv_year, year)
157 0 : call check_ncerr(istat, subname, 'AMIE year')
158 :
159 0 : if (.not. allocated(month)) then
160 0 : allocate(month(ntimes), stat=ier)
161 0 : call check_alloc(ier, subname, 'month', ntimes=ntimes)
162 : end if
163 0 : istat = pio_inq_varid(ncid_nh, 'month', idv_mon)
164 0 : call check_ncerr(istat, subname, 'AMIE month id')
165 0 : istat = pio_get_var(ncid_nh, idv_mon, month)
166 0 : call check_ncerr(istat, subname, 'AMIE month')
167 0 : if (.not. allocated(day)) then
168 0 : allocate(day(ntimes), stat=ier)
169 0 : call check_alloc(ier, subname, 'day', ntimes=ntimes)
170 : end if
171 0 : istat = pio_inq_varid(ncid_nh, 'day', idv_day)
172 0 : call check_ncerr(istat, subname, 'AMIE day id')
173 0 : istat = pio_get_var(ncid_nh, idv_day, day)
174 0 : call check_ncerr(istat, subname, 'AMIE day')
175 :
176 0 : if (.not. allocated(jday)) then
177 0 : allocate(jday(ntimes), stat=ier)
178 0 : call check_alloc(ier, subname, 'jday', ntimes=ntimes)
179 : end if
180 0 : istat = pio_inq_varid(ncid_nh, 'jday', idv_jday)
181 0 : call check_ncerr(istat, subname, 'AMIE jday id')
182 0 : istat = pio_get_var(ncid_nh, idv_jday, jday)
183 0 : call check_ncerr(istat, subname, 'AMIE jday')
184 : !
185 : ! Allocate 1-d fields:
186 0 : if (.not. allocated(amie_nh_ut)) then
187 0 : allocate(amie_nh_ut(ntimes), stat=ier)
188 0 : call check_alloc(ier, subname, 'amie_nh_ut', ntimes=ntimes)
189 : end if
190 0 : if (.not. allocated(cusplat_nh_input)) then
191 0 : allocate(cusplat_nh_input(ntimes), stat=ier)
192 0 : call check_alloc(ier, subname, 'cusplat_nh_input', ntimes=ntimes)
193 : end if
194 0 : if (.not. allocated(cuspmlt_nh_input)) then
195 0 : allocate(cuspmlt_nh_input(ntimes), stat=ier)
196 0 : call check_alloc(ier, subname, 'cuspmlt_nh_input', ntimes=ntimes)
197 : end if
198 0 : if (.not. allocated(hpi_nh_input)) then
199 0 : allocate(hpi_nh_input(ntimes), stat=ier)
200 0 : call check_alloc(ier, subname, 'hpi_nh_input', ntimes=ntimes)
201 : end if
202 0 : if (.not. allocated(pcp_nh_input)) then
203 0 : allocate(pcp_nh_input(ntimes), stat=ier)
204 0 : call check_alloc(ier, subname, 'pcp_nh_input', ntimes=ntimes)
205 : end if
206 : !
207 : ! Get ut
208 0 : istat = pio_inq_varid(ncid_nh, 'ut', idv_ut)
209 0 : call check_ncerr(istat, subname, 'AMIE ut id')
210 0 : istat = pio_get_var(ncid_nh, idv_ut, amie_nh_ut)
211 0 : call check_ncerr(istat, subname, 'AMIE ut')
212 : !
213 : ! Get HPI
214 0 : istat = pio_inq_varid(ncid_nh, 'hpi', idv_hpi)
215 0 : call check_ncerr(istat, subname, 'AMIE hpi id')
216 0 : istat = pio_get_var(ncid_nh, idv_hpi, hpi_nh_input)
217 0 : call check_ncerr(istat, subname, 'AMIE hpi')
218 : !
219 : ! Get PCP
220 0 : istat = pio_inq_varid(ncid_nh, 'pcp', idv_pcp)
221 0 : call check_ncerr(istat, subname, 'AMIE pcp id')
222 0 : istat = pio_get_var(ncid_nh, idv_pcp, pcp_nh_input)
223 0 : call check_ncerr(istat, subname, 'AMIE pcp')
224 : !
225 : ! Get cusplat
226 0 : istat = pio_inq_varid(ncid_nh, 'cusplat', idv_cusplat)
227 0 : call check_ncerr(istat, subname, 'AMIE cusplat id')
228 0 : istat = pio_get_var(ncid_nh, idv_cusplat, cusplat_nh_input)
229 0 : call check_ncerr(istat, subname, 'AMIE cusplat')
230 : !
231 : ! Get cuspmlt
232 0 : istat = pio_inq_varid(ncid_nh, 'cuspmlt', idv_cuspmlt)
233 0 : call check_ncerr(istat, subname, 'AMIE cuspmlt id')
234 0 : istat = pio_get_var(ncid_nh, idv_cuspmlt, cuspmlt_nh_input)
235 0 : call check_ncerr(istat, subname, 'AMIE cuspmlt')
236 : !
237 : ! Allocate 2-d fields:
238 0 : if (.not. allocated(pot_nh_amie)) then
239 0 : allocate(pot_nh_amie(lonp1, latp1), stat=ier)
240 0 : call check_alloc(ier, subname, 'pot_nh_amie', lonp1=lonp1, latp1=latp1)
241 : end if
242 0 : if (.not. allocated(ekv_nh_amie)) then
243 0 : allocate(ekv_nh_amie(lonp1, latp1), stat=ier)
244 0 : call check_alloc(ier, subname, 'ekv_nh_amie', lonp1=lonp1, latp1=latp1)
245 : end if
246 0 : if (.not. allocated(efx_nh_amie)) then
247 0 : allocate(efx_nh_amie(lonp1, latp1), stat=ier)
248 0 : call check_alloc(ier, subname, 'efx_nh_amie', lonp1=lonp1, latp1=latp1)
249 : end if
250 : !
251 : ! Allocate 3-d fields:
252 0 : if (.not. allocated(pot_nh_input)) then
253 0 : allocate(pot_nh_input(lonp1, latp1, 2), stat=ier)
254 : call check_alloc(ier, subname, 'pot_nh_input', &
255 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
256 : end if
257 0 : if (.not. allocated(ekv_nh_input)) then
258 0 : allocate(ekv_nh_input(lonp1, latp1, 2), stat=ier)
259 : call check_alloc(ier, subname, 'ekv_nh_input', &
260 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
261 : end if
262 0 : if (.not. allocated(efx_nh_input)) then
263 0 : allocate(efx_nh_input(lonp1, latp1, 2), stat=ier)
264 : call check_alloc(ier, subname, 'efx_nh_input', &
265 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
266 : end if
267 0 : end subroutine rdamie_nh
268 :
269 : !-----------------------------------------------------------------------
270 0 : subroutine rdamie_sh(amiesh)
271 : !
272 : ! Read AMIE data for the southern hemisphere from amiesh
273 : !
274 :
275 : ! Dummy argument
276 : character(len=*), intent(in) :: amiesh
277 : ! Local variables:
278 : integer :: istat, ntimes, ndims, nvars, ngatts, ier
279 : integer :: idunlim
280 : integer :: id_lon, id_lat, id_time
281 : integer :: idv_year, idv_mon, idv_day, idv_jday
282 : integer :: idv_ut
283 : integer :: idv_cusplat, idv_cuspmlt, idv_hpi, idv_pcp
284 : character(len=*), parameter :: subname = 'rdamie_sh'
285 : !
286 0 : if (masterproc) then
287 0 : write(iulog, "(/, 72('-'))")
288 0 : write(iulog, "(a, ': read AMIE data for southern hemisphere:')") subname
289 : end if
290 : !
291 : ! Open netcdf file:
292 0 : call cam_pio_openfile(ncid_sh, amiesh, pio_nowrite)
293 : !
294 : ! Get AMIE grid dimension:
295 0 : istat = pio_inq_dimid(ncid_sh, 'lon', id_lon)
296 0 : istat = pio_inquire_dimension(ncid_sh, id_lon, len=lonp1)
297 0 : call check_ncerr(istat, subname, 'AMIE longitude dimension')
298 :
299 0 : istat = pio_inq_dimid(ncid_sh, 'lat', id_lat)
300 0 : istat = pio_inquire_dimension(ncid_sh, id_lat, len=latp1)
301 0 : call check_ncerr(istat, subname, 'AMIE latitude dimension')
302 :
303 0 : call time_coord_sh%initialize( amiesh, set_weights=.false. )
304 :
305 : !
306 : ! Get time dimension:
307 0 : istat = pio_inquire(ncid_sh, unlimiteddimid=id_time)
308 0 : istat = pio_inquire_dimension(ncid_sh, id_time, len=ntimes)
309 0 : call check_ncerr(istat, subname, 'AMIE time dimension')
310 : !
311 : ! Search for requested AMIE output fields
312 0 : istat = pio_inquire(ncid_sh, ndims, nvars, ngatts, idunlim)
313 : !
314 : ! Get 1-D AMIE fields (ntimes)
315 0 : if (.not. allocated(year)) then
316 0 : allocate(year(ntimes), stat=ier)
317 0 : call check_alloc(ier, subname, 'year', ntimes=ntimes)
318 : end if
319 0 : istat = pio_inq_varid(ncid_sh, 'year', idv_year)
320 0 : call check_ncerr(istat, subname, 'AMIE year id')
321 0 : istat = pio_get_var(ncid_sh, idv_year, year)
322 0 : call check_ncerr(istat, subname, 'AMIE year')
323 :
324 0 : if (.not. allocated(month)) then
325 0 : allocate(month(ntimes), stat=ier)
326 0 : call check_alloc(ier, subname, 'month', ntimes=ntimes)
327 : end if
328 0 : istat = pio_inq_varid(ncid_sh, 'month', idv_mon)
329 0 : call check_ncerr(istat, subname, 'AMIE month id')
330 0 : istat = pio_get_var(ncid_sh, idv_mon, month)
331 0 : call check_ncerr(istat, subname, 'AMIE month')
332 0 : if (.not. allocated(day)) then
333 0 : allocate(day(ntimes), stat=ier)
334 0 : call check_alloc(ier, subname, 'day', ntimes=ntimes)
335 : end if
336 0 : istat = pio_inq_varid(ncid_sh, 'day', idv_day)
337 0 : call check_ncerr(istat, subname, 'AMIE day id')
338 0 : istat = pio_get_var(ncid_sh, idv_day, day)
339 0 : call check_ncerr(istat, subname, 'AMIE day')
340 :
341 0 : if (.not. allocated(jday)) then
342 0 : allocate(jday(ntimes), stat=ier)
343 0 : call check_alloc(ier, subname, 'jday', ntimes=ntimes)
344 : end if
345 0 : istat = pio_inq_varid(ncid_sh, 'jday', idv_jday)
346 0 : call check_ncerr(istat, subname, 'AMIE jday id')
347 0 : istat = pio_get_var(ncid_sh, idv_jday, jday)
348 0 : call check_ncerr(istat, subname, 'AMIE jday')
349 : !
350 : ! Allocate 1-d fields:
351 0 : if (.not. allocated(amie_sh_ut)) then
352 0 : allocate(amie_sh_ut(ntimes), stat=ier)
353 0 : call check_alloc(ier, subname, 'amie_sh_ut', ntimes=ntimes)
354 : end if
355 0 : if (.not. allocated(cusplat_sh_input)) then
356 0 : allocate(cusplat_sh_input(ntimes), stat=ier)
357 0 : call check_alloc(ier, subname, 'cusplat_sh_input', ntimes=ntimes)
358 : end if
359 0 : if (.not. allocated(cuspmlt_sh_input)) then
360 0 : allocate(cuspmlt_sh_input(ntimes), stat=ier)
361 0 : call check_alloc(ier, subname, 'cuspmlt_sh_input', ntimes=ntimes)
362 : end if
363 0 : if (.not. allocated(hpi_sh_input)) then
364 0 : allocate(hpi_sh_input(ntimes), stat=ier)
365 0 : call check_alloc(ier, subname, 'hpi_sh_input', ntimes=ntimes)
366 : end if
367 0 : if (.not. allocated(pcp_sh_input)) then
368 0 : allocate(pcp_sh_input(ntimes), stat=ier)
369 0 : call check_alloc(ier, subname, 'pcp_sh_input', ntimes=ntimes)
370 : end if
371 : !
372 : ! Get ut
373 0 : istat = pio_inq_varid(ncid_sh, 'ut', idv_ut)
374 0 : call check_ncerr(istat, subname, 'AMIE ut id')
375 0 : istat = pio_get_var(ncid_sh, idv_ut, amie_sh_ut)
376 0 : call check_ncerr(istat, subname, 'AMIE ut')
377 : !
378 : ! Get HPI
379 0 : istat = pio_inq_varid(ncid_sh, 'hpi', idv_hpi)
380 0 : call check_ncerr(istat, subname, 'AMIE hpi id')
381 0 : istat = pio_get_var(ncid_sh, idv_hpi, hpi_sh_input)
382 0 : call check_ncerr(istat, subname, 'AMIE hpi')
383 : !
384 : ! Get PCP
385 0 : istat = pio_inq_varid(ncid_sh, 'pcp', idv_pcp)
386 0 : call check_ncerr(istat, subname, 'AMIE pcp id')
387 0 : istat = pio_get_var(ncid_sh, idv_pcp, pcp_sh_input)
388 0 : call check_ncerr(istat, subname, 'AMIE pcp')
389 : !
390 : ! Get cusplat
391 0 : istat = pio_inq_varid(ncid_sh, 'cusplat', idv_cusplat)
392 0 : call check_ncerr(istat, subname, 'AMIE cusplat id')
393 0 : istat = pio_get_var(ncid_sh, idv_cusplat, cusplat_sh_input)
394 0 : call check_ncerr(istat, subname, 'AMIE cusplat')
395 : !
396 : ! Get cuspmlt
397 0 : istat = pio_inq_varid(ncid_sh, 'cuspmlt', idv_cuspmlt)
398 0 : call check_ncerr(istat, subname, 'AMIE cuspmlt id')
399 0 : istat = pio_get_var(ncid_sh, idv_cuspmlt, cuspmlt_sh_input)
400 0 : call check_ncerr(istat, subname, 'AMIE cuspmlt')
401 : !
402 : ! Allocate 2-d fields:
403 0 : if (.not. allocated(pot_sh_amie)) then
404 0 : allocate(pot_sh_amie(lonp1, latp1), stat=ier)
405 0 : call check_alloc(ier, subname, 'pot_sh_amie', lonp1=lonp1, latp1=latp1)
406 : end if
407 0 : if (.not. allocated(ekv_sh_amie)) then
408 0 : allocate(ekv_sh_amie(lonp1, latp1), stat=ier)
409 0 : call check_alloc(ier, subname, 'ekv_sh_amie', lonp1=lonp1, latp1=latp1)
410 : end if
411 0 : if (.not. allocated(efx_sh_amie)) then
412 0 : allocate(efx_sh_amie(lonp1, latp1), stat=ier)
413 0 : call check_alloc(ier, subname, 'efx_sh_amie', lonp1=lonp1, latp1=latp1)
414 : end if
415 : !
416 : ! Allocate 3-d fields:
417 0 : if (.not. allocated(pot_sh_input)) then
418 0 : allocate(pot_sh_input(lonp1, latp1, 2), stat=ier)
419 : call check_alloc(ier, subname, 'pot_sh_input', &
420 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
421 : end if
422 0 : if (.not. allocated(ekv_sh_input)) then
423 0 : allocate(ekv_sh_input(lonp1, latp1, 2), stat=ier)
424 : call check_alloc(ier, subname, 'ekv_sh_input', &
425 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
426 : end if
427 0 : if (.not. allocated(efx_sh_input)) then
428 0 : allocate(efx_sh_input(lonp1, latp1, 2), stat=ier)
429 : call check_alloc(ier, subname, 'efx_sh_input', &
430 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
431 : end if
432 0 : end subroutine rdamie_sh
433 :
434 : !-----------------------------------------------------------------------
435 0 : subroutine update_3d_fields( ncid, offset, kount, pot_3d,ekv_3d,efx_3d )
436 :
437 : type(file_desc_t), intent(in) :: ncid
438 : integer, intent(in) :: offset(:)
439 : integer, intent(in) :: kount(:)
440 : real(r8),intent(out) :: pot_3d(:,:,:)
441 : real(r8),intent(out) :: ekv_3d(:,:,:)
442 : real(r8),intent(out) :: efx_3d(:,:,:)
443 :
444 :
445 : integer :: istat
446 : integer :: idv_pot, idv_ekv, idv_efx
447 : character(len=*), parameter :: subname = 'update_3d_fields'
448 :
449 : !
450 : ! Get 3-D fields (lon,lat,ntimes)
451 : !
452 : ! electric potential
453 0 : istat = pio_inq_varid(ncid, 'pot', idv_pot)
454 0 : call check_ncerr(istat, subname, 'AMIE pot id')
455 0 : istat = pio_get_var(ncid, idv_pot, offset, kount, pot_3d)
456 0 : call check_ncerr(istat, subname, 'AMIE pot')
457 : !
458 : ! mean energy
459 0 : istat = pio_inq_varid(ncid, 'ekv', idv_ekv)
460 0 : call check_ncerr(istat, subname, 'AMIE ekv id')
461 0 : istat = pio_get_var(ncid, idv_ekv, offset, kount, ekv_3d)
462 0 : call check_ncerr(istat, subname, 'AMIE ekv')
463 : !
464 : ! energy flux
465 0 : istat = pio_inq_varid(ncid, 'efx', idv_efx)
466 0 : call check_ncerr(istat, subname, 'AMIE efx id')
467 0 : istat = pio_get_var(ncid, idv_efx, offset, kount, efx_3d)
468 0 : call check_ncerr(istat, subname, 'AMIE efx')
469 :
470 0 : end subroutine update_3d_fields
471 :
472 : !-----------------------------------------------------------------------
473 0 : subroutine getamie(iyear, imo, iday, iutsec, sunlon, iprint, &
474 0 : iamie, phihm, amie_efxm, amie_kevm, crad)
475 : use cam_history_support, only: fillvalue
476 : use rgrd_mod, only: rgrd2
477 :
478 : !
479 : ! Read AMIE outputs from amie_ncfile file, returning electric potential,
480 : ! auroral mean energy and energy flux at current date and time,
481 : ! and the data is linearly interpolated to the model time
482 : ! gl - 12/07/2002
483 : !
484 : !
485 : ! Args:
486 :
487 : integer, intent(in) :: iyear
488 : integer, intent(in) :: imo
489 : integer, intent(in) :: iday
490 : real(r8), intent(in) :: sunlon
491 : integer, intent(in) :: iutsec
492 : integer, intent(in) :: iprint
493 : integer, intent(out) :: iamie
494 : real(r8), intent(out) :: phihm(nmlonp1,nmlat)
495 : real(r8), intent(out) :: amie_efxm(nmlonp1,nmlat) ! on geomag grid
496 : real(r8), intent(out) :: amie_kevm(nmlonp1,nmlat) ! on geomag grid
497 : real(r8), intent(out) :: crad(2)
498 : !
499 : !
500 : ! Local:
501 0 : real(r8) :: potm(lonp1,jmxm)
502 0 : real(r8) :: efxm(lonp1,jmxm), ekvm(lonp1,jmxm)
503 0 : real(r8) :: alat(jmxm), alon(lonp1)
504 0 : real(r8) :: alatm(jmxm), alonm(lonp1)
505 : integer :: ier, lw, liw, intpol(2)
506 0 : integer, allocatable :: iw(:)
507 0 : real(r8), allocatable :: w(:)
508 : integer :: i, j
509 : integer :: nn, iset1, iset2, m, mp1, n
510 : integer :: iboxcar
511 : real(r8) :: f1, f2
512 : real(r8) :: del, xmlt, dmlat, dlatm, dlonm, dmltm, rot
513 : integer :: offset(3), kount(3)
514 : character(len=*), parameter :: subname = 'getamie'
515 :
516 0 : phihm = fillvalue
517 0 : amie_efxm = fillvalue
518 0 : amie_kevm = fillvalue
519 0 : crad = fillvalue
520 :
521 0 : if (iprint > 0 .and. masterproc) then
522 0 : write(iulog,"(/,72('-'))")
523 0 : write(iulog,"(a,':')") subname
524 : write(iulog,"(a,i4,', iday = ',i3,', iutsec = ',i10)") &
525 0 : 'Initial requested iyear= ', iyear, iday, iutsec
526 : end if
527 :
528 0 : nn = size(amie_sh_ut)
529 : !
530 : ! Check times:
531 : !
532 :
533 0 : iamie = 1 - time_coord_sh%times_check()
534 0 : check_loop: do while( iamie/=1 )
535 0 : if (iamie==2) then
536 0 : if (masterproc) then
537 : write(iulog, "(a,': Model date prior to AMIE first date:',3I5)") &
538 0 : subname, year(1), month(1), day(1)
539 : end if
540 0 : return
541 : end if
542 :
543 0 : if (iamie==0) then
544 0 : if (masterproc) then
545 : write(iulog, "(a,': Model date beyond the AMIE last Data:',3I5)") &
546 0 : subname, year(nn), month(nn), day(nn)
547 : end if
548 :
549 0 : if (file_ndx<num_files) then
550 0 : file_ndx = file_ndx+1
551 0 : call close_files()
552 0 : call open_files()
553 0 : iamie = 1 - time_coord_sh%times_check()
554 : else
555 : return
556 : end if
557 : end if
558 : end do check_loop
559 :
560 : ! get SH AMIE data
561 0 : pot_sh_amie(:,:) = 0._r8
562 0 : ekv_sh_amie(:,:) = 0._r8
563 0 : efx_sh_amie(:,:) = 0._r8
564 0 : cusplat_sh_amie = 0._r8
565 0 : cuspmlt_sh_amie = 0._r8
566 0 : hpi_sh_amie = 0._r8
567 0 : pcp_sh_amie = 0._r8
568 : !
569 :
570 0 : iboxcar = 0
571 :
572 0 : call time_coord_sh%advance()
573 :
574 0 : iset1 = time_coord_sh%indxs(1)
575 0 : iset2 = time_coord_sh%indxs(2)
576 :
577 0 : f1 = time_coord_sh%wghts(1)
578 0 : f2 = time_coord_sh%wghts(2)
579 :
580 0 : cusplat_sh_amie = (f1*cusplat_sh_input(iset1) + &
581 0 : f2*cusplat_sh_input(iset2))
582 0 : cuspmlt_sh_amie = (f1*cuspmlt_sh_input(iset1) + &
583 0 : f2*cuspmlt_sh_input(iset2))
584 0 : hpi_sh_amie = (f1*hpi_sh_input(iset1) + f2*hpi_sh_input(iset2))
585 0 : pcp_sh_amie = (f1*pcp_sh_input(iset1) + f2*pcp_sh_input(iset2))
586 :
587 0 : offset = (/1,1,iset1/)
588 0 : kount = (/lonp1,latp1,2/)
589 :
590 0 : call update_3d_fields( ncid_sh, offset, kount, pot_sh_input,ekv_sh_input,efx_sh_input )
591 0 : if (iboxcar == 0) then
592 0 : pot_sh_amie(:,:) = (f1*pot_sh_input(:,:,1) + &
593 0 : f2*pot_sh_input(:,:,2))
594 0 : ekv_sh_amie(:,:) = (f1*ekv_sh_input(:,:,1) + &
595 0 : f2*ekv_sh_input(:,:,2))
596 0 : efx_sh_amie(:,:) = (f1*efx_sh_input(:,:,1) + &
597 0 : f2*efx_sh_input(:,:,2))
598 : else
599 : call boxcar_ave(pot_sh_input,pot_sh_amie,lonp1,latp1, &
600 0 : nn,iset1,iboxcar)
601 : call boxcar_ave(efx_sh_input,efx_sh_amie,lonp1,latp1, &
602 0 : nn,iset1,iboxcar)
603 : call boxcar_ave(ekv_sh_input,ekv_sh_amie,lonp1,latp1, &
604 0 : nn,iset1,iboxcar)
605 : end if
606 : !
607 : ! get NH AMIE data
608 0 : pot_nh_amie(:,:) = 0._r8
609 0 : ekv_nh_amie(:,:) = 0._r8
610 0 : efx_nh_amie(:,:) = 0._r8
611 0 : cusplat_nh_amie = 0._r8
612 0 : cuspmlt_nh_amie = 0._r8
613 0 : hpi_nh_amie = 0._r8
614 0 : pcp_nh_amie = 0._r8
615 :
616 0 : iboxcar = 0
617 :
618 0 : call time_coord_nh%advance()
619 :
620 0 : iset1 = time_coord_nh%indxs(1)
621 0 : iset2 = time_coord_nh%indxs(2)
622 :
623 0 : f1 = time_coord_nh%wghts(1)
624 0 : f2 = time_coord_nh%wghts(2)
625 :
626 0 : cusplat_nh_amie = (f1*cusplat_nh_input(iset1) + &
627 0 : f2*cusplat_nh_input(iset2))
628 0 : cuspmlt_nh_amie = (f1*cuspmlt_nh_input(iset1) + &
629 0 : f2*cuspmlt_nh_input(iset2))
630 0 : hpi_nh_amie = (f1*hpi_nh_input(iset1) + f2*hpi_nh_input(iset2))
631 0 : pcp_nh_amie = (f1*pcp_nh_input(iset1) + f2*pcp_nh_input(iset2))
632 :
633 0 : offset = (/1,1,iset1/)
634 0 : kount = (/lonp1,latp1,2/)
635 :
636 0 : call update_3d_fields( ncid_nh, offset, kount, pot_nh_input,ekv_nh_input,efx_nh_input )
637 :
638 0 : if (iboxcar == 0) then
639 0 : pot_nh_amie(:,:) = (f1*pot_nh_input(:,:,1) + &
640 0 : f2*pot_nh_input(:,:,2))
641 0 : ekv_nh_amie(:,:) = (f1*ekv_nh_input(:,:,1) + &
642 0 : f2*ekv_nh_input(:,:,2))
643 0 : efx_nh_amie(:,:) = (f1*efx_nh_input(:,:,1) + &
644 0 : f2*efx_nh_input(:,:,2))
645 : else
646 : call boxcar_ave(pot_nh_input,pot_nh_amie,lonp1,latp1, &
647 0 : nn,iset1,iboxcar)
648 : call boxcar_ave(efx_nh_input,efx_nh_amie,lonp1,latp1, &
649 0 : nn,iset1,iboxcar)
650 : call boxcar_ave(ekv_nh_input,ekv_nh_amie,lonp1,latp1, &
651 0 : nn,iset1,iboxcar)
652 : end if
653 : !
654 : ! The OLTMAX latitude also defines the co-latitude theta0, which in
655 : ! turn determines crit1(+2.5deg) and crit2(-12.5deg) which are used
656 : ! in TIE-GCM as the boundaries of the polar cap and the region of
657 : ! influence of the high-lat potential versus the low-lat dynamo potential
658 : ! Define this latitude to be between 70 and 77.5 degrees
659 : !
660 0 : if (cusplat_sh_amie > 75.0_r8) then
661 0 : cusplat_sh_amie = 75.0_r8
662 0 : cuspmlt_sh_amie = 11._r8
663 : end if
664 0 : if (cusplat_sh_amie < 60.0_r8) then
665 0 : cusplat_sh_amie = 60.0_r8
666 0 : cuspmlt_sh_amie = 11._r8
667 : end if
668 0 : if (cusplat_nh_amie > 75.0_r8) then
669 0 : cusplat_nh_amie = 75.0_r8
670 0 : cuspmlt_nh_amie = 11._r8
671 : end if
672 0 : if (cusplat_nh_amie < 60.0_r8) then
673 0 : cusplat_nh_amie = 60.0_r8
674 0 : cuspmlt_nh_amie = 11._r8
675 : end if
676 : ! cusplat_nh_amie = amin1(65.0,cusplat_nh_amie)
677 0 : if (cuspmlt_sh_amie > 12.5_r8) cuspmlt_sh_amie = 12.5_r8
678 0 : if (cuspmlt_sh_amie < 11.0_r8) cuspmlt_sh_amie = 11.0_r8
679 0 : if (cuspmlt_nh_amie > 12.5_r8) cuspmlt_nh_amie = 12.5_r8
680 0 : if (cuspmlt_nh_amie < 11.0_r8) cuspmlt_nh_amie = 11.0_r8
681 0 : crad(1) = (90._r8-cusplat_sh_amie)*pi/180._r8
682 0 : crad(2) = (90._r8-cusplat_nh_amie)*pi/180._r8
683 :
684 0 : active_task: if ( mytid<ntask ) then
685 :
686 : ! mlongitude starts from 180 degree
687 0 : rot = sunlon*rtd
688 0 : if(rot < 0) then
689 0 : rot = rot + 360._r8 ! 0 to 360 degrees
690 : end if
691 0 : rot = rot / 15._r8 ! convert from degree to hrs
692 :
693 0 : dmltm = 24._r8 / real(lonmx, kind=r8)
694 0 : do i = 1, lonp1
695 0 : xmlt = (real(i-1, kind=r8) * dmltm) - rot + 24._r8
696 0 : xmlt = MOD(xmlt, 24._r8)
697 0 : m = int(xmlt/dmltm + 1.01_r8)
698 0 : mp1 = m + 1
699 0 : if (mp1 > lonp1) mp1 = 2
700 0 : del = xmlt - (m-1)*dmltm
701 : ! Initialize arrays around equator
702 0 : do j = latp1+1, ithmx
703 0 : potm(i,j) = 0._r8
704 0 : potm(i,jmxm+1-j) = 0._r8
705 0 : ekvm(i,j) = (1._r8-del)*ekv_sh_amie(m,latp1) + &
706 0 : del*ekv_sh_amie(mp1,latp1)
707 0 : ekvm(i,jmxm+1-j) = (1._r8-del)*ekv_nh_amie(m,latp1) + &
708 0 : del*ekv_nh_amie(mp1,latp1)
709 0 : efxm(i,j) = 0._r8
710 0 : efxm(i,jmxm+1-j) = 0._r8
711 : end do
712 : ! Put in AMIE arrays from pole to latp1
713 0 : do j = 1, latp1
714 0 : potm(i,j) = (1._r8-del)*pot_sh_amie(m,j) + &
715 0 : del*pot_sh_amie(mp1,j)
716 0 : potm(i,jmxm+1-j) = (1._r8-del)*pot_nh_amie(m,j) + &
717 0 : del*pot_nh_amie(mp1,j)
718 0 : ekvm(i,j) = (1._r8-del)*ekv_sh_amie(m,j) + &
719 0 : del*ekv_sh_amie(mp1,j)
720 0 : ekvm(i,jmxm+1-j) = (1._r8-del)*ekv_nh_amie(m,j) + &
721 0 : del*ekv_nh_amie(mp1,j)
722 0 : efxm(i,j) = (1._r8-del)*efx_sh_amie(m,j) + &
723 0 : del*efx_sh_amie(mp1,j)
724 0 : efxm(i,jmxm+1-j) = (1._r8-del)*efx_nh_amie(m,j) + &
725 0 : del*efx_nh_amie(mp1,j)
726 : end do
727 :
728 : end do
729 :
730 : ! Set up coeffs to go between EPOTM(IMXMP,JMNH) and TIEPOT(IMAXM,JMAXMH)
731 :
732 : ! **** SET GRID SPACING DLATM, DLONM
733 : ! DMLAT=lat spacing in degrees of AMIE apex grid
734 0 : dmlat = 180._r8 / real(jmxm-1, kind=r8)
735 0 : dlatm = dmlat * dtr
736 0 : dlonm = 2._r8 * pi / real(lonmx, kind=r8)
737 0 : dmltm = 24._r8 / real(lonmx, kind=r8)
738 : ! ****
739 : ! **** SET ARRAY YLATM (LATITUDE VALUES FOR GEOMAGNETIC GRID
740 : ! ****
741 0 : alatm(1) = -pi / 2._r8
742 0 : alat(1) = -90._r8
743 0 : alatm(jmxm) = pi / 2._r8
744 0 : alat(jmxm) = 90._r8
745 0 : do i = 2, ithmx
746 0 : alat(i) = alat(i-1)+dlatm*rtd
747 0 : alat(jmxm+1-i) = alat(jmxm+2-i)-dlatm*rtd
748 0 : alatm(i) = alatm(i-1)+dlatm
749 0 : alatm(jmxm+1-i) = alatm(jmxm+2-i)-dlatm
750 : end do
751 0 : alon(1) = -pi*rtd
752 0 : alonm(1) = -pi
753 0 : do i = 2, lonp1
754 0 : alon(i) = alon(i-1) + dlonm*rtd
755 0 : alonm(i) = alonm(i-1) + dlonm
756 : end do
757 :
758 : ! ylatm and ylonm are arrays of latitudes and longitudes of the
759 : ! distorted magnetic grids in radian - from consdyn.h
760 : ! Convert from apex magnetic grid to distorted magnetic grid
761 : !
762 : ! Allocate workspace for regrid routine rgrd_mod:
763 0 : lw = nmlonp1+nmlat+2*nmlonp1
764 : if (.not. allocated(w)) then
765 0 : allocate(w(lw), stat=ier)
766 0 : call check_alloc(ier, 'getamie', 'w', lw=lw)
767 : end if
768 0 : liw = nmlonp1 + nmlat
769 : if (.not. allocated(iw)) then
770 0 : allocate(iw(liw), stat=ier)
771 0 : call check_alloc(ier, 'getamie', 'iw', lw=liw)
772 : end if
773 0 : intpol(:) = 1 ! linear (not cubic) interp in both dimensions
774 0 : if (alatm(1) > ylatm(1)) then
775 0 : alatm(1) = ylatm(1)
776 : end if
777 0 : if (alatm(jmxm) < ylatm(nmlat)) then
778 0 : alatm(jmxm) = ylatm(nmlat)
779 : end if
780 0 : if (alonm(1) > ylonm(1)) then
781 0 : alonm(1) = ylonm(1)
782 : end if
783 0 : if (alonm(lonp1) < ylonm(nmlonp1)) then
784 0 : alonm(lonp1) = ylonm(nmlonp1)
785 : end if
786 :
787 : ! ylatm from -pi/2 to pi/2, and ylonm from -pi to pi
788 : call rgrd2(lonp1, jmxm, alonm, alatm, potm, nmlonp1, nmlat, &
789 0 : ylonm, ylatm, phihm, intpol, w, lw, iw, liw, ier)
790 : call rgrd2(lonp1, jmxm, alonm, alatm, ekvm, nmlonp1, nmlat, &
791 0 : ylonm, ylatm, amie_kevm, intpol, w, lw, iw, liw, ier)
792 : call rgrd2(lonp1, jmxm, alonm, alatm, efxm, nmlonp1, nmlat, &
793 0 : ylonm, ylatm, amie_efxm, intpol, w, lw, iw, liw, ier)
794 :
795 0 : if (iprint > 0 .and. masterproc) then
796 0 : write(iulog, *) subname, ': Max, min amie_efxm = ', &
797 0 : maxval(amie_efxm), minval(amie_efxm)
798 0 : write(iulog, "(a,': AMIE data interpolated to date and time')") subname
799 0 : write(iulog,"(a,': iyear,imo,iday,iutsec = ',3i6,i10)") subname, &
800 0 : iyear, imo, iday, iutsec
801 : write(iulog,"(2a,i6,2F9.5,3I6,f10.4)") &
802 0 : subname, ': AMIE iset1 f1,f2,year,mon,day,ut = ', iset1, &
803 0 : f1, f2, year(iset1), month(iset1), day(iset1), amie_nh_ut(iset1)
804 0 : write(iulog,*) subname, ': max,min phihm= ', maxval(phihm), minval(phihm)
805 : end if
806 : end if active_task
807 :
808 0 : end subroutine getamie
809 :
810 : !-----------------------------------------------------------------------
811 0 : subroutine close_files
812 :
813 0 : deallocate( year,month,day )
814 0 : deallocate( cusplat_nh_input, cuspmlt_nh_input, hpi_nh_input, &
815 0 : pcp_nh_input, amie_nh_ut, &
816 0 : cusplat_sh_input, cuspmlt_sh_input, hpi_sh_input, &
817 0 : pcp_sh_input, amie_sh_ut )
818 :
819 0 : call cam_pio_closefile(ncid_nh)
820 0 : call cam_pio_closefile(ncid_sh)
821 :
822 :
823 0 : end subroutine close_files
824 : !-----------------------------------------------------------------------
825 0 : subroutine open_files()
826 :
827 0 : call rdamie_nh(amienh_files(file_ndx))
828 0 : call rdamie_sh(amiesh_files(file_ndx))
829 :
830 0 : end subroutine open_files
831 :
832 : end module amie_module
|
