module module_mozart_lib implicit none !--------------------------------------------------------------- ! include files !--------------------------------------------------------------- include 'netcdf.inc' !--------------------------------------------------------------- ! public procedures !--------------------------------------------------------------- public :: init_mozart_lib public :: bc_interpolate4d public :: ic_interpolate4d public :: exit_mozart_lib public :: get_moz_time_ndx !--------------------------------------------------------------- ! private procedures !--------------------------------------------------------------- private :: handle_error private :: lotim !--------------------------------------------------------------- ! public variables !--------------------------------------------------------------- integer, public :: nlon, nlat, nlev ! data dimension length integer, public :: year, month, day integer, public :: hour, minute, second integer, public :: ntime_m real, public, allocatable :: moz_times(:) character(len=9) :: moz_var_suffix = '_VMR_inst' character(len=16) :: src_lon_dim_name = 'lon' character(len=16) :: src_lon_var_name = 'lon' character(len=16) :: src_lat_dim_name = 'lat' character(len=16) :: src_lat_var_name = 'lat' character(len=16) :: src_lev_dim_name = 'lev' character(len=16) :: src_lev_var_name = 'lev' character(len=16) :: src_hybrid_a_name = 'hyam' character(len=16) :: src_hybrid_b_name = 'hybm' character(len=16) :: src_srf_press_name = 'PS' character(len=16) :: src_ref_press_name = 'P0' !--------------------------------------------------------------- ! variables for reading netCDF data !--------------------------------------------------------------- integer, private :: ncid = 0 ! netCDF file ID integer, private :: nstt(4), ncnt(4) ! start index and counter array integer, private :: start_index = 1 ! monthly mean data integer, private :: ntime ! data dimension length integer, private :: nscalar, nstring ! data dimension length integer, private, allocatable :: moz_date(:) integer, private, allocatable :: moz_datesec(:) !--------------------------------------------------------------- ! variables used by interpolation !--------------------------------------------------------------- ! real, private, parameter :: ps0 = 1.e5 ! Pa real, private, parameter :: THREE60 = 360. integer, private, allocatable :: ix(:,:,:) integer, private, allocatable :: jy(:,:,:) ! index used by interpolation real, private :: ps0 = 1.e5 ! Pa real, private, allocatable :: ax(:,:,:) ! weight coef. all domain real, private, allocatable :: by(:,:,:) ! weight coef. all domain real, private, allocatable :: ps_moz(:,:) real, private, allocatable :: ps_mozi(:,:) real, private, allocatable :: hyam(:) real, private, allocatable :: hybm(:) character(len=16) :: dbg_species = ' ' logical, private :: reorder_lats type t_type integer :: ndx integer :: ncid_lo, ncid_hi integer :: lo_moz_ndx, hi_moz_ndx integer :: lo_buf_ndx, hi_buf_ndx integer :: gap_date, gap_secs real :: dels real, pointer :: conc(:,:,:,:) real, pointer :: ps(:,:,:) logical :: in_gap logical :: t_interp end type t_type type(t_type), private :: time_type contains subroutine init_mozart_lib( moz_dir, moz_fn, x2d, y2d, wrf_date, & wrf_datesec, nx, ny, nspec ) !--------------------------------------------------------------- ! initialize netcdf file !--------------------------------------------------------------- use mo_calendar, only : diffdat use utils, only : wrf2mz_map !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: wrf_date integer, intent(in) :: wrf_datesec integer, intent(in) :: nx, ny integer, intent(in) :: nspec character(*), intent(in) :: moz_dir character(*), intent(inout) :: moz_fn real, intent(in) :: x2d(nx,ny) real, intent(in) :: y2d(nx,ny) !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: status, tstat integer :: dimid, varid integer :: i, j, jtmp, k, n integer :: jbeg, jend, jinc integer :: vid real :: model_x, model_y, model_dx real :: wrk real, allocatable :: x1d(:), y1d(:) character(len=128) :: message character(len=128) :: filenm character(len=32) :: spcnam logical :: found logical :: monotonic_moz_x filenm = trim( moz_dir ) // adjustl( moz_fn ) !--------------------------------------------------------------- ! open the initial mozart netCDF file !--------------------------------------------------------------- status = nf_open( trim(filenm), nf_nowrite, ncid ) if( status /= nf_noerr ) then message = 'init_mozart_lib: failed to open ' // trim(filenm) call handle_error( message,status ) end if write(*,*) 'init_mozart_lib: opened ',trim(filenm) !--------------------------------------------------------------- ! read times !--------------------------------------------------------------- call read_moz_times( ncid ) time_type%ndx = lotim( wrf_date, wrf_datesec, moz_date, moz_datesec, ntime_m ) if( time_type%ndx == 0 ) then status = nf_close( ncid ) call next_flnm( moz_fn, .false. ) filenm = trim( moz_dir ) // adjustl( moz_fn ) !--------------------------------------------------------------- ! open the initial mozart netCDF file !--------------------------------------------------------------- status = nf_open( trim(filenm), nf_nowrite, ncid ) if( status /= nf_noerr ) then message = 'init_mozart_lib: failed to open ' // trim(filenm) call handle_error( message,status ) end if write(*,*) 'init_mozart_lib: opened ',trim(filenm) call read_moz_times( ncid ) time_type%ndx = lotim( wrf_date, wrf_datesec, moz_date, moz_datesec, ntime_m ) if( time_type%ndx == 0 ) then write(*,*) 'init_mozart_lib: time ',wrf_date,' ',wrf_datesec write(*,*) ' is before ' write(*,*) ' ',moz_date(1),' ',moz_datesec(1) write(*,*) ' the first time in mozart file ',trim(moz_fn) stop end if end if !--------------------------------------------------------------- ! check mozart variables !--------------------------------------------------------------- call chk_moz_vars( nspec, filenm ) !--------------------------------------------------------------- ! get the longitude dimension length of the MOZART data !--------------------------------------------------------------- message = 'init_mozart_lib: failed to get ' // trim(src_lon_dim_name) // ' id' status = nf_inq_dimid( ncid, trim(src_lon_dim_name), dimid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to get ' // trim(src_lon_dim_name) // ' dimenstion size' status = nf_inq_dimlen( ncid, dimid, nlon ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------- ! get the latitude dimension length of the MOZART data !--------------------------------------------------------------- message = 'init_mozart_lib: failed to get ' // trim(src_lat_dim_name) // ' id' status = nf_inq_dimid( ncid, trim(src_lat_dim_name), dimid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to get ' // trim(src_lat_dim_name) // ' dimenstion size' status = nf_inq_dimlen( ncid, dimid, nlat ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------- ! get the vertical dimension length of the MOZART data !--------------------------------------------------------------- message = 'init_mozart_lib: failed to get ' // trim(src_lev_dim_name) // ' id' status = nf_inq_dimid( ncid, trim(src_lev_dim_name), dimid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to get ' // trim(src_lev_dim_name) // ' dimenstion size' status = nf_inq_dimlen( ncid, dimid, nlev ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------- ! read vertical coordinates !--------------------------------------------------------------- status = nf_inq_varid( ncid, trim(src_ref_press_name), varid ) if( status == nf_noerr ) then message = 'init_mozart_lib: failed to read ' // trim(src_ref_press_name) // ' value' status = nf_get_var_real( ncid, varid, wrk ) if( status /= nf_noerr ) call handle_error( message,status ) ps0 = wrk endif allocate( hyam(nlev), hybm(nlev), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate hyam, hybm; error = ',status stop 'Allocation error' end if message = 'init_mozart_lib: failed to get ' // trim(src_hybrid_a_name) // ' id' status = nf_inq_varid( ncid, trim(src_hybrid_a_name), varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to read ' // trim(src_hybrid_a_name) // ' array' status = nf_get_var_real( ncid, varid, hyam ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to get ' // trim(src_hybrid_b_name) // ' id' status = nf_inq_varid( ncid, trim(src_hybrid_b_name), varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to read ' // trim(src_hybrid_b_name) // ' array' status = nf_get_var_real( ncid, varid, hybm ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------- ! allocate longitude,latitude arrays !--------------------------------------------------------------- allocate( x1d(nlon), y1d(nlat), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate x1d, y1d; error = ',status stop 'Allocation error' end if !--------------------------------------------------------------- ! read longitudes !--------------------------------------------------------------- message = 'init_mozart_lib: failed to get ' // trim(src_lon_var_name) // ' id' status = nf_inq_varid( ncid, trim(src_lon_var_name), varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to read ' // trim(src_lon_var_name) // ' array' status = nf_get_var_real( ncid, varid, x1d ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------- ! check mozart longitudes for monotonicity !--------------------------------------------------------------- monotonic_moz_x = x1d(nlon) > x1d(1) if( .not. monotonic_moz_x ) then do i = 2,nlon if( x1d(i) < x1d(i-1) ) then x1d(i:nlon) = x1d(i:nlon) + THREE60 exit endif end do endif !--------------------------------------------------------------- ! read latitudes !--------------------------------------------------------------- message = 'init_mozart_lib: failed to get ' // trim(src_lat_var_name) // ' id' status = nf_inq_varid( ncid, trim(src_lat_var_name), varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'init_mozart_lib: failed to read ' // trim(src_lat_var_name) // ' array' status = nf_get_var_real( ncid, varid, y1d ) if( status /= nf_noerr ) call handle_error( message,status ) !--------------------------------------------------------------------- ! are lats monotonically increasing? !--------------------------------------------------------------------- reorder_lats = all( y1d(2:nlat) < y1d(1:nlat-1) ) !--------------------------------------------------------------- ! allocate memory space to store interpolation coefficients !--------------------------------------------------------------- tstat = 0 allocate( ax(nx,ny,0:1), by(nx,ny,0:1), stat=status ) tstat = tstat + status allocate( ix(nx,ny,0:1), jy(nx,ny,0:1), stat=status ) tstat = tstat + status if( tstat /= 0 ) then write(*,*) 'allocate for ax ... jy failed; error = ',tstat end if !--------------------------------------------------------------- ! horizontal interpolation coefs. !--------------------------------------------------------------- !--------------------------------------------------------------- ! all domain !--------------------------------------------------------------- if( .not. reorder_lats ) then jbeg = 1 ; jend = nlat ; jinc = 1 else jbeg = nlat ; jend = 1 ; jinc = -1 endif do j = 1, ny do i = 1, nx model_x = THREE60 + x2d(i,j) if( monotonic_moz_x ) then model_x = mod( model_x,THREE60 ) endif if( model_x >= x1d(nlon) ) then ix(i,j,0) = nlon else if( model_x < x1d(1) ) then ix(i,j,0) = nlon else do n = 1, nlon if( model_x < x1d(n) ) then ix(i,j,0) = min( nlon-1, max(n-1,1) ) exit end if end do end if ix(i,j,1) = mod( ix(i,j,0),nlon ) + 1 model_dx = x1d(ix(i,j,1)) - x1d(ix(i,j,0)) if( model_dx < 0. ) then model_dx = THREE60 + model_dx end if ax(i,j,0) = (model_x - x1d(ix(i,j,0)))/model_dx ax(i,j,1) = 1.0 - ax(i,j,0) model_y = y2d(i,j) do n = jbeg,jend,jinc if( model_y < y1d(n) ) then exit end if end do if( .not. reorder_lats ) then jy(i,j,0) = min( nlat-1, max(n-1,1) ) jy(i,j,1) = jy(i,j,0) + 1 else jy(i,j,1) = min( nlat-1, max(n,1) ) jy(i,j,0) = jy(i,j,1) + 1 endif by(i,j,0) = (model_y - y1d(jy(i,j,0)))/(y1d(jy(i,j,1)) - y1d(jy(i,j,0))) by(i,j,1) = 1.0 - by(i,j,0) end do end do !--------------------------------------------------------------- ! release memory !--------------------------------------------------------------- deallocate( x1d, y1d ) !--------------------------------------------------------------- ! allocate memory space for reading !--------------------------------------------------------------- allocate( time_type%conc(nlon,nlat,nlev,2), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate time_type%conc; error = ',status stop end if !--------------------------------------------------------------- ! allocate mozart surface pressure arrays !--------------------------------------------------------------- allocate( ps_moz(nlon,nlat), ps_mozi(nx,ny), & time_type%ps(nlon,nlat,2), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate ps_moz,ps_mozi; error = ',status stop end if !--------------------------------------------------------------- ! setup time_type !--------------------------------------------------------------- time_type%lo_moz_ndx = 0 time_type%hi_moz_ndx = 0 time_type%lo_buf_ndx = 1 time_type%hi_buf_ndx = 2 time_type%ncid_lo = 0 time_type%ncid_hi = 0 time_type%in_gap = .false. time_type%t_interp = .false. if( time_type%ndx > 0 ) then time_type%lo_moz_ndx = time_type%ndx time_type%ncid_lo = ncid time_type%ncid_hi = ncid endif write(*,*) 'finished init_mozart_lib' end subroutine init_mozart_lib subroutine read_moz_times( ncid ) !--------------------------------------------------------------- ! read times from current mozart input file !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: ncid !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: varid integer :: status character(len=128) :: message message = 'read_moz_times: failed to get time id' status = nf_inq_dimid( ncid, 'time', varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'read_moz_times: failed to get time size' status = nf_inq_dimlen( ncid, varid, ntime_m ) if( status /= nf_noerr ) call handle_error( message,status ) if( allocated( moz_date ) ) then deallocate( moz_date, moz_datesec ) end if allocate( moz_date(ntime_m), moz_datesec(ntime_m), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate date, datesec; error = ',status stop 'Allocation error' end if message = 'read_moz_times: failed to get date id' status = nf_inq_varid( ncid, 'date', varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'read_moz_times: failed to read date array' status = nf_get_var_int( ncid, varid, moz_date ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'read_moz_times: failed to get datesec id' status = nf_inq_varid( ncid, 'datesec', varid ) if( status /= nf_noerr ) call handle_error( message,status ) message = 'read_moz_times: failed to read datesec array' status = nf_get_var_int( ncid, varid, moz_datesec ) if( status /= nf_noerr ) call handle_error( message,status ) end subroutine read_moz_times subroutine get_moz_time_ndx( moz_dir, moz_fn, wrf_date, wrf_datesec, nspec ) !--------------------------------------------------------------- ! get mozart time index for wrf_date, wrf_datesec !--------------------------------------------------------------- use mo_calendar, only : diffdat !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: wrf_date integer, intent(in) :: wrf_datesec integer, intent(in) :: nspec character(len=*), intent(in) :: moz_dir character(len=*), intent(inout) :: moz_fn !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: i, n integer :: status character(len=128) :: filenm character(len=128) :: message logical :: found write(*,*) ' ' write(*,*) 'get_moz_time_ndx; moz_dir,moz_fn = ',trim(moz_dir),trim(moz_fn) write(*,*) 'get_moz_time_ndx; wrf_date,wrf_datesec,ntime_m = ',wrf_date,wrf_datesec,ntime_m n = lotim( wrf_date, wrf_datesec, moz_date, moz_datesec, ntime_m ) if( n > 0 ) then !--------------------------------------------------------------- ! wrf time in present mozart dataset !--------------------------------------------------------------- write(*,*) 'get_moz_time_ndx; moz_tndx = ',n time_type%lo_moz_ndx = n time_type%in_gap = .false. if( time_type%t_interp ) then time_type%hi_moz_ndx = n + 1 endif if( time_type%ncid_hi /= time_type%ncid_lo ) then status = nf_close( time_type%ncid_lo ) if( status /= 0 ) then filenm = trim( moz_dir ) // adjustl( moz_fn ) write(*,*) 'get_moz_time_ndx: failed to close ',trim(filenm),' ; error = ',status stop end if time_type%ncid_lo = time_type%ncid_hi endif time_type%ncid_hi = time_type%ncid_lo else if( n < 0 ) then !--------------------------------------------------------------- ! wrf time after present mozart dataset !--------------------------------------------------------------- time_type%ncid_lo = ncid time_type%lo_moz_ndx = ntime_m call next_flnm( moz_fn, .true. ) filenm = trim( moz_dir ) // adjustl( moz_fn ) !--------------------------------------------------------------- ! open the input netCDF file !--------------------------------------------------------------- message = 'get_moz_time_ndx: failed to open ' // trim(filenm) status = nf_open( trim(filenm), nf_nowrite, ncid ) if( status /= nf_noerr ) call handle_error( message,status ) write(*,*) 'get_moz_time_ndx: opened ',trim(filenm) !--------------------------------------------------------------- ! check mozart variables !--------------------------------------------------------------- call chk_moz_vars( nspec, filenm ) time_type%gap_date = moz_date(ntime_m) time_type%gap_secs = moz_datesec(ntime_m) call read_moz_times( ncid ) n = lotim( wrf_date, wrf_datesec, moz_date, moz_datesec, ntime_m ) time_type%ndx = n if( n > 0 ) then write(*,*) 'get_moz_time_ndx; moz_tndx = ',n status = nf_close( time_type%ncid_lo ) if( status /= 0 ) then filenm = trim( moz_dir ) // adjustl( moz_fn ) write(*,*) 'get_moz_time_ndx: failed to close ',trim(filenm),' ; error = ',status stop end if time_type%in_gap = .false. time_type%ncid_lo = ncid time_type%ncid_hi = ncid time_type%lo_moz_ndx = n else if( n == 0 ) then time_type%in_gap = .true. time_type%hi_moz_ndx = 1 time_type%ncid_hi = ncid time_type%dels = time_type%dels/diffdat( time_type%gap_date, time_type%gap_secs, moz_date(1), moz_datesec(1) ) time_type%t_interp = .true. else write(*,*) 'get_moz_time_ndx: failed to find ',wrf_date,' : ',wrf_datesec write(*,*) ' in file ',trim(filenm) stop end if else if( time_type%in_gap ) then time_type%dels = time_type%dels/diffdat( time_type%gap_date, time_type%gap_secs, moz_date(1), moz_datesec(1) ) end if end subroutine get_moz_time_ndx subroutine next_flnm( filenm, incr ) !--------------------------------------------------------------- ! increment mozart filename !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- character(len=*), intent(inout) :: filenm logical, intent(in) :: incr !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: il, iu, ios integer :: file_number character(len=5) :: numa logical :: found found = .false. iu = scan( trim(filenm), '0123456789', back=.true. ) do il = iu-1,1,-1 if( index( '0123456789', filenm(il:il) ) == 0 ) then found = .true. exit end if end do if( .not. found ) then write(*,*) 'next_filenm: mozart file ',trim(filenm),' is not proper format' stop end if il = il + 1 write(*,*) ' ' write(*,*) 'next_flnm; trying to increment file ',trim(filenm) write(*,*) 'next_flnm; il, iu = ',il,iu read(filenm(il:iu),*,iostat=ios) file_number if( ios /= 0 ) then write(*,*) 'next_filenm: failed to read ',filenm(il:iu),' ; error = ',ios stop end if write(*,*) 'next_flnm; file_number = ',file_number if( incr ) then file_number = file_number + 1 else file_number = file_number - 1 endif write(numa,'(i5)') file_number+10000 filenm(il:iu) = numa(2:5) write(*,*) 'next_flnm; new file = ',trim(filenm) end subroutine next_flnm subroutine read_mozart_ps( bndx, tndx, ncid ) !--------------------------------------------------------------- ! read mozart surface pressure !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: bndx ! buffer index integer, intent(in) :: tndx ! time index integer, intent(in) :: ncid ! netcdf file index !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: status integer :: varid character(len=128) :: message status = nf_inq_varid( ncid, trim(src_srf_press_name), varid ) if( status /= nf_noerr ) then message = 'read_mozart_ps: failed to get ' // trim(src_srf_press_name) // ' id' call handle_error( message,status ) end if nstt(1:3) = (/ 1, 1, tndx /) ncnt(1:3) = (/ nlon, nlat, 1 /) status = nf_get_vara_real( ncid, varid, nstt(1:3), ncnt(1:3), time_type%ps(:,:,bndx) ) if( status /= nf_noerr ) then message = 'read_mozart_ps: failed to read ' // trim(src_srf_press_name) // ' array' call handle_error( message,status ) end if end subroutine read_mozart_ps subroutine tinterp_mozart_ps( dels ) !--------------------------------------------------------------- ! time interpolation of mozart surface pressure !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- real, intent(in) :: dels ! linear interp factor !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: j integer :: lo_ndx, hi_ndx real :: delsm1 if( time_type%t_interp ) then delsm1 = 1. - dels !--------------------------------------------------------------- ! interpolate mozart surface pressure to wrf grid !--------------------------------------------------------------- do j = 1,nlat ps_moz(:,j) = time_type%ps(:,j,1) * dels & + time_type%ps(:,j,2) * delsm1 end do else do j = 1,nlat ps_moz(:,j) = time_type%ps(:,j,1) end do endif end subroutine tinterp_mozart_ps subroutine hinterp_mozart_ps( nx, ny ) !--------------------------------------------------------------- ! horizontal interpolation of mozart surface pressure !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: nx ! wrf x dimension integer, intent(in) :: ny ! wrf y dimension !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: i, j !--------------------------------------------------------------- ! interpolate mozart surface pressure to wrf grid !--------------------------------------------------------------- do j = 1,ny do i = 1,nx ps_mozi(i,j) = ps_moz(ix(i,j,0),jy(i,j,0))*ax(i,j,1)*by(i,j,1) & + ps_moz(ix(i,j,0),jy(i,j,1))*ax(i,j,1)*by(i,j,0) & + ps_moz(ix(i,j,1),jy(i,j,0))*ax(i,j,0)*by(i,j,1) & + ps_moz(ix(i,j,1),jy(i,j,1))*ax(i,j,0)*by(i,j,0) end do end do end subroutine hinterp_mozart_ps subroutine bc_interpolate4d( ndx, wrfxs, wrfxe, wrfys, wrfye, & ps_wrf, znu, p_top, nx, ny, & nz, nw ) !--------------------------------------------------------------- ! interpolate four-dimensional field ... !--------------------------------------------------------------- use utils, only : wrf2mz_map !--------------------------------------------------------------- ! input arguments !--------------------------------------------------------------- integer, intent(in) :: ndx integer, intent(in) :: nx, ny, nz, nw ! dimensions real, intent(in) :: p_top ! wrf top reference pressure (Pa) real, intent(in) :: ps_wrf(nx,ny) ! wrf surface pressure (Pa) real, intent(in) :: znu(nz) ! sigma coordinates real, intent(out) :: wrfxs(ny,nz,nw) ! wrfchem vmr(ppm) real, intent(out) :: wrfxe(ny,nz,nw) ! wrfchem vmr(ppm) real, intent(out) :: wrfys(nx,nz,nw) ! wrfchem vmr(ppm) real, intent(out) :: wrfye(nx,nz,nw) ! wrfchem vmr(ppm) !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: status, varid real, allocatable :: wrk(:) real, allocatable :: wrk1(:) real, allocatable :: p_moz(:) real, allocatable :: p_wrf(:) real :: mozval(nlon,nlat,nlev) character(len=20) :: mozspn integer :: i, j, k, n if( wrf2mz_map(ndx)%moz_cnt > 0 ) then allocate( wrk(nlev), wrk1(nz), p_moz(nlev), p_wrf(nz), stat=status ) if( status /= 0 ) then write(*,*) 'bc_interpolated4d: failed to allocate wrk ... p_wrf; error = ',status stop end if !--------------------------------------------------------------- ! read mozart pressure !--------------------------------------------------------------- if( ndx == 1 ) then call read_mozart_ps( time_type%lo_buf_ndx, time_type%lo_moz_ndx, time_type%ncid_lo ) if( time_type%t_interp ) then call read_mozart_ps( time_type%hi_buf_ndx, time_type%hi_moz_ndx, time_type%ncid_hi ) endif call tinterp_mozart_ps( time_type%dels ) call hinterp_mozart_ps( nx, ny ) endif !--------------------------------------------------------------- ! read mozart species !--------------------------------------------------------------- call read_moz_species( ndx, time_type%lo_buf_ndx, time_type%lo_moz_ndx, time_type%ncid_lo ) if( time_type%t_interp ) then call read_moz_species( ndx, time_type%hi_buf_ndx, time_type%hi_moz_ndx, time_type%ncid_hi ) endif call tinterp_moz_species( time_type%dels, mozval ) !--------------------------------------------------------------- ! horizontally interpolate species at boundaries !--------------------------------------------------------------- ! west !--------------------------------------------------------------- do j = 1,ny wrk(:) = mozval(ix(1,j,0),jy(1,j,0),:)*ax(1,j,1)*by(1,j,1) & + mozval(ix(1,j,0),jy(1,j,1),:)*ax(1,j,1)*by(1,j,0) & + mozval(ix(1,j,1),jy(1,j,0),:)*ax(1,j,0)*by(1,j,1) & + mozval(ix(1,j,1),jy(1,j,1),:)*ax(1,j,0)*by(1,j,0) p_moz(:) = ps_mozi(1,j)*hybm(:) + ps0*hyam(:) p_wrf(:) = ps_wrf(1,j)*znu(nz:1:-1) + (1. - znu(nz:1:-1))*p_top call vinterp( p_moz, p_wrf, wrk, wrk1, nz, nlev ) do n = 1,nw wrfxs(j,:,n) = wrk1(nz:1:-1) end do end do !--------------------------------------------------------------- ! east !--------------------------------------------------------------- do j = 1, ny wrk(:) = mozval(ix(nx,j,0),jy(nx,j,0),:)*ax(nx,j,1)*by(nx,j,1) & + mozval(ix(nx,j,0),jy(nx,j,1),:)*ax(nx,j,1)*by(nx,j,0) & + mozval(ix(nx,j,1),jy(nx,j,0),:)*ax(nx,j,0)*by(nx,j,1) & + mozval(ix(nx,j,1),jy(nx,j,1),:)*ax(nx,j,0)*by(nx,j,0) p_moz(:) = ps_mozi(nx,j)*hybm(:) + ps0*hyam(:) p_wrf(:) = ps_wrf(nx,j)*znu(nz:1:-1) + (1. - znu(nz:1:-1))*p_top call vinterp( p_moz, p_wrf, wrk, wrk1, nz, nlev ) do n = 1,nw wrfxe(j,:,n) = wrk1(nz:1:-1) end do end do !--------------------------------------------------------------- ! north !--------------------------------------------------------------- do i = 1,nx wrk(:) = mozval(ix(i,ny,0),jy(i,ny,0),:)*ax(i,ny,1)*by(i,ny,1) & + mozval(ix(i,ny,0),jy(i,ny,1),:)*ax(i,ny,1)*by(i,ny,0) & + mozval(ix(i,ny,1),jy(i,ny,0),:)*ax(i,ny,0)*by(i,ny,1) & + mozval(ix(i,ny,1),jy(i,ny,1),:)*ax(i,ny,0)*by(i,ny,0) p_moz(:) = ps_mozi(i,ny)*hybm(:) + ps0*hyam(:) p_wrf(:) = ps_wrf(i,ny)*znu(nz:1:-1) + (1. - znu(nz:1:-1))*p_top call vinterp( p_moz, p_wrf, wrk, wrk1, nz, nlev ) do n = 1,nw wrfye(i,:,n) = wrk1(nz:1:-1) end do end do !--------------------------------------------------------------- ! south !--------------------------------------------------------------- do i = 1,nx wrk(:) = mozval(ix(i,1,0),jy(i,1,0),:)*ax(i,1,1)*by(i,1,1) & + mozval(ix(i,1,0),jy(i,1,1),:)*ax(i,1,1)*by(i,1,0) & + mozval(ix(i,1,1),jy(i,1,0),:)*ax(i,1,0)*by(i,1,1) & + mozval(ix(i,1,1),jy(i,1,1),:)*ax(i,1,0)*by(i,1,0) p_moz(:) = ps_mozi(i,1)*hybm(:) + ps0*hyam(:) p_wrf(:) = ps_wrf(i,1)*znu(nz:1:-1) + (1. - znu(nz:1:-1))*p_top call vinterp( p_moz, p_wrf, wrk, wrk1, nz, nlev ) do n = 1,nw wrfys(i,:,n) = wrk1(nz:1:-1) end do end do deallocate( wrk, wrk1, p_moz, p_wrf ) else do n = 1,nw wrfxs(:,:,n) = wrf2mz_map(ndx)%wrf_conc wrfxe(:,:,n) = wrf2mz_map(ndx)%wrf_conc wrfys(:,:,n) = wrf2mz_map(ndx)%wrf_conc wrfye(:,:,n) = wrf2mz_map(ndx)%wrf_conc end do end if if( trim(wrf2mz_map(ndx)%wrf_name) == trim(dbg_species) ) then write(*,*) ' ' write(*,*) 'mozart ',trim(dbg_species),' east bc values' write(*,'(10(1p,g15.7))') wrfxe(:,1,1) end if end subroutine bc_interpolate4d subroutine ic_interpolate4d( ndx, conc, ps_wrf, znu, p_top, & nx, ny, nz ) !--------------------------------------------------------------- ! interpolate four-dimensional field ... !--------------------------------------------------------------- use utils, only : wrf2mz_map !--------------------------------------------------------------- ! input arguments !--------------------------------------------------------------- integer, intent(in) :: ndx integer, intent(in) :: nx, ny, nz ! dimensions real, intent(in) :: p_top ! wrf top reference pressure (Pa) real, intent(in) :: ps_wrf(nx,ny) ! wrf surface pressure (Pa) real, intent(in) :: znu(nz) ! sigma coordinates real, intent(out) :: conc(nx,ny,nz) ! wrfchem vmr(ppm) !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: i, j, k, n integer :: status, varid real :: mozval(nlon,nlat,nlev) real, allocatable :: wrk(:) real, allocatable :: wrk1(:) real, allocatable :: p_moz(:) real, allocatable :: p_wrf(:) character(len=20) :: mozspn if( wrf2mz_map(ndx)%moz_cnt > 0 ) then allocate( wrk(nlev), wrk1(nz), p_moz(nlev), p_wrf(nz), stat=status ) if( status /= 0 ) then write(*,*) 'failed to allocate wrk ... p_wrf; error = ',status stop end if !--------------------------------------------------------------- ! read mozart pressure !--------------------------------------------------------------- if( ndx == 1 ) then call read_mozart_ps( time_type%lo_buf_ndx, time_type%lo_moz_ndx, time_type%ncid_lo ) if( time_type%t_interp ) then call read_mozart_ps( time_type%hi_buf_ndx, time_type%hi_moz_ndx, time_type%ncid_hi ) endif call tinterp_mozart_ps( time_type%dels ) call hinterp_mozart_ps( nx, ny ) endif !--------------------------------------------------------------- ! read mozart species !--------------------------------------------------------------- call read_moz_species( ndx, time_type%lo_buf_ndx, time_type%lo_moz_ndx, time_type%ncid_lo ) if( time_type%t_interp ) then call read_moz_species( ndx, time_type%hi_buf_ndx, time_type%hi_moz_ndx, time_type%ncid_hi ) endif call tinterp_moz_species( time_type%dels, mozval ) !--------------------------------------------------------------- ! horizontally interpolate species !--------------------------------------------------------------- do j = 1,ny do i = 1,nx wrk(:) = mozval(ix(i,j,0),jy(i,j,0),:)*ax(i,j,1)*by(i,j,1) & + mozval(ix(i,j,0),jy(i,j,1),:)*ax(i,j,1)*by(i,j,0) & + mozval(ix(i,j,1),jy(i,j,0),:)*ax(i,j,0)*by(i,j,1) & + mozval(ix(i,j,1),jy(i,j,1),:)*ax(i,j,0)*by(i,j,0) p_moz(:) = ps_mozi(i,j)*hybm(:) + ps0*hyam(:) p_wrf(:) = ps_wrf(i,j)*znu(nz:1:-1) + (1. - znu(nz:1:-1))*p_top call vinterp( p_moz, p_wrf, wrk, wrk1, nz, nlev ) conc(i,j,:) = wrk1(nz:1:-1) if( trim(wrf2mz_map(ndx)%wrf_name) == trim(dbg_species) ) then if( i == 91 .and. j == 93 ) then write(*,*) 'interpolation points' write(*,*) 'ix = ',ix(i,j,:) write(*,*) 'jy = ',jy(i,j,:) write(*,*) 'interpolation wghts' write(*,*) 'ax = ',ax(i,j,:) write(*,*) 'by = ',by(i,j,:) write(*,*) 'mozart ',trim(dbg_species),' values @ i,j =',ix(i,j,0),jy(i,j,0) write(*,'(10(1p,g15.7))') mozval(ix(i,j,0),jy(i,j,0),:) write(*,*) 'horiz interpolated mozart ',trim(dbg_species),' values' write(*,'(10(1p,g15.7))') wrk write(*,*) 'vert interpolated mozart ',trim(dbg_species),' values' write(*,'(10(1p,g15.7))') wrk1 end if end if end do end do if( trim(wrf2mz_map(ndx)%wrf_name) == trim(dbg_species) ) then write(*,*) ' ' write(*,*) 'mozart ',trim(dbg_species),' ic values' write(*,'(10(1p,g15.7))') conc(nx,:,1) end if deallocate( wrk, wrk1, p_moz, p_wrf ) else conc(:,:,:) = wrf2mz_map(ndx)%wrf_conc end if end subroutine ic_interpolate4d subroutine read_moz_species( sndx, bndx, tndx, ncid ) !--------------------------------------------------------------- ! read mozart species !--------------------------------------------------------------- use utils, only : wrf2mz_map !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: sndx ! species index integer, intent(in) :: bndx ! buffer index integer, intent(in) :: tndx ! time index integer, intent(in) :: ncid ! netcdf file index !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: n integer :: nspec integer :: status integer :: varid real :: accum(nlon,nlat,nlev) real :: tmpval(nlon,nlat,nlev) character(len=20) :: mozspn nspec = wrf2mz_map(sndx)%moz_cnt if( nspec > 0 ) then nstt(1:4) = (/ 1, 1, 1, tndx /) ncnt(1:4) = (/ nlon, nlat, nlev, 1 /) accum(:,:,:) = 0. do n = 1,nspec if( wrf2mz_map(sndx)%moz_ext(n) ) then mozspn = trim(wrf2mz_map(sndx)%moz_names(n)) // moz_var_suffix else mozspn = trim(wrf2mz_map(sndx)%moz_names(n)) end if status = nf_inq_varid( ncid, mozspn, varid ) if( status /= nf_noerr ) then write(*,*) 'read_moz_species: failed to get id of ',mozspn call handle_error( ' ',status ) end if status = nf_get_vara_real( ncid, varid, nstt, ncnt, tmpval ) if( status /= nf_noerr ) then write(*,*) 'read_moz_species: failed to read ',mozspn call handle_error( ' ',status ) end if if( wrf2mz_map(sndx)%moz_wght(n) == 1. ) then accum(:,:,:) = accum(:,:,:) + tmpval(:,:,:) else accum(:,:,:) = accum(:,:,:) + wrf2mz_map(sndx)%moz_wght(n)*tmpval(:,:,:) end if end do accum(:,:,:) = accum(:,:,:)*wrf2mz_map(sndx)%wrf_wght end if time_type%conc(:,:,:,bndx) = accum(:,:,:) end subroutine read_moz_species subroutine tinterp_moz_species( dels, mozval ) !--------------------------------------------------------------- ! time interpolation of mozart species concentration !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- real, intent(in) :: dels ! linear interp factor real, intent(out) :: mozval(nlon,nlat,nlev) !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: j, k integer :: lo_ndx, hi_ndx real :: delsm1 lo_ndx = time_type%lo_buf_ndx if( time_type%t_interp ) then hi_ndx = time_type%hi_buf_ndx delsm1 = 1. - dels do k = 1,nlev do j = 1,nlat mozval(:,j,k) = time_type%conc(:,j,k,lo_ndx) * dels & + time_type%conc(:,j,k,hi_ndx) * delsm1 end do end do else do k = 1,nlev do j = 1,nlat mozval(:,j,k) = time_type%conc(:,j,k,lo_ndx) end do end do endif end subroutine tinterp_moz_species subroutine handle_error( message, status ) !--------------------------------------------------------------- ! handle errors produced by calling netCDF functions !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments : !--------------------------------------------------------------- integer, intent(in) :: status character(len=*), intent(in) :: message !--------------------------------------------------------------- ! output the error information from processing NETcdf file !--------------------------------------------------------------- write(*,*) ' ' write(*,*) trim(message) write(*,*) nf_strerror( status ) !--------------------------------------------------------------- ! exit from the bconLib !--------------------------------------------------------------- call exit_mozart_lib( flag=1 ) end subroutine handle_error subroutine exit_mozart_lib( flag ) !--------------------------------------------------------------- ! exit from module_mozart_lib !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, optional, intent(in) :: flag !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: status !--------------------------------------------------------------- ! deallocate !--------------------------------------------------------------- if( allocated(ix) ) deallocate( ix ) if( allocated(jy) ) deallocate( jy ) if( allocated(ax) ) deallocate( ax ) if( allocated(by) ) deallocate( by ) if( associated(time_type%conc) ) deallocate( time_type%conc ) if( associated(time_type%ps) ) deallocate( time_type%ps ) !----------------------------------------------------------------------- ! close netCDF file !----------------------------------------------------------------------- if( time_type%ncid_lo /= 0 ) then status = nf_close( time_type%ncid_lo ) if( time_type%ncid_hi /= 0 .and. & (time_type%ncid_hi /= time_type%ncid_lo) ) then status = nf_close( time_type%ncid_hi ) end if end if !----------------------------------------------------------------------- ! output information !----------------------------------------------------------------------- if( present(flag) ) then select case( flag ) case( 1 ); write(*,*) 'fail to process netCDF file...' case default; write(*,*) 'unknown error(s) occurred ...' endselect stop ' in module_mozart_lib ...' else write(*,*) 'successfully exited from module_mozart_lib ...' end if end subroutine exit_mozart_lib subroutine vinterp( p_moz, p_wrf, src, interp, nz, nlev ) !----------------------------------------------------------------------- ! ... vertically interpolate input data !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: nz, nlev real, intent(in) :: p_moz(nlev) real, intent(in) :: p_wrf(nz) real, intent(in) :: src(nlev) real, intent(out) :: interp(nz) !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: i integer :: k, kl, ku real :: delp, pinterp level_loop : & do k = 1,nz pinterp = p_wrf(k) if( pinterp <= p_moz(1) ) then interp(k) = src(1) else if( pinterp > p_moz(nlev) ) then interp(k) = src(nlev) else do ku = 2,nlev if( pinterp <= p_moz(ku) ) then kl = ku - 1 delp = log( pinterp/p_moz(kl) ) & / log( p_moz(ku)/p_moz(kl) ) interp(k) = src(kl) + delp * (src(ku) - src(kl)) exit end if end do end if end do level_loop end subroutine vinterp integer function lotim( cdate, csec, date, datesec, ntim ) !----------------------------------------------------------------------- ! ... return the index of the time sample that is the lower ! bound of the interval that contains the input date. if ! (cdate,csec) is earlier than the first time sample then 0 is ! returned. if (cdate,csec) is later than the last time sample then ! -index is returned. if (cdate,csec) is equal to the date of a ! dynamics time sample then that index is returned. !----------------------------------------------------------------------- use mo_calendar, only : diffdat implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: cdate ! date in yyyymmdd integer, intent(in) :: csec ! seconds relative to date integer, intent(in) :: ntim ! number times integer, intent(in) :: date(ntim) integer, intent(in) :: datesec(ntim) !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: n real :: dtime !----------------------------------------------------------------------- ! ... find latest date that is earlier than or equal to (date,sec) !----------------------------------------------------------------------- do n = 1,ntim dtime = diffdat( cdate, csec, date(n), datesec(n) ) if( dtime > 0. ) then lotim = n - 1 if( lotim > 0 ) then time_type%dels = dtime/diffdat( date(lotim), datesec(lotim), date(n), datesec(n) ) time_type%t_interp = time_type%dels /= 1. else time_type%dels = dtime endif exit endif end do if( n > ntim ) then if( dtime == 0. ) then lotim = ntim time_type%dels = 1. time_type%t_interp = .false. else lotim = -ntim time_type%dels = dtime endif endif end function lotim subroutine chk_moz_vars( nspec, filenm ) !--------------------------------------------------------------- ! check wrf to mozart variable mapping !--------------------------------------------------------------- use utils, only : wrf2mz_map !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- integer, intent(in) :: nspec character(len=*), intent(in) :: filenm !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- integer :: i, n integer :: vid integer :: status character(len=32) :: spcnam do n = 1,nspec write(*,*) 'checking wrf variable ',trim(wrf2mz_map(n)%wrf_name) do i = 1,wrf2mz_map(n)%moz_cnt spcnam = ' ' if( wrf2mz_map(n)%moz_ext(i) ) then spcnam = trim(wrf2mz_map(n)%moz_names(i)) // trim(moz_var_suffix) else spcnam = trim(wrf2mz_map(n)%moz_names(i)) end if write(*,*) 'len ',trim(wrf2mz_map(n)%moz_names(i)),' = ',len_trim(wrf2mz_map(n)%moz_names(i)) status = nf_inq_varid( ncid, trim(spcnam), vid ) if( status /= nf_noerr ) then write(*,*) 'chk_moz_vars: could not find ',spcnam,' in ',trim(filenm) call handle_error( ' ',status ) end if end do end do end subroutine chk_moz_vars end module module_mozart_lib