Line data Source code
1 : module ltr_module
2 : !
3 : ! Module used to read data from the LFM/LTR 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
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_ltr
26 : public :: getltr
27 :
28 : ! Grid dimension sizes for LTR input data file:
29 : integer :: lonp1,latp1
30 : !
31 : ! Define fields for LTR input data file:
32 : ! electric potential in Volt
33 : ! mean energy in KeV
34 : ! energy flux in W/m^2
35 : ! Time interpolated LTR outputs with suffix _ltr
36 : !
37 : real(r8),allocatable,dimension(:,:,:) :: & ! (lonp1,latp1,ntimes)
38 : pot_input, ekv_input, efx_input
39 : real(r8),allocatable,dimension(:,:) :: & ! (lonp1,latp1)
40 : pot_ltr, ekv_ltr, efx_ltr
41 : integer, allocatable,dimension(:) :: & ! (ntimes)
42 : year,month,day,jday
43 : real(r8), allocatable,dimension(:) :: & ! (ntimes)
44 : hpi_input, pcp_input, ltr_ut
45 : real(r8) :: hpi_ltr, pcp_ltr
46 : !
47 : type(file_desc_t) :: ncid
48 :
49 : character(len=cl), allocatable :: ltr_files(:)
50 : integer :: num_files, file_ndx
51 :
52 : type(time_coordinate) :: time_coord
53 :
54 : contains
55 :
56 : !-----------------------------------------------------------------------
57 0 : subroutine init_ltr(ltr_list)
58 :
59 : character(len=*),intent(in) :: ltr_list(:)
60 :
61 : integer :: n, nfiles
62 :
63 0 : nfiles = size(ltr_list)
64 0 : num_files = 0
65 :
66 0 : count_files: do n = 1,nfiles
67 0 : if (len_trim(ltr_list(n))<1 .or. &
68 0 : trim(ltr_list(n))=='NONE') then
69 : exit count_files
70 : else
71 0 : num_files = num_files + 1
72 : end if
73 : end do count_files
74 :
75 0 : allocate(ltr_files(num_files))
76 0 : ltr_files(:num_files) = ltr_list(:num_files)
77 0 : file_ndx = 1
78 0 : call open_files()
79 :
80 0 : end subroutine init_ltr
81 :
82 : !-----------------------------------------------------------------------
83 0 : subroutine rdltr(ltrfile)
84 : !
85 : ! Read LTR data
86 : !
87 : character(len=*), intent(in) :: ltrfile
88 : ! Local:
89 : integer :: istat, ntimes, ndims, nvars, ngatts, ier
90 : integer :: idunlim
91 : integer :: id_lon, id_lat, id_time
92 : integer :: idv_year, idv_mon, idv_day, idv_jday
93 : integer :: idv_ut, idv_hpi, idv_pcp
94 : character(len=*), parameter :: subname = 'rdltr'
95 :
96 : !
97 : !
98 0 : if (masterproc) then
99 0 : write(iulog, "(/, 72('-'))")
100 0 : write(iulog, "(a, ': read LTR data:')") subname
101 : end if
102 : !
103 : ! Open netcdf file:
104 0 : call cam_pio_openfile(ncid, ltrfile, pio_nowrite)
105 : !
106 : ! Get LTR grid dimension:
107 0 : istat = pio_inq_dimid(ncid, 'lon', id_lon)
108 0 : istat = pio_inquire_dimension(ncid, id_lon, len=lonp1)
109 0 : call check_ncerr(istat, subname, 'LTR longitude dimension')
110 :
111 0 : istat = pio_inq_dimid(ncid, 'lat', id_lat)
112 0 : istat = pio_inquire_dimension(ncid, id_lat, len=latp1)
113 0 : call check_ncerr(istat, subname, 'LTR latitude dimension')
114 :
115 0 : call time_coord%initialize( ltrfile, set_weights=.false. )
116 : !
117 : ! Get time dimension:
118 0 : istat = pio_inquire(ncid, unlimiteddimid=id_time)
119 0 : istat = pio_inquire_dimension(ncid, id_time, len=ntimes)
120 0 : call check_ncerr(istat, subname, 'LTR time dimension')
121 : !
122 : ! Search for requested LTR output fields
123 0 : istat = pio_inquire(ncid,ndims,nvars,ngatts,idunlim)
124 : !
125 : ! Get 1-D LTR fields (ntimes)
126 0 : if (.not. allocated(year)) then
127 0 : allocate(year(ntimes), stat=ier)
128 0 : call check_alloc(ier, subname, 'year', ntimes=ntimes)
129 : end if
130 0 : istat = pio_inq_varid(ncid, 'year', idv_year)
131 0 : call check_ncerr(istat, subname, 'LTR year id')
132 0 : istat = pio_get_var(ncid, idv_year, year)
133 0 : call check_ncerr(istat, subname, 'LTR year')
134 :
135 0 : if (.not. allocated(month)) then
136 0 : allocate(month(ntimes), stat=ier)
137 0 : call check_alloc(ier, subname, 'month', ntimes=ntimes)
138 : end if
139 0 : istat = pio_inq_varid(ncid, 'month', idv_mon)
140 0 : call check_ncerr(istat, subname, 'LTR month id')
141 0 : istat = pio_get_var(ncid, idv_mon, month)
142 0 : call check_ncerr(istat, subname, 'LTR month')
143 0 : if (.not. allocated(day)) then
144 0 : allocate(day(ntimes), stat=ier)
145 0 : call check_alloc(ier, subname, 'day', ntimes=ntimes)
146 : end if
147 0 : istat = pio_inq_varid(ncid, 'day', idv_day)
148 0 : call check_ncerr(istat, subname, 'LTR day id')
149 0 : istat = pio_get_var(ncid, idv_day, day)
150 0 : call check_ncerr(istat, subname, 'LTR day')
151 :
152 0 : if (.not. allocated(jday)) then
153 0 : allocate(jday(ntimes), stat=ier)
154 0 : call check_alloc(ier, subname, 'jday', ntimes=ntimes)
155 : end if
156 0 : istat = pio_inq_varid(ncid, 'jday', idv_jday)
157 0 : call check_ncerr(istat, subname, 'LTR jday id')
158 0 : istat = pio_get_var(ncid, idv_jday, jday)
159 0 : call check_ncerr(istat, subname, 'LTR jday')
160 : !
161 : ! Allocate 1-d fields:
162 0 : if (.not. allocated(ltr_ut)) then
163 0 : allocate(ltr_ut(ntimes), stat=ier)
164 0 : call check_alloc(ier, subname, 'ltr_ut', ntimes=ntimes)
165 : end if
166 0 : if (.not. allocated(hpi_input)) then
167 0 : allocate(hpi_input(ntimes), stat=ier)
168 0 : call check_alloc(ier, subname, 'hpi_input', ntimes=ntimes)
169 : end if
170 0 : if (.not. allocated(pcp_input)) then
171 0 : allocate(pcp_input(ntimes), stat=ier)
172 0 : call check_alloc(ier, subname, 'pcp_input', ntimes=ntimes)
173 : end if
174 : !
175 : ! Get ut
176 0 : istat = pio_inq_varid(ncid, 'ut', idv_ut)
177 0 : call check_ncerr(istat, subname, 'LTR ut id')
178 0 : istat = pio_get_var(ncid, idv_ut, ltr_ut)
179 0 : call check_ncerr(istat, subname, 'LTR ut')
180 : !
181 : ! Get HPI
182 0 : istat = pio_inq_varid(ncid, 'hpiN', idv_hpi)
183 0 : call check_ncerr(istat, subname, 'LTR hpi id')
184 0 : istat = pio_get_var(ncid, idv_hpi, hpi_input)
185 0 : call check_ncerr(istat, subname, 'LTR hpi')
186 : !
187 : ! Get PCP
188 0 : istat = pio_inq_varid(ncid, 'pcpN', idv_pcp)
189 0 : call check_ncerr(istat, subname, 'LTR pcp id')
190 0 : istat = pio_get_var(ncid, idv_pcp, pcp_input)
191 0 : call check_ncerr(istat, subname, 'LTR pcp')
192 : !
193 : ! Allocate 2-d fields:
194 0 : if (.not. allocated(pot_ltr)) then
195 0 : allocate(pot_ltr(lonp1, latp1), stat=ier)
196 0 : call check_alloc(ier, subname, 'pot_ltr', lonp1=lonp1, latp1=latp1)
197 : end if
198 0 : if (.not. allocated(ekv_ltr)) then
199 0 : allocate(ekv_ltr(lonp1, latp1), stat=ier)
200 0 : call check_alloc(ier, subname, 'ekv_ltr', lonp1=lonp1, latp1=latp1)
201 : end if
202 0 : if (.not. allocated(efx_ltr)) then
203 0 : allocate(efx_ltr(lonp1, latp1), stat=ier)
204 0 : call check_alloc(ier, subname, 'efx_ltr', lonp1=lonp1, latp1=latp1)
205 : end if
206 : !
207 : ! Allocate 3-d fields:
208 0 : if (.not. allocated(pot_input)) then
209 0 : allocate(pot_input(lonp1, latp1, 2), stat=ier)
210 : call check_alloc(ier, subname, 'pot_input', &
211 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
212 : end if
213 0 : if (.not. allocated(ekv_input)) then
214 0 : allocate(ekv_input(lonp1, latp1, 2), stat=ier)
215 : call check_alloc(ier, subname, 'ekv_input', &
216 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
217 : end if
218 0 : if (.not. allocated(efx_input)) then
219 0 : allocate(efx_input(lonp1, latp1, 2), stat=ier)
220 : call check_alloc(ier, subname, 'efx_input', &
221 0 : lonp1=lonp1, latp1=latp1, ntimes=ntimes)
222 : end if
223 0 : end subroutine rdltr
224 :
225 : !-----------------------------------------------------------------------
226 0 : subroutine update_3d_fields( ncid, offset, kount, pot_3d,ekv_3d,efx_3d )
227 :
228 : type(file_desc_t), intent(in) :: ncid
229 : integer, intent(in) :: offset(:)
230 : integer, intent(in) :: kount(:)
231 : real(r8),intent(out) :: pot_3d(:,:,:)
232 : real(r8),intent(out) :: ekv_3d(:,:,:)
233 : real(r8),intent(out) :: efx_3d(:,:,:)
234 :
235 :
236 : integer :: istat
237 : integer :: idv_pot,idv_ekv, idv_efx
238 : character(len=*), parameter :: subname = 'update_3d_fields'
239 :
240 : !
241 : ! Get 3-D fields (lon,lat,ntimes)
242 : !
243 : ! electric potential
244 0 : istat = pio_inq_varid(ncid, 'pot', idv_pot)
245 0 : call check_ncerr(istat, subname, 'LTR pot id')
246 0 : istat = pio_get_var(ncid, idv_pot, offset, kount, pot_3d)
247 0 : call check_ncerr(istat, subname, 'LTR pot')
248 : !
249 : ! mean energy
250 0 : istat = pio_inq_varid(ncid, 'ekv', idv_ekv)
251 0 : call check_ncerr(istat, subname, 'LTR ekv id')
252 0 : istat = pio_get_var(ncid, idv_ekv, offset, kount, ekv_3d)
253 0 : call check_ncerr(istat, subname, 'LTR ekv')
254 : !
255 : ! energy flux
256 0 : istat = pio_inq_varid(ncid, 'efx', idv_efx)
257 0 : call check_ncerr(istat, subname, 'LTR efx id')
258 0 : istat = pio_get_var(ncid, idv_efx, offset, kount, efx_3d)
259 0 : call check_ncerr(istat, subname, 'LTR efx')
260 :
261 0 : end subroutine update_3d_fields
262 :
263 : !-----------------------------------------------------------------------
264 0 : subroutine getltr(iyear, imo, iday, iutsec, sunlon, iprint, &
265 0 : iltr, phihm, ltr_efxm, ltr_kevm)
266 : use cam_history_support, only: fillvalue
267 : use rgrd_mod, only: rgrd2
268 : !
269 : ! Read LTR outputs from ltr_ncfile file, returning electric potential,
270 : ! auroral mean energy and energy flux at current date and time,
271 : ! and the data is linearly interpolated to the model time
272 : !
273 : !
274 : ! Args:
275 :
276 : integer, intent(in) :: iyear
277 : integer, intent(in) :: imo
278 : integer, intent(in) :: iday
279 : real(r8), intent(in) :: sunlon
280 : integer, intent(in) :: iutsec
281 : integer, intent(in) :: iprint
282 : integer, intent(out) :: iltr
283 : real(r8), intent(out) :: phihm(nmlonp1,nmlat)
284 : real(r8), intent(out) :: ltr_efxm(nmlonp1,nmlat) ! on geomag grid
285 : real(r8), intent(out) :: ltr_kevm(nmlonp1,nmlat) ! on geomag grid
286 : !
287 : !
288 : ! Local:
289 0 : real(r8) :: potm(lonp1,latp1)
290 0 : real(r8) :: efxm(lonp1,latp1), ekvm(lonp1,latp1)
291 0 : real(r8) :: alat(latp1), alon(lonp1)
292 0 : real(r8) :: alatm(latp1), alonm(lonp1)
293 : integer :: ier, lw, liw, intpol(2)
294 0 : integer, allocatable :: iw(:)
295 0 : real(r8), allocatable :: w(:)
296 : integer :: i, j, ithmx
297 : integer :: nn, iset2, iset1, m, mp1, n
298 : real(r8) :: f1, f2
299 : real(r8) :: del, xmlt, dmlat, dlatm, dlonm, dmltm, rot
300 : integer :: offset(3), kount(3)
301 : character(len=*), parameter :: subname = 'getltr'
302 :
303 0 : phihm = fillvalue
304 0 : ltr_efxm = fillvalue
305 0 : ltr_kevm = fillvalue
306 :
307 0 : if (iprint > 0 .and. masterproc) then
308 0 : write(iulog,"(/,72('-'))")
309 0 : write(iulog,"(a,':')") subname
310 : write(iulog,"(a,i4,', iday = ',i3,', iutsec = ',i10)") &
311 0 : 'Initial requested iyear= ', iyear, iday, iutsec
312 : end if
313 :
314 0 : nn = size(ltr_ut)
315 :
316 : !
317 : ! Check times:
318 : !
319 0 : iltr = 1 - time_coord%times_check()
320 :
321 0 : check_loop: do while( iltr/=1 )
322 :
323 0 : if (masterproc) write(iulog,*) 'file_ndx = ',file_ndx
324 :
325 0 : if (iltr==2) then
326 0 : if (masterproc) then
327 : write(iulog, "(a,': Model date prior to LTR first date:',3I5)") &
328 0 : subname, year(1), month(1), day(1)
329 : end if
330 0 : return
331 : endif
332 :
333 0 : if (iltr==0) then
334 0 : if (masterproc) then
335 : write(iulog, "(a,': Model date beyond the LTR last Data:',3I5)") &
336 0 : subname, year(nn), month(nn), day(nn)
337 : end if
338 :
339 0 : if (file_ndx<num_files) then
340 0 : file_ndx = file_ndx+1
341 0 : call close_files()
342 0 : call open_files()
343 0 : iltr = 1 - time_coord%times_check()
344 : else
345 : return
346 : end if
347 :
348 : endif
349 : end do check_loop
350 :
351 : !
352 : ! get LTR data
353 0 : pot_ltr(:,:) = 0._r8
354 0 : ekv_ltr(:,:) = 0._r8
355 0 : efx_ltr(:,:) = 0._r8
356 0 : hpi_ltr = 0._r8
357 0 : pcp_ltr = 0._r8
358 :
359 0 : call time_coord%advance()
360 :
361 0 : iset1 = time_coord%indxs(1)
362 0 : iset2 = time_coord%indxs(2)
363 :
364 0 : f1 = time_coord%wghts(1)
365 0 : f2 = time_coord%wghts(2)
366 :
367 0 : hpi_ltr = (f1*hpi_input(iset1) + f2*hpi_input(iset2))
368 0 : pcp_ltr = (f1*pcp_input(iset1) + f2*pcp_input(iset2))
369 :
370 0 : offset = (/1,1,iset1/)
371 0 : kount = (/lonp1,latp1,2/)
372 0 : call update_3d_fields( ncid, offset, kount, pot_input,ekv_input,efx_input )
373 0 : pot_ltr(:,:) = (f1*pot_input(:,:,1) + f2*pot_input(:,:,2))
374 0 : ekv_ltr(:,:) = (f1*ekv_input(:,:,1) + f2*ekv_input(:,:,2))
375 0 : efx_ltr(:,:) = (f1*efx_input(:,:,1) + f2*efx_input(:,:,2))
376 :
377 0 : active_task: if ( mytid<ntask ) then
378 :
379 : ! mlongitude starts from 180 degree
380 0 : rot = sunlon*rtd
381 0 : if(rot < 0) then
382 0 : rot = rot + 360._r8 ! 0 to 360 degrees
383 : end if
384 0 : rot = rot / 15._r8 ! convert from degree to hrs
385 :
386 0 : dmltm = 24._r8 / real(lonp1, kind=r8)
387 :
388 0 : do i = 1, lonp1
389 0 : xmlt = (real(i-1, kind=r8) * dmltm) - rot + 24._r8
390 0 : xmlt = MOD(xmlt, 24._r8)
391 0 : m = int(xmlt/dmltm + 1.001_r8)
392 0 : mp1 = m + 1
393 0 : if (mp1 > lonp1) mp1 = 2
394 0 : del = xmlt - (m-1)*dmltm
395 : ! Put in LTR arrays from south pole to north pole
396 0 : do j=1,latp1
397 0 : potm(i,j) = (1._r8-del)*pot_ltr(m,j) + &
398 0 : del*pot_ltr(mp1,j)
399 0 : ekvm(i,j) = (1._r8-del)*ekv_ltr(m,j) + &
400 0 : del*ekv_ltr(mp1,j)
401 0 : if (ekvm(i,j) == 0._r8) ekvm(i,j)=1._r8
402 0 : efxm(i,j) = (1._r8-del)*efx_ltr(m,j) + &
403 0 : del*efx_ltr(mp1,j)
404 : end do
405 :
406 : end do
407 :
408 : ! Set up coeffs to go between EPOTM(IMXMP,JMNH) and TIEPOT(IMAXM,JMAXMH)
409 :
410 : ! **** SET GRID SPACING DLATM, DLONG, DLONM
411 : ! DMLAT=lat spacing in degrees of LTR apex grid
412 0 : dmlat = 180._r8 / real(latp1-1, kind=r8)
413 0 : dlatm = dmlat * dtr
414 0 : dlonm = 2._r8 * pi / real(lonp1, kind=r8)
415 0 : dmltm = 24._r8 / real(lonp1, kind=r8)
416 : ! ****
417 : ! **** SET ARRAY YLATM (LATITUDE VALUES FOR GEOMAGNETIC GRID
418 : ! ****
419 0 : alatm(1) = -pi / 2._r8
420 0 : alat(1) = -90._r8
421 0 : alatm(latp1) = pi / 2._r8
422 0 : alat(latp1) = 90._r8
423 0 : ithmx = (latp1+1)/2
424 0 : do i = 2, ithmx
425 0 : alat(i) = alat(i-1)+dlatm*rtd
426 0 : alat(latp1+1-i) = alat(latp1+2-i)-dlatm*rtd
427 0 : alatm(i) = alatm(i-1)+dlatm
428 0 : alatm(latp1+1-i) = alatm(latp1+2-i)-dlatm
429 : end do
430 0 : alon(1) = -pi*rtd
431 0 : alonm(1) = -pi
432 0 : do i = 2, lonp1
433 0 : alon(i) = alon(i-1) + dlonm*rtd
434 0 : alonm(i) = alonm(i-1) + dlonm
435 : end do
436 :
437 : ! ylatm and ylonm are arrays of latitudes and longitudes of the
438 : ! distorted magnetic grids in radian - from consdyn.h
439 : ! Convert from apex magnetic grid to distorted magnetic grid
440 : !
441 : ! Allocate workspace for regrid routine rgrd_mod:
442 0 : lw = nmlonp1+nmlat+2*nmlonp1
443 : if (.not. allocated(w)) then
444 0 : allocate(w(lw), stat=ier)
445 0 : call check_alloc(ier, 'getltr', 'w', lw=lw)
446 : end if
447 0 : liw = nmlonp1 + nmlat
448 : if (.not. allocated(iw)) then
449 0 : allocate(iw(liw), stat=ier)
450 0 : call check_alloc(ier, 'getltr', 'iw', lw=liw)
451 : end if
452 0 : intpol(:) = 1 ! linear (not cubic) interp in both dimensions
453 0 : if (alatm(1) > ylatm(1)) then
454 0 : alatm(1) = ylatm(1)
455 : end if
456 0 : if (alatm(latp1) < ylatm(nmlat)) then
457 0 : alatm(latp1) = ylatm(nmlat)
458 : end if
459 0 : if (alonm(1) > ylonm(1)) then
460 0 : alonm(1) = ylonm(1)
461 : end if
462 0 : if (alonm(lonp1) < ylonm(nmlonp1)) then
463 0 : alonm(lonp1) = ylonm(nmlonp1)
464 : end if
465 :
466 : ! ylatm from -pi/2 to pi/2, and ylonm from -pi to pi
467 : call rgrd2(lonp1, latp1, alonm, alatm, potm, nmlonp1, nmlat, &
468 0 : ylonm, ylatm, phihm, intpol, w, lw, iw, liw, ier)
469 : call rgrd2(lonp1, latp1, alonm, alatm, ekvm, nmlonp1, nmlat, &
470 0 : ylonm, ylatm, ltr_kevm, intpol, w, lw, iw, liw, ier)
471 : call rgrd2(lonp1, latp1, alonm, alatm, efxm, nmlonp1, nmlat, &
472 0 : ylonm, ylatm, ltr_efxm, intpol, w, lw, iw, liw, ier)
473 :
474 0 : if (iprint > 0 .and. masterproc) then
475 0 : write(iulog, *) subname, ': Max, min ltr_efxm = ', &
476 0 : maxval(ltr_efxm), minval(ltr_efxm)
477 0 : write(iulog, "('getltr: LTR data interpolated to date and time')")
478 : write(iulog,"('getltr: iyear,imo,iday,iutsec = ',3i6,i10)") &
479 0 : iyear,imo,iday,iutsec
480 : write(iulog,"('getltr: LTR iset1 f1,f2,year,mon,day,ut = ', &
481 : i6,2F9.5,3I6,f10.4)") &
482 0 : iset1,f1,f2,year(iset1),month(iset1),day(iset1),ltr_ut(iset1)
483 0 : write(iulog,*)'getltr: max,min phihm= ', maxval(phihm),minval(phihm)
484 : end if
485 :
486 : end if active_task
487 :
488 0 : end subroutine getltr
489 : !-------------------------------------------------------------------
490 :
491 0 : subroutine close_files
492 :
493 0 : deallocate( year,month,day )
494 0 : deallocate( hpi_input, pcp_input, ltr_ut )
495 :
496 0 : call cam_pio_closefile(ncid)
497 :
498 0 : end subroutine close_files
499 : !-----------------------------------------------------------------------
500 0 : subroutine open_files()
501 :
502 0 : call rdltr(ltr_files(file_ndx))
503 :
504 0 : end subroutine open_files
505 :
506 : end module ltr_module
|
