module grid ! интерфейс ! --------------------------------------------------------------------------------- use const, only : miss_v, pi, r_earth, deg2rad, day2sec, len_default use datetime, only : date_type use netcdf use netcdf_kit, only : nc_errhand implicit none public private :: spatial_grid_init private :: temporal_grid_init private :: vertical_grid_init ! переменные ! --------------------------------------------------------------------------------- ! широтно-долготная сетка: integer :: nlon, nlat real :: dlon, dlat, lon_ref, lat_ref integer :: i0 !< Longitude of the Western edge of regional domain, number of grid integer :: i1 !< Longitude of the Eastern edge of regional domain, number of grid integer :: j0 !< Latitude of the Southern edge of regional domain, number of grid integer :: j1 !< Latitude of the Northern edge of regional domain, number of grid integer :: ni, nj integer :: ni_nc, nj_nc, i0_nc, i1_nc, j0_nc, j1_nc integer :: nlon_nc, nlat_nc real :: dlon_nc, dlat_nc, lon_ref_nc, lat_ref_nc real :: lon0_nc, lon1_nc, lat0_nc, lat1_nc integer :: ncell_tot, ncell_mask integer, allocatable :: id(:,:), id_to_i(:), id_to_j(:) integer, allocatable :: ich(:,:), ich_to_i(:), ich_to_j(:) real, allocatable :: lon(:), lat(:), area(:,:) ! SC1*SC(J) integer, allocatable :: mask(:,:) real, allocatable :: mask_nc(:,:) integer :: nid, nich real, parameter :: dlon_default = 0.5, dlat_default = 0.5 ! вертикальная сетка: integer, parameter :: ms = 4 !< номер уровня по высоте, соотвествующего поверхности почвы integer, parameter :: ml = 27 !< номер уровня по высоте, соотвествующего нижней границе почвы real, parameter :: z(ml) = (/-1.,-2/3.,-1/3.,0.,1.,2.,4.,8.,15.,25.,35.,45.,55.,65.,75.,85.,95.,105.,115.,125.,135.,145.,155.,200.,300.,500.,1000./) real :: dz(ml) ! сетка по времени: integer :: dt !< Timestep, sec integer :: ntime !< Number of timesteps integer :: nt_nc integer :: dt_nc real(8) :: jday0_nc character(len_default) :: DL character(len_default) :: DI type(date_type) :: date_c, date_fst, date_lst integer :: year_min, year_max ! текущие индексы в циклах: integer :: tt integer :: ii integer :: jj integer :: nn contains ! внешние процедуры ! --------------------------------------------------------------------------------- subroutine grid_init() ! --------------------------------------- call spatial_grid_init() call temporal_grid_init() call vertical_grid_init() end subroutine ! внутренние процедуры ! --------------------------------------------------------------------------------- subroutine spatial_grid_init() ! --------------------------------------- !< @brief широтно-долготная сетка use config, only : environment_data_type, & & lsm_datafile, & & nc_namespace, & & spatial_grid_mode, & & spatial_grid_res, & & spatial_sample_mode, & & id_user => id, & & ich_user => ich, & & point, & & polygon use paths, only : path_inmcm_data integer :: i, j, n integer :: ncid, dimid, varid, ios real, allocatable :: work1r(:), work2r(:,:) integer, allocatable :: work2i(:,:) real :: fill_v if (environment_data_type == 'lsm_offline') then call nc_errhand( nf90_open(trim(lsm_datafile), nf90_nowrite, ncid) ) call nc_errhand( nf90_inq_dimid(ncid, 'lon', dimid) ) call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = ni_nc) ) call nc_errhand( nf90_inq_varid(ncid, 'lon', varid) ) allocate(work1r(ni_nc)) call nc_errhand( nf90_get_var(ncid, varid, work1r(:), (/1/), (/ni_nc/)) ) if (ni_nc > 1) then dlon_nc = work1r(2) - work1r(1) else dlon_nc = dlon_default endif if (dlon_nc < 0.) stop 'check failed : dlon_nc < 0' lon0_nc = work1r(1) lon1_nc = work1r(ni_nc) deallocate(work1r) call nc_errhand( nf90_inq_dimid(ncid, 'lat', dimid) ) call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = nj_nc) ) call nc_errhand( nf90_inq_varid(ncid, 'lat', varid) ) allocate(work1r(nj_nc)) call nc_errhand( nf90_get_var(ncid, varid, work1r(:), (/1/), (/nj_nc/)) ) if (nj_nc > 1) then dlat_nc = work1r(2) - work1r(1) else dlat_nc = dlat_default endif if (dlat_nc < 0.) stop 'check failed : dlat_nc < 0' lat0_nc = work1r(1) lat1_nc = work1r(nj_nc) deallocate(work1r) nlon_nc = 360/dlon_nc nlat_nc = 180/dlat_nc lon_ref_nc = -180. + 0.5*dlon_nc lat_ref_nc = -90. + 0.5*dlat_nc i0_nc = nint((lon0_nc - lon_ref_nc)/dlon_nc) + 1 i1_nc = nint((lon1_nc - lon_ref_nc)/dlon_nc) + 1 j0_nc = nint((lat0_nc - lat_ref_nc)/dlat_nc) + 1 j1_nc = nint((lat1_nc - lat_ref_nc)/dlat_nc) + 1 allocate(mask_nc(i0_nc:i1_nc,j0_nc:j1_nc)) ios = nf90_inq_varid(ncid, nc_namespace%mask, varid) if (ios == 0) then call nc_errhand( nf90_get_var(ncid, varid, mask_nc(:,:), (/1,1/), (/ni_nc,nj_nc/)) ) else call nc_errhand( nf90_inq_varid(ncid, nc_namespace%ta, varid) ) call nc_errhand( nf90_get_att(ncid, varid, '_FillValue', fill_v) ) allocate(work2r(i0_nc:i1_nc,j0_nc:j1_nc)) call nc_errhand( nf90_get_var(ncid, varid, work2r(:,:), (/1,1/), (/ni_nc,nj_nc/)) ) do i = i0_nc, i1_nc do j = j0_nc, j1_nc if (work2r(i,j) == fill_v) then mask_nc(i,j) = 0 else mask_nc(i,j) = 1 endif enddo enddo endif call nc_errhand( nf90_close(ncid) ) else ni_nc = miss_v nj_nc = miss_v dlon_nc = miss_v dlat_nc = miss_v lon0_nc = miss_v lon1_nc = miss_v lat0_nc = miss_v lat1_nc = miss_v nlon_nc = miss_v nlat_nc = miss_v lon_ref_nc = miss_v lat_ref_nc = miss_v i0_nc = miss_v i1_nc = miss_v j0_nc = miss_v j1_nc = miss_v endif select case(spatial_grid_mode) case('none') dlon = miss_v dlat = miss_v case('auto') dlon = dlon_nc dlat = dlon_nc case('manual') dlon = spatial_grid_res(1) dlat = spatial_grid_res(2) if (environment_data_type == 'lsm_offline') then if (dlon /= dlon_nc .or. dlat /= dlat_nc) then print*, 'user grid :', dlon, dlat print*, 'netcdf grid :', dlon_nc, dlat_nc stop "check failed : interpolation is not supported yet" endif endif end select if (spatial_grid_mode /= 'none' .and. (dlon /= 0.5 .or. dlat /= 0.5)) then stop "check failed : now only 05x05 grid is supported (because of INMCM external data)" endif select case(spatial_grid_mode) case('none') nlon = 1 nlat = 1 lon_ref = miss_v lat_ref = miss_v case default nlon = 360/dlon nlat = 180/dlat lon_ref = -180. + 0.5*dlon lat_ref = -90. + 0.5*dlat end select allocate(id(nlon,nlat)) nid = nlon*nlat allocate(id_to_i(nid)) allocate(id_to_j(nid)) n = 0 do i = 1, nlon do j = nlat, 1, -1 ! порядок обхода с СЗ на юг (как в QGIS) n = n + 1 id(i,j) = n id_to_i(n) = i id_to_j(n) = j enddo enddo if (dlon == 0.5 .and. dlat == 0.5) then allocate(work2i(nlon+2,nlat+1)) open(1, file=path_inmcm_data//'OLIJAN', status='old') read(1,'(80i1)') ((work2i(i,j),i=1,nlon+2),j=1,nlat+1) close(1) allocate(ich(nlon,nlat)) nich = count(work2i(2:nlon+1,2:nlat+1) == 1) allocate(ich_to_i(nich)) allocate(ich_to_j(nich)) ich(:,:) = 0 n = 0 do j = 1, nlat do i = 1, nlon ! порядок обхода с ЮЗ на восток (как в INMCM) if (work2i(i+1,j+1) == 1) then n = n + 1 ich(i,j) = n ich_to_i(n) = i ich_to_j(n) = j endif enddo enddo deallocate(work2i) else allocate(ich, source = id) allocate(ich_to_i, source = id_to_i) allocate(ich_to_j, source = id_to_j) ich = id ich_to_i = id_to_i ich_to_j = id_to_j print*, 'warning : grid resolution is not suitable for OLIJAN' if (spatial_sample_mode == 'cell_ich') then stop "check failed : can not realize spatial_sample_mode = 'cell_ich'" endif endif select case(spatial_grid_mode) case('none') i0 = 1 i1 = 1 j0 = 1 j1 = 1 allocate(lon(i0:i1)) allocate(lat(j0:j1)) allocate(area(i0:i1,j0:j1)) lon = miss_v lat = miss_v area = miss_v case default select case(spatial_sample_mode) case('id') i0 = id_to_i(id_user) j0 = id_to_j(id_user) i1 = i0 j1 = j0 case('ich') i0 = ich_to_i(ich_user) j0 = ich_to_j(ich_user) i1 = i0 j1 = j0 case('point') i0 = nint((point(1) - lon_ref)/dlon) + 1 j0 = nint((point(2) - lat_ref)/dlat) + 1 i1 = i0 j1 = j0 case('polygon') i0 = nint((polygon(1) - lon_ref)/dlon) + 1 i1 = nint((polygon(2) - lon_ref)/dlon) + 1 j0 = nint((polygon(3) - lat_ref)/dlat) + 1 j1 = nint((polygon(4) - lat_ref)/dlat) + 1 case('all') i0 = i0_nc i1 = i1_nc j0 = j0_nc j1 = j1_nc end select allocate(lon(nlon)) do i = 1, nlon lon(i) = lon_ref + (i-1)*dlon enddo allocate(lat(nlat)) do j = 1, nlat lat(j) = lat_ref + (j-1)*dlat enddo allocate(area(nlon,nlat)) do j = 1, nlat do i = 1, nlon area(i,j) = r_earth*abs(dlat)*deg2rad * r_earth*cos(lat(j)*deg2rad)*abs(dlon)*deg2rad enddo enddo end select allocate(mask(i0:i1,j0:j1)) select case(environment_data_type) case('lsm_offline') mask(:,:) = mask_nc(i0:i1,j0:j1) case default mask(:,:) = 1. end select ni = i1-i0+1 nj = j1-j0+1 ncell_tot = ni*nj ncell_mask = count(mask == 1) end subroutine subroutine temporal_grid_init() ! --------------------------------------- !< @brief сетка по времени use config, only : environment_data_type, & & lsm_datafile, & & dt_mode, & & dt_user => dt, & & datetime_init_mode, & & timestamp_fst => datetime_init, & & ntime_mode, & & ntime_user => ntime, & & timestamp_lst => datetime_last, & & UTC_user => UTC, & & datetime_init_1, & & datetime_last_1, & & datetime_init_2, & & datetime_last_2, & & datetime_init_3, & & datetime_last_3, & & datetime_init_4, & & datetime_last_4, & & datetime_init_5, & & datetime_last_5, & & station_name integer :: ncid, dimid, varid real(8), allocatable :: work1r8(:) select case(environment_data_type) case('lsm_offline') call nc_errhand( nf90_open(trim(lsm_datafile), nf90_nowrite, ncid) ) call nc_errhand( nf90_inq_dimid(ncid, 'time', dimid) ) call nc_errhand( nf90_inquire_dimension(ncid, dimid, len = nt_nc) ) call nc_errhand( nf90_inq_varid(ncid, 'time', varid) ) allocate(work1r8(2)) call nc_errhand( nf90_get_var(ncid, varid, work1r8(:), (/1/), (/2/)) ) dt_nc = nint((work1r8(2) - work1r8(1)) * day2sec) jday0_nc = work1r8(1) deallocate(work1r8) call nc_errhand( nf90_close(ncid) ) case('station') select case(station_name) case('DAO4') DI = datetime_init_1 DL = datetime_last_1 case('DAO3') DI = datetime_init_2 DL = datetime_last_2 case('TRGK') DI = datetime_init_3 DL = datetime_last_3 case('VLDMR') DI = datetime_init_4 DL = datetime_last_4 case('Rostov') DI = datetime_init_5 DL = datetime_last_5 end select case default nt_nc = miss_v dt_nc = miss_v jday0_nc = miss_v end select select case(dt_mode) case('auto') dt = dt_nc case('manual') dt = dt_user if (environment_data_type == 'lsm_offline') then if (dt /= dt_nc) then print*, 'user dt :', dt print*, 'netcdf dt :', dt_nc stop "check failed : interpolation is not supported yet" endif endif end select select case(datetime_init_mode) case('auto') select case(environment_data_type) case('lsm_offline') call date_fst%init(jday = jday0_nc, UTC = UTC_user) case('station') call date_fst%init(timestamp = DI, UTC = UTC_user) case default stop "check failed : this option can not work with these settings" end select case('manual') call date_fst%init(timestamp = timestamp_fst, UTC = UTC_user) end select select case(ntime_mode) case('auto') select case(environment_data_type) case('lsm_offline') ntime = nt_nc call date_lst%init(jday = date_fst%jday + ntime*dble(dt)/day2sec) case('station') call date_lst%init(timestamp = DL, UTC = UTC_user) ntime = nint((date_lst%jday - date_fst%jday) * day2sec/dt) case default stop "check failed : this option can not work with these settings" end select case('ntime') ntime = ntime_user call date_lst%init(jday = date_fst%jday + ntime*dble(dt)/day2sec) case('datetime_last') call date_lst%init(timestamp = timestamp_lst) ntime = nint((date_lst%jday - date_fst%jday) * day2sec/dt) end select date_c = date_fst year_min = date_fst%y year_max = date_lst%y end subroutine subroutine vertical_grid_init() ! --------------------------------------- !< @brief вертикальная сетка integer :: k dz(:) = miss_v do k = ms+1, ml-1 dz(k) = 0.5*(z(k+1) - z(k-1)) enddo end subroutine end module