Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Provids electron fluxes read from input data set
3 : !--------------------------------------------------------------------------------
4 : module mee_fluxes
5 : use shr_kind_mod, only : r8 => shr_kind_r8, cl=> shr_kind_cl
6 : use spmd_utils, only : masterproc
7 : use cam_logfile, only : iulog
8 : use cam_abortutils, only : endrun
9 : use input_data_utils, only : time_coordinate
10 : use pio, only : file_desc_t, var_desc_t, pio_get_var, pio_inq_varid
11 : use pio, only : PIO_NOWRITE, pio_inq_dimid, pio_inq_dimlen
12 : use infnan, only: isnan
13 :
14 : implicit none
15 :
16 : private
17 : public :: mee_fluxes_readnl
18 : public :: mee_fluxes_init
19 : public :: mee_fluxes_final
20 : public :: mee_fluxes_adv ! read and time interpolate fluxes
21 : public :: mee_fluxes_extract ! interpolate flux to column L-shell
22 : public :: mee_fluxes_active ! true when input flux file is specified
23 : public :: mee_fluxes_denergy ! energy bin widths
24 : public :: mee_fluxes_energy ! center of each energy bin
25 : public :: mee_fluxes_nenergy ! number of energy bins
26 :
27 : real(r8),protected, pointer :: mee_fluxes_denergy(:) => null()
28 : real(r8),protected, pointer :: mee_fluxes_energy(:) => null()
29 : integer, protected :: mee_fluxes_nenergy
30 : logical, protected :: mee_fluxes_active = .false.
31 :
32 : real(r8), allocatable :: lshell(:)
33 : real(r8), allocatable :: indata(:,:,:)
34 : real(r8), allocatable :: influx(:,:)
35 : logical , allocatable :: valflx(:,:)
36 :
37 : character(len=cl) :: mee_fluxes_filepath = 'NONE'
38 : logical :: mee_fluxes_fillin = .false.
39 :
40 : type(time_coordinate) :: time_coord
41 : integer :: nlshells
42 :
43 : type(file_desc_t) :: file_id
44 : type(var_desc_t) :: flux_var_id
45 :
46 : contains
47 :
48 : !-----------------------------------------------------------------------------
49 : ! read namelist options
50 : !-----------------------------------------------------------------------------
51 0 : subroutine mee_fluxes_readnl(nlfile)
52 :
53 : use namelist_utils, only: find_group_name
54 : use spmd_utils, only: mpicom, mpi_character, mpi_logical, masterprocid
55 :
56 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
57 :
58 : ! Local variables
59 : integer :: unitn, ierr
60 : character(len=*), parameter :: subname = 'mee_fluxes_readnl'
61 :
62 : namelist /mee_fluxes_opts/ mee_fluxes_filepath, mee_fluxes_fillin
63 :
64 : ! Read namelist
65 0 : if (masterproc) then
66 0 : open( newunit=unitn, file=trim(nlfile), status='old' )
67 0 : call find_group_name(unitn, 'mee_fluxes_opts', status=ierr)
68 0 : if (ierr == 0) then
69 0 : read(unitn, mee_fluxes_opts, iostat=ierr)
70 0 : if (ierr /= 0) then
71 0 : call endrun(subname // ':: ERROR reading namelist')
72 : end if
73 : end if
74 0 : close(unitn)
75 : end if
76 :
77 : ! Broadcast namelist variables
78 0 : call mpi_bcast(mee_fluxes_filepath, len(mee_fluxes_filepath), mpi_character, masterprocid, mpicom, ierr)
79 0 : call mpi_bcast(mee_fluxes_fillin, 1, mpi_logical, masterprocid, mpicom, ierr)
80 :
81 0 : mee_fluxes_active = mee_fluxes_filepath /= 'NONE'
82 :
83 0 : if ( masterproc ) then
84 0 : if ( mee_fluxes_active ) then
85 0 : write(iulog,*) subname//':: Input electron fluxes filepath: '//trim(mee_fluxes_filepath)
86 0 : write(iulog,*) subname//':: Fill in missing fluxes with vdk-derived fluxes: ', mee_fluxes_fillin
87 : else
88 0 : write(iulog,*) subname//':: Electron fluxes are not prescribed'
89 : end if
90 : end if
91 :
92 0 : end subroutine mee_fluxes_readnl
93 :
94 : !-----------------------------------------------------------------------------
95 : ! intialize -- allocate memory and read coordinate data
96 : !-----------------------------------------------------------------------------
97 0 : subroutine mee_fluxes_init()
98 : use cam_pio_utils, only : cam_pio_openfile
99 : use ioFileMod, only : getfil
100 :
101 : character(len=cl) :: filen
102 : integer :: ierr, dimid, varid
103 0 : real(r8), allocatable :: logdelta(:)
104 : character(len=*), parameter :: subname = 'mee_fluxes_init: '
105 :
106 0 : if (.not.mee_fluxes_active) return
107 :
108 0 : call time_coord%initialize( mee_fluxes_filepath, force_time_interp=.true. )
109 :
110 0 : call getfil( mee_fluxes_filepath, filen, 0 )
111 0 : call cam_pio_openfile( file_id, filen, PIO_NOWRITE )
112 :
113 0 : ierr = pio_inq_dimid(file_id, 'energy', dimid)
114 0 : ierr = pio_inq_dimlen(file_id, dimid, mee_fluxes_nenergy )
115 :
116 0 : ierr = pio_inq_dimid(file_id, 'lshell', dimid)
117 0 : ierr = pio_inq_dimlen(file_id, dimid, nlshells )
118 :
119 0 : ierr = pio_inq_varid(file_id, 'RBSP_flux_scaled', flux_var_id)
120 :
121 0 : allocate( indata( mee_fluxes_nenergy, nlshells, 2 ), stat=ierr)
122 0 : if (ierr/=0) call endrun(subname//'not able to allocate indata')
123 0 : allocate( influx( mee_fluxes_nenergy, nlshells ), stat=ierr )
124 0 : if (ierr/=0) call endrun(subname//'not able to allocate influx')
125 0 : allocate( valflx( mee_fluxes_nenergy, nlshells ), stat=ierr )
126 0 : if (ierr/=0) call endrun(subname//'not able to allocate valflx')
127 0 : allocate( mee_fluxes_energy( mee_fluxes_nenergy ), stat=ierr )
128 0 : if (ierr/=0) call endrun(subname//'not able to allocate mee_fluxes_energy')
129 0 : allocate( mee_fluxes_denergy( mee_fluxes_nenergy ), stat=ierr )
130 0 : if (ierr/=0) call endrun(subname//'not able to allocate mee_fluxes_denergy')
131 0 : allocate( logdelta( mee_fluxes_nenergy ), stat=ierr )
132 0 : if (ierr/=0) call endrun(subname//'not able to allocate logdelta')
133 0 : allocate( lshell( nlshells ), stat=ierr )
134 0 : if (ierr/=0) call endrun(subname//'not able to allocate lshell')
135 :
136 :
137 0 : ierr = pio_inq_varid(file_id, 'energy', varid)
138 0 : ierr = pio_get_var( file_id, varid, mee_fluxes_energy)
139 :
140 0 : ierr = pio_inq_varid(file_id, 'lshell', varid)
141 0 : ierr = pio_get_var( file_id, varid, lshell)
142 :
143 0 : logdelta(2:) = log(mee_fluxes_energy(2:mee_fluxes_nenergy))-log(mee_fluxes_energy(1:mee_fluxes_nenergy-1))
144 0 : logdelta(1) = logdelta(2)
145 0 : mee_fluxes_denergy(:) = exp( log(mee_fluxes_energy(:)) + 0.5_r8*logdelta(:) ) &
146 0 : - exp( log(mee_fluxes_energy(:)) - 0.5_r8*logdelta(:) )
147 :
148 0 : deallocate(logdelta)
149 :
150 0 : call read_fluxes()
151 :
152 0 : end subroutine mee_fluxes_init
153 :
154 : !-----------------------------------------------------------------------------
155 : !-----------------------------------------------------------------------------
156 0 : subroutine mee_fluxes_final()
157 0 : use cam_pio_utils, only : cam_pio_closefile
158 :
159 0 : if (.not.mee_fluxes_active) return
160 :
161 0 : call cam_pio_closefile(file_id)
162 :
163 0 : deallocate(indata)
164 0 : deallocate(influx)
165 0 : deallocate(valflx)
166 0 : deallocate(lshell)
167 :
168 0 : deallocate(mee_fluxes_energy)
169 0 : deallocate(mee_fluxes_denergy)
170 :
171 0 : nullify(mee_fluxes_energy)
172 : nullify(mee_fluxes_denergy)
173 :
174 0 : end subroutine mee_fluxes_final
175 :
176 : !-----------------------------------------------------------------------------
177 : ! time interpolate the input fluxes
178 : ! reads data as needed
179 : !-----------------------------------------------------------------------------
180 0 : subroutine mee_fluxes_adv
181 :
182 0 : if (.not.mee_fluxes_active) return
183 :
184 0 : call time_coord%advance()
185 :
186 0 : if ( time_coord%read_more() ) then
187 0 : call read_fluxes( )
188 : endif
189 :
190 0 : influx(:,:) = 0._r8
191 :
192 0 : valflx(:,:) = (.not.isnan(indata(:,:,1))) .and. (.not.isnan(indata(:,:,2)))
193 0 : where ( valflx(:,:) )
194 0 : influx(:,:) = time_coord%wghts(1)*indata(:,:,1) + time_coord%wghts(2)*indata(:,:,2)
195 : end where
196 :
197 0 : if (any(isnan(influx))) then
198 0 : call endrun('mee_fluxes_adv -- influx has NaNs')
199 : end if
200 :
201 0 : end subroutine mee_fluxes_adv
202 :
203 : !-----------------------------------------------------------------------------
204 : ! linear interpolate fluxes in L-shell where the fluxes are valid
205 : !-----------------------------------------------------------------------------
206 0 : subroutine mee_fluxes_extract( l_shell, fluxes, valid )
207 :
208 : real(r8), intent(in) :: l_shell
209 : real(r8), intent(out) :: fluxes(mee_fluxes_nenergy)
210 : logical, intent(out) :: valid(mee_fluxes_nenergy)
211 :
212 : integer :: i, ndx1, ndx2
213 : logical :: found
214 : real(r8) :: wght1,wght2
215 :
216 0 : valid(:) = .not. mee_fluxes_fillin
217 0 : fluxes(:) = 0._r8
218 :
219 0 : if (.not.mee_fluxes_active) return
220 :
221 0 : found = .false.
222 :
223 0 : findloop: do i = 1,nlshells-1
224 0 : if ( l_shell>=lshell(i) .and. l_shell<=lshell(i+1) ) then
225 0 : ndx1=i
226 0 : ndx2=i+1
227 0 : wght2 = (l_shell-lshell(ndx1))/(lshell(ndx2)-lshell(ndx1))
228 0 : wght1 = 1._r8 - wght2
229 0 : found = .true.
230 0 : exit findloop
231 : endif
232 : end do findloop
233 :
234 0 : if (found) then
235 0 : if (mee_fluxes_fillin) then
236 0 : valid(:) = valflx(:,ndx1) .and. valflx(:,ndx2)
237 : end if
238 0 : where( valid(:) )
239 0 : fluxes(:) = wght1*influx(:,ndx1) + wght2*influx(:,ndx2)
240 : end where
241 : end if
242 :
243 : end subroutine mee_fluxes_extract
244 :
245 : !-----------------------------------------------------------------------------
246 : !-----------------------------------------------------------------------------
247 0 : subroutine read_fluxes()
248 :
249 : ! local vars
250 : integer :: ierr, cnt(4), start(4)
251 :
252 0 : cnt = (/1, mee_fluxes_nenergy, nlshells, 2/)
253 :
254 : ! use the 50 percentile level data ( index 3 )
255 0 : start = (/3, 1, 1, time_coord%indxs(1)/)
256 :
257 : ! float RBSP_flux_scaled(time, lshell, energy, percentiles) ;
258 0 : ierr = pio_get_var( file_id, flux_var_id, start, cnt, indata )
259 :
260 0 : end subroutine read_fluxes
261 :
262 : end module mee_fluxes
|
