program gen_spectra use da_control use da_spectral use da_tracing implicit none integer, parameter :: iunit = 11 real, parameter :: tolerance = 1.0e-6 character(len=filename_len) :: filename character*10 :: date, cvar, ctime integer :: time, member ! Loop counters. integer :: i, j, k, n, m ! Loop counters. integer :: jts, jte ! start, end indices. integer :: nt ! Number of times to process. integer :: ne ! Number of ensemble members. integer :: ni, nj, nk ! Dimensions read in. integer :: max_wavenumber ! Smallest scale required (ni/2 - 1). integer :: index_k, index_start, index_end ! Array indexer. integer :: inc ! Jump between elements of vector in array. integer :: lenr ! FFT array dimension (at least inc*(n-1)+1). integer :: lensav ! wsave dimension (n+int(log(real(ni)))+4). integer :: lenwrk ! Dimension of work array. integer :: ier ! FFT error flag. integer :: count ! Counter integer :: index_m, index_j ! Array indices. integer :: alp_size ! Ass. Leg. Pol. array size. integer :: cv_size ! Complex control variable array size. logical :: first_time ! True if first time through loop. logical :: testing_spectral ! Test spectral transform. real :: pi_over_180 ! pi / 180 real :: diff_rms ! RMS error measure. real, allocatable :: ps_prime(:,:) ! Surface pressure perturbation. real, allocatable :: t_prime(:,:,:) ! Temperature perturbation. real, allocatable :: psi_prime(:,:,:) ! Streamfunction perturbation. real, allocatable :: chi_prime(:,:,:) ! Velocity Potential perturbation. real, allocatable :: rh_prime(:,:,:) ! Relative Humidity Perturbation. real, allocatable :: height(:,:,:) ! Geopotential height. real, allocatable :: latitude(:,:) ! Latitude (radians) real, allocatable :: field(:,:) ! Gridpoint field to be decomposed. real, allocatable :: field_out(:,:) ! Output test field. real, allocatable :: lat(:) ! Latitude (radians). real, allocatable :: sinlat(:) ! sine(latitude). real, allocatable :: coslat(:) ! cosine(latitude). real, allocatable :: int_wgts(:) ! Legendre integration weights. real, allocatable :: lon(:) ! Longitude (radians). real, allocatable :: sinlon(:) ! sine(longitude). real, allocatable :: coslon(:) ! cosine(longitude). real, allocatable :: wsave(:) ! Prime values for fast Fourier transform. real, allocatable :: alp(:) ! Associated Legendre Polynomial. real, allocatable :: power(:) ! Power spectrum (n). complex, allocatable:: cv(:) ! Control variable vector. complex, allocatable:: cv_out(:) ! Control variable vector copy for tests. if (trace_use) call da_trace_init if (trace_use) call da_trace_entry("gen_spectra") pi = 4.0 * atan(1.0) pi_over_180 = pi / 180.0 gaussian_lats = .false. testing_spectral = .true. nt = 1 ne = 1 !--------------------------------------------------------------------------------------------- ! gen_statistics2: Read standard fields, calculate regression coefficients, and create and ! output "control variables". !--------------------------------------------------------------------------------------------- ! Loop over times, ensemble members. !--------------------------------------------------------------------------------------------- write(6,'(a)')' [1] Read in perturbation fields.' !--------------------------------------------------------------------------------------------- open (iunit, file='/data2/hcshin/youn/DIFF63/GHMD_DIFF.2004030212_old', form='unformatted') read(iunit)date, ni, nj, nk write(6,'(a,a10)')' Date = ', date write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk allocate( ps_prime(1:ni,1:nj) ) allocate( t_prime(1:ni,1:nj,1:nk) ) allocate( psi_prime(1:ni,1:nj,1:nk) ) allocate( chi_prime(1:ni,1:nj,1:nk) ) allocate( rh_prime(1:ni,1:nj,1:nk) ) allocate( height(1:ni,1:nj,1:nk) ) allocate( latitude(1:ni,1:nj) ) read(iunit)ps_prime read(iunit)t_prime read(iunit)psi_prime read(iunit)chi_prime read(iunit)rh_prime read(iunit)height read(iunit)latitude ! Not currently used. filename = "psi"//"date"//"member"//"mode" open (iunit, file=filename, form='unformatted') write(iunit)ni, nj write(iunit)gaussian_lats write(iunit)psi_prime(1:ni,1:nj,1) close(iunit) ! Calculate regression coefficients. ! Calculate control variables. ! Output vp(i,j,k) = vp^-1 vs(i,j,k). !--------------------------------------------------------------------------------------------- ! gen_statistics3: Calculate Vertical Eigenvectors, Eigenvalues for a given variable. !--------------------------------------------------------------------------------------------- ! Loop-over/distribute control variables. ! Loop over times, ensemble members. ! Read control variables ! Calculate/accumulate Bv(j,k1,k2). ! End loop over times, ensemble members. ! Calculate E, L (v,j,k,k). ! Loop over times, ensemble members: ! Output v(i,j,m) = vv^-1 vp(i,j,k) !--------------------------------------------------------------------------------------------- ! gen_statistics4: Calculate Power spectra for given 2D field. !--------------------------------------------------------------------------------------------- ! Loop over/distribute variables, vertical modes: ! Loop over times, ensemble members: cvar = 'psi' first_time = .true. do time = 1, nt write(ctime,'(i4)')time do member = 1, ne ! write(cmember,'(i4)')member ! Read 2D fields. filename = trim(cvar)//"date"//"member"//"mode" open (iunit, file=filename, form='unformatted') read(iunit)ni, nj read(iunit)gaussian_lats if ( first_time ) then ! Initialize allocate( field(1:ni,1:nj) ) read(iunit)field !--------------------------------------------------------------------------------------------- write(6,'(a)')' Initialize spectral transforms.' !--------------------------------------------------------------------------------------------- inc = 1 lenr = inc * (ni - 1 ) + 1 lensav = ni + int(log(real(ni))) + 4 lenwrk = ni allocate( lat(1:nj) ) allocate( sinlat(1:nj) ) allocate( coslat(1:nj) ) allocate( int_wgts(1:nj) ) allocate( lon(1:ni) ) allocate( sinlon(1:ni) ) allocate( coslon(1:ni) ) allocate( wsave(1:lensav) ) max_wavenumber = ni / 2 - 1 allocate ( power(0:max_wavenumber) ) power(:) = 0.0 alp_size = ( nj + 1 ) * ( max_wavenumber + 1 ) * ( max_wavenumber + 2 ) / 4 allocate( alp( 1:alp_size) ) call da_initialize_h( ni, nj, max_wavenumber, lensav, alp_size, & wsave, lon, sinlon, coslon, lat, sinlat, coslat, & int_wgts, alp ) cv_size = ( max_wavenumber + 1 ) * ( max_wavenumber + 2 ) / 2 allocate( cv( 1:cv_size) ) ! Test horizontal transforms: if ( testing_spectral ) then call da_test_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, & alp_size, cv_size, alp, wsave, int_wgts, field ) end if first_time = .false. else ! Just read field read(iunit)field end if !--------------------------------------------------------------------------------------------- write(6,'(a)')' Perform gridpoint to spectral decomposition.' !--------------------------------------------------------------------------------------------- call da_vv_to_v_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, & alp_size, cv_size, alp, wsave, int_wgts, field, cv ) !--------------------------------------------------------------------------------------------- write(6,'(a)')' Calculate power spectra.' !--------------------------------------------------------------------------------------------- call da_calc_power( max_wavenumber, cv_size, cv, power ) do n = 0, max_wavenumber write(6,'(i4,1pe15.5)')n, power(n) end do end do ! End loop over ensemble members. end do ! End loop over times. if (trace_use) call da_trace_exit("gen_spectra") if (trace_use) call da_trace_report end program gen_spectra