diff --git a/var/gen_be_v3/README.gen_be_v3 b/var/gen_be_v3/README.gen_be_v3 index 1d35f08654..9121c16331 100644 --- a/var/gen_be_v3/README.gen_be_v3 +++ b/var/gen_be_v3/README.gen_be_v3 @@ -1,11 +1,25 @@ *** gen_be_v3.F90 is a user-contribute program and is provided as it is. *** -*** gen_be_v3.F90 is for univariate processing only. *** *** gen_be_v3.F90 is for bin_type=5 (domain-average BE statistics) only. *** 1. To compile: - See the sample compile_cheyenne/compile_casper and the info in gen_be_v3.F90 + Option 1: use lapack lib compiled in WRFDA + (1) need to compile WRFDA first, + cd your_WRFDA_dir + ./clean -a + ./configure wrfda + ./compile all_wrfvar + (2) cd your_WRFDA_dir/var/gen_be_v3 + edit makefile + make + + Option 2: use general lapack lib (see the desciption in gen_be_v3.F90) + + Note about precision: + For write_ep = .true. application, the code should be compiled with real-4. + For gen_be applications, especially cv_options=5, real-4 works, but real-8 is recommended if memory is not an issue. + Hopefully the gen_be_v3 code will be updated in the future to better handle the variable declaration. 2. To run: @@ -49,6 +63,7 @@ be_method = 'NMC' varnames = 'T', 'U', 'V', 'RH', 'PSFC' ! for cv_options=7. ! can also be 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW', 'QGRAUP', 'W' + !varnames = 'PSI', 'CHI_U', 'T_U', 'RH', 'PSFC_U' ! for cv_options=5 do_pert_calc = .true. do_slen_calc = .true. slen_opt = 2 ! 1: curve-fitting method (extremely slow, need to run with OPENMP on) @@ -60,6 +75,7 @@ ! 2: read from pert1 file when need it ! This option writes out additional vertical-mode-projected fields when do_eof_transform=.true. ! Set pert1_read_opt=2 for large number of cases and large domain sizes that encounter memory insufficiency. + ! Note: when unbalanced variables are requested in varnames, pert1_read_opt=2 is used regardless of the namelist setting. / &ens_nml / @@ -68,11 +84,10 @@ Each variable is in separate output file. - For cv_options=7 purpose, util/combine_be_cv7.f90 can be used to combine be_[T/U/V/RH/PSFC].dat into be.dat that WRFDA can read. - For cv_options=7 and cloud_cv_options=2 purpose, util/combine_be_cv7_ccv2.f90 can be used to combine be_[T/U/V/RH/PSFC/QCLOUD/QRAIN/QICE/QSNOW/QGRAUP/W].dat into be.dat that WRFDA can read. - pert1.* files are intermediate output that can be removed after the successful execution. + For cv_options=5 purpose, util/combine_be_cv5.f90 can be used to combine be_[PSI/CHI_U/T_U/RH/PSFC_U].dat and regcoeff[1/2/3].dat into be.dat that WRFDA can read. + + pert* [psi.*] files are intermediate output that can be removed after the successful execution. - pertv_* files are intermediate pert1_read_opt=2 output that can be removed after the successful execution. diff --git a/var/gen_be_v3/compile_casper b/var/gen_be_v3/compile_casper deleted file mode 100755 index 1c698bfa9b..0000000000 --- a/var/gen_be_v3/compile_casper +++ /dev/null @@ -1,44 +0,0 @@ -#!/bin/csh -#set echo - -set MACHINE = `hostname -s` - -if ( `echo $MACHINE|cut -c1-6` == casper ) then #NCAR casper DAV - - # WRFIO_LIB = WRF/external/io_netcdf/libwrfio_nf.a - set WRFIO_LIB = /gpfs/u/home/hclin/extlib/intel/libwrfio_nf.a - if ( ! $?WRFIO_LIB ) then - echo "environment variable WRFIO_LIB not set" - exit 1 - endif - - # LAPACK library for EOF decomposition - set LAPACK_LIB = /glade/u/apps/opt/intel/2017u1/mkl/lib/intel64_lin/libmkl_lapack95_lp64.a - if ( ! $?LAPACK_LIB ) then - echo "environment variable LAPACK_LIB not set" - exit 1 - endif - - module purge - module load ncarenv - module load intel - module load ncarcompilers - module load netcdf - module load mkl - module load openmpi - module list - -set echo - - #MPI/OMP - cpp -P -traditional -DDM_PARALLEL gen_be_v3.F90 > be.f90 - mpif90 be.f90 -o gen_be_v3_mpp.exe -qopenmp -convert big_endian ${WRFIO_LIB} -L${NETCDF}/lib -lnetcdf -lnetcdff -lm -I${NETCDF}/include ${LAPACK_LIB} - - #serial/OMP - cpp -P -traditional gen_be_v3.F90 > be.f90 - ifort be.f90 -o gen_be_v3.exe -qopenmp -convert big_endian ${WRFIO_LIB} -L${NETCDF}/lib -lnetcdf -lnetcdff -lm -I${NETCDF}/include ${LAPACK_LIB} - -unset echo - - \rm -f be.f90 -endif diff --git a/var/gen_be_v3/compile_cheyenne b/var/gen_be_v3/compile_cheyenne deleted file mode 100755 index bed06ab65d..0000000000 --- a/var/gen_be_v3/compile_cheyenne +++ /dev/null @@ -1,44 +0,0 @@ -#!/bin/csh -#set echo - -set MACHINE = `hostname -s` - -if ( `echo $MACHINE|cut -c1-8` == cheyenne ) then #NCAR HPC - - # WRFIO_LIB = WRF/external/io_netcdf/libwrfio_nf.a - set WRFIO_LIB = /gpfs/u/home/hclin/extlib/intel/libwrfio_nf.a - if ( ! $?WRFIO_LIB ) then - echo "environment variable WRFIO_LIB not set" - exit 1 - endif - - # LAPACK library for EOF decomposition - set LAPACK_LIB = /glade/u/apps/opt/intel/2017u1/compilers_and_libraries/linux/mkl/lib/intel64_lin/libmkl_lapack95_lp64.a - if ( ! $?LAPACK_LIB ) then - echo "environment variable LAPACK_LIB not set" - exit 1 - endif - - module purge - module load ncarenv - module load intel - module load ncarcompilers - module load netcdf - module load mkl - module load mpt - module list - -set echo - - #MPI/OMP - cpp -P -traditional -DDM_PARALLEL gen_be_v3.F90 > be.f90 - mpif90 be.f90 -o gen_be_v3_mpp.exe -qopenmp -convert big_endian ${WRFIO_LIB} -L${NETCDF}/lib -lnetcdf -lnetcdff -lm -I${NETCDF}/include ${LAPACK_LIB} - - #serial/OMP - cpp -P -traditional gen_be_v3.F90 > be.f90 - ifort be.f90 -o gen_be_v3.exe -qopenmp -convert big_endian ${WRFIO_LIB} -L${NETCDF}/lib -lnetcdf -lnetcdff -lm -I${NETCDF}/include ${LAPACK_LIB} - -unset echo - - \rm -f be.f90 -endif diff --git a/var/gen_be_v3/gen_be_v3.F90 b/var/gen_be_v3/gen_be_v3.F90 index 95f0a7b83f..b81db79090 100644 --- a/var/gen_be_v3/gen_be_v3.F90 +++ b/var/gen_be_v3/gen_be_v3.F90 @@ -49,14 +49,16 @@ program gen_be_v3 real, parameter :: cp = 7.0*gas_constant/2.0 real, parameter :: kappa = gas_constant/cp real, parameter :: t_kelvin = 273.15 + real, parameter :: t_triple = 273.16 ! triple point of water real, parameter :: gas_constant_v = 461.6 real, parameter :: rd_over_rv = gas_constant/gas_constant_v real, parameter :: rd_over_rv1 = 1.0 - rd_over_rv namelist /gen_be_nml/ nc_list_file, be_method, varnames, outnames, & do_pert_calc, do_eof_transform, do_slen_calc, slen_opt, stride, yr_cutoff, & - verbose, level_start, level_end, aux_file, write_be_dat_r8, pert1_read_opt - namelist /ens_nml/ nens, write_ep, write_ens_stdv, ep_format, nproc_out + verbose, level_start, level_end, aux_file, write_be_dat_r8, pert1_read_opt, & + vertloc_opt, es_ice + namelist /ens_nml/ nens, write_ep, write_ens_stdv, ep_format, nproc_out, time_lag_ens character(len=filename_len) :: nc_list_file ! the text file that contains a list of netcdf files to process character (len=3) :: be_method ! "NMC" or "ENS" @@ -70,6 +72,7 @@ program gen_be_v3 integer :: stride ! for slen_opt=1 to skip every stride grid point real :: yr_cutoff(nvar_max) ! for slen_opt=1 logical :: verbose + logical :: es_ice ! whether to consider ice effect or not for saturation vapor pressure integer :: nens ! set in namelist for ENS method. hard-coded to 2 for NMC method. logical :: write_ep ! if writing out ensemble perturbations logical :: write_ens_stdv ! if writing out stdv of ensemble perturbations @@ -78,11 +81,14 @@ program gen_be_v3 ! 2: all members in one file (full-domain) (real*4) ! 3: all members in one file (sub-domain), need to specify nproc_out (real*4) integer :: nproc_out ! number of processors to decompose for ep_format=3 + logical :: time_lag_ens ! if calculating statistics for time lagged ensembles integer :: level_start, level_end ! for do_slen_calc=true. Can be set to be other than 1 and nvert for quick debugging/testing integer :: pert1_read_opt ! how to access pert1 data for compute_bv_sl ! 1: read and store all cases in memory at once ! 2: read from pert1 file when need it + integer :: vertloc_opt ! 0: no vertical localization + ! 1: with vertical localization using log-P algorithm logical :: write_pert0 logical :: write_pert1 logical :: write_be_dat_r8 ! if writing out be.dat in real*8 @@ -94,14 +100,31 @@ program gen_be_v3 integer :: nvar, nvar_avail integer :: mp_physics real :: ds + character(len=VarNameLen) :: vartmp + integer :: nvar_all - logical, allocatable :: read_it(:) + integer, allocatable :: istart_ens(:), iend_ens(:) + integer, allocatable :: istart_case(:), iend_case(:) + integer :: ens_istart, ens_iend + integer :: case_istart, case_iend + + logical, allocatable :: do_this_var(:) integer, allocatable :: var_dim(:) + real, allocatable :: vertloc_rho(:,:) + real, allocatable :: vertloc_kcutoff(:) + logical :: read_t logical :: read_qv logical :: read_prs + logical :: read_psfc logical :: read_ght = .false. + logical :: read_u + logical :: read_v + logical :: calc_psi + logical :: calc_chi + logical :: calc_regcoeff + integer :: indx_psi integer :: nk_2d = 1 integer :: sl_unit = 77 @@ -113,7 +136,7 @@ program gen_be_v3 integer, allocatable :: bin(:,:,:) integer, allocatable :: bin2d(:,:) - integer :: ii, i, j, n, k, ierr + integer :: ii, i, j, n, k, ierr, iv integer :: icase, icode integer :: istat logical :: isfile @@ -144,9 +167,15 @@ program gen_be_v3 real(r_double), allocatable :: r8tmp2d2(:,:) real(r_double), allocatable :: r8tmp3d(:,:,:) + ! for slen_opt=2 real, allocatable :: mapfac_x(:,:) real, allocatable :: mapfac_y(:,:) + ! for psi/chi calculation + real, allocatable :: mapfac_m(:,:) + real, allocatable :: mapfac_u(:,:) + real, allocatable :: mapfac_v(:,:) + real :: mean_scale real :: spike_tolerance = 1.5 integer :: n_smth_sl = 2 @@ -189,9 +218,12 @@ program gen_be_v3 ep_format = 2 nproc_out = 1 verbose = .false. + es_ice = .false. aux_file = 'none' write_be_dat_r8 = .true. pert1_read_opt = 1 + vertloc_opt = 0 + time_lag_ens = .false. ! initialize internal options write_pert0 = .false. @@ -218,7 +250,7 @@ program gen_be_v3 if ( myproc == root ) write(stdout,*) 'num_procs = ', num_procs if ( myproc == root ) write(stdout,*) 'nvar = ', nvar if ( nvar > 0 ) then - allocate (read_it(nvar)) + allocate (do_this_var(nvar)) allocate (var_dim(nvar)) else call error_handler(-1, 'varnames not set in namelist.gen_be_v3') @@ -231,6 +263,22 @@ program gen_be_v3 end if end do + ! convert varnames to be in upper case + do iv = 1, nvar + vartmp = varnames(iv) + do i = 1, len_trim(vartmp) + icode = ichar(vartmp(i:i)) + if (icode>=97 .and. icode<=122) then + vartmp(i:i) = char(icode - 97 + 65) + end if + end do + varnames(iv) = vartmp + ! allow ps_u in user namelist, but set it to psfc_u internally + if ( trim(varnames(iv)) == 'PS_U' ) then + varnames(iv) = 'PSFC_U' + end if + end do + ! convert be_method to be in upper case do i = 1, len_trim(be_method) icode = ichar(be_method(i:i)) @@ -242,6 +290,8 @@ program gen_be_v3 if ( trim(be_method) == 'NMC' ) then ! hard-coded nens=2 for NMC application nens = 2 + ! make sure time_lag_ens is false for NMC method + time_lag_ens = .false. else if ( trim(be_method) == 'ENS' ) then if ( nens == 0 ) then call error_handler(-1, & @@ -260,13 +310,23 @@ program gen_be_v3 ! as well as dimensions (ni, nj, nk), ds, and mp_physics option. call get_info(nc_list_file) - nvar_avail = count(read_it(:)) + nvar_avail = count(do_this_var(:)) if ( nvar_avail == 0 ) then call error_handler(-1, 'no valid variables found from the varnames in namelist.gen_be_v3') end if if ( myproc == root ) then write(stdout,*) 'ncase, nens = ', ncase, nens + if ( trim(be_method) == 'ENS' ) then + if ( num_procs > nens ) then + call error_handler(-1, 'num_procs must be smaller or equal to nens') + end if + end if + if ( trim(be_method) == 'NMC' ) then + if ( num_procs > ncase ) then + call error_handler(-1, 'num_procs must be smaller or equal to ncase') + end if + end if if ( verbose ) then if ( ncase > 0 .and. nens > 0 ) then do i = 1, ncase @@ -280,6 +340,47 @@ program gen_be_v3 if ( myproc == root ) write(stdout, '(a,3i5,f10.2)') ' ni, nj, nk, ds = ', ni, nj, nk, ds + ni1 = ni + 1 + nj1 = nj + 1 + + ! divide nens among available processors + allocate (istart_ens(0:num_procs-1)) + allocate (iend_ens (0:num_procs-1)) + do i = 0, num_procs - 1 + call para_range(1, nens, num_procs, i, istart_ens(i), iend_ens(i)) + end do + + ! divide ncase among available processors + allocate (istart_case(0:num_procs-1)) + allocate (iend_case (0:num_procs-1)) + do i = 0, num_procs - 1 + call para_range(1, ncase, num_procs, i, istart_case(i), iend_case(i)) + end do + + ! loop indices for distributing pert calulcation among processors + if ( trim(be_method) == 'ENS' ) then + case_istart = 1 + case_iend = ncase + ens_istart = istart_ens(myproc) + ens_iend = iend_ens(myproc) + else if ( trim(be_method) == 'NMC' ) then + case_istart = istart_case(myproc) + case_iend = iend_case(myproc) + ens_istart = 1 + ens_iend = 2 + end if + + if ( calc_psi .or. calc_chi ) then + allocate(mapfac_m(ni,nj)) + allocate(mapfac_u(ni1,nj)) + allocate(mapfac_v(ni,nj1)) + if ( trim(aux_file) == 'none' ) then + ! read MAPFAC_M, MAPFAC_U and MAPFAC_V from file filenames(1,1) + aux_file = filenames(1,1) + end if + call read_fixed_fields_mapfac(trim(aux_file)) + end if + if ( do_slen_calc .and. slen_opt==2 ) then allocate(mapfac_x(ni,nj)) allocate(mapfac_y(ni,nj)) @@ -290,13 +391,31 @@ program gen_be_v3 call read_fixed_fields(trim(aux_file)) end if + if ( vertloc_opt > 0 ) then + allocate(vertloc_rho(nk,nk)) + allocate(vertloc_kcutoff(nk)) + if ( trim(aux_file) == 'none' ) then + ! read vertical level info from file filenames(1,1) + aux_file = filenames(1,1) + end if + call get_vertloc(trim(aux_file), nk, vertloc_rho, vertloc_kcutoff) + if ( myproc == root ) then + write(stdout,*) + do k = 1,nk + write(stdout, '(a,i4,2x,f4.1)') ' vertloc: k, kcutoff = ', k, vertloc_kcutoff(k) + write(stdout, '(a,i4,100(f6.3))') ' vertloc: k, rho(k:nk) = ', k, vertloc_rho(k,k:nk) + end do + write(stdout,*) + end if + end if + if ( do_pert_calc ) then call compute_pert end if call ext_ncd_ioexit(ierr) - if ( do_slen_calc ) then + if ( calc_regcoeff .or. do_slen_calc ) then ! hard-code bin_type to be 5 ! set num_bins, num_bins2d, bin, bin2d with bin_type=5 values allocate(bin(ni,nj,nk)) @@ -307,7 +426,19 @@ program gen_be_v3 bin(:,:,k) = k end do bin2d(:,:) = 1 + end if + if ( calc_regcoeff ) then +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif + call compute_regcoeff_unbalanced + ! force pert1_read_opt=2 in order to read in unbalanced fields that are in separate files + pert1_read_opt = 2 + !if ( myproc == root ) write(stdout, '(a,3i5,f10.2)') ' --- Using pert1_read_opt=2 for unbalanced fields ---' + end if + + if ( do_slen_calc ) then allocate( be_data(nvar) ) do i = 1, nvar allocate( be_data(i)%e_vec(1:nk,1:nk) ) @@ -324,10 +455,16 @@ program gen_be_v3 end if if ( do_slen_calc ) then +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif call compute_bv_sl end if if ( do_slen_calc ) then +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif allocate(sl_smth(nk)) if ( write_be_dat_r8 ) then allocate(r8tmp1d(1:nk)) @@ -336,7 +473,7 @@ program gen_be_v3 allocate(r8tmp3d(1:nk,1:nk,1:num_bins2d)) end if do i = 1, nvar - if ( .not. read_it(i) ) cycle + if ( .not. do_this_var(i) ) cycle if ( myproc /= MOD((i-1), num_procs) ) cycle ! deal with zero scale_length @@ -455,11 +592,27 @@ program gen_be_v3 if ( allocated(mapfac_y) ) deallocate (mapfac_y) end if + if ( calc_psi .or. calc_chi ) then + if ( allocated(mapfac_m) ) deallocate (mapfac_m) + if ( allocated(mapfac_u) ) deallocate (mapfac_u) + if ( allocated(mapfac_v) ) deallocate (mapfac_v) + end if + if ( allocated(filenames) ) deallocate(filenames) if ( allocated(filedates) ) deallocate(filedates) - if ( allocated(read_it) ) deallocate(read_it) + if ( allocated(do_this_var) ) deallocate(do_this_var) if ( allocated(var_dim) ) deallocate(var_dim) + if ( vertloc_opt > 0 ) then + if ( allocated(vertloc_rho) ) deallocate(vertloc_rho) + if ( allocated(vertloc_kcutoff) ) deallocate(vertloc_kcutoff) + end if + + deallocate (istart_ens) + deallocate (iend_ens ) + deallocate (istart_case) + deallocate (iend_case ) + if ( myproc == root ) write(stdout,'(a)')' All Done!' #ifdef DM_PARALLEL @@ -524,7 +677,11 @@ subroutine get_info(nc_list_file) close(iunit) allocate (file_init_dates(nfile)) allocate (file_valid_dates(nfile)) - allocate (ilist_file(nens,nfile)) + if ( .not. time_lag_ens ) then + allocate (ilist_file(nens,nfile)) + else + allocate (ilist_file(nfile,nfile)) + end if allocate (valid(nfile)) file_init_dates(:) = '0000-00-00_00:00:00' file_valid_dates(:) = '0000-00-00_00:00:00' @@ -579,30 +736,75 @@ subroutine get_info(nc_list_file) read_t = .false. read_qv = .false. read_prs = .false. + read_psfc= .false. + read_u = .false. + read_v = .false. + calc_psi = .false. + calc_chi = .false. + calc_regcoeff = .false. + indx_psi = -1 do i = 1, nvar call ext_ncd_get_var_info (fid, trim(varnames(i)), ndim, ordering, staggering, & start_index, end_index, wrftype, ierr) var_dim(i) = ndim - if ( ierr == 0 ) then - read_it(i) = .true. + if ( ierr == 0 ) then ! direct variables + do_this_var(i) = .true. if ( trim(varnames(i)) == 'T' ) then read_prs = .true. end if - else + else ! possible derived variables if ( trim(varnames(i)) == 'RH' ) then + do_this_var(i) = .true. read_t = .true. read_qv = .true. read_prs = .true. - read_it(i) = .true. var_dim(i) = 3 + else if ( trim(varnames(i)) == 'PSI' ) then + do_this_var(i) = .true. + read_u = .true. + read_v = .true. + var_dim(i) = 3 + calc_psi = .true. + indx_psi = i + else if ( trim(varnames(i)) == 'CHI' ) then + do_this_var(i) = .true. + read_u = .true. + read_v = .true. + var_dim(i) = 3 + calc_chi = .true. + else if ( trim(varnames(i)) == 'CHI_U' ) then + do_this_var(i) = .true. + read_u = .true. + read_v = .true. + var_dim(i) = 3 + calc_psi = .true. + calc_chi = .true. + calc_regcoeff = .true. + else if ( trim(varnames(i)) == 'T_U' ) then + do_this_var(i) = .true. + read_prs = .true. + read_t = .true. + read_u = .true. + read_v = .true. + var_dim(i) = 3 + calc_psi = .true. + calc_regcoeff = .true. + else if ( trim(varnames(i)) == 'PSFC_U' ) then + do_this_var(i) = .true. + read_psfc = .true. + read_u = .true. + read_v = .true. + var_dim(i) = 2 + calc_psi = .true. + calc_regcoeff = .true. else - read_it(i) = .false. + do_this_var(i) = .false. if ( myproc == root ) write(stdout,'(3a)') & - ' Warning: ', trim(varnames(i)), ' not available in the input file, skipping it' + ' Warning: ', trim(varnames(i)), ' not available in the input file or not implemented, skipping it' end if - end if - end do - end if + end if ! ext_ncd_get_var_info check + end do ! nvar loop + end if ! one-time check call ext_ncd_ioclose(fid, ierr) @@ -631,26 +833,43 @@ subroutine get_info(nc_list_file) if ( icount == 1 ) then iens = 1 icase = 1 + ilist_file(iens,icase) = ifile else if ( icount > 1 ) then - find_match_loop_ens: do ii = ifile-1, 1, -1 - if ( valid(ii) ) then - if ( file_init_dates(ifile) == file_init_dates(ii) ) then - if ( file_valid_dates(ifile) == file_valid_dates(ii) ) then - iens = iens + 1 - ilist_file(iens,icase) = ifile - if ( ii == 1 ) ilist_file(iens-1,icase) = ii + if ( .not. time_lag_ens ) then + find_match_loop_ens: do ii = ifile-1, 1, -1 + if ( valid(ii) ) then + if ( file_init_dates(ifile) == file_init_dates(ii) ) then + if ( file_valid_dates(ifile) == file_valid_dates(ii) ) then + iens = iens + 1 + ilist_file(iens,icase) = ifile + exit find_match_loop_ens + end if + else + icase = icase + 1 + iens = 1 + ilist_file(iens,icase) = ifile exit find_match_loop_ens - end if - else - icase = icase + 1 - end if + end if ! if same init_date + end if ! valid file + end do find_match_loop_ens + else ! time_lag_ens + ! only check_valid_date for time_lag_ens + if ( file_valid_dates(ifile) == file_valid_dates(ilist_file(iens,icase)) ) then + iens = iens + 1 + ilist_file(iens,icase) = ifile end if - end do find_match_loop_ens + end if ! time_lag_ens end if end if end do file_loop1 ncase = icase + if ( time_lag_ens ) then + nens = iens + if ( myproc == root ) then + write(stdout,'(a,i4,a)') ' --- Resetting nens = ', nens, ' for time_lag_ens ---' + end if + end if if ( ncase > 0 .and. nens > 0 ) then if ( .not. allocated(filenames) ) allocate(filenames(1:nens,1:ncase)) @@ -753,6 +972,239 @@ subroutine read_fixed_fields(input_file) end subroutine read_fixed_fields +subroutine read_fixed_fields_mapfac(input_file) + + implicit none + + integer :: fid, ierr + character(len=*),intent(in) :: input_file + character(len=DateStrLen) :: DateStr + real(r_single), allocatable :: xfield2d(:,:) + real(r_single), allocatable :: xfield2d_u(:,:) + real(r_single), allocatable :: xfield2d_v(:,:) + character(len=80), dimension(3) :: dimnames + character(len=4) :: staggering = ' N/A' !dummy + character(len=3) :: ordering + integer, dimension(4) :: start_index, end_index + integer :: ndim, wrftype + character(len=512) :: message + character(len=VarNameLen) :: varname + + !input_file = trim(filenames(1,1)) + + call ext_ncd_open_for_read(trim(input_file), 0, 0, "", fid, ierr) + write(message,'(a,i3,a,a)') ' Proc ', myproc, ' opening ', trim(input_file) + call error_handler(ierr, trim(message)) + + call ext_ncd_get_next_time(fid, DateStr, ierr) + call error_handler(ierr, 'ext_ncd_get_next_time') + + if ( .not. allocated(xfield2d) ) allocate (xfield2d(ni, nj)) + + if ( .not. allocated(mapfac_m) ) allocate (mapfac_m(ni, nj)) + varname = 'MAPFAC_M' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield2d(:,:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + mapfac_m(:,:) = xfield2d(:,:) + deallocate (xfield2d) + + if ( .not. allocated(xfield2d_u) ) allocate (xfield2d_u(ni1, nj)) + + if ( .not. allocated(mapfac_u) ) allocate (mapfac_u(ni1, nj)) + varname = 'MAPFAC_U' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield2d_u(:,:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + mapfac_u(:,:) = xfield2d_u(:,:) + deallocate (xfield2d_u) + + if ( .not. allocated(xfield2d_v) ) allocate (xfield2d_v(ni, nj1)) + + if ( .not. allocated(mapfac_v) ) allocate (mapfac_v(ni, nj1)) + varname = 'MAPFAC_V' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield2d_v(:,:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + mapfac_v(:,:) = xfield2d_v(:,:) + deallocate (xfield2d_v) + + call ext_ncd_ioclose(fid, ierr) + +end subroutine read_fixed_fields_mapfac + +subroutine get_vertloc(input_file, kz, rho, kcutoff) + + implicit none + + character(len=*),intent(in) :: input_file + integer, intent(in) :: kz + real, intent(inout) :: rho(kz,kz) + real, intent(inout) :: kcutoff(kz) + + integer :: fid, ierr + character(len=DateStrLen) :: DateStr + character(len=80), dimension(3) :: dimnames + character(len=4) :: staggering = ' N/A' !dummy + character(len=3) :: ordering + integer, dimension(4) :: start_index, end_index + integer :: ndim, wrftype + character(len=512) :: message + character(len=VarNameLen) :: varname + real(r_single), allocatable :: xfield(:) + + integer :: k, k1, i + real :: kscale, kscale_invsq, kdist + real :: cutoff + real :: p_top, base_pres + real, allocatable :: znw(:) + real, allocatable :: p_w(:) ! pressure at full levels + real, allocatable :: ln_p_w(:) + real, allocatable :: dlnp(:,:) + real, allocatable :: kcutoff_tmp(:) + logical :: smooth_kcutoff + + call ext_ncd_open_for_read(trim(input_file), 0, 0, "", fid, ierr) + write(message,'(a,i3,a,a)') ' Proc ', myproc, ' opening ', trim(input_file) + call error_handler(ierr, trim(message)) + + call ext_ncd_get_next_time(fid, DateStr, ierr) + call error_handler(ierr, 'ext_ncd_get_next_time') + + if ( .not. allocated(xfield) ) allocate (xfield(kz+1)) + varname = 'P_TOP' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield(1), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + p_top = xfield(1) + + varname = 'P00' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield(1), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + base_pres = xfield(1) + + if ( .not. allocated(znw) ) allocate(znw(kz+1)) + varname = 'ZNW' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, trim(varname), & + xfield(:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + znw(:) = xfield(:) + + call ext_ncd_ioclose(fid, ierr) + if ( allocated(xfield) ) deallocate (xfield) + + allocate(p_w(1:kz+1)) + allocate(ln_p_w(1:kz+1)) + allocate(dlnp(1:kz,1:kz)) + + ! empirical settings + cutoff = 1.0/2.71828 !1/e + + do k = 1, kz+1 + p_w(k) = znw(k)*(base_pres - p_top) + p_top + ln_p_w(k) = log(max(p_w(k),0.0001)) + end do + kcutoff(:) = 1.0 ! initialize + do k = 1, kz + do k1 = k+1, kz + dlnp(k,k1) = abs(ln_p_w(k)-ln_p_w(k1)) + if ( dlnp(k,k1) <= cutoff ) then + kcutoff(k) = k1-k+1 + !if ( 0.5*(p_w(k)+p_w(k+1)) <= 15000.0 ) then ! 15000.0Pa + ! kcutoff(k) = min(kcutoff(k), 1.0) + !end if + end if + end do + end do + + ! smoothing + !smooth_kcutoff = .false. + smooth_kcutoff = .true. + if ( smooth_kcutoff ) then + allocate(kcutoff_tmp(1:kz)) + kcutoff_tmp(:) = kcutoff(:) + do i = 1, 2 ! two passes + do k = 2, kz-1 + kcutoff_tmp(k) = kcutoff(k)+ 0.25 * ( kcutoff(k-1) + kcutoff(k+1) -2.0*kcutoff(k) ) + end do + kcutoff(:) = kcutoff_tmp(:) + end do + deallocate(kcutoff_tmp) + end if + + deallocate(dlnp) + deallocate(ln_p_w) + deallocate(p_w) + + ! specify probability densities: + rho(:,:) = 1.0 ! initialize + do k = 1, kz + kscale = kcutoff(k) + kscale_invsq = 1.0 / ( kscale * kscale ) + do k1 = k, kz + kdist = k1 - k + rho(k,k1) = exp ( -1.0 * real(kdist * kdist) * kscale_invsq ) + rho(k1,k) = rho(k,k1) + end do + end do + +end subroutine get_vertloc + subroutine compute_pert implicit none @@ -775,21 +1227,25 @@ subroutine compute_pert character(len=DateStrLen) :: DateStr character(len=512) :: message - integer, allocatable :: istart_ens(:), iend_ens(:) - integer, allocatable :: istart_case(:), iend_case(:) - integer, allocatable :: istart_var(:), iend_var(:) integer, allocatable :: iproc_var(:) - integer :: ens_istart, ens_iend - integer :: case_istart, case_iend integer :: ivs, ive, ie_indx integer :: my_nvar, dest_proc + integer :: this_mpi_real real(r_single), allocatable :: xfield(:,:,:) real(r_single), allocatable :: xfield_u(:,:,:) real(r_single), allocatable :: xfield_v(:,:,:) real(r_single), allocatable :: xfield_w(:,:,:) + real, allocatable :: vor(:,:,:) ! vorticity + real, allocatable :: div(:,:,:) ! divergence + real, allocatable :: psi_s(:,:,:) ! stream function + real, allocatable :: psi(:,:,:) ! stream function on mass grid + real, allocatable :: chi(:,:,:) ! velocity potential + real, allocatable :: ustag(:,:,:) ! u on staggered grid + real, allocatable :: vstag(:,:,:) ! v on staggered grid + real, allocatable :: psfc(:,:) ! surface pressure real, allocatable :: prs(:,:,:) ! prs on half levels real, allocatable :: prs_w(:,:,:) ! prs on full levels real, allocatable :: ght(:,:,:) @@ -799,8 +1255,8 @@ subroutine compute_pert real, allocatable :: mub(:,:) real, allocatable :: znw(:) real :: p_top - real :: tc, es, qs - logical :: got_prs, got_th, got_qv + real :: tk, es, qs + logical :: got_prs, got_th, got_qv, got_u, got_v, got_psfc real, allocatable :: time_mean(:,:,:,:) real(r_double), allocatable :: xsum_loc(:,:,:,:), xsum(:,:,:,:) @@ -835,45 +1291,28 @@ subroutine compute_pert nk1 = nk + 1 ijk = ni * nj * nk - ! divide nens among available processors - allocate (istart_ens(0:num_procs-1)) - allocate (iend_ens (0:num_procs-1)) - do i = 0, num_procs - 1 - call para_range(1, nens, num_procs, i, istart_ens(i), iend_ens(i)) - end do - - ! divide ncase among available processors - allocate (istart_case(0:num_procs-1)) - allocate (iend_case (0:num_procs-1)) - do i = 0, num_procs - 1 - call para_range(1, ncase, num_procs, i, istart_case(i), iend_case(i)) - end do - - if ( trim(be_method) == 'ENS' ) then - case_istart = 1 - case_iend = ncase - ens_istart = istart_ens(myproc) - ens_iend = iend_ens(myproc) - else if ( trim(be_method) == 'NMC' ) then - case_istart = istart_case(myproc) - case_iend = iend_case(myproc) - ens_istart = 1 - ens_iend = nens - end if - write(stdout,'(a,i4,a,i4,a,i4)') & ' Proc ', myproc, ' will read ens files ', ens_istart, ' - ', ens_iend write(stdout,'(a,i4,a,i4,a,i4)') & ' Proc ', myproc, ' will read case files ', case_istart, ' - ', case_iend if ( ens_iend >= ens_istart .and. case_iend >= case_istart ) then - allocate(xdata(nvar,ens_istart:ens_iend,case_istart:case_iend), stat=ierr) + if ( calc_psi .and. indx_psi < 0 ) then + ! when psi is not in the var list but is needed for unbalanced fields + ! allocate index nvar+1 for psi + allocate(xdata(nvar+1,ens_istart:ens_iend,case_istart:case_iend), stat=ierr) + indx_psi = nvar + 1 + nvar_all = nvar + 1 + else + allocate(xdata(nvar,ens_istart:ens_iend,case_istart:case_iend), stat=ierr) + nvar_all = nvar + end if call error_handler(ierr, 'allocating xdata in compute_pert') end if do ic = case_istart, case_iend do ie = ens_istart, ens_iend - do iv = 1, nvar + do iv = 1, nvar_all allocate(xdata(iv,ie,ic)%value(ni,nj,nk), stat=ierr) call error_handler(ierr, 'allocating xdata%value in compute_pert') xdata(iv,ie,ic)%value(:,:,:) = 0.0 @@ -892,6 +1331,9 @@ subroutine compute_pert got_qv = .false. got_th = .false. + got_u = .false. + got_v = .false. + got_psfc = .false. call ext_ncd_open_for_read(trim(input_file), 0, 0, "", fid, ierr) write(message,'(a,i3,a,a)') ' Proc ', myproc, ' opening ', trim(input_file) @@ -941,6 +1383,26 @@ subroutine compute_pert end if ! p, pb end if ! read_prs + if ( read_psfc ) then + if ( .not. allocated(xfield) ) allocate (xfield (ni, nj, nk)) + varname = 'PSFC' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, varname, & + xfield(:,:,1), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + if ( .not. allocated(psfc) ) allocate (psfc(ni, nj)) + psfc(:,:) = xfield(:,:,1) + got_psfc = .true. + end if + if ( read_ght ) then if ( .not. allocated(ght) ) allocate (ght(ni, nj, nk)) if ( .not. allocated(xfield_w) ) allocate (xfield_w(ni, nj ,nk1)) @@ -1086,7 +1548,7 @@ subroutine compute_pert start_index, end_index, wrftype, ierr) call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) call ext_ncd_read_field(fid, DateStr, trim(varname), & - xfield_w(1,1,:), wrftype, & + xfield_w(1,1,:), wrftype, & 0, 0, 0, ordering, & staggering, dimnames, & !dummy start_index, end_index, & !dom @@ -1116,19 +1578,93 @@ subroutine compute_pert do k = nk, 1, -1 do j = 1, nj do i = 1, ni - prs(i,j,k) = 0.5*(prs_w(i,j,k)+prs_w(i,j,k+1)) + prs(i,j,k) = 0.5*(prs_w(i,j,k)+prs_w(i,j,k+1)) + end do + end do + end do + deallocate (prs_w) + deallocate (mu) + deallocate (mub) + deallocate (znw) + end if ! read_prs, got_prs + + if ( read_u ) then + if ( .not. allocated(xfield_u) ) allocate (xfield_u(ni1,nj, nk)) + if ( .not. allocated(ustag) ) allocate (ustag(ni1, nj, nk)) + varname = 'U' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, varname, & + xfield_u(:,:,:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + ustag(:,:,:) = xfield_u(:,:,:) + got_u = .true. + end if + + if ( read_v ) then + if ( .not. allocated(xfield_v) ) allocate (xfield_v(ni, nj1,nk)) + if ( .not. allocated(vstag) ) allocate (vstag(ni, nj1, nk)) + varname = 'V' + call ext_ncd_get_var_info (fid, trim(varname), ndim, ordering, staggering, & + start_index, end_index, wrftype, ierr) + call error_handler(ierr, 'ext_ncd_get_var_info '//trim(varname)) + call ext_ncd_read_field(fid, DateStr, varname, & + xfield_v(:,:,:), wrftype, & + 0, 0, 0, ordering, & + staggering, dimnames, & !dummy + start_index, end_index, & !dom + start_index, end_index, & !mem + start_index, end_index, & !pat + ierr ) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + vstag(:,:,:) = xfield_v(:,:,:) + got_v = .true. + end if + + if ( calc_psi .and. got_u .and. got_v ) then + if ( .not. allocated(vor) ) allocate (vor(ni+1, nj+1, nk)) + if ( .not. allocated(psi_s) ) allocate (psi_s(ni+1, nj+1, nk)) + if ( .not. allocated(psi) ) allocate (psi(ni, nj, nk)) + ! Calculate vorticity (in center of mass grid on WRF's Arakawa C-grid) + call da_uv_to_vor_c(ni, nj, nk, ds, & + mapfac_m, mapfac_u, mapfac_v, ustag, vstag, vor) + ! Calculate streamfunction + ! Assumes vor converted to Del**2 psi + call da_del2a_to_a(ni+1, nj+1, nk, ds, vor, psi_s) + + ! interpolate psi to mass points + do k = 1, nk + do j = 1, nj + do i = 1, ni + psi(i,j,k) = 0.25 * ( psi_s(i,j,k) + psi_s(i+1,j,k) + & + psi_s(i,j+1,k) + psi_s(i+1,j+1,k) ) end do end do end do - deallocate (prs_w) - deallocate (mu) - deallocate (mub) - deallocate (znw) - end if ! read_prs, got_prs + if ( allocated(psi_s) ) deallocate (psi_s) + end if ! calc_psi + + if ( calc_chi .and. got_u .and. got_v ) then + if ( .not. allocated(div) ) allocate (div(ni, nj, nk)) + if ( .not. allocated(chi) ) allocate (chi(ni, nj, nk)) + ! Calculate divergence (at mass pts. on WRF's Arakawa C-grid) + call da_uv_to_div_c(ni, nj, nk, ds, & + mapfac_m, mapfac_u, mapfac_v, ustag, vstag, div) + ! Calculate velocity potential + ! Assumes div converted to Del**2 chi + call da_del2a_to_a(ni, nj, nk, ds, div, chi) + end if ! calc_chi var_loop1: do iv = 1, nvar - if ( .not. read_it(iv) ) cycle var_loop1 + if ( .not. do_this_var(iv) ) cycle var_loop1 varname = trim(varnames(iv)) call ext_ncd_get_var_info (fid, varname, ndim, ordering, staggering, & @@ -1138,9 +1674,10 @@ subroutine compute_pert !write(stdout, '(a,i3,a,a8,a,a)') ' Proc ', myproc, ' Reading ', & ! trim(varname), ' from ', trim(input_file) - if ( varname == 'PSFC' ) then + if ( varname == 'PSFC' .or. varname == 'PSFC_U' ) then if ( .not. allocated(xfield) ) allocate (xfield (ni, nj, nk)) - call ext_ncd_read_field(fid, DateStr, varname, & + if ( .not. got_psfc ) then + call ext_ncd_read_field(fid, DateStr, varname, & xfield(:,:,1), wrftype, & 0, 0, 0, ordering, & staggering, dimnames, & !dummy @@ -1148,11 +1685,15 @@ subroutine compute_pert start_index, end_index, & !mem start_index, end_index, & !pat ierr ) - call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + else + xfield(:,:,1) = psfc(:,:) + end if xdata(iv,ie,ic)%value(:,:,1) = xfield(:,:,1) else if ( varname == 'U' ) then if ( .not. allocated(xfield_u) ) allocate (xfield_u(ni1,nj, nk)) - call ext_ncd_read_field(fid, DateStr, varname, & + if ( .not. got_u ) then + call ext_ncd_read_field(fid, DateStr, varname, & xfield_u(:,:,:), wrftype, & 0, 0, 0, ordering, & staggering, dimnames, & !dummy @@ -1160,7 +1701,10 @@ subroutine compute_pert start_index, end_index, & !mem start_index, end_index, & !pat ierr ) - call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + else + xfield_u(:,:,:) = ustag(:,:,:) + end if do k = 1, nk do j = 1, nj do i = 1, ni @@ -1171,7 +1715,8 @@ subroutine compute_pert end do else if ( varname == 'V' ) then if ( .not. allocated(xfield_v) ) allocate (xfield_v(ni, nj1,nk)) - call ext_ncd_read_field(fid, DateStr, varname, & + if ( .not. got_v ) then + call ext_ncd_read_field(fid, DateStr, varname, & xfield_v(:,:,:), wrftype, & 0, 0, 0, ordering, & staggering, dimnames, & !dummy @@ -1179,7 +1724,10 @@ subroutine compute_pert start_index, end_index, & !mem start_index, end_index, & !pat ierr ) - call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + call error_handler(ierr, 'ext_ncd_read_field '//trim(varname)) + else + xfield_v(:,:,:) = vstag(:,:,:) + end if do k = 1, nk do j = 1, nj do i = 1, ni @@ -1211,9 +1759,11 @@ subroutine compute_pert do k = 1, nk do j = 1, nj do i = 1, ni - tc = (theta(i,j,k)*(prs(i,j,k)/p00)**kappa) - t_kelvin - es = 611.2 * exp( 17.67 * tc / (tc + 243.5) ) - qs = rd_over_rv * es /( prs(i,j,k) - rd_over_rv1 * es) + !tc = (theta(i,j,k)*(prs(i,j,k)/p00)**kappa) - t_kelvin + !es = 611.2 * exp( 17.67 * tc / (tc + 243.5) ) + !qs = rd_over_rv * es /( prs(i,j,k) - rd_over_rv1 * es) + tk = theta(i,j,k)*(prs(i,j,k)/p00)**kappa + call calc_es_qs(tk, prs(i,j,k), es, qs, es_ice) ! get RH is ratio, not percentage xdata(iv,ie,ic)%value(i,j,k) = (q(i,j,k)/(1.0+q(i,j,k))) / qs end do @@ -1238,7 +1788,7 @@ subroutine compute_pert xdata(iv,ie,ic)%value(:,:,:) = & xfield(:,:,:) / ( 1.0 + xfield(:,:,:) ) end if ! got_qv - else if ( varname == 'T' ) then + else if ( varname == 'T' .or. varname == 'T_U' ) then if ( got_th ) then ! potential temperature to temperature xdata(iv,ie,ic)%value(:,:,:) = & @@ -1258,6 +1808,10 @@ subroutine compute_pert xdata(iv,ie,ic)%value(:,:,:) = & (t00+xfield(:,:,:))*(prs(:,:,:)/p00)**kappa end if ! got_th + else if ( varname == 'PSI' ) then + xdata(iv,ie,ic)%value(:,:,:) = psi(:,:,:) + else if ( varname == 'CHI' .or. varname == 'CHI_U' ) then + xdata(iv,ie,ic)%value(:,:,:) = chi(:,:,:) else if ( .not. allocated(xfield) ) allocate (xfield (ni, nj, nk)) call ext_ncd_read_field(fid, DateStr, varname, & @@ -1279,6 +1833,10 @@ subroutine compute_pert end do var_loop1 ! nvar loop + if ( indx_psi > nvar ) then + xdata(indx_psi,ie,ic)%value(:,:,:) = psi(:,:,:) + end if + call ext_ncd_ioclose(fid, ierr) end do ens_loop1 ! ensemble member loop @@ -1291,13 +1849,13 @@ subroutine compute_pert if ( remove_time_mean ) then ! allocate for summing up the cases each processor sees if ( .not. allocated(xsum_loc) ) then - allocate(xsum_loc(ni,nj,nk,nvar), stat=ierr) + allocate(xsum_loc(ni,nj,nk,nvar_all), stat=ierr) call error_handler(ierr, 'allocating xsum_loc') xsum_loc = 0.0 end if end if - do iv = 1, nvar + do iv = 1, nvar_all ! store forecast difference in index 1 of ie xdata(iv,1,ic)%value(:,:,:) = xdata(iv,1,ic)%value(:,:,:) - xdata(iv,2,ic)%value(:,:,:) if ( remove_time_mean ) then @@ -1329,12 +1887,14 @@ subroutine compute_pert if ( .not. allocated(ens_sum_loc) ) allocate(ens_sum_loc(ni,nj,nk)) if ( .not. allocated(ens_sum) ) allocate(ens_sum(ni,nj,nk)) - if ( .not. allocated(ens_mean) ) allocate(ens_mean(ni,nj,nk,nvar)) + if ( .not. allocated(ens_mean) ) allocate(ens_mean(ni,nj,nk,nvar_all)) if ( myproc == root ) write(stdout,'(a)') ' ====== Computing ensemble mean ======' - var_loop2: do iv = 1, nvar - if ( .not. read_it(iv) ) cycle var_loop2 + var_loop2: do iv = 1, nvar_all + if ( iv <= nvar ) then + if ( .not. do_this_var(iv) ) cycle var_loop2 + end if ens_sum_loc(:,:,:) = 0.0 ! initialize for each variable ens_sum (:,:,:) = 0.0 ! initialize for each variable do n = istart_ens(myproc),iend_ens(myproc) @@ -1365,8 +1925,10 @@ subroutine compute_pert end if if ( myproc == root ) write(stdout,'(a)') ' ====== Computing ensemble perturbations ======' - do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + do iv = 1, nvar_all + if ( iv <= nvar ) then + if ( .not. do_this_var(iv) ) cycle + end if do ie = ens_istart, ens_iend ! each proc loops over a subset of ens if ( remove_ens_mean ) then xdata(iv,ie,ic)%value(:,:,:) = xdata(iv,ie,ic)%value(:,:,:) - ens_mean(:,:,:,iv) @@ -1378,7 +1940,7 @@ subroutine compute_pert if ( myproc == root ) write(stdout,'(a)') ' ====== Computing ensemble stdv ======' allocate(ens_stdv(ni,nj,nk,nvar)) do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + if ( .not. do_this_var(iv) ) cycle ens_sum_loc(:,:,:) = 0.0 ! initialize for each variable ens_sum (:,:,:) = 0.0 ! initialize for each variable do n = istart_ens(myproc),iend_ens(myproc) @@ -1401,6 +1963,21 @@ subroutine compute_pert deallocate(ens_sum) end if + if ( calc_regcoeff ) then + ! write out psi + do ie = ens_istart, ens_iend ! each proc loops over a subset of ens + write(ce,'(i3.3)') ie + output_file = 'psi.'//filedates(ie,ic)//'.e'//trim(ce) + !write(stdout,'(a,i3,a,a)') ' Proc', myproc, ' writing to ', trim(output_file) + open (ounit, file = output_file, form='unformatted') + write(ounit) 1, ni, nj, nk + write(ounit) 3 + write(ounit) filedates(ie,ic), 'PSI ' + write(ounit) xdata(indx_psi,ie,ic)%value(:,:,:) + close(ounit) + end do + end if + if ( write_pert1 .or. do_slen_calc ) then do ie = ens_istart, ens_iend ! each proc loops over a subset of ens write(ce,'(i3.3)') ie @@ -1426,7 +2003,7 @@ subroutine compute_pert allocate(r8tmp(ni,nj,nk)) end if do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + if ( .not. do_this_var(iv) ) cycle do ie = ens_istart, ens_iend ! each proc loops over a subset of ens write(ce,'(i3.3)') ie !output_file = trim(varnames(iv))//'.e'//trim(ce) @@ -1498,9 +2075,17 @@ subroutine compute_pert end if allocate (locbuf(ni, nj, nk, nens)) +#ifdef DM_PARALLEL + if ( kind(locbuf) == 4 ) then + this_mpi_real = mpi_real + else + this_mpi_real = mpi_real8 + end if +#endif + my_nvar = 0 var_loop3: do iv = 1, nvar - if ( .not. read_it(iv) ) cycle var_loop3 + if ( .not. do_this_var(iv) ) cycle var_loop3 ! initialize to zero for each variable if ( myproc == root ) globuf(:,:,:,:) = 0.0 @@ -1512,7 +2097,7 @@ subroutine compute_pert #ifdef DM_PARALLEL call mpi_reduce(locbuf,globuf,ijk*nens, & - mpi_real,mpi_sum,root,mpi_comm_world,ierr) + this_mpi_real,mpi_sum,root,mpi_comm_world,ierr) if ( ierr /= 0 ) then write(stdout, '(a, i3)') 'Error mpi_reduce on proc', myproc call mpi_abort(mpi_comm_world,1,ierr) @@ -1612,6 +2197,13 @@ subroutine compute_pert end do case_loop1 ! case loop + if ( allocated(vor) ) deallocate (vor) + if ( allocated(div) ) deallocate (div) + if ( allocated(psi) ) deallocate (psi) + if ( allocated(chi) ) deallocate (chi) + if ( allocated(ustag) ) deallocate (ustag) + if ( allocated(vstag) ) deallocate (vstag) + if ( allocated(psfc) ) deallocate (psfc) if ( allocated(prs) ) deallocate (prs) if ( allocated(ght) ) deallocate (ght) if ( allocated(theta) ) deallocate (theta) @@ -1623,22 +2215,22 @@ subroutine compute_pert if ( trim(be_method) == 'NMC' ) then if ( remove_time_mean ) then - if ( .not. allocated(time_mean) ) allocate(time_mean(ni,nj,nk,nvar)) + if ( .not. allocated(time_mean) ) allocate(time_mean(ni,nj,nk,nvar_all)) #ifdef DM_PARALLEL - if ( .not. allocated(xsum) ) allocate(xsum(ni,nj,nk,nvar)) - do iv = 1, nvar + if ( .not. allocated(xsum) ) allocate(xsum(ni,nj,nk,nvar_all)) + do iv = 1, nvar_all call mpi_allreduce(xsum_loc(:,:,:,iv),xsum(:,:,:,iv), ijk, & mpi_real8, mpi_sum, mpi_comm_world,ierr) time_mean(:,:,:,iv) = xsum(:,:,:,iv)/real(ncase) end do deallocate(xsum) #else - do iv = 1, nvar + do iv = 1, nvar_all time_mean(:,:,:,iv) = xsum_loc(:,:,:,iv)/real(ncase) end do #endif do ic = case_istart, case_iend - do iv = 1, nvar + do iv = 1, nvar_all xdata(iv,1,ic)%value(:,:,:) = xdata(iv,1,ic)%value(:,:,:) - time_mean(:,:,:,iv) end do end do @@ -1646,6 +2238,20 @@ subroutine compute_pert deallocate(xsum_loc) end if ! remove_time_mean + if ( calc_regcoeff ) then + ! write out psi + do ic = case_istart, case_iend + output_file = 'psi.'//filedates(1,ic)//'.e001' + !write(stdout,'(a,i3,a,a)') ' Proc', myproc, ' writing to ', trim(output_file) + open (ounit, file = output_file, form='unformatted') + write(ounit) 1, ni, nj, nk + write(ounit) 3 + write(ounit) filedates(1,ic), 'PSI ' + write(ounit) xdata(indx_psi,1,ic)%value(:,:,:) + close(ounit) + end do + end if + if ( write_pert1 .or. do_slen_calc ) then do ic = case_istart, case_iend output_file = 'pert1.'//filedates(1,ic)//'.e001' @@ -1666,15 +2272,12 @@ subroutine compute_pert end if end if ! be_method=NMC - deallocate (istart_ens) - deallocate (iend_ens ) - deallocate (istart_case) - deallocate (iend_case ) - do ic = case_istart, case_iend do ie = ens_istart, ens_iend - do iv = 1, nvar - deallocate(xdata(iv,ie,ic)%value) + do iv = 1, nvar_all + if ( allocated(xdata(iv,ie,ic)%value) ) then + deallocate(xdata(iv,ie,ic)%value) + end if end do end do end do @@ -1682,16 +2285,509 @@ subroutine compute_pert end subroutine compute_pert +subroutine compute_regcoeff_unbalanced +! adapted from var/gen_be/gen_be_stage2.f90 and var/gen_be/gen_be_stage2a.f90 + + implicit none + + real, parameter :: variance_threshold = 1e-6 ! Percentage of variance discarded. + + character(len=filename_len) :: input_file, output_file + character*3 :: ce ! Member index -> character + integer :: i, j, k, member, k2, k3, m ! Loop counters + integer :: b ! Bin marker + integer :: mmax ! Maximum mode (after variance truncation) + real :: coeffa, coeffb ! Accumulating mean coefficients + real :: total_variance ! Total variance of matrix + real :: cumul_variance ! Cumulative variance of matrix + real :: summ ! Summation dummy + + real, allocatable :: psi(:,:,:) ! psi + real, allocatable :: chi(:,:,:) ! chi + real, allocatable :: temp(:,:,:) ! Temperature + real, allocatable :: ps(:,:,:) ! Surface pressure + + integer, allocatable:: bin_pts(:) ! Number of points in bin (3D fields) + integer, allocatable:: bin_pts2d(:) ! Number of points in bin (2D fields) + real, allocatable :: covar1(:) ! Covariance between input fields + real, allocatable :: covar2(:,:) ! Covariance between input fields + real, allocatable :: covar3(:,:,:) ! Covariance between input fields + real, allocatable :: var1(:) ! Autocovariance of field + real, allocatable :: var2(:,:,:) ! Autocovariance of field + real, allocatable :: var2_inv(:,:,:) ! Inverse Autocovariance of field + + real, allocatable :: work(:,:) ! EOF work array + real*8, allocatable :: evec(:,:) ! Gridpoint eigenvectors + real*8, allocatable :: eval(:) ! Gridpoint sqrt(eigenvalues) + real, allocatable :: LamInvET(:,:) ! ET/sqrt(Eigenvalue) + real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient + real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient + real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient + + integer :: ic, ie, iv + integer :: istart_member4cov, iend_member4cov ! index for covariance loop + integer :: istart_member4bal, iend_member4bal ! index for unbalanced loop + integer :: ounit + integer :: this_mpi_real + integer :: proc_send + + logical :: got_var2_inv + + if ( myproc == root ) then + write(stdout,'(a)')' ====== Computing regcoeff and unbalanced fields ======' + end if + + ounit = 24 + + allocate( psi(1:ni,1:nj,1:nk) ) + allocate( bin_pts(1:num_bins) ) + allocate( bin_pts2d(1:num_bins2d) ) + + if ( trim(be_method) == 'NMC' ) then + istart_member4cov = 1 + iend_member4cov = 1 + istart_member4bal = 1 + iend_member4bal = 1 + else if ( trim(be_method) == 'ENS' ) then + istart_member4cov = 1 + iend_member4cov = nens + istart_member4bal = istart_ens(myproc) + iend_member4bal = iend_ens(myproc) + end if + + allocate( regcoeff1(1:num_bins) ) + allocate( regcoeff2(1:nk,1:num_bins2d) ) + allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) ) + regcoeff1 = 0.0 + regcoeff2 = 0.0 + regcoeff3 = 0.0 + +#ifdef DM_PARALLEL + if ( kind(regcoeff1) == 4 ) then + this_mpi_real = mpi_real + else + this_mpi_real = mpi_real8 + end if +#endif + + ! nvar distributed among processors + var_loop_regcoef: do iv = 1, nvar + + if ( .not. do_this_var(iv) ) cycle var_loop_regcoef + if ( myproc /= MOD((iv-1), num_procs) ) cycle var_loop_regcoef + if ( trim(varnames(iv)) /= 'CHI_U' .and. & + trim(varnames(iv)) /= 'PSFC_U' .and. & + trim(varnames(iv)) /= 'T_U' ) cycle var_loop_regcoef + + !write(stdout,'(a,i3,2a)') ' Proc', myproc, ' Computing regcoeff for variable ', trim(varnames(iv)) + + bin_pts(:) = 0 + bin_pts2d(:) = 0 + + if ( trim(varnames(iv)) == 'CHI_U' ) then + if ( .not. allocated(chi) ) allocate( chi(1:ni,1:nj,1:nk) ) + if ( .not. allocated(covar1) ) allocate( covar1(1:num_bins) ) + if ( .not. allocated(var1) ) allocate( var1(1:num_bins) ) + chi(:,:,:) = 0.0 + covar1(:) = 0.0 + var1(:) = 0.0 + end if + + if ( trim(varnames(iv)) == 'PSFC_U' .or. trim(varnames(iv)) == 'T_U' ) then + if ( .not. allocated(var2) ) allocate( var2(1:nk,1:nk,1:num_bins2d) ) + var2(:,:,:) = 0.0 + if ( trim(varnames(iv)) == 'PSFC_U' ) then + ! still allocate 1:nk for ps even though only 1 will be used + ! this is needed to make get_pert1_data work for ps + if ( .not. allocated(ps) ) allocate( ps(1:ni,1:nj,1:nk) ) + if ( .not. allocated(covar2) ) allocate( covar2(1:nk,1:num_bins2d) ) + ps(:,:,:) = 0.0 + covar2(:,:) = 0.0 + end if + if ( trim(varnames(iv)) == 'T_U' ) then + if ( .not. allocated(temp) ) allocate( temp(1:ni,1:nj,1:nk) ) + if ( .not. allocated(covar3) ) allocate( covar3(1:nk,1:nk,1:num_bins2d) ) + temp(:,:,:) = 0.0 + covar3(:,:,:) = 0.0 + end if + end if + + write(stdout,'(a,i3,2a)') ' Proc', myproc, ' Computing psi covariances for variable ', trim(varnames(iv)) + + do ic = 1, ncase + do ie = istart_member4cov, iend_member4cov + !write(stdout,'(5a,i4)')' Processing data for date ', filedates(ie,ic), ', variable ', trim(varnames(iv)), & + ! ', member ', ie + + ! psi is needed for all variables + write(ce,'(i3.3)') ie + input_file = 'psi.'//filedates(ie,ic)//'.e'//trim(ce) + call get_pert1_data(trim(input_file), 'PSI', ni, nj, nk, psi) + + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + if ( trim(varnames(iv)) == 'CHI_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, chi) + ! Calculate psi/chi covariances: + do k = 1, nk + do j = 1, nj + do i = 1, ni + b = bin(i,j,k) + bin_pts(b) = bin_pts(b) + 1 + coeffa = 1.0 / real(bin_pts(b)) + coeffb = real(bin_pts(b)-1) * coeffa + covar1(b) = coeffb * covar1(b) + coeffa * psi(i,j,k) * chi(i,j,k) + var1(b) = coeffb * var1(b) + coeffa * psi(i,j,k) * psi(i,j,k) + end do + end do + end do + end if ! chi + + if ( trim(varnames(iv)) == 'PSFC_U' .or. trim(varnames(iv)) == 'T_U' ) then + + if ( trim(varnames(iv)) == 'PSFC_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, ps) + end if ! ps + if ( trim(varnames(iv)) == 'T_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, temp) + end if ! t + + ! Calculate psi/ps, and psi/T, and psi/psi covariances: + do j = 1, nj + do i = 1, ni + b = bin2d(i,j) + bin_pts2d(b) = bin_pts2d(b) + 1 + coeffa = 1.0 / real(bin_pts2d(b)) + coeffb = real(bin_pts2d(b)-1) * coeffa + do k = 1, nk + + if ( trim(varnames(iv)) == 'PSFC_U' ) then + ! psi/ps: + covar2(k,b) = coeffb * covar2(k,b) + coeffa * ps(i,j,1) * psi(i,j,k) + end if + + if ( trim(varnames(iv)) == 'T_U' ) then + ! psi/T: + do k2 = 1, nk + covar3(k,k2,b) = coeffb * covar3(k,k2,b) + & + coeffa * temp(i,j,k) * psi(i,j,k2) + end do + end if + + ! psi/psi (symmetric): + do k2 = 1, k + var2(k,k2,b) = coeffb * var2(k,k2,b) + & + coeffa * psi(i,j,k) * psi(i,j,k2) + end do + + end do ! k loop + end do ! i loop + end do ! j loop + + end if ! ps_u or t_u + + end do ! ens member loop + + end do ! ncase loop + + got_var2_inv = .false. ! initialize + + if ( trim(varnames(iv)) == 'PSFC_U' .or. trim(varnames(iv)) == 'T_U' ) then + + ! Fill in psi/psi covariance by symmetry: + do b = 1, num_bins2d + do k = 1, nk + do k2 = k+1, nk ! Symmetry. + var2(k,k2,b) = var2(k2,k,b) + end do + end do + end do + + if ( .not. got_var2_inv ) then + + write(stdout,'(a,i3,2a)') ' Proc', myproc, ' Computing eigenvectors, eigenvalues and inverse for psi/psi covariance' + + if ( .not. allocated(work) ) allocate( work(1:nk,1:nk) ) + if ( .not. allocated(evec) ) allocate( evec(1:nk,1:nk) ) + if ( .not. allocated(eval) ) allocate( eval(1:nk) ) + if ( .not. allocated(LamInvET) ) allocate( LamInvET(1:nk,1:nk) ) + if ( .not. allocated(var2_inv) ) allocate( var2_inv(1:nk,1:nk,1:num_bins2d) ) + + do b = 1, num_bins2d + LamInvET(:,:) = 0.0 + work(1:nk,1:nk) = var2(1:nk,1:nk,b) + call da_eof_decomposition( nk, work, evec, eval ) + + ! Truncate eigenvalues to ensure inverse is not dominated by rounding error: + summ = 0.0 + do m = 1, nk + summ = summ + eval(m) + end do + total_variance = summ + + cumul_variance = 0.0 + mmax = nk + do m = 1, nk + cumul_variance = cumul_variance + eval(m) / total_variance + if ( cumul_variance > 1.0 - variance_threshold ) then + mmax = m - 1 + exit + end if + end do + !write(6,'(2(a,i6),2(a,1pe11.5))') ' Bin = ', b, ', truncation = ', mmax, & + ! ', Total Variance = ', total_variance, & + ! ', Condition number = ', eval(1) / eval(nk-1) + + ! Lam{-1} . E^T: + do k = 1, nk + do m = 1, mmax + LamInvET(m,k) = evec(k,m) / eval(m) + end do + end do + + ! ^{-1} = E . Lam{-1} . E^T: + do k = 1, nk + do k2 = 1, k + summ = 0.0 + do m = 1, nk + summ = summ + evec(k,m) * LamInvET(m,k2) + end do + var2_inv(k,k2,b) = summ + end do + end do + + do k = 1, nk + do k2 = k+1, nk ! Symmetry. + var2_inv(k,k2,b) = var2_inv(k2,k,b) + end do + end do + end do + got_var2_inv = .true. + end if ! got_var2_inv + end if ! ps_u or t_u + + write(stdout,'(a,i3,2a)') ' Proc', myproc, ' Computing regression coefficients for variable ', trim(varnames(iv)) + + if ( trim(varnames(iv)) == 'CHI_U' ) then + ! psi/chi: + do b = 1, num_bins + regcoeff1(b) = covar1(b) / var1(b) + end do + ! Output regression coefficients + output_file = 'regcoeff1.dat' + open (ounit, file=trim(output_file), form='unformatted', status='unknown') + write(ounit)ni, nj, nk + write(ounit)num_bins, num_bins2d + write(ounit)regcoeff1 + close(ounit) + end if ! chi_u + + if ( trim(varnames(iv)) == 'PSFC_U' ) then + ! psi/ps: + do b = 1, num_bins2d + do k = 1, nk + summ = 0.0 + do k2 = 1, nk + summ = summ + covar2(k2,b) * var2_inv(k2,k,b) + end do + regcoeff2(k,b) = summ + end do + end do + ! Output regression coefficients + output_file = 'regcoeff2.dat' + open (ounit, file=trim(output_file), form='unformatted', status='unknown') + write(ounit)ni, nj, nk + write(ounit)num_bins, num_bins2d + write(ounit)regcoeff2 + close(ounit) + end if ! ps_u + + if ( trim(varnames(iv)) == 'T_U' ) then + ! psi/T: + do b = 1, num_bins2d + do k = 1, nk + do k2 = 1, nk + summ = 0.0 + do k3 = 1, nk + summ = summ + covar3(k,k3,b) * var2_inv(k3,k2,b) + end do + regcoeff3(k,k2,b) = summ + end do + end do + end do + ! Output regression coefficients + output_file = 'regcoeff3.dat' + open (ounit, file=trim(output_file), form='unformatted', status='unknown') + write(ounit)ni, nj, nk + write(ounit)num_bins, num_bins2d + write(ounit)regcoeff3 + close(ounit) + end if ! t_u + + if ( trim(varnames(iv)) == 'PSFC_U' .or. trim(varnames(iv)) == 'T_U' ) then + if ( allocated(var2_inv) ) deallocate( var2_inv ) + if ( allocated(LamInvET) ) deallocate( LamInvET ) + if ( allocated(eval) ) deallocate( eval ) + if ( allocated(evec) ) deallocate( evec ) + if ( allocated(work) ) deallocate( work ) + if ( allocated(var2) ) deallocate( var2 ) + end if + + if ( trim(varnames(iv)) == 'CHI_U' ) then + deallocate( covar1 ) + deallocate( var1 ) + end if + if ( trim(varnames(iv)) == 'PSFC_U' ) then + deallocate( covar2 ) + end if + if ( trim(varnames(iv)) == 'T_U' ) then + deallocate( covar3 ) + end if + + end do var_loop_regcoef + +#ifdef DM_PARALLEL + do iv = 1, nvar + if ( .not. do_this_var(iv) ) cycle + proc_send = MOD((iv-1), num_procs) + if ( trim(varnames(iv)) == 'CHI_U' ) then + call mpi_bcast(regcoeff1(:), num_bins, this_mpi_real, proc_send, mpi_comm_world, ierr ) + end if + if ( trim(varnames(iv)) == 'PSFC_U' ) then + call mpi_bcast(regcoeff2(:,:), nk*num_bins2d, this_mpi_real, proc_send, mpi_comm_world, ierr ) + end if + if ( trim(varnames(iv)) == 'T_U' ) then + call mpi_bcast(regcoeff3(:,:,:), nk*nk*num_bins2d, this_mpi_real, proc_send, mpi_comm_world, ierr ) + end if + end do +#endif + + var_loop_unbal: do iv = 1, nvar + + if ( .not. do_this_var(iv) ) cycle var_loop_unbal + if ( trim(varnames(iv)) /= 'CHI_U' .and. & + trim(varnames(iv)) /= 'PSFC_U' .and. & + trim(varnames(iv)) /= 'T_U' ) cycle var_loop_unbal + + if ( myproc == root ) write(stdout,'(2a)') ' Computing unbalanced field for variable ', trim(varnames(iv)) + + if ( trim(varnames(iv)) == 'CHI_U' ) then + if ( .not. allocated(chi) ) allocate( chi(1:ni,1:nj,1:nk) ) + chi(:,:,:) = 0.0 + end if + if ( trim(varnames(iv)) == 'PSFC_U' ) then + if ( .not. allocated(ps) ) allocate( ps(1:ni,1:nj,1:nk) ) + ps(:,:,:) = 0.0 + end if + if ( trim(varnames(iv)) == 'T_U' ) then + if ( .not. allocated(temp) ) allocate( temp(1:ni,1:nj,1:nk) ) + temp(:,:,:) = 0.0 + end if + + ! ncase/nens distributed among processors + do ic = case_istart, case_iend + do ie = istart_member4bal, iend_member4bal + + !write(stdout,'(a,i3,4a)') ' Proc', myproc, ' Computing unbalanced field for variable ', trim(varnames(iv)), & + ! ' date ', filedates(ie,ic) + + ! psi is needed for all variables + write(ce,'(i3.3)') ie + input_file = 'psi.'//filedates(ie,ic)//'.e'//trim(ce) + call get_pert1_data(trim(input_file), 'PSI', ni, nj, nk, psi) + + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + output_file = 'pert1_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) + + if ( trim(varnames(iv)) == 'CHI_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, chi) + do k = 1, nk + do j = 1, nj + do i = 1, ni + b = bin(i,j,k) + chi(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k) + end do + end do + end do + open (ounit, file=output_file, form='unformatted') + write(ounit) 1, ni, nj, nk + write(ounit) 3 + write(ounit) filedates(ie,ic), varnames(iv) + write(ounit) chi(:,:,:) + close(ounit) + end if ! chi_u + + if ( trim(varnames(iv)) == 'PSFC_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, ps) + do j = 1, nj + do i = 1, ni + b = bin2d(i,j) + ps(i,j,1) = ps(i,j,1) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk)) + end do + end do + open (ounit, file=output_file, form='unformatted') + write(ounit) 1, ni, nj, nk + write(ounit) 2 + write(ounit) filedates(ie,ic), varnames(iv) + write(ounit) ps(:,:,1) + close(ounit) + end if ! ps_u + + if ( trim(varnames(iv)) == 'T_U' ) then + call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, temp) + do j = 1, nj + do i = 1, ni + b = bin2d(i,j) + do k = 1, nk + temp(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk)) + end do + end do + end do + open (ounit, file=output_file, form='unformatted') + write(ounit) 1, ni, nj, nk + write(ounit) 3 + write(ounit) filedates(ie,ic), varnames(iv) + write(ounit) temp(:,:,:) + close(ounit) + end if ! t_u + + end do ! ens member loop + end do ! ncase loop + +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif + + if ( trim(varnames(iv)) == 'CHI_U' ) then + if ( allocated(chi) ) deallocate( chi ) + end if + if ( trim(varnames(iv)) == 'PSFC_U' ) then + if ( allocated(ps) ) deallocate( ps ) + end if + if ( trim(varnames(iv)) == 'T_U' ) then + if ( allocated(temp) ) deallocate( temp ) + end if + + end do var_loop_unbal + + deallocate( regcoeff1 ) + deallocate( regcoeff2 ) + deallocate( regcoeff3 ) + deallocate( psi ) + deallocate( bin_pts ) + deallocate( bin_pts2d ) + +end subroutine compute_regcoeff_unbalanced + subroutine compute_bv_sl - !use da_lapack, only : dsyev implicit none character(len=filename_len) :: filename ! Input filename. character(len=filename_len) :: output_file integer :: i, j, k, k1, k2, b ! Loop counters. integer :: ic, ie, iv, m - integer :: istart_member, iend_member + integer :: istart_member, iend_member ! loop index for nens + integer :: istart_member_p, iend_member_p ! loop index for projecting to modes real :: inv_nij ! 1 / (ni*nj). real :: mean_field ! Mean field. real :: coeffa, coeffb ! Accumulating mean coefficients. @@ -1722,6 +2818,7 @@ subroutine compute_bv_sl integer :: nn integer :: nvar_read, ni_read, nj_read, nk_read, nkk integer :: k_start, k_end + integer :: this_mpi_real integer, allocatable:: nr(:), icount(:) real, allocatable :: cov(:,:) real, allocatable :: ml(:), sl(:) @@ -1737,12 +2834,24 @@ subroutine compute_bv_sl inv_nij = 1.0 / real(ni*nj) allocate( field(1:ni,1:nj,1:nk) ) +#ifdef DM_PARALLEL + if ( kind(field) == 4 ) then + this_mpi_real = mpi_real + else + this_mpi_real = mpi_real8 + end if +#endif + if ( trim(be_method) == 'NMC' ) then - istart_member = 1 - iend_member = 1 + istart_member = 1 + iend_member = 1 + istart_member_p = 1 + iend_member_p = 1 else if ( trim(be_method) == 'ENS' ) then - istart_member = 1 - iend_member = nens + istart_member = 1 + iend_member = nens + istart_member_p = istart_ens(myproc) + iend_member_p = iend_ens(myproc) end if if ( pert1_read_opt == 1 ) then @@ -1781,7 +2890,7 @@ subroutine compute_bv_sl if ( do_eof_transform ) then - if ( myproc == root ) write(stdout,'(4a)')'====== Computing vertical error eigenvalues, eigenvectors ======' + if ( myproc == root ) write(stdout,'(4a)')' ====== Computing vertical error eigenvalues, eigenvectors ======' allocate( bv(1:nk,1:nk,1:num_bins2d) ) allocate( bin_pts2d(1:num_bins2d) ) @@ -1797,10 +2906,11 @@ subroutine compute_bv_sl allocate( fieldv(1:ni,1:nj,1:nk) ) end if - var_loop: do iv = 1, nvar + ! nvar distributed among processors + var_loop_bv: do iv = 1, nvar - if ( .not. read_it(iv) ) cycle var_loop - if ( myproc /= MOD((iv-1), num_procs) ) cycle var_loop + if ( .not. do_this_var(iv) ) cycle var_loop_bv + if ( myproc /= MOD((iv-1), num_procs) ) cycle var_loop_bv write(stdout,'(a,i3,2a)') ' Proc', myproc, ' Processing vertical error stats for variable ', trim(varnames(iv)) @@ -1825,7 +2935,13 @@ subroutine compute_bv_sl else if ( pert1_read_opt == 2 ) then write(ce,'(i3.3)') ie - input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + if ( trim(varnames(iv)) == 'CHI_U' .or. & + trim(varnames(iv)) == 'T_U' .or. & + trim(varnames(iv)) == 'PSFC_U' ) then + input_file = 'pert1_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) + else + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + end if call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, field) @@ -1858,6 +2974,12 @@ subroutine compute_bv_sl do b = 1, num_bins2d do k1 = 1, nkk do k2 = k1+1, nkk ! Symmetry. + if ( vertloc_opt > 0 ) then + bv(k2,k1,b) = bv(k2,k1,b) * vertloc_rho(k2,k1) + end if + !if ( k1 > 40 ) then ! arbitrarily remove negative values near top + ! bv(k2,k1,b) = max(0.0, bv(k2,k1,b)) + !end if bv(k1,k2,b) = bv(k2,k1,b) end do end do @@ -1931,10 +3053,40 @@ subroutine compute_bv_sl end do end do + ! Output eigenvectors, eigenvalues + filename = 'EV_'//trim(varnames(iv))//'.dat' + open (ounit, file = filename, form='unformatted') + write(ounit)varnames(iv) + write(ounit)ni, nj, nkk + write(ounit)evec(1:nj,1:nkk,1:nkk) + write(ounit)eval(1:nj,1:nkk) + close(ounit) + end do var_loop_bv + +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif + + var_loop_proj: do iv = 1, nvar + + if ( .not. do_this_var(iv) ) cycle var_loop_proj + if ( var_dim(iv) == 3 ) then + + ! raed eigenvectors, eigenvalues + filename = 'EV_'//trim(varnames(iv))//'.dat' + open (iunit, file = filename, form='unformatted') + read(iunit) !varnames(iv) + read(iunit) !ni, nj, nkk + read(iunit) evec(1:nj,1:nk,1:nk) + read(iunit) eval(1:nj,1:nk) + close(iunit) + if ( myproc == root ) write(stdout,'(a)')' Projecting fields onto vertical modes' - do ic = 1, ncase - do ie = istart_member, iend_member + + ! ncase/nens distributed among processors + do ic = case_istart, case_iend + do ie = istart_member_p, iend_member_p !write(stdout,'(5a,i4)')' Date = ', filedates(ie,ic), ', variable ', trim(varnames(iv)), & ! ', member ', ie ! Project fields onto vertical modes @@ -1951,7 +3103,13 @@ subroutine compute_bv_sl xdata(iv,ie,ic)%value(:,:,:) = field(:,:,:) else if ( pert1_read_opt == 2 ) then write(ce,'(i3.3)') ie - input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + if ( trim(varnames(iv)) == 'CHI_U' .or. & + trim(varnames(iv)) == 'T_U' .or. & + trim(varnames(iv)) == 'PSFC_U' ) then + input_file = 'pert1_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) + else + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + end if call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, field) do m = 1, nk do j = 1, nj @@ -1971,22 +3129,26 @@ subroutine compute_bv_sl close(ounit) end if ! pert1_read_opt end do ! End loop over members. - end do - end if - end do var_loop + end do ! case loop + end if ! var_dim=3 + end do var_loop_proj + +#ifdef DM_PARALLEL + call mpi_barrier(mpi_comm_world,ierr) +#endif if ( pert1_read_opt == 1 ) then #ifdef DM_PARALLEL ijk = ni*nj*nk do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + if ( .not. do_this_var(iv) ) cycle proc_send = MOD((iv-1), num_procs) do ic = 1, ncase do ie = istart_member, iend_member if ( myproc == MOD((iv-1), num_procs) ) then field(:,:,:) = xdata(iv,ie,ic)%value(:,:,:) end if - call mpi_bcast(field(:,:,:), ijk, mpi_real, proc_send, mpi_comm_world, ierr ) + call mpi_bcast(field(:,:,:), ijk, this_mpi_real, proc_send, mpi_comm_world, ierr ) if ( myproc /= MOD((iv-1), num_procs) ) then xdata(iv,ie,ic)%value(:,:,:) = field(:,:,:) end if @@ -2029,7 +3191,7 @@ subroutine compute_bv_sl allocate(ml(nk)) do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + if ( .not. do_this_var(iv) ) cycle sl(:) = 0.0 ml(:) = 0.0 if ( var_dim(iv) == 2 ) then @@ -2061,7 +3223,13 @@ subroutine compute_bv_sl if ( do_eof_transform .and. (var_dim(iv) == 3) ) then input_file = 'pertv_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) else - input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + if ( trim(varnames(iv)) == 'CHI_U' .or. & + trim(varnames(iv)) == 'T_U' .or. & + trim(varnames(iv)) == 'PSFC_U' ) then + input_file = 'pert1_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) + else + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + end if ! unbalanced variables end if call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, field) field2d(:,:) = field(:,:,k) @@ -2095,7 +3263,7 @@ subroutine compute_bv_sl allocate( hl(1:num_bins2d) ) allocate( hlt(1:num_bins2d) ) do iv = 1, nvar - if ( .not. read_it(iv) ) cycle + if ( .not. do_this_var(iv) ) cycle sl(:) = 0.0 if ( var_dim(iv) == 2 ) then k_start = 1 @@ -2119,7 +3287,13 @@ subroutine compute_bv_sl if ( do_eof_transform .and. (var_dim(iv) == 3) ) then input_file = 'pertv_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) else - input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + if ( trim(varnames(iv)) == 'CHI_U' .or. & + trim(varnames(iv)) == 'T_U' .or. & + trim(varnames(iv)) == 'PSFC_U' ) then + input_file = 'pert1_'//trim(varnames(iv))//'.'//filedates(ie,ic)//'.e'//trim(ce) + else + input_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) + end if ! unbalanced end if call get_pert1_data(trim(input_file), trim(varnames(iv)), ni, nj, nk, field) field2d(:,:) = field(:,:,k) @@ -2135,12 +3309,12 @@ subroutine compute_bv_sl end if sl(k) = hl(1) !hcl-num_bins2d=1 sl(k) = sl(k) / ds ! convert to grid unit - write(stdout,'(a,a,i3,f10.3)') 'ScaleLength: ', varnames(iv), k, sl(k)*ds*0.001 + write(stdout,'(a,a,i3,f10.3,a)') 'ScaleLength: ', varnames(iv), k, sl(k)*ds*0.001, ' km' be_data(iv)%scale_length(k) = sl(k) end do k_loop_opt2 #ifdef DM_PARALLEL call mpi_allreduce(sl(:),sl_g(:), nk, & - mpi_real, mpi_sum, mpi_comm_world, ierr) + this_mpi_real, mpi_sum, mpi_comm_world, ierr) #else sl_g(:) = sl(:) #endif @@ -2462,6 +3636,9 @@ subroutine da_eof_decomposition (kz, bx, e, l) ! BE field. !--------------------------------------------------------------------------- +#ifdef USE_WRFDA + use da_lapack, only : dsyev +#endif implicit none integer, intent(in) :: kz ! Dimension of error matrix. @@ -2714,6 +3891,8 @@ subroutine get_pert1_data(pert_file, vname, nx, ny, nz, vdata) character(len=DateStrLen) :: DateStr_tmp character(len=VarNameLen) :: var_tmp + iunit = 21 + !pert_file = 'pert1.'//filedates(ie,ic)//'.e'//trim(ce) inquire(file=trim(pert_file), exist=isfile) if ( .not. isfile ) then @@ -2752,6 +3931,71 @@ subroutine get_pert1_data(pert_file, vname, nx, ny, nz, vdata) end subroutine get_pert1_data +subroutine calc_es_qs (t, p, es, qs, wrt_ice) + +! Purpose: calculate saturation vapor pressure and saturation specific humidity +! given temperature and pressure + + implicit none + + real, intent(in) :: t, p + real, intent(out) :: es, qs + logical, intent(in), optional :: wrt_ice + + real, parameter :: es_alpha = 611.2 + real, parameter :: es_beta = 17.67 + real, parameter :: es_gamma = 243.5 + real, parameter :: t_c_ref1 = 0.0 ! C + real, parameter :: t_c_ref2 = -20.0 ! C + real, parameter :: a0 = 6.107799961 + real, parameter :: a1 = 4.436518521e-01 + real, parameter :: a2 = 1.428945805e-02 + real, parameter :: a3 = 2.650648471e-04 + real, parameter :: a4 = 3.031240396e-06 + real, parameter :: a5 = 2.034080948e-08 + real, parameter :: a6 = 6.136820929e-11 + real, parameter :: c1 = 9.09718 + real, parameter :: c2 = 3.56654 + real, parameter :: c3 = 0.876793 + real, parameter :: c4 = 6.1071 + real :: t_c, t1_c ! T in degree C + logical :: ice + + ice = .false. + if ( present(wrt_ice) ) then + if ( wrt_ice ) ice = .true. + end if + + t_c = t - t_kelvin + t1_c = t - t_triple + + ! Calculate saturation vapor pressure es + + if ( .not. ice ) then ! over water only + + es = es_alpha * exp( es_beta * t_c / ( t_c + es_gamma ) ) + + else ! consider ice-water and ice effects + + if ( t1_c > t_c_ref1 ) then ! vapor pressure over water + es = es_alpha * exp( es_beta * t_c / ( t_c + es_gamma ) ) + else if ( (t1_c <= t_c_ref1) .and. (t1_c >= t_c_ref2) ) then ! vapor pressure over water below 0C + es = a0 + t1_c * (a1 + t1_c * (a2 + t1_c * (a3 + t1_c * (a4 + t1_c * (a5 + t1_c * a6))))) + es = es * 100.0 ! to Pa + else ! vapor pressure over ice + es = 10.0 ** ( -c1 * (t_triple / t - 1.0) - c2 * alog10(t_triple / t) + & + c3 * (1.0 - t / t_triple) + alog10(c4)) + es = es * 100.0 ! to Pa + end if + + end if + + ! Calculate saturation specific humidity qs + + qs = rd_over_rv * es / ( p - rd_over_rv1 * es ) + +end subroutine calc_es_qs + end program gen_be_v3 subroutine para_range(n1, n2, nprocs, myrank, ista, iend) @@ -3068,3 +4312,1797 @@ subroutine dsyev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) end subroutine dsyev #endif +subroutine da_uv_to_div_c( dim1, dim2, dim3, ds, & + mapfac_m, mapfac_u, mapfac_v, & + u, v, div ) + +!------------------------------------------------------------------------------ +! PURPOSE: Calculate divergence on a co-ordinate surface, given an input +! wind field on an Arakawa C-grid. +! +! NOTE: No boundary conditions required on the WRF Arakawa C-grid as +! divergence (mass) points are all within the outer u/v pts. +! +! +! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker +! d U d V +! Div = m^2 *[---(---) + ---(---) ] +! dx m dy M +!------------------------------------------------------------------------------ + + implicit none + + integer, intent(in):: dim1, dim2, dim3 ! Dimensions + real, intent(in) :: ds ! Resolution + real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts + real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points + real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points + real, intent(in) :: u(1:dim1+1,1:dim2,1:dim3) ! v wind + real, intent(in) :: v(1:dim1,1:dim2+1,1:dim3) ! v wind + real, intent(out) :: div(1:dim1,1:dim2,1:dim3) ! Divergence + + integer :: i, j, k ! Loop counters + real :: ds_inv ! 1/ds + real :: coeff(1:dim1,1:dim2) ! Coefficient + real :: um(1:dim1+1,1:dim2,1:dim3) ! u-wind copy + real :: vm(1:dim1,1:dim2+1,1:dim3) ! v-wind copy + +!------------------------------------------------------------------------------ +! [1.0] Initialise: +!------------------------------------------------------------------------------ + + div(:,:,:) = 0.0 + + ds_inv = 1.0 / ds + coeff(:,:) = ds_inv + +!------------------------------------------------------------------------------ +! [2] Calculate scaled u/v field: +!------------------------------------------------------------------------------ + + do k = 1, dim3 + do j = 1, dim2 + do i = 1, dim1+1 + um(i,j,k) = u(i,j,k) / mapfac_u(i,j) + end do + end do + end do + + do k = 1, dim3 + do j = 1, dim2+1 + do i = 1, dim1 + if (mapfac_v(i,j) > 0.000001) then + vm(i,j,k) = v(i,j,k) / mapfac_v(i,j) + else + vm(i,j,k)=0.0 + end if + end do + end do + end do + +!------------------------------------------------------------------------------ +! [3] Calculate divergence field: +!------------------------------------------------------------------------------ + + do k = 1, dim3 + do j = 1, dim2 + do i = 1, dim1 + div(i,j,k) = coeff(i,j) * ( um(i+1,j,k) - um(i,j,k) + vm(i,j+1,k) - vm(i,j,k) ) + end do + end do + end do + +end subroutine da_uv_to_div_c + +subroutine da_uv_to_vor_c( dim1, dim2, dim3, ds, & + mapfac_m, mapfac_u, mapfac_v, & + u, v, vor ) + +!------------------------------------------------------------------------------ +! PURPOSE: Calculate vorticity on a co-ordinate surface, given an input +! wind field on an Arakawa C-grid. +! +! NOTE: Zero vorticity boundary conditions. +! +! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker +! d V d U +! Vor = m^2 *[---(---) - ---(---) ] +! dx m dy M +!------------------------------------------------------------------------------ + + implicit none + + integer, intent(in):: dim1, dim2, dim3 ! Dimensions + real, intent(in) :: ds ! Resolution + real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts + real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points + real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points + real, intent(in) :: u(1:dim1+1,1:dim2,1:dim3) ! v wind + real, intent(in) :: v(1:dim1,1:dim2+1,1:dim3) ! v wind + real, intent(out) :: vor(1:dim1+1,1:dim2+1,1:dim3)! Vorticity + + integer :: i, j, k ! Loop counters + real :: ds_inv ! 1/ds + real :: coeff(1:dim1,1:dim2) ! Coefficient + real :: um(1:dim1+1,1:dim2,1:dim3) ! u-wind copy + real :: vm(1:dim1,1:dim2+1,1:dim3) ! v-wind copy + +!------------------------------------------------------------------------------ +! [1.0] Initialise: +!------------------------------------------------------------------------------ + + vor(:,:,:) = 0.0 + + ds_inv = 1.0 / ds + coeff(:,:) = ds_inv + +!------------------------------------------------------------------------------ +! [2] Calculate scaled u/v field: +!------------------------------------------------------------------------------ + + do k = 1, dim3 + do j = 1, dim2 + do i = 1, dim1+1 + um(i,j,k) = u(i,j,k) / mapfac_u(i,j) + end do + end do + end do + + do k = 1, dim3 + do j = 1, dim2+1 + do i = 1, dim1 + if (mapfac_v(i,j) > 0.000001) then + vm(i,j,k) = v(i,j,k) / mapfac_v(i,j) + else + vm(i,j,k)=0.0 + end if + end do + end do + end do + +!------------------------------------------------------------------------------ +! [4] Calculate vorticity field: +!------------------------------------------------------------------------------ + + do k = 1, dim3 + do j = 2, dim2 + do i = 2, dim1 + vor(i,j,k) = coeff(i,j) * ( vm(i,j,k) - vm(i-1,j,k) - um(i,j,k) + um(i,j-1,k) ) + end do + end do + end do + +! Boundary values (extrapolation): +! Note - not used in Del**2 calculation if overwritten with bcs there). +! vor(1,1:dim2+1) = 2.0 * vor(2,1:dim2+1) - vor(3,1:dim2+1) ! West +! vor(dim1+1,1:dim2+1) = 2.0 * vor(dim1,1:dim2+1) - vor(dim1-1,1:dim2+1) ! East +! vor(1:dim1+1,1) = 2.0 * vor(1:dim1+1,2) - vor(1:dim1+1,3) ! South +! vor(1:dim1+1,dim2+1) = 2.0 * vor(1:dim1+1,dim2) - vor(1:dim1+1,dim2-1) ! North + +! Boundary values (zero gradient): +! Note - not used in Del**2 calculation if overwritten with bcs there). + vor(1,2:dim2,:) = vor(2,2:dim2,:) ! West + vor(dim1+1,2:dim2,:) = vor(dim1,2:dim2,:) ! East + vor(1:dim1+1,1,:) = vor(1:dim1+1,2,:) ! South + vor(1:dim1+1,dim2+1,:) = vor(1:dim1+1,dim2,:) ! North + +end subroutine da_uv_to_vor_c + +subroutine da_del2a_to_a( dim1, dim2, dim3, ds, del2a, a ) +! adapted from var/gen_be/gen_be_stage0_wrf.f90 + + implicit none + + integer, intent(in):: dim1, dim2, dim3 + real, intent(in) :: ds ! grid distance + real, intent(in) :: del2a(1:dim1,1:dim2,1:dim3) ! Del**2 a + real, intent(out) :: a(1:dim1,1:dim2,1:dim3) ! Field a + + integer, parameter :: num_fft_factors = 10 + integer, parameter :: nrange = 1000 + real, parameter :: pi = 3.1415926 + + integer :: n1, n2 ! Padded dimensions (n=dim-1+pad) + integer :: fft_method ! 1=Cosine, 2=Sine transform + integer :: ifax1(1:num_fft_factors) ! FFT factors + integer :: ifax2(1:num_fft_factors) ! FFT factors + real, allocatable :: trigs1(:) ! FFT trig functions + real, allocatable :: trigs2(:) ! FFT trig functions + real, allocatable :: fft_coeffs(:,:) ! FFT coefficients + + integer :: fft_pad1, fft_pad2 ! Range to search for efficient FFT + logical :: found_magic ! True if 2**p 3**p 5**r dimension found + integer :: fft_factors(1:num_fft_factors) ! FFT factors + real :: const ! Multiplicative constant + real :: coeff_nx ! Multiplicative constant + real :: coeff_ny ! Multiplicative constant + real :: cos_coeff_nx ! Multiplicative constant + real :: cos_coeff_ny ! Multiplicative constant + + integer :: i, j, k ! Loop counters + integer :: ij ! 1D array counter + integer :: isign ! -1=Grid>spec, 1=Spec>Grid + integer :: inc ! Stride between data points + integer :: jump ! Increment between start of data vectors + integer :: lot ! Number of data vectors + integer :: n ! n+1 is the length of the data + integer :: work_area ! Dimension of workspace + !real :: a2d(1:n1+1,1:n2+1) ! 2D data array + !real :: a1d(1:(n1+1)*(n2+1)) ! 1D data array + real, allocatable :: a2d(:,:) ! 2D data array + real, allocatable :: a1d(:) ! 1D data array + +! Ensure efficient FFT dimensions by padding if necessary: + n1 = dim1 - 1 + do n = n1, n1 + nrange + call da_find_fft_factors( n, found_magic, fft_factors ) + if ( found_magic .and. mod(n,2) == 0 ) then ! Even magic number found. + fft_pad1 = n - n1 + ifax1 = fft_factors + exit + end if + end do + n1 = n1 + fft_pad1 + + n2 = dim2 - 1 + do n = n2, n2 + nrange + call da_find_fft_factors( n, found_magic, fft_factors ) + if ( found_magic .and. mod(n,2) == 0 ) then ! Even magic number found. + fft_pad2 = n - n2 + ifax2 = fft_factors + exit + end if + end do + n2 = n2 + fft_pad2 + + allocate( trigs1(1:3*n1) ) + allocate( trigs2(1:3*n2) ) + allocate( fft_coeffs(1:n1+1,1:n2+1) ) + + const = -0.5 * ds * ds + coeff_nx = pi / real(n1) + coeff_ny = pi / real(n2) + +! Calculate spectral Del**2 coefficients for C-grid (all pts. except i=j=1): + fft_coeffs(1,1) = 0.0 ! Not used? + do j = 2, n2+1 + cos_coeff_ny = cos(coeff_ny * real(j - 1)) + do i = 1, n1+1 + cos_coeff_nx = cos(coeff_nx * real(i - 1)) + fft_coeffs(i,j) = const / ( 2.0 - cos_coeff_nx - cos_coeff_ny) + end do + end do + j = 1 + cos_coeff_ny = cos(coeff_ny * real(j - 1)) + do i = 2, n1+1 + cos_coeff_nx = cos(coeff_nx * real(i - 1)) + fft_coeffs(i,j) = const / ( 2.0 - cos_coeff_nx - cos_coeff_ny) + end do + + call da_find_fft_trig_funcs( n1, trigs1 ) + call da_find_fft_trig_funcs( n2, trigs2 ) + + fft_method = 2 + + work_area = ( n1 + 1 ) * ( n2 + 1 ) + + allocate ( a2d(1:n1+1,1:n2+1) ) + allocate ( a1d(1:(n1+1)*(n2+1)) ) + +do k = 1, dim3 + +! Fill 2D array structure + do j = 1, dim2 + do i = 1, dim1 + a2d(i,j) = del2a(i,j,k) + end do + +! Fill pad zone (and force b.c.s to satisfy solution type): + if ( fft_method == 1 ) then ! Cosine transform. + a2d(1,j) = a2d(2,j) + do i = dim1, n1+1 + a2d(i,j) = a2d(dim1-1,j) + end do + else if ( fft_method == 2 ) then ! Sine transform: + a2d(1,j) = 0.0 + do i = dim1, n1+1 + a2d(i,j) = 0.0 + end do + end if + end do + + if ( fft_method == 1 ) then ! Cosine transform. + do i = 1, n1+1 + a2d(i,1) = a2d(i,2) + do j = dim2, n2+1 + a2d(i,j) = a2d(i,dim2-1) + end do + end do + else if ( fft_method == 2 ) then ! Sine transform: + do i = 1, n1+1 + a2d(i,1) = 0.0 + do j = dim2, n2+1 + a2d(i,j) = 0.0 + end do + end do + end if + +! Transfer to data array: + do j = 1, n2+1 + do i = 1, n1+1 + ij = (j-1) * (n1+1) + i + a1d(ij) = a2d(i,j) + end do + end do + +!------------------------------------------------------------------------------ +! Perform double fast sine/cosine transform to get spectral del2a: +!------------------------------------------------------------------------------ + + isign = -1 ! Grid to spectral + +! 1st dimension: + inc = 1 ! Stride between data points. + jump = n1+1! Increment between start of data vectors. + lot = n2+1 ! Number of data vectors. + n = n1 ! n+1 is the length of the data. + if ( fft_method == 1 ) then + call fft551( isign, inc, jump, lot, n, & + ifax1, trigs1, a1d, work_area ) + else if ( fft_method == 2 ) then + call fft661( isign, inc, jump, lot, n, & + ifax1, trigs1, a1d, work_area ) + end if + +! 2nd dimension: + inc = n1+1 ! Stride between data points. + jump = 1 ! Increment between start of data vectors. + lot = n1+1 ! Number of data vectors. + n = n2 ! n+1 is the length of the data. + + if ( fft_method == 1 ) then + call fft551( isign, inc, jump, lot, n, & + ifax2, trigs2, a1d, work_area ) + else if ( fft_method == 2 ) then + call fft661( isign, inc, jump, lot, n, & + ifax2, trigs2, a1d, work_area ) + end if + +!------------------------------------------------------------------------------ +! Perform conversion from del2a to a in spectral space: +!------------------------------------------------------------------------------ + +! Note fft_coeffs(1,1)=0 so a(k=0,l=0) is also 0. + do j = 1, n2+1 + do i = 1, n1+1 + ij = (j-1) * (n1+1) + i + a1d(ij) = fft_coeffs(i,j) * a1d(ij) + end do + end do + +!------------------------------------------------------------------------------ +! Perform double fast sine/cosine transform to get gridpoint a: +!------------------------------------------------------------------------------ + + isign = 1 ! Spectral to grid. + +! 1st dimension: + inc = 1 ! Stride between data points. + jump = n1+1! Increment between start of data vectors. + lot = n2+1 ! Number of data vectors. + n = n1 ! n+1 is the length of the data. + + if ( fft_method == 1 ) then + call fft551( isign, inc, jump, lot, n, & + ifax1, trigs1, a1d, work_area ) + else if ( fft_method == 2 ) then + call fft661( isign, inc, jump, lot, n, & + ifax1, trigs1, a1d, work_area ) + end if + +! 2nd dimension: + inc = n1+1 ! Stride between data points. + jump = 1 ! Increment between start of data vectors. + lot = n1+1 ! Number of data vectors. + n = n2 ! n+1 is the length of the data. + + if ( fft_method == 1 ) then + call fft551( isign, inc, jump, lot, n, & + ifax2, trigs2, a1d, work_area ) + else if ( fft_method == 2 ) then + call fft661( isign, inc, jump, lot, n, & + ifax2, trigs2, a1d, work_area ) + end if + +! Transfer grid-point chi to 2D-array (throwing away pad): + do j = 1, dim2 + do i = 1, dim1 + ij = (j-1) * (n1+1) + i + a(i,j,k) = a1d(ij) + end do + end do + +end do ! k loop + + deallocate (a2d) + deallocate (a1d) + deallocate (fft_coeffs) + deallocate (trigs2) + deallocate (trigs1) + +contains + +subroutine da_find_fft_factors(n, n_ok, fft_factors) +! taken from var/da/da_tools/da_find_fft_factors.inc + + !--------------------------------------------------------------------------- + ! Purpose: Calculates prime factors of input number. + !--------------------------------------------------------------------------- + + implicit none + + integer, intent(in) :: n + logical, intent(out) :: n_ok + integer, intent(out) :: fft_factors(:) + + integer :: i, k, l + integer :: nfax, nu, ifac + integer :: jfax(num_fft_factors) + integer :: lfax(7) + + data lfax /6,8,5,4,3,2,1/ + + ! in da_control + !if (trace_use) call da_trace_entry("da_find_fft_factors") + + !--------------------------------------------------------------------------- + ! [1.0] Find factors of vector size (8,6,5,4,3,2; only one 8 allowed): + !--------------------------------------------------------------------------- + + n_ok = .false. + fft_factors(:) = 0 + + ! look for sixes first, store factors in descending order + nu=n + ifac=6 + k=0 + l=1 + +20 continue + + if (mod(nu,ifac).ne.0) goto 30 + + ! 6 is a factor: + k=k+1 + jfax(k)=ifac + if (ifac.ne.8) goto 25 + if (k.eq.1) goto 25 + jfax(1)=8 + jfax(k)=6 + +25 continue + nu=nu/ifac + if (nu.eq.1) goto 50 + if (ifac.ne.8) goto 20 + +30 continue + l=l+1 + ifac=lfax(l) + if (ifac .gt. 1) goto 20 + + ! illegal factors: + ! write (unit=message(1),fmt='(a,i4,a)') 'n = ', n, ' contains illegal factors.' + ! call da_warning(__file__,__line__,message(1:1)) + + goto 9 + + ! now reverse order of factors +50 continue + nfax=k + fft_factors(1)=nfax + do i=1,nfax + fft_factors(nfax+2-i)=jfax(i) + end do + + n_ok = .true. + +9 continue + + !if (trace_use) call da_trace_exit("da_find_fft_factors") + +end subroutine da_find_fft_factors + +subroutine da_find_fft_trig_funcs(n, trig_functs) +! taken from var/da/da_tools/da_find_fft_trig_funcs.inc + + !--------------------------------------------------------------------------- + ! Purpose: Set up constants required for Fourier, sine and cosine transforms + !--------------------------------------------------------------------------- + + implicit none + + integer, intent(in) :: n + real, intent(out) :: trig_functs(:) + + integer :: k, nil, nhl + real :: del, angle + + ! in da_control + ! if (trace_use) call da_trace_entry("da_find_fft_trig_funcs") + + !--------------------------------------------------------------------------- + ! [1.0] Trig functions for real periodic transform: + !--------------------------------------------------------------------------- + + trig_functs(:) = 0.0 + + del=4.0*(pi/2.0)/float(n) + nil=0 + nhl=(n/2)-1 + + do k=nil,nhl + angle=float(k)*del + trig_functs(2*k+1)=cos(angle) + trig_functs(2*k+2)=sin(angle) + end do + + ! [1.1] extra trig functions for cosine transform: + + del=0.5*del + do k=1,nhl + angle=float(k)*del + trig_functs(2*n+k)=sin(angle) + end do + + ! [1.2] extra trig functions for shifted cosine transform: + + del=0.5*del + do k=1,n + angle=float(k)*del + trig_functs(n+k)=sin(angle) + end do + + !if (trace_use) call da_trace_exit("da_find_fft_trig_funcs") + +end subroutine da_find_fft_trig_funcs + +end subroutine da_del2a_to_a + +!DECK FFT551 +! SUBROUTINE 'FFT551' - MULTIPLE FAST COSINE TRANSFORM +! +! AUTHOR: CLIVE TEMPERTON, MAY 1988 +! [ALL-FORTRAN VERSION: C.T., OCTOBER 1995] +! +! COSINE TRANSFORM OF LENGTH N IS CONVERTED TO +! REAL PERIODIC TRANSFORM OF LENGTH N BY PRE- AND POST- +! PROCESSING. REAL PERIODIC TRANSFORM IS PERFORMED BY +! PRUNING REDUNDANT OPERATIONS FROM COMPLEX TRANSFORM. +! +! SEE FOR EXAMPLE PAUL SWARZTRAUBER, "SYMMETRIC FFT'S", +! MATH. COMP. 47 (1986), 323-346. +! +! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA +! WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64) +! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES +! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N +! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' +! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) +! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR +! N+1 IS THE LENGTH OF THE DATA VECTORS +! (WHICH INCLUDE NONZERO VALUES AT BOTH ENDPOINTS) +! LOT IS THE NUMBER OF DATA VECTORS +! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT +! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL +! +! ORDERING OF COEFFICIENTS: Z(0) , Z(1) , Z(2) , ... , Z(N) +! +! ORDERING OF DATA: X(0) , X(1) , X(2) , ... , X(N) +! +! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS +! IN PARALLEL +! +! N MUST BE COMPOSED OF FACTORS 2,3 & 5 AND MUST BE EVEN +! +! DEFINITION OF TRANSFORMS: +! ------------------------- +! +! ISIGN=+1: X(I)=SUM(J=0,...,N)(E(J)*Z(J)*COS(I*J*PI/N)) +! WHERE E(J)=0.5 FOR J=0,N --- ELSE E(J)=1 +! +! ISIGN=-1: Z(J)=(2/N)*SUM(I=0,...,N)(E(I)*X(I)*COS(I*J*PI/N)) +! +! N.B. FFT551 has an unusual definition of the FFTs, +! such that the the coeff of wave0 is NOT the mean. +! +!--------------------------------------------------------------------- +Subroutine FFT551 & ! in + ( ISIGN, & ! in + INC, & ! in + JUMP, & ! in + LOT, & ! in + N, & ! in + IFAX, & ! in + TRIGS, & ! in + A, & ! inout + IDIM ) ! in + +! Code Description: ORIGINAL CODE F77 IS HARDLY TOUCHED !!! + + Integer , intent (in) :: ISIGN ! Switch forward (-1) or inverse (+1) + Integer , intent (in) :: INC ! increment within each data + ! vector (e.g. INC=1 for + ! consecutively stored data) + Integer , intent (in) :: Jump ! increment between start of + ! data vectors + Integer , intent (in) :: LOT ! Number of data vectors + Integer , intent (in) :: N ! N+1 is the length of the data + + Integer , intent (in) :: IFAX(10) ! previously prepared list of + ! factors of N + + Real , intent (in) :: TRIGS(3*N) ! previously prepared list of + ! trigonometric function values + Real , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array ! vectors (which include zeros + ! at the endpoints) + Integer , intent (in) :: IDIM ! dimension workspace + + Real :: WORK(IDIM) ! size (n+1)*min(lot,VectorLength) + Integer :: NFAX,NX,NH + Integer :: NBLOX,NVEX,NB + Integer :: K, IC, J, LA, IGO, JA,JB,IA,IB + Integer :: IFAC,IERR,ISTART + + Real :: CO,S, t1,t2,si,scale, vectorlength + +CHARACTER (LEN=*), PARAMETER :: RoutineName = "Var_FFT551" + + VectorLength = LOT + NFAX=IFAX(1) + NX=N+1 + NH=N/2 + NBLOX=1+(LOT-1)/VectorLength + NVEX=LOT-(NBLOX-1)*VectorLength + ISTART=1 +! + DO 200 NB=1,NBLOX +! +! PREPROCESSING +! ------------- + IA=ISTART + IB=IA+NH*INC + IC=IA+N*INC + JA=1 + JB=NH+1 + IF (MOD(NFAX,2).EQ.1) THEN +!DIR$ IVDEP + DO 105 J=1,NVEX + T1=0.5*(A(IA)+A(IC)) + T2=0.5*(A(IA)-A(IC)) + A(IA)=T1 + A(IC)=T2 + IA=IA+JUMP + IC=IC+JUMP + 105 CONTINUE + ELSE +!DIR$ IVDEP + DO 110 J=1,NVEX + WORK(JA)=0.5*(A(IA)+A(IC)) + WORK(JB)=A(IB) + A(IC)=0.5*(A(IA)-A(IC)) + IA=IA+JUMP + IB=IB+JUMP + IC=IC+JUMP + JA=JA+NX + JB=JB+NX + 110 CONTINUE + ENDIF +! + DO 130 K=1,NH-1 + JA=K+1 + JB=N+1-K + IA=ISTART+K*INC + IB=ISTART+(JB-1)*INC + IC=ISTART+N*INC + SI=TRIGS(2*N+K) + CO=TRIGS(2*N+NH-K) + IF (MOD(NFAX,2).EQ.1) THEN +!DIR$ IVDEP + DO 115 J=1,NVEX + T1 = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB)) + T2 = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB)) + A(IC) = A(IC) + CO*(A(IA)-A(IB)) + A(IA) = T1 + A(IB) = T2 + IA=IA+JUMP + IB=IB+JUMP + IC=IC+JUMP + 115 CONTINUE + ELSE +!DIR$ IVDEP + DO 120 J=1,NVEX + WORK(JA) = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB)) + WORK(JB) = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB)) + A(IC) = A(IC) + CO*(A(IA)-A(IB)) + IA=IA+JUMP + IB=IB+JUMP + IC=IC+JUMP + JA=JA+NX + JB=JB+NX + 120 CONTINUE + ENDIF + 130 CONTINUE +! +! PERIODIC FOURIER ANALYSIS +! ------------------------- + IA=ISTART + LA=N + IGO=1-2*MOD(NFAX,2) +! + DO 140 K=1,NFAX + IFAC=IFAX(NFAX+2-K) + LA=LA/IFAC + IERR=-1 + IF (IGO.EQ.+1) THEN + CALL qpassm(WORK,WORK(IFAC*LA+1),A(IA),A(IA+LA*INC), & + TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) + ELSE IF (IGO.EQ.-1) THEN + CALL qpassm(A(IA),A(IA+IFAC*LA*INC),WORK,WORK(LA+1), & + TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) + ENDIF + IF (IERR.NE.0) GO TO 500 + IGO=-IGO + 140 CONTINUE +! +! POSTPROCESSING +! -------------- + SCALE=2.0 + IF (ISIGN.EQ.+1) SCALE = FLOAT(N) + S=1.0 + IF (ISIGN.EQ.-1) S = 2.0/FLOAT(N) + JA=ISTART + JB=JA+N*INC + IA=1 + IB=N +!DIR$ IVDEP + DO 150 J=1,NVEX + A(JA)=SCALE*WORK(IA) + A(JA+INC)=S*A(JB) + A(JB)=SCALE*WORK(IB) + IA=IA+NX + IB=IB+NX + JA=JA+JUMP + JB=JB+JUMP + 150 CONTINUE +! + DO 170 K=2,N-2,2 + JA=ISTART+K*INC + IA=K +!DIR$ IVDEP + DO 160 J=1,NVEX + A(JA)=SCALE*WORK(IA) + A(JA+INC)=-SCALE*WORK(IA+1)+A(JA-INC) + IA=IA+NX + JA=JA+JUMP + 160 CONTINUE + 170 CONTINUE +! + ISTART=ISTART+NVEX*JUMP + NVEX=VectorLength + 200 CONTINUE + GO TO 570 +! +! ERROR MESSAGES +! -------------- + 500 CONTINUE + GO TO (510,530,550) IERR + 510 CONTINUE + WRITE(UNIT=0,FMT='(A,I4,A)') & + 'VECTOR LENGTH =',NVEX,', GREATER THAN VectorLength' + GO TO 570 + 530 CONTINUE + WRITE(UNIT=0,FMT='(A,I3,A)') & + 'FACTOR =',IFAC,', NOT CATERED FOR' + GO TO 570 + 550 CONTINUE + WRITE(UNIT=0,FMT='(A,I3,A)') & + 'FACTOR =',IFAC,', ONLY CATERED FOR IF LA*IFAC=N' + 570 CONTINUE + + RETURN + END SUBROUTINE FFT551 + + +Subroutine FFT661 & ! in + ( ISIGN, & ! in + INC, & ! in + JUMP, & ! in + LOT, & ! in + N, & ! in + IFAX, & ! in + TRIGS, & ! in + A, & ! inout + DIM ) ! in +! +! +! Description: +! MULTIPLE FAST SINE TRANSFORM +! (Originally called FFT661, then Var_SinTrans) +! author: clive temperton, may 1988 +! (slightly modified for all-fortran version) +! +! Sine transform of length n is converted to +! Real periodic transform of length n by pre- and post- +! processing. Real periodic transform is performed by +! pruning redundant operations from complex transform. +! +! see for example paul swarztrauber, "symmetric fft's", +! math. comp. 47 (1986), 323-346. +! +! Method: +! +! ordering of coefficients: z(0) , z(1) , z(2) , ... , z(n) +! ordering of data: x(0) , x(1) , x(2) , ... , x(n) +! +! vectorization is achieved on cray by doing the transforms +! in parallel +! +! N must be composed of factors 2,3 & 5 and must be even +! +! definition of transforms: +! ------------------------- +! +! isign=+1: x(i)=sum(j=1,...,n-1)(z(j)*sin(i*j*pi/n)) +! +! isign=-1: z(j)=(2/n)*sum(i=1,...,n-1)(x(i)*sin(i*j*pi/n)) +! +! Current Code Owner: Andrew Lorenc +! +! History: +! Version Date Comment +! ------- ---- ------- +! 0.1 14/12/93 Original code. Phil Andrews +! 0.2 16/09/94 Small Modifications for the +! incorporation in the VAR project. HB +! 1.1 21/04/95 placed under control. JB +! 1.2 01/06/95 Tracing added. JB +! +! Code Description: +! NB BECAUSE OF THE TRICKY NESTED LOOPS +! ORIGINAL CODE F77 IS HARDLY TOUCHED !!! + +Implicit none + +! Subroutine arguments + Integer , intent (in) :: ISIGN ! Switch forward (-1) or inverse (+1) + Integer , intent (in) :: INC ! increment within each data + ! vector (e.g. INC=1 for + ! consecutively stored data) + Integer , intent (in) :: Jump ! increment between start of + ! data vectors + Integer , intent (in) :: LOT ! Number of data vectors + Integer , intent (in) :: N ! N+1 is the length of the data + ! vectors (which include zeros + ! at the endpoints) + Integer , intent (in) :: DIM ! dimension workspace + Integer , intent (in) :: IFAX(10) ! previously prepared list of + ! factors of N + + Real , intent (in) :: TRIGS(3*N) ! previously prepared list of + ! trigonometric function values + Real , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array + + ! No descriptions given + Integer :: NFAX,NX,NH + Integer :: NBLOX,NVEX,NB + Integer :: K,JA,JB,IA,IB,IGO,LA,J + Integer :: IFAC,IERR,ISTART + + Real :: SI,T1,T2,SCALE, vectorlength + Real :: WORK(DIM) ! size (n+1)*min(lot,VectorLength) + + VectorLength = LOT + NFAX=IFAX(1) + NX=N+1 + NH=N/2 + NBLOX=1+(LOT-1)/VectorLength + NVEX=LOT-(NBLOX-1)*VectorLength + ISTART=1 +! + DO 200 NB=1,NBLOX +! +! PREPROCESSING +! ------------- + DO 120 K=1,NH-1 + JA=K+1 + JB=N+1-K + IA=ISTART+K*INC + IB=ISTART+(JB-1)*INC + SI=TRIGS(2*N+K) + IF (MOD(NFAX,2).EQ.0) THEN +!DIR$ IVDEP + DO 110 J=1,NVEX + WORK(JA) = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB)) + WORK(JB) = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB)) + IA=IA+JUMP + IB=IB+JUMP + JA=JA+NX + JB=JB+NX + 110 CONTINUE + ELSE +!DIR$ IVDEP + DO 115 J=1,NVEX + T1 = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB)) + T2 = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB)) + A(IA) = T1 + A(IB) = T2 + IA=IA+JUMP + IB=IB+JUMP + 115 CONTINUE + ENDIF + 120 CONTINUE + + JA=1 + JB=NH+1 + IA=ISTART + IB=ISTART+NH*INC + IF (MOD(NFAX,2).EQ.0) THEN +!DIR$ IVDEP + DO 130 J=1,NVEX + WORK(JA)=0.0 + WORK(JB)=2.0*A(IB) + IB=IB+JUMP + JA=JA+NX + JB=JB+NX + 130 CONTINUE + IGO = +1 + ELSE +!DIR$ IVDEP + DO 135 J=1,NVEX + A(IA)=0.0 + A(IB)=2.0*A(IB) + IA=IA+JUMP + IB=IB+JUMP + 135 CONTINUE + IGO = -1 + ENDIF +! +! PERIODIC FOURIER ANALYSIS +! ------------------------- + IA=ISTART + LA=N +! + DO 140 K=1,NFAX + IFAC=IFAX(NFAX+2-K) + LA=LA/IFAC + IERR=-1 + IF (IGO.EQ.+1) THEN + CALL qpassm(WORK,WORK(IFAC*LA+1),A(IA),A(LA*INC+IA), & + TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) + ELSE IF (IGO.EQ.-1) THEN + CALL qpassm(A(IA),A(IFAC*LA*INC+IA),WORK,WORK(LA+1), & + TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) + ENDIF + IF (IERR.NE.0) GO TO 500 + IGO=-IGO + 140 CONTINUE +! +! POSTPROCESSING +! -------------- + SCALE=2.0 + IF (ISIGN.EQ.+1) SCALE = FLOAT(N) + JA=ISTART + JB=JA+N*INC + IA=1 +!DIR$ IVDEP + DO 150 J=1,NVEX + A(JA)=0.0 + A(JA+INC)=0.5*SCALE*WORK(IA) + A(JB)=0.0 + IA=IA+NX + JA=JA+JUMP + JB=JB+JUMP + 150 CONTINUE +! + DO 170 K=2,N-2,2 + JA=ISTART+K*INC + IA=K +!DIR$ IVDEP + DO 160 J=1,NVEX + A(JA)=-SCALE*WORK(IA+1) + A(JA+INC)=SCALE*WORK(IA)+A(JA-INC) + IA=IA+NX + JA=JA+JUMP + 160 CONTINUE + 170 CONTINUE +! + ISTART=ISTART+NVEX*JUMP + NVEX=VectorLength + 200 CONTINUE + + Go To 570 +! +! ERROR MESSAGES +! -------------- + 500 CONTINUE + GO TO (510,530,550) IERR + 510 CONTINUE + WRITE(UNIT=0,FMT='(A,I5,A)') 'NVEX=', NVEX ,' GREATER THAN VectorLength' + GO TO 570 + 530 CONTINUE + WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, 'NOT CATERED FOR' + GO TO 570 + 550 CONTINUE + WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, ' ONLY CATERED FOR IF LA*IFAC=N' + 570 CONTINUE + + + RETURN + + +End subroutine FFT661 + + +!C SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART!C OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE +!C +!C A IS FIRST REAL INPUT VECTOR +!C EQUIVALENCE B(1) WITH A(IFAC*LA*INC1+1) +!C C IS FIRST REAL OUTPUT VECTOR +!C EQUIVALENCE D(1) WITH C(LA*INC2+1) +!C TRIGS IS A PRECALCULATED LIST OF SINES & COSINES +!C INC1 IS THE ADDRESSING INCREMENT FOR A +!C INC2 IS THE ADDRESSING INCREMENT FOR C +!C INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A +!C INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C +!C LOT IS THE NUMBER OF VECTORS +!C N IS THE LENGTH OF THE VECTORS +!C IFAC IS THE CURRENT FACTOR OF N +!C LA = N/(PRODUCT OF FACTORS USED SO FAR) +!C IERR IS AN ERROR INDICATOR: +!C 0 - PASS COMPLETED WITHOUT ERROR +!C 1 - LOT GREATER THAN VectorLength +!C 2 - IFAC NOT CATERED FOR +!C 3 - IFAC ONLY CATERED FOR IF LA=N/IFAC +!C +!C----------------------------------------------------------------------- +!C + SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA,IERR) + + INTEGER :: inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr + + REAL :: a, b, c, d, trigs + DIMENSION A(*),B(*),C(*),D(*),TRIGS(N) +! Local named constants + CHARACTER (LEN=*), PARAMETER :: RoutineName = "QPASSM" +! + REAL :: SIN36, SIN72, QRT5, SIN60 + DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/, & + QRT5/0.559016994374947/,SIN60/0.866025403784437/ + + REAL :: s1, s2, s3, s4, s5 + REAL :: sin45, zsin36, zsin72, zqrt5, zsin60, zsin45, z + REAL :: a0, a1, a2, a3, a4, a5, a6, a10, a11, a20, a21 + REAL :: b0, b1, b2, b3, b4, b5, b6, b10, b11, b20, b21 +! REAL :: c0, c1, c2, c3, c4, c5 + REAL :: c1, c2, c3, c4, c5 + INTEGER :: i, ijk, l, k, kb, m, iink, jink, ijump, kstop + INTEGER :: ibad, igo, ia, ie, je, ibase, jbase, ja, jb, j, ic + INTEGER :: if, jf, kf, ib, jc, kc, id, jd, kd, ke, ig, ih + INTEGER :: vectorlength +! +!- End of header --------------------------------------------------------------- + + + M=N/IFAC + IINK=LA*INC1 + JINK=LA*INC2 + IJUMP=(IFAC-1)*IINK + KSTOP=(N-IFAC)/(2*IFAC) +! + IBAD=1 + VectorLength = lot + IF (LOT.GT.VectorLength) GO TO 910 + IBASE=0 + JBASE=0 + IGO=IFAC-1 + IF (IGO.EQ.7) IGO=6 + IBAD=2 + IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910 + GO TO (200,300,400,500,600,800),IGO +! +! CODING FOR FACTOR 2 +! ------------------- + 200 CONTINUE + IA=1 + IB=IA+IINK + JA=1 + JB=JA+(2*M-LA)*INC2 +! + IF (LA.EQ.M) GO TO 290 +! + DO 220 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 210 IJK=1,LOT + C(JA+J)=A(IA+I)+A(IB+I) + C(JB+J)=A(IA+I)-A(IB+I) + I=I+INC3 + J=J+INC4 + 210 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 220 CONTINUE + JA=JA+JINK + JINK=2*JINK + JB=JB-JINK + IBASE=IBASE+IJUMP + IJUMP=2*IJUMP+IINK + IF (JA.EQ.JB) GO TO 260 + DO 250 K=LA,KSTOP,LA + KB=K+K + C1=TRIGS(KB+1) + S1=TRIGS(KB+2) + JBASE=0 + DO 240 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 230 IJK=1,LOT + C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I)) + C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I)) + D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I) + D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I) + I=I+INC3 + J=J+INC4 + 230 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 240 CONTINUE + IBASE=IBASE+IJUMP + JA=JA+JINK + JB=JB-JINK + 250 CONTINUE + IF (JA.GT.JB) GO TO 900 + 260 CONTINUE + JBASE=0 + DO 280 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 270 IJK=1,LOT + C(JA+J)=A(IA+I) + D(JA+J)=-A(IB+I) + I=I+INC3 + J=J+INC4 + 270 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 280 CONTINUE + GO TO 900 +! + 290 CONTINUE + Z=1.0/FLOAT(N) + DO 294 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 292 IJK=1,LOT + C(JA+J)=Z*(A(IA+I)+A(IB+I)) + C(JB+J)=Z*(A(IA+I)-A(IB+I)) + I=I+INC3 + J=J+INC4 + 292 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 294 CONTINUE + GO TO 900 +! +! CODING FOR FACTOR 3 +! ------------------- + 300 CONTINUE + IA=1 + IB=IA+IINK + IC=IB+IINK + JA=1 + JB=JA+(2*M-LA)*INC2 + JC=JB +! + IF (LA.EQ.M) GO TO 390 +! + DO 320 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 310 IJK=1,LOT + C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) + C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I)) + D(JB+J)=SIN60*(A(IC+I)-A(IB+I)) + I=I+INC3 + J=J+INC4 + 310 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 320 CONTINUE + JA=JA+JINK + JINK=2*JINK + JB=JB+JINK + JC=JC-JINK + IBASE=IBASE+IJUMP + IJUMP=2*IJUMP+IINK + IF (JA.EQ.JC) GO TO 360 + DO 350 K=LA,KSTOP,LA + KB=K+K + KC=KB+KB + C1=TRIGS(KB+1) + S1=TRIGS(KB+2) + C2=TRIGS(KC+1) + S2=TRIGS(KC+2) + JBASE=0 + DO 340 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 330 IJK=1,LOT + A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I)) + B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I)) + A2=A(IA+I)-0.5*A1 + B2=B(IA+I)-0.5*B1 + A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I))) + B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I))) + C(JA+J)=A(IA+I)+A1 + D(JA+J)=B(IA+I)+B1 + C(JB+J)=A2+B3 + D(JB+J)=B2-A3 + C(JC+J)=A2-B3 + D(JC+J)=-(B2+A3) + I=I+INC3 + J=J+INC4 + 330 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 340 CONTINUE + IBASE=IBASE+IJUMP + JA=JA+JINK + JB=JB+JINK + JC=JC-JINK + 350 CONTINUE + IF (JA.GT.JC) GO TO 900 + 360 CONTINUE + JBASE=0 + DO 380 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 370 IJK=1,LOT + C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I)) + D(JA+J)=-SIN60*(A(IB+I)+A(IC+I)) + C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I)) + I=I+INC3 + J=J+INC4 + 370 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 380 CONTINUE + GO TO 900 +! + 390 CONTINUE + Z=1.0/FLOAT(N) + ZSIN60=Z*SIN60 + DO 394 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 392 IJK=1,LOT + C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I))) + C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I))) + D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I)) + I=I+INC3 + J=J+INC4 + 392 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 394 CONTINUE + GO TO 900 +! +! CODING FOR FACTOR 4 +! ------------------- + 400 CONTINUE + IA=1 + IB=IA+IINK + IC=IB+IINK + ID=IC+IINK + JA=1 + JB=JA+(2*M-LA)*INC2 + JC=JB+2*M*INC2 + JD=JB +! + IF (LA.EQ.M) GO TO 490 +! + DO 420 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 410 IJK=1,LOT + C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) + C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) + C(JB+J)=A(IA+I)-A(IC+I) + D(JB+J)=A(ID+I)-A(IB+I) + I=I+INC3 + J=J+INC4 + 410 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 420 CONTINUE + JA=JA+JINK + JINK=2*JINK + JB=JB+JINK + JC=JC-JINK + JD=JD-JINK + IBASE=IBASE+IJUMP + IJUMP=2*IJUMP+IINK + IF (JB.EQ.JC) GO TO 460 + DO 450 K=LA,KSTOP,LA + KB=K+K + KC=KB+KB + KD=KC+KB + C1=TRIGS(KB+1) + S1=TRIGS(KB+2) + C2=TRIGS(KC+1) + S2=TRIGS(KC+2) + C3=TRIGS(KD+1) + S3=TRIGS(KD+2) + JBASE=0 + DO 440 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 430 IJK=1,LOT + A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I)) + A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I)) + A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I)) + A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I)) + B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I)) + B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I)) + B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I)) + B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I)) + C(JA+J)=A0+A1 + C(JC+J)=A0-A1 + D(JA+J)=B0+B1 + D(JC+J)=B1-B0 + C(JB+J)=A2+B3 + C(JD+J)=A2-B3 + D(JB+J)=B2-A3 + D(JD+J)=-(B2+A3) + I=I+INC3 + J=J+INC4 + 430 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 440 CONTINUE + IBASE=IBASE+IJUMP + JA=JA+JINK + JB=JB+JINK + JC=JC-JINK + JD=JD-JINK + 450 CONTINUE + IF (JB.GT.JC) GO TO 900 + 460 CONTINUE + SIN45=SQRT(0.5) + JBASE=0 + DO 480 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 470 IJK=1,LOT + C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I)) + C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I)) + D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) + D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) + I=I+INC3 + J=J+INC4 + 470 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 480 CONTINUE + GO TO 900 +! + 490 CONTINUE + Z=1.0/FLOAT(N) + DO 494 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 492 IJK=1,LOT + C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))) + C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) + C(JB+J)=Z*(A(IA+I)-A(IC+I)) + D(JB+J)=Z*(A(ID+I)-A(IB+I)) + I=I+INC3 + J=J+INC4 + 492 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 494 CONTINUE + GO TO 900 +! +! CODING FOR FACTOR 5 +! ------------------- + 500 CONTINUE + IA=1 + IB=IA+IINK + IC=IB+IINK + ID=IC+IINK + IE=ID+IINK + JA=1 + JB=JA+(2*M-LA)*INC2 + JC=JB+2*M*INC2 + JD=JC + JE=JB +! + IF (LA.EQ.M) GO TO 590 +! + DO 520 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 510 IJK=1,LOT + A1=A(IB+I)+A(IE+I) + A3=A(IB+I)-A(IE+I) + A2=A(IC+I)+A(ID+I) + A4=A(IC+I)-A(ID+I) + A5=A(IA+I)-0.25*(A1+A2) + A6=QRT5*(A1-A2) + C(JA+J)=A(IA+I)+(A1+A2) + C(JB+J)=A5+A6 + C(JC+J)=A5-A6 + D(JB+J)=-SIN72*A3-SIN36*A4 + D(JC+J)=-SIN36*A3+SIN72*A4 + I=I+INC3 + J=J+INC4 + 510 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 520 CONTINUE + JA=JA+JINK + JINK=2*JINK + JB=JB+JINK + JC=JC+JINK + JD=JD-JINK + JE=JE-JINK + IBASE=IBASE+IJUMP + IJUMP=2*IJUMP+IINK + IF (JB.EQ.JD) GO TO 560 + DO 550 K=LA,KSTOP,LA + KB=K+K + KC=KB+KB + KD=KC+KB + KE=KD+KB + C1=TRIGS(KB+1) + S1=TRIGS(KB+2) + C2=TRIGS(KC+1) + S2=TRIGS(KC+2) + C3=TRIGS(KD+1) + S3=TRIGS(KD+2) + C4=TRIGS(KE+1) + S4=TRIGS(KE+2) + JBASE=0 + DO 540 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 530 IJK=1,LOT + A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I)) + A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I)) + A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I)) + A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I)) + B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I)) + B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I)) + B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I)) + B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I)) + A5=A(IA+I)-0.25*(A1+A2) + A6=QRT5*(A1-A2) + B5=B(IA+I)-0.25*(B1+B2) + B6=QRT5*(B1-B2) + A10=A5+A6 + A20=A5-A6 + B10=B5+B6 + B20=B5-B6 + A11=SIN72*B3+SIN36*B4 + A21=SIN36*B3-SIN72*B4 + B11=SIN72*A3+SIN36*A4 + B21=SIN36*A3-SIN72*A4 + C(JA+J)=A(IA+I)+(A1+A2) + C(JB+J)=A10+A11 + C(JE+J)=A10-A11 + C(JC+J)=A20+A21 + C(JD+J)=A20-A21 + D(JA+J)=B(IA+I)+(B1+B2) + D(JB+J)=B10-B11 + D(JE+J)=-(B10+B11) + D(JC+J)=B20-B21 + D(JD+J)=-(B20+B21) + I=I+INC3 + J=J+INC4 + 530 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 540 CONTINUE + IBASE=IBASE+IJUMP + JA=JA+JINK + JB=JB+JINK + JC=JC+JINK + JD=JD-JINK + JE=JE-JINK + 550 CONTINUE + IF (JB.GT.JD) GO TO 900 + 560 CONTINUE + JBASE=0 + DO 580 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 570 IJK=1,LOT + A1=A(IB+I)+A(IE+I) + A3=A(IB+I)-A(IE+I) + A2=A(IC+I)+A(ID+I) + A4=A(IC+I)-A(ID+I) + A5=A(IA+I)+0.25*(A3-A4) + A6=QRT5*(A3+A4) + C(JA+J)=A5+A6 + C(JB+J)=A5-A6 + C(JC+J)=A(IA+I)-(A3-A4) + D(JA+J)=-SIN36*A1-SIN72*A2 + D(JB+J)=-SIN72*A1+SIN36*A2 + I=I+INC3 + J=J+INC4 + 570 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 580 CONTINUE + GO TO 900 +! + 590 CONTINUE + Z=1.0/FLOAT(N) + ZQRT5=Z*QRT5 + ZSIN36=Z*SIN36 + ZSIN72=Z*SIN72 + DO 594 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 592 IJK=1,LOT + A1=A(IB+I)+A(IE+I) + A3=A(IB+I)-A(IE+I) + A2=A(IC+I)+A(ID+I) + A4=A(IC+I)-A(ID+I) + A5=Z*(A(IA+I)-0.25*(A1+A2)) + A6=ZQRT5*(A1-A2) + C(JA+J)=Z*(A(IA+I)+(A1+A2)) + C(JB+J)=A5+A6 + C(JC+J)=A5-A6 + D(JB+J)=-ZSIN72*A3-ZSIN36*A4 + D(JC+J)=-ZSIN36*A3+ZSIN72*A4 + I=I+INC3 + J=J+INC4 + 592 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 594 CONTINUE + GO TO 900 +! +! CODING FOR FACTOR 6 +! ------------------- + 600 CONTINUE + IA=1 + IB=IA+IINK + IC=IB+IINK + ID=IC+IINK + IE=ID+IINK + IF=IE+IINK + JA=1 + JB=JA+(2*M-LA)*INC2 + JC=JB+2*M*INC2 + JD=JC+2*M*INC2 + JE=JC + JF=JB +! + IF (LA.EQ.M) GO TO 690 +! + DO 620 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 610 IJK=1,LOT + A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) + C(JA+J)=(A(IA+I)+A(ID+I))+A11 + C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11) + D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) + A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) + C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11 + D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) + C(JD+J)=(A(IA+I)-A(ID+I))+A11 + I=I+INC3 + J=J+INC4 + 610 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 620 CONTINUE + JA=JA+JINK + JINK=2*JINK + JB=JB+JINK + JC=JC+JINK + JD=JD-JINK + JE=JE-JINK + JF=JF-JINK + IBASE=IBASE+IJUMP + IJUMP=2*IJUMP+IINK + IF (JC.EQ.JD) GO TO 660 + DO 650 K=LA,KSTOP,LA + KB=K+K + KC=KB+KB + KD=KC+KB + KE=KD+KB + KF=KE+KB + C1=TRIGS(KB+1) + S1=TRIGS(KB+2) + C2=TRIGS(KC+1) + S2=TRIGS(KC+2) + C3=TRIGS(KD+1) + S3=TRIGS(KD+2) + C4=TRIGS(KE+1) + S4=TRIGS(KE+2) + C5=TRIGS(KF+1) + S5=TRIGS(KF+2) + JBASE=0 + DO 640 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 630 IJK=1,LOT + A1=C1*A(IB+I)+S1*B(IB+I) + B1=C1*B(IB+I)-S1*A(IB+I) + A2=C2*A(IC+I)+S2*B(IC+I) + B2=C2*B(IC+I)-S2*A(IC+I) + A3=C3*A(ID+I)+S3*B(ID+I) + B3=C3*B(ID+I)-S3*A(ID+I) + A4=C4*A(IE+I)+S4*B(IE+I) + B4=C4*B(IE+I)-S4*A(IE+I) + A5=C5*A(IF+I)+S5*B(IF+I) + B5=C5*B(IF+I)-S5*A(IF+I) + A11=(A2+A5)+(A1+A4) + A20=(A(IA+I)+A3)-0.5*A11 + A21=SIN60*((A2+A5)-(A1+A4)) + B11=(B2+B5)+(B1+B4) + B20=(B(IA+I)+B3)-0.5*B11 + B21=SIN60*((B2+B5)-(B1+B4)) + C(JA+J)=(A(IA+I)+A3)+A11 + D(JA+J)=(B(IA+I)+B3)+B11 + C(JC+J)=A20-B21 + D(JC+J)=A21+B20 + C(JE+J)=A20+B21 + D(JE+J)=A21-B20 + A11=(A2-A5)+(A4-A1) + A20=(A(IA+I)-A3)-0.5*A11 + A21=SIN60*((A4-A1)-(A2-A5)) + B11=(B5-B2)-(B4-B1) + B20=(B3-B(IA+I))-0.5*B11 + B21=SIN60*((B5-B2)+(B4-B1)) + C(JB+J)=A20-B21 + D(JB+J)=A21-B20 + C(JD+J)=A11+(A(IA+I)-A3) + D(JD+J)=B11+(B3-B(IA+I)) + C(JF+J)=A20+B21 + D(JF+J)=A21+B20 + I=I+INC3 + J=J+INC4 + 630 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 640 CONTINUE + IBASE=IBASE+IJUMP + JA=JA+JINK + JB=JB+JINK + JC=JC+JINK + JD=JD-JINK + JE=JE-JINK + JF=JF-JINK + 650 CONTINUE + IF (JC.GT.JD) GO TO 900 + 660 CONTINUE + JBASE=0 + DO 680 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 670 IJK=1,LOT + C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I)) + D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I)) + C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I)) + D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I)) + C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I)) + D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I)) + I=I+INC3 + J=J+INC4 + 670 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 680 CONTINUE + GO TO 900 +! + 690 CONTINUE + Z=1.0/FLOAT(N) + ZSIN60=Z*SIN60 + DO 694 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 692 IJK=1,LOT + A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) + C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11) + C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11) + D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) + A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) + C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11) + D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) + C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11) + I=I+INC3 + J=J+INC4 + 692 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 694 CONTINUE + GO TO 900 +! +! CODING FOR FACTOR 8 +! ------------------- + 800 CONTINUE + IBAD=3 + IF (LA.NE.M) GO TO 910 + IA=1 + IB=IA+IINK + IC=IB+IINK + ID=IC+IINK + IE=ID+IINK + IF=IE+IINK + IG=IF+IINK + IH=IG+IINK + JA=1 + JB=JA+LA*INC2 + JC=JB+2*M*INC2 + JD=JC+2*M*INC2 + JE=JD+2*M*INC2 + Z=1.0/FLOAT(N) + ZSIN45=Z*SQRT(0.5) +! + DO 820 L=1,LA + I=IBASE + J=JBASE +!DIR$ IVDEP + DO 810 IJK=1,LOT + C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+ & + ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) + C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))- & + ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) + C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I))) + D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I))) + C(JB+J)=Z*(A(IA+I)-A(IE+I))+ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) + C(JD+J)=Z*(A(IA+I)-A(IE+I))-ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) + D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))+Z*(A(IG+I)-A(IC+I)) + D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))-Z*(A(IG+I)-A(IC+I)) + I=I+INC3 + J=J+INC4 + 810 CONTINUE + IBASE=IBASE+INC1 + JBASE=JBASE+INC2 + 820 CONTINUE +! + + 900 CONTINUE + IBAD=0 + 910 CONTINUE + IERR=IBAD + + RETURN + END SUBROUTINE QPASSM + diff --git a/var/gen_be_v3/makefile b/var/gen_be_v3/makefile new file mode 100644 index 0000000000..2a19536f49 --- /dev/null +++ b/var/gen_be_v3/makefile @@ -0,0 +1,42 @@ +#!/bin/sh + +#serial +FC = gfortran +CPPFLAG = -DUSE_WRFDA + +#mpp +#FC = mpif90 +#CPPFLAG = -DUSE_WRFDA -DDM_PARALLEL + +#openmp (optional, used only for slen_opt=1) +OMP_FLAG = #-fopenmp +OMP_LIB = #-L/usr/local/lib -lgfortran -lgomp + +#for write_ep = .true. application, the code should be compiled with real-4 +#for gen_be applications, especially cv_options=5, real-4 works, but real-8 is recommended if memory is not an issue +PROMOTION = #-fdefault-real-8 +FCBASEOPTS = -ffree-form -fconvert=big-endian #-ffree-line-length-none +FCDEBUG = #-g -O0 -fbacktrace -ggdb -fcheck=bounds,do,mem,pointer -ffpe-trap=invalid,zero,overflow +FFLAGS = $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS) $(OMP_FLAG) +CPP = cpp -P -traditional + +WRF_SRC_ROOT_DIR = ../.. +LIBS = -L$(WRF_SRC_ROOT_DIR)/external/io_netcdf -lwrfio_nf -L$(WRF_SRC_ROOT_DIR)/var/build -lwrfvar -L$(NETCDF)/lib -lnetcdff -lnetcdf +INCS = -I$(WRF_SRC_ROOT_DIR)/var/build + +OBJS = gen_be_v3.o + +gen_be_v3: ${OBJS} + ${FC} -o gen_be_v3.exe ${FFLAGS} ${OBJS} ${LIBS} $(OMP_LIB) + +.SUFFIXES : .F90 .f .o + +.F90.f : + $(RM) $@ + $(CPP) $(CPPFLAG) $*.F90 > $@ + +.f.o : + ${FC} ${FFLAGS} ${INCS} -c $*.f + +clean: + rm -f *.o *.exe diff --git a/var/gen_be_v3/util/combine_be_cv5.f90 b/var/gen_be_v3/util/combine_be_cv5.f90 new file mode 100644 index 0000000000..503fedbb5c --- /dev/null +++ b/var/gen_be_v3/util/combine_be_cv5.f90 @@ -0,0 +1,296 @@ +program be_readwrite + +! gfortran -o combine_be_cv5.exe -fconvert=big-endian combine_be_cv5.f90 + +! input: be_[PSI|CHI_U|T_U|RH|PSFC_U].dat +! regcoeff[1|2|3].dat +! output: be.dat (in real(8) and big_endian) + + implicit none + + integer :: ni, nj, nk ! dimensions read in + integer :: bin_type ! type of bin to average over + integer :: num_bins ! number of 3D bins + integer :: num_bins2d ! number of 2D bins + integer, allocatable :: bin(:,:,:) ! bin assigned to each 3D point + integer, allocatable :: bin2d(:,:) ! bin assigned to each 2D point + +!for be_[PSI|CHI_U|T_U|RH|PSFC_U].dat in real(4) when gen_be_v3 is run with write_be_dat_r8=.false. + !real, allocatable :: e_vec(:,:) ! domain-averaged eigenvectors + !real, allocatable :: e_val(:) ! domain-averaged eigenvalues + !real, allocatable :: e_vec_loc(:,:,:) ! latitudinally varying eigenvectors + !real, allocatable :: e_val_loc(:,:) ! latitudinally varying eigenvalues + !real, allocatable :: scale_length(:) ! scale length for regional application + !real, allocatable :: regcoeff1(:) ! psi/chi regcoeff + !real, allocatable :: regcoeff2(:,:) ! psi/ps regcoeff + !real, allocatable :: regcoeff3(:,:,:) ! psi/t regcoeff +!for be_[PSI|CHI_U|T_U|RH|PSFC_U].dat in real(8) when gen_be_v3 is run with write_be_dat_r8=.true.(default) + real*8, allocatable :: e_vec(:,:) ! domain-averaged eigenvectors + real*8, allocatable :: e_val(:) ! domain-averaged eigenvalues + real*8, allocatable :: e_vec_loc(:,:,:) ! latitudinally varying eigenvectors + real*8, allocatable :: e_val_loc(:,:) ! latitudinally varying eigenvalues + real*8, allocatable :: scale_length(:) ! scale length for regional application + real*8, allocatable :: regcoeff1(:) ! psi/chi regcoeff + real*8, allocatable :: regcoeff2(:,:) ! psi/ps regcoeff + real*8, allocatable :: regcoeff3(:,:,:) ! psi/t regcoeff + + real*8, allocatable :: e_vec_r8(:,:) ! domain-averaged eigenvectors + real*8, allocatable :: e_val_r8(:) ! domain-averaged eigenvalues + real*8, allocatable :: e_vec_loc_r8(:,:,:) ! latitudinally varying eigenvectors + real*8, allocatable :: e_val_loc_r8(:,:) ! latitudinally varying eigenvalues + real*8, allocatable :: scale_length_r8(:) ! scale length for regional application + + real*8, allocatable :: regcoeff1_r8(:) ! psi/chi regcoeff + real*8, allocatable :: regcoeff2_r8(:,:) ! psi/ps regcoeff + real*8, allocatable :: regcoeff3_r8(:,:,:) ! psi/t regcoeff + + character(len=256) :: infile ! input filename + character(len=256) :: outfile ! output filename + logical :: isfile + integer :: iunit, ounit + + integer, parameter :: nvar = 5 + character(len=10) :: varnames(nvar) = (/ "PSI ", & + "CHI_U ", & + "T_U ", & + "RH ", & + "PSFC_U " /) + character(len=10) :: varout(nvar) = (/ "psi ", & + "chi_u ", & + "t_u ", & + "rh ", & + "ps_u " /) + character(len=10) :: this_var + integer :: var_dim + integer :: i + integer :: nk_2d = 1 + + real*8 :: lat_min, lat_max, binwidth_lat + real*8 :: hgt_min, hgt_max, binwidth_hgt + + lat_min = 0.0 + lat_max = 0.0 + binwidth_lat = 0.0 + hgt_min = 0.0 + hgt_max = 0.0 + binwidth_hgt = 0.0 + + iunit = 81 + ounit = 82 + + outfile = 'be.dat' + open (ounit, file=trim(outfile), form='unformatted', status='replace') + + ! read and write dimension and bin info + + infile = 'be_'//trim(varnames(1))//'.dat' + inquire(file=trim(infile), exist=isfile) + if ( .not. isfile ) then + write(*,*) 'STOP: ', trim(infile), ' does not exist.' + stop + end if + open (iunit, file=trim(infile), form='unformatted', status='old') + write(*,*) 'Reading from ', trim(infile), ' for dimension and bin info' + + read(iunit) ni, nj, nk + + allocate( bin(1:ni,1:nj,1:nk) ) + allocate( bin2d(1:ni,1:nj) ) + + read(iunit) bin_type + read(iunit) num_bins, num_bins2d + read(iunit) bin + read(iunit) bin2d + close(iunit) + + write(ounit) ni, nj, nk + write(ounit) bin_type + write(ounit) lat_min, lat_max, binwidth_lat + write(ounit) hgt_min, hgt_max, binwidth_hgt + write(ounit) num_bins, num_bins2d + write(ounit) bin + write(ounit) bin2d + + ! read and write regcoeff info + + allocate( regcoeff1(1:num_bins) ) + allocate( regcoeff2(1:nk,1:num_bins2d) ) + allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) ) + allocate( regcoeff1_r8(1:num_bins) ) + allocate( regcoeff2_r8(1:nk,1:num_bins2d) ) + allocate( regcoeff3_r8(1:nk,1:nk,1:num_bins2d) ) + + infile = 'regcoeff1.dat' + inquire(file=trim(infile), exist=isfile) + if ( .not. isfile ) then + write(*,*) 'STOP: ', trim(infile), ' does not exist.' + stop + end if + open (iunit, file=trim(infile), form='unformatted', status='old') + write(*,*) 'Reading from ', trim(infile) + read(iunit) ni, nj, nk + read(iunit) num_bins, num_bins2d + read(iunit) regcoeff1 + regcoeff1_r8(:) = regcoeff1(:) + write(ounit) regcoeff1_r8 + close(iunit) + + infile = 'regcoeff2.dat' + inquire(file=trim(infile), exist=isfile) + if ( .not. isfile ) then + write(*,*) 'STOP: ', trim(infile), ' does not exist.' + stop + end if + open (iunit, file=trim(infile), form='unformatted', status='old') + write(*,*) 'Reading from ', trim(infile) + read(iunit) ni, nj, nk + read(iunit) num_bins, num_bins2d + read(iunit) regcoeff2 + regcoeff2_r8(:,:) = regcoeff2(:,:) + write(ounit) regcoeff2_r8 + close(iunit) + + infile = 'regcoeff3.dat' + inquire(file=trim(infile), exist=isfile) + if ( .not. isfile ) then + write(*,*) 'STOP: ', trim(infile), ' does not exist.' + stop + end if + open (iunit, file=trim(infile), form='unformatted', status='old') + write(*,*) 'Reading from ', trim(infile) + read(iunit) ni, nj, nk + read(iunit) num_bins, num_bins2d + read(iunit) regcoeff3 + regcoeff3_r8(:,:,:) = regcoeff3(:,:,:) + write(ounit) regcoeff3_r8 + close(iunit) + + ! done writing dimension, bin and regcoeff info + + deallocate( regcoeff1 ) + deallocate( regcoeff2 ) + deallocate( regcoeff3 ) + deallocate( regcoeff1_r8 ) + deallocate( regcoeff2_r8 ) + deallocate( regcoeff3_r8 ) + + allocate( e_vec(1:nk,1:nk) ) + allocate( e_val(1:nk) ) + allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) ) + allocate( e_val_loc(1:nk,1:num_bins2d) ) + allocate( e_vec_r8(1:nk,1:nk) ) + allocate( e_val_r8(1:nk) ) + allocate( e_vec_loc_r8(1:nk,1:nk,1:num_bins2d) ) + allocate( e_val_loc_r8(1:nk,1:num_bins2d) ) + + ! loop over individual files to gather and write out e_vec/e_val info + + do i = 1, nvar + + infile = 'be_'//trim(varnames(i))//'.dat' + inquire(file=trim(infile), exist=isfile) + if ( .not. isfile ) then + write(*,*) 'STOP: ', trim(infile), ' does not exist.' + stop + end if + open (iunit, file=trim(infile), form='unformatted', status='old') + write(*,*) 'Reading from ', trim(infile) + + read(iunit) ni, nj, nk + + read(iunit) bin_type + read(iunit) num_bins, num_bins2d + read(iunit) bin + read(iunit) bin2d + + read(iunit) this_var + write(ounit) varout(i) + + read(iunit) var_dim + if ( var_dim == 3 ) then + read(iunit) e_vec + read(iunit) e_val + read(iunit) e_vec_loc + read(iunit) e_val_loc + e_vec_r8 = e_vec + e_val_r8 = e_val + e_vec_loc_r8 = e_vec_loc + e_val_loc_r8 = e_val_loc + write(ounit) nk, num_bins2d + write(ounit) e_vec_r8 + write(ounit) e_val_r8 + write(ounit) e_vec_loc_r8 + write(ounit) e_val_loc_r8 + else if ( var_dim == 2 ) then + read(iunit) e_vec(1:1,1:1) + read(iunit) e_val(1:1) + read(iunit) e_vec_loc(1:1,1:1,:) + read(iunit) e_val_loc(1:1,:) + e_vec_r8 = e_vec + e_val_r8 = e_val + e_vec_loc_r8 = e_vec_loc + e_val_loc_r8 = e_val_loc + write(ounit) nk_2d, num_bins2d + write(ounit) e_vec_r8(1:1,1:1) + write(ounit) e_val_r8(1:1) + write(ounit) e_vec_loc_r8(1:1,1:1,:) + write(ounit) e_val_loc_r8(1:1,:) + end if + close(iunit) + end do ! var loop for e_val/e_vec + + ! read individual files again to gather and write out lengthscales + + allocate( scale_length(1:nk) ) + allocate( scale_length_r8(1:nk) ) + + do i = 1, nvar + infile = 'be_'//trim(varnames(i))//'.dat' + open (iunit, file=trim(infile), form='unformatted', status='old') + !write(*,*) 'Reading from ', trim(infile) + read(iunit) ni, nj, nk + read(iunit) bin_type + read(iunit) num_bins, num_bins2d + read(iunit) bin + read(iunit) bin2d + + read(iunit) this_var + write(ounit) varout(i) + + read(iunit) var_dim + if ( var_dim == 3 ) then + read(iunit) e_vec + read(iunit) e_val + read(iunit) e_vec_loc + read(iunit) e_val_loc + read(iunit) scale_length + scale_length_r8 = scale_length + write(ounit) scale_length_r8 + else if ( var_dim == 2 ) then + read(iunit) e_vec(1:1,1:1) + read(iunit) e_val(1:1) + read(iunit) e_vec_loc(1:1,1:1,:) + read(iunit) e_val_loc(1:1,:) + read(iunit) scale_length(1:1) + scale_length_r8 = scale_length + write(ounit) scale_length_r8(1:1) + end if + close(iunit) + end do + close(ounit) + + deallocate( e_vec ) + deallocate( e_val ) + deallocate( e_vec_loc ) + deallocate( e_val_loc ) + deallocate( e_vec_r8 ) + deallocate( e_val_r8 ) + deallocate( e_vec_loc_r8 ) + deallocate( e_val_loc_r8 ) + + deallocate (scale_length) + deallocate (scale_length_r8) + + write(*,*) 'Done writing be.dat for cv_options=5' + +end program be_readwrite diff --git a/var/gen_be_v3/util/makefile b/var/gen_be_v3/util/makefile new file mode 100644 index 0000000000..37eb30dc94 --- /dev/null +++ b/var/gen_be_v3/util/makefile @@ -0,0 +1,26 @@ +#!/bin/sh + +FC = gfortran +FFLAGS = -fconvert=big-endian + +all: combine_be_cv5.exe combine_be_cv7.exe combine_be_cv7_ccv2.exe + +combine_be_cv5.exe: combine_be_cv5.o + rm -f $@ + ${FC} -o $@ $(FFLAGS) combine_be_cv5.o + +combine_be_cv7.exe: combine_be_cv7.o + rm -f $@ + ${FC} -o $@ $(FFLAGS) combine_be_cv7.o + +combine_be_cv7_ccv2.exe: combine_be_cv7_ccv2.o + rm -f $@ + ${FC} -o $@ $(FFLAGS) combine_be_cv7_ccv2.o + +.SUFFIXES : .f90 .o + +.f90.o : + $(FC) $(FFLAGS) -c $*.f90 + +clean: + rm -f *.o *.exe