Line data Source code
1 : module mo_snoe
2 : !----------------------------------------------------------------------
3 : ! An empirical model of nitric oxide (NO) in the lower thermosphere
4 : ! (100 - 150 km altitude), based on measurements from the Student
5 : ! Nitric Oxide Explorer (SNOE). Model uses empirical orthogonal functions
6 : ! (EOFs) derived from the SNOE dataset to describe spatial variability
7 : ! in NO. Model NO is the sum of a mean distribution and EOFs multiplied
8 : ! by coefficients based on geophysical parameters. These geophysical
9 : ! parameters are day of year, Kp magnetic index and F10.7 solar uv index.
10 : !
11 : ! Model is utilized by calling subroutine snoe_3d(), which returns
12 : ! a 3-D distribution of NO on geographic coordinates (lon/lat are
13 : ! supplied by user). Altitude is fixed to SNOE grid (every 3.33 km).
14 : !
15 : !
16 : ! DRM 4/03
17 : !----------------------------------------------------------------------
18 :
19 : use shr_kind_mod, only : r8 => shr_kind_r8
20 : use ppgrid, only : pcols
21 : use physconst, only : pi
22 : use mo_constants, only : r2d, d2r
23 : use cam_abortutils, only : endrun
24 : use cam_logfile, only : iulog
25 :
26 : implicit none
27 :
28 : private
29 :
30 : public :: snoe_inti, snoe_timestep_init, set_no_ubc
31 : public :: ndx_no
32 :
33 : save
34 :
35 : integer, parameter :: nmodes = 3
36 : real(r8), parameter :: delz = 104.6_r8 ! delta z (km)
37 : real(r8), parameter :: re = 6471._r8 ! radial distance from ED at 100 km altitude (km)
38 : real(r8), parameter :: delx = -399.1_r8 ! delta x (km)
39 : real(r8), parameter :: dely = -286.1_r8 ! delta y (km)
40 : real(r8), parameter :: twopi = 2._r8*pi, pid2=0.5_r8*pi
41 : real(r8), parameter :: theta_n = d2r*(90._r8 - 78.98_r8) ! co-latitude of CD North Pole (radians)
42 : real(r8), parameter :: phi_n = d2r*289.1_r8 ! longitude of CD North Pole (radians)
43 :
44 : integer :: nmlat
45 : integer :: nlev
46 : integer :: ndx_no
47 :
48 :
49 : !----------------------------------------------------------------------
50 : ! ... snoe mean and eof data
51 : !----------------------------------------------------------------------
52 : real(r8), allocatable :: mlat(:,:) ! magnetic latitudes corresponding to geo latitudes (radians)
53 : real(r8), allocatable :: mlat2d(:) ! snoe latitudes (degrees)
54 : real(r8), allocatable :: lev(:) ! snoe levels (km)
55 : real(r8), allocatable :: no_mean(:,:) ! mean no
56 : real(r8), allocatable :: eofs(:,:,:) ! empirical orthogonal ftns
57 : real(r8), allocatable :: snoe_no(:,:,:) ! snoe no interpolated to model lons (molecules/cm^3)
58 :
59 : real(r8) :: cthetan
60 : real(r8) :: sthetan
61 :
62 : contains
63 :
64 0 : subroutine snoe_inti(snoe_ubc_file)
65 : !----------------------------------------------------------------------
66 : ! ... initialize snoe
67 : !----------------------------------------------------------------------
68 :
69 : use ppgrid, only : begchunk, endchunk
70 : use constituents, only : cnst_get_ind, cnst_fixed_ubc
71 :
72 : !----------------------------------------------------------------------
73 : ! ... dummy arguments
74 : !----------------------------------------------------------------------
75 : character(len=*), intent(in) :: snoe_ubc_file
76 :
77 : !----------------------------------------------------------------------
78 : ! ... local variables
79 : !----------------------------------------------------------------------
80 : integer :: astat
81 :
82 : !----------------------------------------------------------------------
83 : ! ... do we have no with a fixed ubc ?
84 : !----------------------------------------------------------------------
85 0 : call cnst_get_ind( 'NO', ndx_no, abort=.false. )
86 0 : if( ndx_no > 0 ) then
87 0 : if( .not. cnst_fixed_ubc(ndx_no) ) then
88 0 : return
89 : end if
90 : else
91 : return
92 : end if
93 :
94 :
95 :
96 0 : cthetan = cos( theta_n )
97 0 : sthetan = sin( theta_n )
98 : !----------------------------------------------------------------------
99 : ! ... read snoe netcdf file
100 : !----------------------------------------------------------------------
101 0 : call snoe_rdeof(snoe_ubc_file)
102 : allocate( snoe_no(pcols,nlev,begchunk:endchunk), &
103 0 : mlat(pcols,begchunk:endchunk),stat=astat )
104 0 : if( astat /= 0 ) then
105 0 : write(iulog,*) 'snoe_inti: failed to allocate snoe_no,mlat; error = ',astat
106 0 : call endrun
107 : end if
108 :
109 : !----------------------------------------------------------------------
110 : ! ... get lon/lat transformed to magnetic coords
111 : !----------------------------------------------------------------------
112 0 : call geo2mag
113 :
114 : end subroutine snoe_inti
115 :
116 0 : subroutine snoe_rdeof(snoe_ubc_file)
117 : !----------------------------------------------------------------------
118 : ! ... read in eofs from netcdf file
119 : !----------------------------------------------------------------------
120 :
121 : use ioFileMod, only : getfil
122 : use cam_pio_utils, only : cam_pio_openfile
123 : use pio, only : file_desc_t, pio_get_var, pio_closefile, &
124 : pio_nowrite, pio_inq_varid, pio_inq_dimid, pio_inq_dimlen
125 : implicit none
126 :
127 : !----------------------------------------------------------------------
128 : ! ... dummy arguments
129 : !----------------------------------------------------------------------
130 : character(len=*), intent(in) :: snoe_ubc_file
131 :
132 : !----------------------------------------------------------------------
133 : ! ... local variables
134 : !----------------------------------------------------------------------
135 : integer :: istat, ierr
136 : integer :: dim_id, var_id
137 : type(file_desc_t) :: ncid
138 : character(len=256) :: locfn
139 :
140 : !----------------------------------------------------------------------
141 :
142 : #ifdef SNOE_DIAGS
143 : write(iulog,*) 'snoe_rdeof: entered routine'
144 : #endif
145 :
146 : !----------------------------------------------------------------------
147 : ! ... open the netcdf file
148 : !----------------------------------------------------------------------
149 0 : call getfil(snoe_ubc_file, locfn, 0)
150 0 : call cam_pio_openfile( ncid, trim(locfn), PIO_NOWRITE)
151 :
152 : !----------------------------------------------------------------------
153 : ! ... read the snoe dimensions
154 : !----------------------------------------------------------------------
155 0 : ierr = pio_inq_dimid( ncid, 'lat', dim_id )
156 0 : ierr = pio_inq_dimlen( ncid, dim_id, nmlat )
157 0 : ierr = pio_inq_dimid( ncid, 'lev', dim_id )
158 0 : ierr = pio_inq_dimlen( ncid, dim_id, nlev )
159 :
160 : !----------------------------------------------------------------------
161 : ! ... allocate snoe variables
162 : !----------------------------------------------------------------------
163 : allocate( mlat2d(nmlat), lev(nlev), &
164 0 : no_mean(nmlat,nlev), eofs(nmlat,nlev,nmodes), stat=istat )
165 0 : if( istat /= 0 ) then
166 0 : write(iulog,*) 'snoe_rdeof: failed to allocate mlat2d ... eofs; error = ',istat
167 0 : call endrun
168 : end if
169 :
170 : !----------------------------------------------------------------------
171 : ! ... read the snoe variables
172 : !----------------------------------------------------------------------
173 0 : ierr = pio_inq_varid( ncid, 'NO', var_id )
174 0 : ierr = pio_get_var( ncid, var_id, no_mean )
175 :
176 0 : ierr = pio_inq_varid( ncid, 'EOF', var_id )
177 0 : ierr = pio_get_var( ncid, var_id, (/1,1,1/), (/nmlat, nlev, nmodes/), eofs )
178 :
179 0 : ierr = pio_inq_varid( ncid, 'lat', var_id )
180 0 : ierr = pio_get_var( ncid, var_id, mlat2d )
181 0 : mlat2d(:) = d2r*mlat2d(:)
182 :
183 0 : ierr = pio_inq_varid( ncid, 'z', var_id )
184 0 : ierr = pio_get_var( ncid, var_id, lev )
185 :
186 : !----------------------------------------------------------------------
187 : ! ... close the netcdf file
188 : !----------------------------------------------------------------------
189 0 : call pio_closefile( ncid)
190 :
191 0 : end subroutine snoe_rdeof
192 :
193 0 : subroutine geo2mag
194 : !----------------------------------------------------------------------
195 : ! ... converts geographic latitude and longitude pair
196 : ! to eccentric dipole (ED) geomagnetic coordinates
197 : !----------------------------------------------------------------------
198 :
199 0 : use ppgrid, only : begchunk, endchunk
200 : use phys_grid, only : get_ncols_p, get_rlon_all_p, get_rlat_all_p
201 :
202 : implicit none
203 :
204 : !----------------------------------------------------------------------
205 : ! ... local variables
206 : !----------------------------------------------------------------------
207 : integer :: c, i, j
208 : integer :: ncol
209 : real(r8) :: cdlat
210 : real(r8) :: wrk
211 : real(r8) :: smaglat, cmaglat, cmaglon, smaglon
212 : real(r8) :: edlon, sedlon, cedlon, edlat, cedlat
213 : real(r8) :: singlat, cosglat
214 : real(r8) :: rlons(pcols)
215 : real(r8) :: rlats(pcols)
216 :
217 : !----------------------------------------------------------------------
218 : ! ... start with centered dipole (CD) coordinates
219 : !----------------------------------------------------------------------
220 : chunk_loop : &
221 0 : do c = begchunk,endchunk
222 0 : ncol = get_ncols_p( c )
223 0 : call get_rlon_all_p( c, ncol, rlons )
224 0 : call get_rlat_all_p( c, ncol, rlats )
225 : column_loop : &
226 0 : do i = 1,ncol
227 0 : singlat = sin( rlats(i) )
228 0 : cosglat = cos( rlats(i) )
229 0 : wrk = rlons(i) - phi_n
230 0 : cdlat = acos( singlat * cthetan + cosglat * sthetan * cos(wrk) )
231 0 : smaglat = sin( cdlat )
232 0 : cmaglat = cos( cdlat )
233 0 : cmaglon = -(singlat - cthetan * cmaglat) / (sthetan * smaglat)
234 0 : smaglon = cosglat * sin(wrk) / smaglat
235 : !----------------------------------------------------------------------
236 : ! ... convert CD coords to ED coords- equation (39 a,b) of Fraser-Smith
237 : !----------------------------------------------------------------------
238 0 : sedlon = re * smaglat * smaglon - dely
239 0 : cedlon = re * smaglat * cmaglon - delx
240 0 : edlon = atan2( sedlon, cedlon ) ! ED magnetic long. (degrees)
241 0 : cedlat = (re * cmaglat - delz) * cos(edlon)
242 0 : edlat = atan2( cedlon, cedlat ) ! ED magnetic latitude (degrees)
243 : !----------------------------------------------------------------------
244 : ! ... convert co-latitudes into latitudes
245 : !----------------------------------------------------------------------
246 0 : if( edlat < 0._r8 ) then
247 0 : edlat = edlat + pi
248 : end if
249 0 : mlat(i,c) = pid2 - edlat
250 : end do column_loop
251 : end do chunk_loop
252 :
253 0 : end subroutine geo2mag
254 :
255 0 : subroutine snoe_timestep_init( kp, f107 )
256 : !----------------------------------------------------------------------
257 : ! ... driver routine that provides access to empirical model
258 : ! data returned is three dimensional (on geographic coordinates)
259 : !----------------------------------------------------------------------
260 :
261 0 : use time_manager, only : get_curr_calday
262 : use ppgrid, only : pcols, begchunk, endchunk
263 : use phys_grid, only : get_ncols_p
264 : use spmd_utils, only : masterproc
265 :
266 : implicit none
267 :
268 : !----------------------------------------------------------------------
269 : ! ... dummy arguments
270 : !----------------------------------------------------------------------
271 : real(r8), intent(in) :: kp ! solar activity index
272 : real(r8), intent(in) :: f107 ! solar activity index
273 :
274 : !----------------------------------------------------------------------
275 : ! ... local variables
276 : !----------------------------------------------------------------------
277 : integer :: astat, c
278 : integer :: ncol
279 : real(r8) :: doy ! day of year
280 0 : real(r8), allocatable :: zm(:,:) ! zonal mean nitric oxide distribution (molecules/cm^3)
281 :
282 0 : allocate( zm(nmlat,nlev),stat=astat )
283 0 : if( astat /= 0 ) then
284 0 : write(iulog,*) 'snoe_3d: failed to allocate zm; error = ',astat
285 0 : call endrun
286 : end if
287 0 : doy = aint( get_curr_calday() )
288 : #ifdef SNOE_DIAGS
289 : if( masterproc ) then
290 : write(iulog,*) ' '
291 : write(iulog,*) 'set_snoe_no: doy,kp,f107 = ',doy,kp,f107
292 : write(iulog,*) ' '
293 : end if
294 : #endif
295 : !----------------------------------------------------------------------
296 : ! ... obtain SNOE zonal mean data in geomagnetic coordinates
297 : !----------------------------------------------------------------------
298 0 : call snoe_zm( doy, kp, f107, zm )
299 :
300 : !----------------------------------------------------------------------
301 : ! ... map mean to model longitudes
302 : !----------------------------------------------------------------------
303 0 : do c = begchunk,endchunk
304 0 : ncol = get_ncols_p( c )
305 0 : call snoe_zmto3d( c, ncol, zm, nmlat, nlev )
306 : end do
307 0 : deallocate( zm )
308 :
309 0 : end subroutine snoe_timestep_init
310 :
311 0 : subroutine set_no_ubc( lchunk, ncol, zint, mmr, rho )
312 : !----------------------------------------------------------------------
313 : ! ... set no mixing ratio at the model top interface level
314 : !----------------------------------------------------------------------
315 :
316 0 : use constituents, only : pcnst, cnst_mw
317 :
318 : implicit none
319 :
320 : !----------------------------------------------------------------------
321 : ! ... dummy arguments
322 : !----------------------------------------------------------------------
323 : integer, intent(in) :: lchunk
324 : integer, intent(in) :: ncol
325 : real(r8), intent(in) :: zint(pcols) ! geopot height (km)
326 : real(r8), intent(in) :: rho(pcols) ! total atm density (kg/m^3)
327 : real(r8), intent(inout) :: mmr(pcols,pcnst) ! concentration at top interface (kg/kg)
328 :
329 : !----------------------------------------------------------------------
330 : ! ... local variables
331 : !----------------------------------------------------------------------
332 : real(r8), parameter :: amu_fac = 1.65979e-27_r8 ! kg/amu
333 : real(r8), parameter :: cm3_2_m3 = 1.e6_r8
334 : integer :: astat, i, k, ks, ks1
335 : real(r8) :: zinterp, zfac, mfac
336 :
337 : !----------------------------------------------------------------------
338 : ! ... map to model levels
339 : !----------------------------------------------------------------------
340 0 : if( ndx_no > 0 ) then
341 : column_loop : &
342 0 : do i = 1,ncol
343 0 : zinterp = zint(i)
344 0 : if( zinterp >= lev(1) ) then
345 0 : mmr(i,ndx_no) = snoe_no(i,1,lchunk)
346 0 : else if( zint(i) < lev(nlev) ) then
347 0 : mmr(i,ndx_no) = 0._r8
348 : else
349 0 : do ks = 2,nlev
350 0 : if( zinterp >= lev(ks) ) then
351 0 : ks1 = ks - 1
352 0 : zfac = (zinterp - lev(ks))/(lev(ks1) - lev(ks))
353 0 : mmr(i,ndx_no) = snoe_no(i,ks,lchunk) &
354 0 : + zfac*(snoe_no(i,ks1,lchunk) - snoe_no(i,ks,lchunk))
355 0 : exit
356 : end if
357 : end do
358 : end if
359 : end do column_loop
360 0 : mfac = amu_fac * cm3_2_m3 * cnst_mw(ndx_no)
361 0 : mmr(:ncol,ndx_no) = mmr(:ncol,ndx_no) * mfac /rho(:ncol)
362 : end if
363 :
364 0 : end subroutine set_no_ubc
365 :
366 0 : subroutine snoe_zm( doy, kp, f107, zm )
367 : !----------------------------------------------------------------------
368 : ! ... calculates zonal mean nitric oxide distribution on a given day
369 : ! and solar conditions (represented by the f10.7 and kp indices)
370 : !----------------------------------------------------------------------
371 :
372 : implicit none
373 :
374 : !----------------------------------------------------------------------
375 : ! ... dummy arguments
376 : !----------------------------------------------------------------------
377 : real(r8), intent(in) :: doy
378 : real(r8), intent(in) :: kp
379 : real(r8), intent(in) :: f107
380 : real(r8), intent(out) :: zm(:,:)
381 :
382 : !----------------------------------------------------------------------
383 : ! ... local variables
384 : !----------------------------------------------------------------------
385 : integer :: k
386 : real(r8) :: theta0 ! day number in degrees
387 : real(r8) :: dec ! solar declination angle
388 : real(r8) :: m1, m2, m3 ! coefficients for first 3 eofs
389 :
390 : #ifdef SNOE_DIAGS
391 : write(iulog,*) 'snoe_zm: doy, kp, f107 = ',doy,kp,f107
392 : #endif
393 :
394 : !----------------------------------------------------------------------
395 : ! ... calculate coefficients (m1 to m3) for eofs based on
396 : ! geophysical parametes eof1 - kp
397 : !----------------------------------------------------------------------
398 :
399 0 : m1 = kp * 0.785760_r8 - 1.94262_r8
400 :
401 : !----------------------------------------------------------------------
402 : ! ... eof2 - declination
403 : !----------------------------------------------------------------------
404 0 : theta0 = twopi * (doy - 1._r8) / 365._r8
405 :
406 : dec = .006918_r8 &
407 : - .399912_r8 * cos(theta0) + .070257_r8 * sin(theta0) &
408 : - .006758_r8 * cos(2._r8*theta0) + .000907_r8 * sin(2._r8*theta0) &
409 0 : - .002697_r8 * cos(3._r8*theta0) + .001480_r8 * sin(3._r8*theta0)
410 :
411 0 : dec = dec * r2d
412 :
413 : #ifdef SNOE_DIAGS
414 : write(iulog,*) 'snoe_zm: dec = ',dec
415 : #endif
416 0 : m2 = -.319782_r8 + dec*(.0973109_r8 + dec*(.000489814_r8 - dec*.000103608_r8))
417 :
418 : !----------------------------------------------------------------------
419 : ! ... eof3 - f107
420 : !----------------------------------------------------------------------
421 0 : m3 = log10(f107) * 6.44069_r8 - 13.9832_r8
422 :
423 : #ifdef SNOE_DIAGS
424 : write(iulog,*) 'snoe_zm: m1, m2, m3 = ',m1,m2,m3
425 : #endif
426 :
427 : !----------------------------------------------------------------------
428 : ! ... zonal mean distrib. is sum of mean and eofs
429 : !----------------------------------------------------------------------
430 0 : do k = 1,nlev
431 0 : zm(:,k) = no_mean(:,k) - m1 * eofs(:,k,1) + m2 * eofs(:,k,2) - m3 * eofs(:,k,3)
432 : end do
433 :
434 0 : end subroutine snoe_zm
435 :
436 0 : subroutine snoe_zmto3d( lchunk, ncol, zm, nmlat, nlev )
437 : !----------------------------------------------------------------------
438 : ! ... interpolate zonal mean on magnetic grid
439 : ! to 3d field in geographic coords
440 : !----------------------------------------------------------------------
441 :
442 : implicit none
443 :
444 : !----------------------------------------------------------------------
445 : ! ... dummy arguments
446 : !----------------------------------------------------------------------
447 : integer, intent(in) :: lchunk
448 : integer, intent(in) :: ncol
449 : integer, intent(in) :: nmlat
450 : integer, intent(in) :: nlev
451 : real(r8), intent(in) :: zm(nmlat,nlev) ! zonal mean snoe no concentration
452 :
453 : !----------------------------------------------------------------------
454 : ! ... local variables
455 : !----------------------------------------------------------------------
456 : integer :: k
457 :
458 0 : do k = 1,nlev
459 0 : call interpol( nmlat, mlat2d, zm(1,k), ncol, mlat(1,lchunk), snoe_no(1,k,lchunk) )
460 : end do
461 :
462 0 : end subroutine snoe_zmto3d
463 :
464 0 : subroutine interpol( nin, xin, yin, nout, xout, yout )
465 : !-----------------------------------------------------------------------
466 : ! ... linear interpolation
467 : ! does not extrapolate, but repeats edge values
468 : !-----------------------------------------------------------------------
469 :
470 : implicit none
471 :
472 : !-----------------------------------------------------------------------
473 : ! ... dummy arguments
474 : !-----------------------------------------------------------------------
475 : integer, intent(in) :: nin, nout
476 : real(r8), intent(in) :: xin(nin)
477 : real(r8), intent(in) :: yin(nin)
478 : real(r8), intent(in) :: xout(nout)
479 : real(r8), intent(out) :: yout(nout)
480 :
481 : !-----------------------------------------------------------------------
482 : ! ... local variables
483 : !-----------------------------------------------------------------------
484 : integer :: i, j
485 :
486 0 : do j = 1,nout
487 0 : if( xout(j) < xin(1) ) then
488 0 : yout(j) = yin(1)
489 : else
490 0 : if( xout(j) > xin(nin) ) then
491 0 : yout(j) = yin(nin)
492 : else
493 0 : do i = 1, nin-1
494 0 : if ((xout(j) >= xin(i)) .and. (xout(j) < xin(i+1)) ) then
495 : yout(j) = yin(i) &
496 0 : + (yin(i+1) - yin(i)) * (xout(j) - xin(i)) / (xin(i+1) - xin(i))
497 : end if
498 : end do
499 : end if
500 : end if
501 : end do
502 :
503 0 : end subroutine interpol
504 :
505 : end module mo_snoe
|
