program add_zeros include 'cgcam.h' parameter(n_prec=4) real(n_prec), allocatable :: u(:,:) integer file_size character(20) arg, fname_in, fname_out character( 6) ext logical file_exists, extension_read file_exists = .false. extension_read = .false. n = 1 n_args = command_argument_count() do while( n .le. n_args ) call get_command_argument( n, arg, len_arg, ierr ) if( scan(arg(1:1),'0123456789') .eq. 0 ) then print *, 'Bad input' goto 11 else read(arg(1:len_arg),*) i_ext if( i_ext .le. 0 ) then print *, 'Bad input' goto 11 else extension_read = .true. end if end if n = n + 1 end do goto 12 11 print *, 'USAGE: add_zeros [file_extension]' print *, ' ' stop 12 continue c *** Prompt if the file extension was not read on the command line. if( .not. extension_read ) then print *, 'Enter the vel file extension' read(5,*) i_ext end if write(ext,'(i6.6)') i_ext fname_in = 'vel_short.'//ext inquire( file=fname_in, exist=file_exists ) if( .not. file_exists ) then print *, 'ERROR: file ',trim(fname_in),' does not exist' stop end if c *** Read input parameters call input( .false., ierr ) if( ierr .ne. 0 ) stop Lu = 5 + n_scalar Lb = 5 + n_scalar Lp = Lb + 4 call set_dims( Nx, Ny, Nz, Nx, Ny, Nz, nzp, Lx, Ly, & Lz, Lzp, Nxp1, Nxp2, Nyp1, Nyp2, Nzp1, Nzp2 ) c *** Open vel file and determine size differential. inquire( file=fname_in, size=file_size ) Lz_cut = file_size/(n_prec*Lx*Ly*Lu) L_add = Lz - Lz_cut if( Lz_cut .lt. 1 .or. Lz_cut .ge. Lz ) then print *, 'ERROR: Cannot determine Lz_cut' print *, 'file_size, Lz_cut, Lz', file_size, Lz_cut, Lz stop end if c print *, 'file_size, Lz_cut, Lz, L_add', c & file_size, Lz_cut, Lz, L_add open( unit=2, file=fname_in, form='unformatted', & access='stream', action='read') fname_out = 'vel.'//ext call delete_file(fname_out) open( unit=3, file=fname_out, form='unformatted', & access='stream', action='write') c *** Write new vel file with zeros inserted in the correct locations. allocate( u(Lx,Ly) ) do n=1, Lu do k=1, Lz_cut read(2) u write(3) u end do u = 0.0 do k=1, L_add write(3) u end do end do close(2) close(3) deallocate( u ) end subroutine input( write_params, ierr ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com include 'cgcam.h' character( 8) default_flag character(80) fname logical write_params, file_exists call parse_command_line( fname, ierr ) if( ierr .ne. 0 ) return call set_labels( labels, values, required, fixed, L_params, & n_inputs_r, n_inputs_i, n_inputs, & n_dynpar_r, n_dynpar_i, n_dynpar, n_params, ierr) if( ierr .ne. 0 ) return call read_params( fname, labels, values, found, n_inputs, & mult_restart, ierr ) if( ierr .ne. 0 ) return call assign_params( loc(yL2), loc(Nx), values, & n_inputs_r, n_inputs_i ) c *** Check for bad inputs call check_inputs( ierr ) if( ierr .ne. 0 ) return call check_static_dims( ierr ) if( ierr .ne. 0 ) return c *** Set parameters that depend on the inputs call set_dep_params() c *** Write the input parameters to file params.out if( write_params ) then open(newunit=i_unit,file='params.out') do i=1, n_inputs default_flag = ' ' if( .not. found(i) ) default_flag = 'DEFAULT' write(i_unit,*) labels(i), values(i), default_flag end do close(i_unit) end if return end subroutine assign_params( ipr, ipi, values, n_real, n_int ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com pointer(ipr, rdat) pointer(ipi, idat) real rdat(n_real), values(n_real+n_int) integer idat(n_int) rdat(1:n_real) = values(1:n_real) idat(1:n_int ) = values(n_real+1:n_real+n_int) return end subroutine check_inputs( ierr ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com c *** Error codes: c 0 - normal exit c 6 - required input not found c 7 - static dimensions too small c 9 - inconsistent input parameter values include 'cgcam.h' logical match character(20), dimension(7) :: programs character( 1) ext, newline ierr = 0 newline = achar(10) c *** Parse the i_sponge variable i_spongex1 = 0; i_spongex2 = 0; i_spongey1 = 0; i_spongey2 = 0 i_spongez1 = 0; i_spongez2 = 0 is = i_sponge if( float(is)*0.00001 .ge. 1.0 ) then is = is - 100000*(is/100000) i_spongex1 = 1 end if if( float(is)*0.0001 .ge. 1.0 ) then i_spongex2 = 1 is = is - 10000*(is/10000) end if if( float(is)*0.001 .ge. 1.0 ) then i_spongey1 = 1 is = is - 1000*(is/1000) end if if( float(is)*0.01 .ge. 1.0 ) then i_spongey2 = 1 is = is - 100*(is/100) end if if( float(is)*0.1 .ge. 1.0 ) then i_spongez1 = 1 is = is - 10*(is/10) end if if( float(is)*1.0 .ge. 1.0 ) then i_spongez2 = 1 end if c *** Make sure the static dimensions are sufficient if( Nx .gt. mLx .or. Ny .gt. mLy .or. Nz .gt. mLz ) then print *, 'ERROR: static dimensions are too small' print *, 'mLx, mLy, mLz = ', mLx, mLy, mLz, newline ierr = 7 return end if c *** Check for conflicting wind inputs if( i_read_w .eq. 1 .or. i_read_w .eq. 2 ) then if( i_wind .ne. 0 ) then print *, 'ERROR: input variable i_wind can not be set ', & 'when i_read_w=1,2', newline ierr = 9 return end if i = index_param( 'Uo', labels, n_inputs ) if( found(i) ) then print *, 'ERROR: input variable Uo can not be set when ', & 'i_read_w=1,2', newline ierr = 9 return end if end if if( i_read_w .eq. 2 ) then if( i_grid .ne. 3 ) then print *, 'ERROR: i_read_w=2 currently ONLY supported ', & 'for terrain grids', newline ierr = 9 end if if( i_backg_flux .ne. 1 ) then print *, 'ERROR: i_backg_flux must be removed (=1) for ', & 'i_read_w=2 spacially variable winds!', newline ierr = 9 end if if( i_damp_wind .eq. 1 ) then if( i_force .ne. 4 ) then print *, 'ERROR: i_force 4 (hyperbolic tangent) is the ', & 'only supported forcing', newline ierr = 9 end if if( 4*t_force/dt_ramp .lt. 1 ) then print *, 'ERROR: dt_ramp must less than 4*t_force!', & newline ierr = 9 end if if( t_adjust_wind .lt. (4*t_force + dt_ramp) .and. & int(t_adjust_wind/dt_ramp)*dt_ramp.ne.t_adjust_wind)then print *, 'ERROR: dt_ramp must be factor of t_adjust_wind' & , ' < 4*t_force!', newline ierr = 9 end if if( i_cfl .eq. 1 ) then print *, 'WARNING: Fixed cfl ramping does NOT like ', & 'restarting during ramp', newline else if( int( float(n_skip_v)*dt0/dt_ramp )*dt_ramp .ne. & float(n_skip_v)*dt0 ) then print *, 'ERROR: dt_ramp should be a factor of Volumes ', & 'for restart.' print *, ' * Volumes ', n_skip_v*dt0 print *, ' * dt_ramp ', dt_ramp, newline ierr = 9 end if end if end if c *** If this is a gravity wave problem check inputs for consistency if( i_initial .eq. 4 .or. i_initial .eq. 5 .or. & i_initial .eq. 10 .or. i_initial .eq. 11 .or. & i_initial .eq. 12 ) then i0 = index_param( 'lambda_x', labels, n_inputs ) i1 = index_param( 'Gam' , labels, n_inputs ) i2 = index_param( 'omega' , labels, n_inputs ) i3 = index_param( 'Uo' , labels, n_inputs ) if( .not. found(i0) ) then print *, 'ERROR: lambda_x is a required input for ', & 'a GW problem', newline ierr = 9 return end if if( i_read_w .eq. 1 .and. found(i1) .and. found(i2) ) then print *, 'ERROR: You can not specify both Gam and omega ', & 'when i_read_w=1', newline ierr = 9 return end if if( i_wind .ne. 0 .and. found(i1) .and. found(i2) ) then print *, 'ERROR: You can not specify both Gam and omega ', & 'when i_wind != 0', newline ierr = 9 return end if if( found(i1) .and. found(i2) .and. found(i3) ) then print *, 'ERROR: You can only specify two of Gam, ', & 'omega, Uo for a GW problem', newline ierr = 9 return end if if( (i_wind .ne. 0 .or. i_read_w .eq. 1) .and. & .not.found(i1) .and. .not.found(i2) ) then print *, 'ERROR: You must specify either Gam or omega ', & 'whenn i_wind != 0 or i_read_w=1 for a GW ', & 'problem', newline ierr = 9 return end if if( i_wind .eq. 0 .and. i_read_w .eq. 0 .and. & (found(i1) .and. .not.found(i2) .and. .not.found(i3)) .or. & (.not.found(i1) .and. found(i2) .and. .not.found(i3)) .or. & (.not.found(i1) .and. .not.found(i2) .and. found(i3))) then print *, 'ERROR: You must specify exactly two of Gam, ', & 'omega, Uo for a GW problem', newline ierr = 9 return end if end if c *** Check decomposition parameters. id = index_param( 'i_decomp' , labels, n_inputs ) iny = index_param( 'n_pencil_y', labels, n_inputs ) inz = index_param( 'n_pencil_z', labels, n_inputs ) if( i_decomp .lt. 1 .or. i_decomp .gt. 3 ) then print *, 'ERROR: decomposition ', i_decomp, & ' not supported!',newline ierr = 9 end if if( num_procs .eq. 1 ) goto 11 if( found(iny) .or. found(inz) ) then print *, 'WARNING: n_pencil_y and/or n_pencil_z specified.' if( i_decomp .eq. 1 ) then print *, 'These inputs will be recomputed!' end if print *, '' end if programs = [ character(20) :: "cgcam", "vel2vtk", "tst_decomp", & "compare_vel", "convert_filt", "tst_procs", & "remesh_zones" ] if( i_decomp .ge. 2 ) then match = .false. do n=1, size(programs) if( trim(command_name) .eq. trim(programs(n)) ) match=.true. end do if( .not. match ) then print *, 'WARNING: Pencils not implemented for ', & trim(command_name), ', reverting...', newline i_decomp = 1 values(id) = i_decomp end if end if if( i_decomp .eq. 3 ) then if( i_grid .ne. 1 .and. i_grid .ne. 3 ) then print *, 'ERROR: Block decomposition not implemented for ', & 'i_grid =', i_grid, newline ierr = 9 end if n_count = 1 do i=1, size(k_interface) if( k_interface(i) .lt. 0 .or. k_interface(i) .gt. Nz ) then print *, 'ERROR: Interface plane index out of range.' print *, 'k_interface =', k_interface(i), 'Nz =', Nz ierr = 9 end if if( k_interface(i) .gt. 0 ) then if( int(Nx/(2**i))*(2**i) .ne. Nx ) then print *, 'ERROR: Nx is not divisable at interface ', i ierr = 9 end if if( int(Ny/(2**i))*(2**i) .ne. Ny ) then print *, 'ERROR: Ny is not divisable at interface ', i ierr = 9 end if n_count = n_count + 1 end if end do if( n_count .gt. num_procs ) then print *, 'ERROR: Not enough processors for each interface!' ierr = 9 end if else if( k_interface(1) .ne. 0 ) then print *, 'ERROR: Interface index supplied without block', & ' decomposition.', newline ierr = 9 end if if( ierr .ne. 0 ) return 11 continue c *** Stretching of z_bak not currently allowed if( i_stretch_z_bak .ne. 0 ) then print *, 'ERROR: Stretching is not currently allowed for ', & 'z_bak. SEE set_z_bak.', newline ierr = 9 end if c *** Check atmospheric composition inputs and set default inputs if( i_atmos_var .eq. 1 ) then i_g = index_param( 'grav', labels, n_inputs ) i_R = index_param( 'R_gas', labels, n_inputs ) i_gam = index_param( 'gamma_gas', labels, n_inputs ) required(i_g) = .false. required(i_R) = .false. required(i_gam) = .false. if( found(i_g) .or. found(i_R) .or. found(i_gam) ) then print *, 'WARNING: Ignored inputs for variable ', & 'atmospheric composition:' print *, 'grav, R_gas, gamma_gas', newline end if i = index_param( 'R_earth', labels, n_inputs ) required(i) = .true. end if c *** Perform an initial check for missing required inputs do i=1, n_inputs if( required(i) .and. .not.found(i) ) goto 50 end do c *** Promote selected optional inputs to required based on c *** other read-in parameter values. if( n_les .eq. 1 ) then i = index_param( 'c_smag', labels, n_inputs ) required(i) = .true. end if if( i_cfl .eq. 0 ) then i = index_param( 'dt0', labels, n_inputs ) required(i) = .true. end if if( i_restart .ne. 1 ) then i = index_param( 'i_initial', labels, n_inputs ) required(i) = .true. end if if( i_read_b .eq. 0 .and. i_read_w .ne. 2 ) then i = index_param( 'i_background', labels, n_inputs ) required(i) = .true. end if if( i_read_b .eq. 0 .and. i_background .eq. 1 ) then i = index_param( 'lapse', labels, n_inputs ) required(i) = .true. end if if( i_read_b .eq. 0 .and. i_background .eq. 2 ) then i = index_param( 'lapse_pot', labels, n_inputs ) required(i) = .true. end if if( i_force .gt. 0 ) then i = index_param( 't_force', labels, n_inputs ) required(i) = .true. end if if( i_diss .ne. 0 ) then i = index_param( 'eps_diss', labels, n_inputs ) required(i) = .true. end if if( n_skip_s .ne. 0 ) then i = index_param( 'i_beg_sub', labels, n_inputs ) required(i) = .true. i = index_param( 'i_end_sub', labels, n_inputs ) required(i) = .true. i = index_param( 'j_beg_sub', labels, n_inputs ) required(i) = .true. i = index_param( 'j_end_sub', labels, n_inputs ) required(i) = .true. i = index_param( 'k_beg_sub', labels, n_inputs ) required(i) = .true. i = index_param( 'k_end_sub', labels, n_inputs ) required(i) = .true. end if do n=1, n_ts write(ext,'(i1.1)') n i = index_param( 'i_dir_ts'//ext, labels, n_inputs ) required(i) = .true. i = index_param( 'pos_ts'//ext, labels, n_inputs ) required(i) = .true. i = index_param( 'i_beg_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 1 ) required(i) = .true. i = index_param( 'i_end_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 1 ) required(i) = .true. i = index_param( 'j_beg_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 2 ) required(i) = .true. i = index_param( 'j_end_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 2 ) required(i) = .true. i = index_param( 'k_beg_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 3 ) required(i) = .true. i = index_param( 'k_end_ts'//ext, labels, n_inputs ) if( i_dir_ts(n) .ne. 3 ) required(i) = .true. end do if( i_restart .eq. 1 ) then i = index_param( 'nt_restart', labels, n_inputs ) required(i) = .true. end if if( iubc_z1 .eq. 15 ) then i = index_param( 'i_body_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'i_body_end', labels, n_inputs ) required(i) = .true. end if if( i_spongex1 .ne. 0 ) then i = index_param( 'w_spongex1', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongex1', labels, n_inputs ) required(i) = .true. end if if( i_spongex2 .ne. 0 ) then i = index_param( 'w_spongex2', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongex2', labels, n_inputs ) required(i) = .true. end if if( i_spongey1 .ne. 0 ) then i = index_param( 'w_spongey1', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongey1', labels, n_inputs ) required(i) = .true. end if if( i_spongey2 .ne. 0 ) then i = index_param( 'w_spongey2', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongey2', labels, n_inputs ) required(i) = .true. end if if( i_spongez1 .ne. 0 ) then i = index_param( 'w_spongez1', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongez1', labels, n_inputs ) required(i) = .true. end if if( i_spongez2 .ne. 0 ) then i = index_param( 'w_spongez2', labels, n_inputs ) required(i) = .true. i = index_param( 't_spongez2', labels, n_inputs ) required(i) = .true. end if if( i_grid .eq. 2 ) then i = index_param( 'yl', labels, n_inputs ) required(i) = .true. end if if( n_les .eq. 2 ) then i = index_param( 'i_dm_avg', labels, n_inputs ) required(i) = .true. end if if( n_les .eq. 2 .and. i_dm_avg .eq. 12 ) then i = index_param( 'idm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'idm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'jdm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'jdm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'kdm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'kdm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'n_dm_avg_x2', labels, n_inputs ) required(i) = .true. i = index_param( 'n_dm_avg_y2', labels, n_inputs ) required(i) = .true. end if if( n_les .eq. 2 .and. i_dm_avg .eq. 13 ) then i = index_param( 'idm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'idm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'jdm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'jdm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'kdm_beg', labels, n_inputs ) required(i) = .true. i = index_param( 'kdm_end', labels, n_inputs ) required(i) = .true. i = index_param( 'n_dm_avg_x2', labels, n_inputs ) required(i) = .true. i = index_param( 'n_dm_avg_y2', labels, n_inputs ) required(i) = .true. i = index_param( 'n_dm_avg_z2', labels, n_inputs ) required(i) = .true. end if if( i_damp_wind .ne. 0 ) then i = index_param( 'zc_damp_wind', labels, n_inputs ) required(i) = .true. i = index_param( 'zw_damp_wind', labels, n_inputs ) required(i) = .true. end if if( i_initial .eq. 8 ) then i = index_param( 'duct_cent_z', labels, n_inputs ) required(i) = .true. i = index_param( 'duct_width_z', labels, n_inputs ) required(i) = .true. i = index_param( 'duct_N2', labels, n_inputs ) required(i) = .true. end if if( i_initial .gt. 11 .and. i_initial .le. 12 ) then i = index_param( 'packet_width_x', labels, n_inputs ) required(i) = .true. end if if( i_initial .eq. 12 ) then i = index_param( 'packet_width_y', labels, n_inputs ) required(i) = .true. end if if( i_coriolis .ne. 0 ) then i1 = index_param( 'R_earth', labels, n_inputs ) i2 = index_param( 'omega_earth', labels, n_inputs ) i3 = index_param( 'lat_ref', labels, n_inputs ) i4 = index_param( 'x_ref', labels, n_inputs ) i5 = index_param( 'y_ref', labels, n_inputs ) required(i1) = .true. required(i2) = .true. required(i3) = .true. required(i4) = .true. required(i5) = .true. end if do i=1, size(z_gauss_avg) write(ext,'(i1.1)') i i0 = index_param( 'z_gauss_avg' //ext, labels, n_inputs ) i1 = index_param( 'w_gauss_FWHM'//ext, labels, n_inputs ) if( found(i0) .or. found(i1) ) then required(i0) = .true. required(i1) = .true. end if end do c *** Check again for missing required inputs do i=1, n_inputs if( required(i) .and. .not. found(i) ) goto 50 end do c *** Set default values for [ijk]_sub_end ii = index_param( 'i_end_sub', labels, n_inputs ) ij = index_param( 'j_end_sub', labels, n_inputs ) ik = index_param( 'k_end_sub', labels, n_inputs ) if( .not. found(ii) ) then i_end_sub = Nx+1 values(ii) = Nx+1 end if if( .not. found(ij) ) then j_end_sub = Ny+1 values(ij) = Ny+1 end if if( .not. found(ik) ) then k_end_sub = Nz+1 values(ik) = Nz+1 end if c *** Check for inconsistent inputs. if( To .lt. 0.0 .or. rho_o .lt. 0.0 ) then print*, 'ERROR: To and rho_o must be positive' print *, 'To = ', To print *, 'rho_o = ', rho_o, newline ierr = 9 end if if( iubc_x1 .lt. 0 .or. iubc_x1 .eq. 2 .or. & iubc_x1 .eq. 3 .or. iubc_x1 .eq. 7 .or. & iubc_x1 .ge. 10 ) then print *, 'ERROR: invalid iubc_x1 specified, iubc_x1 = ', & iubc_x1, newline ierr = 9 end if if( iubc_x2 .lt. 0 .or. iubc_x2 .eq. 2 .or. & iubc_x2 .eq. 3 .or. iubc_x2 .eq. 7 .or. & iubc_x2 .ge. 10 ) then print *, 'ERROR: invalid iubc_x2 specified, iubc_x2 = ', & iubc_x2, newline ierr = 9 end if if( (iubc_x1 .eq. 5 .and. iubc_x2 .ne. 5) .or. & (iubc_x2 .eq. 5 .and. iubc_x1 .ne. 5) ) then print *, 'Error: only one iubc_x flag set to 5 ' print *, 'iubc_x1, iubc_x2 = ', iubc_x1, iubc_x2, newline ierr = 9 end if if( (iubc_x1 .eq. 4 .and. iubc_x2 .ne. 4) .or. & (iubc_x2 .eq. 4 .and. iubc_x1 .ne. 4) ) then print *, 'Error: only one iubc_x flag set to 4 ' print *, 'iubc_x1, iubc_x2 = ', iubc_x1, iubc_x2, newline ierr = 9 end if if( iubc_y1 .lt. 0 .or. iubc_y1 .eq. 2 .or. & iubc_y1 .eq. 3 .or. iubc_y1 .eq. 4 .or. & iubc_y1 .eq. 7 .or. iubc_y1 .ge. 10 ) then print *, 'ERROR: invalid iubc_y1 specified, iubc_y1 = ', & iubc_y1, newline ierr = 9 end if if( iubc_y2 .lt. 0 .or. iubc_y2 .eq. 2 .or. & iubc_y2 .eq. 3 .or. iubc_y2 .eq. 4 .or. & iubc_y2 .eq. 7 .or. iubc_y2 .ge. 10 ) then print *, 'ERROR: invalid iubc_y2 specified, iubc_y2 = ', & iubc_y2, newline ierr = 9 end if if( (iubc_y1 .eq. 5 .and. iubc_y2 .ne. 5) .or. & (iubc_y2 .eq. 5 .and. iubc_y1 .ne. 5) ) then print *, 'Error: only one iubc_y flag set to 5 ' print *, 'iubc_y1, iubc_y2 = ', iubc_y1, iubc_y2, newline ierr = 9 end if select case( iubc_z1 ) case( :-1, 4, 10:14, 16:22, 24: ) print *, 'ERROR: invalid iubc_z1 specified, iubc_z1 = ', & iubc_z1, newline ierr = 9 end select select case( iubc_z2 ) case( :-1, 10:14, 16:22, 24: ) print *, 'ERROR: invalid iubc_z2 specified, iubc_z2 = ', & iubc_z2, newline ierr = 9 end select if( itbc_x1 .lt. 0 .or. itbc_x1 .eq. 2 .or. & itbc_x1 .eq. 3 .or. itbc_x1 .eq. 7 .or. & itbc_x1 .ge. 10 ) then print *, 'ERROR: invalid itbc_x1 specified, itbc_x1 = ', & itbc_x1, newline ierr = 9 end if if( itbc_x2 .lt. 0 .or. itbc_x2 .eq. 2 .or. & itbc_x2 .eq. 3 .or. itbc_x2 .eq. 7 .or. & itbc_x2 .ge. 10 ) then print *, 'ERROR: invalid itbc_x2 specified, itbc_x2 = ', & itbc_x2, newline ierr = 9 end if if( (itbc_x1 .eq. 5 .and. itbc_x2 .ne. 5) .or. & (itbc_x2 .eq. 5 .and. itbc_x1 .ne. 5) ) then print *, 'Error: only one itbc_x flag set to 5 ' print *, 'itbc_x1, itbc_x2 = ', itbc_x1, itbc_x2, newline ierr = 9 end if if( (itbc_x1 .eq. 4 .and. itbc_x2 .ne. 4) .or. & (itbc_x2 .eq. 4 .and. itbc_x1 .ne. 4) ) then print *, 'Error: only one itbc_x flag set to 4 ' print *, 'itbc_x1, itbc_x2 = ', itbc_x1, itbc_x2, newline ierr = 9 end if if( itbc_y1 .lt. 0 .or. itbc_y1 .eq. 2 .or. & itbc_y1 .eq. 3 .or. itbc_y1 .eq. 4 .or. & itbc_y1 .eq. 7 .or. itbc_y1 .ge. 10 ) then print *, 'ERROR: invalid itbc_y1 specified, itbc_y1 = ', & itbc_y1, newline ierr = 9 end if if( itbc_y2 .lt. 0 .or. itbc_y2 .eq. 2 .or. & itbc_y2 .eq. 3 .or. itbc_y2 .eq. 4 .or. & itbc_y2 .eq. 7 .or. itbc_y2 .ge. 10 ) then print *, 'ERROR: invalid itbc_y2 specified, itbc_y2 = ', & itbc_y2, newline ierr = 9 end if if( (itbc_y1 .eq. 5 .and. itbc_y2 .ne. 5) .or. & (itbc_y2 .eq. 5 .and. itbc_y1 .ne. 5) ) then print *, 'Error: only one itbc_y flag set to 5 ' print *, 'itbc_y1, itbc_y2 = ', itbc_y1, itbc_y2, newline ierr = 9 end if if( (itbc_y1 .eq. 4 .and. itbc_y2 .ne. 4) .or. & (itbc_y2 .eq. 4 .and. itbc_y1 .ne. 4) ) then print *, 'Error: only one itbc_y flag set to 4 ' print *, 'itbc_y1, itbc_y2 = ', itbc_y1, itbc_y2, newline ierr = 9 end if if( (iubc_x1 .eq. 5 .and. itbc_x1 .ne. 5) .or. & (itbc_x1 .eq. 5 .and. iubc_x1 .ne. 5) ) then print *, 'Error: onlx one iubc_x itbc_x flag set to 5 ' print *, 'iubc_x1, itbc_x1 = ', iubc_x1, itbc_x1, newline ierr = 9 end if if( (iubc_x1 .eq. 4 .and. itbc_x1 .ne. 4) .or. & (itbc_x1 .eq. 4 .and. iubc_x1 .ne. 4) ) then print *, 'Error: onlx one iubc_x itbc_x flag set to 4 ' print *, 'iubc_x1, itbc_x1 = ', iubc_x1, itbc_x1, newline ierr = 9 end if if( (iubc_y1 .eq. 5 .and. itbc_y1 .ne. 5) .or. & (itbc_y1 .eq. 5 .and. iubc_y1 .ne. 5) ) then print *, 'Error: only one iubc_y itbc_y flag set to 5 ' print *, 'iubc_y1, itbc_y1 = ', iubc_y1, itbc_y1, newline ierr = 9 end if if( (iubc_y1 .eq. 4 .and. itbc_y1 .ne. 4) .or. & (itbc_y1 .eq. 4 .and. iubc_y1 .ne. 4) ) then print *, 'Error: only one iubc_y itbc_y flag set to 4 ' print *, 'iubc_y1, itbc_y1 = ', iubc_y1, itbc_y1, newline ierr = 9 end if if( (iubc_z1 .eq. 15 .and. itbc_z1 .ne. 15) .or. & (itbc_z1 .eq. 15 .and. iubc_z1 .ne. 15) ) then print *, 'Error: only one iubc_z itbc_z flag set to 15 ' print *, 'iubc_z1, itbc_z1 = ', iubc_z1, itbc_z1, newline ierr = 9 end if if( i_beg_sub .lt. 0 .or. i_beg_sub .gt. Nx+2 ) then print *, 'bad value for parameter i_beg_sub' print *, 'i_beg_sub = ', i_beg_sub, newline ierr = 9 end if if( i_end_sub .lt. 0 .or. i_end_sub .gt. Nx+2 .or. & i_end_sub .lt. i_beg_sub ) then print *, 'bad value for parameter i_end_sub' print *, 'i_end_sub = ', i_end_sub, newline ierr = 9 end if if( j_beg_sub .lt. 0 .or. j_beg_sub .gt. Ny+2 ) then print *, 'bad value for parameter j_beg_sub' print *, 'j_beg_sub = ', j_beg_sub, newline ierr = 9 end if if( j_end_sub .lt. 0 .or. j_end_sub .gt. Ny+2 .or. & j_end_sub .lt. j_beg_sub ) then print *, 'bad value for parameter j_end_sub' print *, 'j_end_sub = ', j_end_sub, newline ierr = 9 end if if( k_beg_sub .lt. 0 .or. k_beg_sub .gt. Nz+2 ) then print *, 'bad value for parameter k_beg_sub' print *, 'k_beg_sub = ', k_beg_sub, newline ierr = 9 end if if( k_end_sub .lt. 0 .or. k_end_sub .gt. Nz+2 .or. & k_end_sub .lt. k_beg_sub ) then print *, 'bad value for parameter k_end_sub' print *, 'k_end_sub = ', k_end_sub, newline ierr = 9 end if do n=1, n_ts if( i_dir_ts(n) .ne. 1 .and. & ( i_beg_ts(n) .lt. 1 .or. i_beg_ts(n) .gt. Nx+2 ) ) then print *, 'bad value for parameter i_beg_ts(',n,')' print *, 'i_beg_ts = ', i_beg_ts(n), newline ierr = 9 end if if( i_dir_ts(n) .ne. 1 .and. & ( i_end_ts(n) .lt. 1 .or. i_end_ts(n) .gt. Nx+2 ) ) then print *, 'bad value for parameter i_end_ts(',n,')' print *, 'i_end_ts = ', i_end_ts(n), newline ierr = 9 end if if( i_dir_ts(n) .ne. 1 .and. & ( j_beg_ts(n) .lt. 1 .or. j_beg_ts(n) .gt. Ny+2 ) ) then print *, 'bad value for parameter j_beg_ts(',n,')' print *, 'j_beg_ts = ', j_beg_ts(n), newline ierr = 9 end if if( i_dir_ts(n) .ne. 1 .and. & ( j_end_ts(n) .lt. 1 .or. j_end_ts(n) .gt. Ny+2 ) ) then print *, 'bad value for parameter j_end_ts(',n,')' print *, 'j_end_ts = ', j_end_ts(n), newline ierr = 9 end if if( i_dir_ts(n) .ne. 1 .and. & ( k_beg_ts(n) .lt. 1 .or. k_beg_ts(n) .gt. Nz+2 ) ) then print *, 'bad value for parameter k_beg_ts(',n,')' print *, 'k_beg_ts = ', k_beg_ts(n), newline ierr = 9 end if if( i_dir_ts(n) .ne. 1 .and. & ( k_end_ts(n) .lt. 1 .or. k_end_ts(n) .gt. Nz+2 ) ) then print *, 'bad value for parameter k_end_ts(',n,')' print *, 'k_end_ts = ', k_end_ts(n), newline ierr = 9 end if end do if( grav .eq. 0.0 .and. i_hydrostat .ne. 0 ) then print *, 'ERROR: grav=0 but i_hydrostat!=0', newline ierr = 9 end if if( i_dm_avg .eq. 12 ) then idm_beg1 = idm_beg-n_dm_avg_x2 idm_end1 = idm_end+n_dm_avg_x2 jdm_beg1 = jdm_beg-n_dm_avg_y2 jdm_end1 = jdm_end+n_dm_avg_y2 if( idm_beg1 .lt. 2 .or. idm_end1 .gt. Nx+1 .or. & jdm_beg1 .lt. 2 .or. jdm_end1 .gt. Ny+1 ) then print *, 'Error: for dynamic_model local average ', & 'stencil extends beyond domain boundary', newline ierr = 9 end if end if if( n_scalar .lt. 0 .or. n_scalar .gt. 3 ) then print *, 'ERROR: n_scalar<0 or n_scalar>3', newline ierr = 9 end if if( n_scalar .gt. 1 ) then print *, 'ERROR: n_scalar>1 is not yet implemented', newline ierr = 9 end if if( n_scalar .ne. 0 .and. i_read_w .eq. 2 ) then print *, 'ERROR: n_scalar not yet implemented for i_read_w=2', & newline ierr = 9 end if if( grav .ne. 0.0 .and. grav_s .ne. 0.0 ) then print *, 'ERROR: only one of grav, grav_s should be nonzero', & newline ierr = 9 end if if( i_damp_wind .eq. 1 .and. i_force .eq. 0 ) then print *, 'ERROR: i_damp_wind=1 but i_force=0', newline ierr = 9 end if i = index_param( 'eps_filt', labels, n_inputs ) if( found(i) .and. i_restart .eq. 0 .and. & trim(command_name) .eq. 'cgcam' ) then print *, '' print *, 'WARNING - Fixed filtering is no longer supported ', & 'so the eps_filt input ' print *, 'will be ignored. See' print *, 'http://boulder.gats-inc.com/cgcam/', & 'source/user/filter.html', newline end if if( n_les .eq. 1 ) then print *, 'ERROR: The fixed-constant Smagorinsky model ', & '(n_les=1) is not implemented.' print *, 'Use the dynamic model (n_les=2) instead.', newline ierr = 9 end if return 50 print *, 'ERROR: Required input ',trim(labels(i)), & ' was not provided', newline ierr = 6 return 60 print *, 'ERROR: Invalid value for parameter ',trim(labels(i)), & ' =', values(i), newline ierr = 9 return return end subroutine check_static_dims( ierr ) include 'cgcam.h' ierr = 0 if( i_decomp .eq. 1 .and. num_procs .gt. Nz ) then ierr = 1 print *, 'ERROR: ncpu exceeds Nz' print *, 'Reduce the number of processes and resubmit' end if if( num_procs .gt. maxp ) then ierr = 1 print *, 'ERROR: ncpu exceeds maxp' print *, 'Edit cgcam.h and recompile' end if if( Nx+2 .gt. mLx ) then ierr = 1 print *, 'ERROR: Lx exceeds mLx' print *, 'Edit cgcam.h and recompile' end if if( Ny+2 .gt. mLy ) then ierr = 1 print *, 'ERROR: Ly exceeds mLy' print *, 'Edit cgcam.h and recompile' end if if( Nz+2 .gt. mLz ) then ierr = 1 print *, 'ERROR: Lz exceeds mLz' print *, 'Edit cgcam.h and recompile' end if return end subroutine parse_command_line( fname, ierr ) c *** Currently the only command line argument to cgcam is the input c *** file name. It can be specified as 'cgcam -f input_file' or simply c *** 'cgcam input_file'. If no input file is specified, we check c *** for the presence of either cgcam.inp or input.dat in that order. include 'cgcam.h' character( *) fname character(80) arg logical input_file_specified, file_exists ierr = 0 input_file_specified = .false. c *** Get and store the name of the main program. call get_command_argument( 0, command_name ) ind = index(command_name,'/',.true.) if( ind .ne. 0 ) then command_name = command_name(ind+1:) end if c *** Now read and parse the list of command arguments. n_args = command_argument_count() n = 1 do while( n .le. n_args ) call get_command_argument( n, arg, len_arg, ierr ) if( arg .eq. '-i' ) then n = n + 1 if( n .le. n_args ) then call get_command_argument( n, fname, len_fname, ierr ) input_file_specified = .true. end if else if( trim(command_name) .eq. 'cgcam' .and. & n .eq. n_args .and. & scan(arg(1:1),'./abcdefghijklmnopqrstuvwxyz'// & 'ABCDEFGHIJKLMNOPQRSTUVWXYZ') .eq. 1 ) then fname = arg len_fname = len_arg input_file_specified = .true. else c goto 10 ! Ignore other arguments for the time being end if n = n + 1 end do goto 12 10 print *, ' ' print *, 'Bad argument. ', & 'Usage: cgcam [-f] [input_file]' ierr = 1 return 12 continue if( input_file_specified ) then inquire(file=fname,exist=file_exists) if(file_exists) then return else print *, 'ERROR: The input file that you specified on the ', & 'command line does not exist' print *, 'The specified file was ',trim(fname) ierr = 1 return end if end if c *** Check for the existence of one of the default input file names. fname = 'cgcam.inp' inquire(file=fname,exist=file_exists) if(file_exists) return fname = 'input.dat' inquire(file=fname,exist=file_exists) if(file_exists) return print *, 'ERROR: Neither default input file cgcam.inp nor', & ' input.dat exist' ierr = 1 return end subroutine read_params( fname, labels, values, found, n_vars, & mult_restart, ierr ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com real values(n_vars) character(* ) fname character(18) labels(n_vars), label_lower, label logical found(n_vars), done found = .false. done = .false. open(newunit=i_unit,file=fname) do while( .not. done ) call read_line( i_unit, label, value, mult, done, ierr ) if( ierr .ne. 0 ) return if( label .eq. 'nt_restart' ) mult_restart = mult do i=1, n_vars call to_lower( labels(i), label_lower, 18 ) if( label .eq. label_lower ) then found( i) = .true. values(i) = value goto 20 end if end do print *, 'Warning: extraneous variable ', trim(label), & ' found in file ', trim(fname) 20 continue end do close(i_unit) return end subroutine set_dep_params( ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com include 'cgcam.h' c *** Set default planes output positions i = index_param( 'n_skip_p', labels, n_inputs ) if( found(i) .and. k_xy_plane(1) .eq. 0 .and. & Nx .gt. 1 .and. Ny .gt. 1 ) then k_xy_plane(1) = (Nz+3)/2 ii = index_param( 'k_xy_plane1', labels, n_inputs ) values(ii) = k_xy_plane(1) end if if( found(i) .and. j_xz_plane(1) .eq. 0 .and. & Nx .gt. 1 .and. Nz .gt. 1 ) then j_xz_plane(1) = (Ny+3)/2 ii = index_param( 'j_xz_plane1', labels, n_inputs ) values(ii) = j_xz_plane(1) end if if( found(i) .and. i_yz_plane(1) .eq. 0 .and. & Ny .gt. 1 .and. Nz .gt. 1 ) then i_yz_plane(1) = (Nx+3)/2 ii = index_param( 'i_yz_plane1', labels, n_inputs ) values(ii) = i_yz_plane(1) end if c *** Determine the number of planes to be written do i=1, size(k_xy_plane) if( k_xy_plane(i) .eq. 0 ) goto 1 end do 1 n_xy_planes = i-1 do i=1, size(j_xz_plane) if( j_xz_plane(i) .eq. 0 ) goto 2 end do 2 n_xz_planes = i-1 do i=1, size(i_yz_plane) if( i_yz_plane(i) .eq. 0 ) goto 3 end do 3 n_yz_planes = i-1 do i=1, size(z_gauss_avg) if( z_gauss_avg(i) .eq. 0 ) goto 4 end do 4 n_xy_averages = i-1 if( n_skip_t .ne. 0 ) then do n=1, 9 if( i_dir_ts(n) .eq. 0 ) goto 5 end do 5 n_ts = n - 1 else n_skip_t = 1000000000 ii = index_param( 'n_skip_t', labels, n_inputs ) values(ii) = n_skip_t end if c *** Determine number of interfaces do i=1, size(k_interface) if( k_interface(i) .eq. 0 ) goto 6 end do 6 n_interfaces = i-1 c *** Set the gravity wave type, with the default being 0 (not a GW problem). i_gw_type = 0 if( i_initial .eq. 4 .or. i_initial .eq. 5 .or. & i_initial .eq. 10 .or. i_initial .eq. 11 .or. & i_initial .eq. 12 ) then i0 = index_param( 'lambda_x', labels, n_inputs ) i1 = index_param( 'Gam' , labels, n_inputs ) i2 = index_param( 'omega' , labels, n_inputs ) i3 = index_param( 'Uo' , labels, n_inputs ) if( found(i1) .and. found(i2) ) then i_gw_type = 1 else if( found(i1) .and. found(i3) ) then i_gw_type = 2 else if( found(i1) .and. i_read_w .eq. 1 ) then i_gw_type = 2 else if( found(i2) .and. found(i3) ) then i_gw_type = 3 end if end if c *** Parse the i_sponge variable i_spongex1 = 0; i_spongex2 = 0; i_spongey1 = 0; i_spongey2 = 0 i_spongez1 = 0; i_spongez2 = 0 is = i_sponge if( float(is)*0.00001 .ge. 1.0 ) then is = is - 100000*(is/100000) i_spongex1 = 1 end if if( float(is)*0.0001 .ge. 1.0 ) then i_spongex2 = 1 is = is - 10000*(is/10000) end if if( float(is)*0.001 .ge. 1.0 ) then i_spongey1 = 1 is = is - 1000*(is/1000) end if if( float(is)*0.01 .ge. 1.0 ) then i_spongey2 = 1 is = is - 100*(is/100) end if if( float(is)*0.1 .ge. 1.0 ) then i_spongez1 = 1 is = is - 10*(is/10) end if if( float(is)*1.0 .ge. 1.0 ) then i_spongez2 = 1 end if c *** Parse the i_beam variable beam_amp = 0.0 ib = i_beam if( float(ib)*0.001 .ge. 1.0 ) then ib = ib - 1000*(ib/1000) beam_amp(1) = 1.0 end if if( float(ib)*0.01 .ge. 1.0 ) then beam_amp(2) = 1.0 ib = ib - 100*(ib/100) end if if( float(ib)*0.1 .ge. 1.0 ) then beam_amp(3) = 1.0 ib = ib - 10*(ib/10) end if if( float(ib)*1.0 .ge. 1.0 ) then beam_amp(4) = 1.0 end if c *** Set the read_background and outside_wind logical variables read_background = .false. if( i_read_b .eq. 1 ) read_background = .true. if( i_read_w .eq. 2 ) then outside_wind = .false. first_cross = .false. end if c *** Set the i_periodic variables i_periodic_y = 0 i_periodic_z = 0 if( iubc_y1 .eq. 5 ) i_periodic_y = 1 if( iubc_z1 .eq. 5 ) i_periodic_z = 1 c *** Set the skip factors to large values if they were not provided if( n_skip_h .le. 0 ) then n_skip_h = 1000000000 ii = index_param( 'n_skip_h', labels, n_inputs ) values(ii) = n_skip_h end if if( n_skip_p .le. 0 ) then n_skip_p = 1000000000 ii = index_param( 'n_skip_p', labels, n_inputs ) values(ii) = n_skip_p end if if( n_skip_v .le. 0 ) then n_skip_v = 1000000000 ii = index_param( 'n_skip_v', labels, n_inputs ) values(ii) = n_skip_v end if if( n_skip_s .le. 0 ) then n_skip_s = 1000000000 ii = index_param( 'n_skip_s', labels, n_inputs ) values(ii) = n_skip_s end if if( n_skip_b .le. 0 ) then n_skip_b = 1000000000 ii = index_param( 'n_skip_b', labels, n_inputs ) values(ii) = n_skip_b end if c *** Set default values for the branch cut location if( iubc_z1 .ne. 15 ) then i_body_beg = 2 ii = index_param( 'i_body_beg', labels, n_inputs ) values(ii) = i_body_beg i_body_end = Nx+1 ii = index_param( 'i_body_end', labels, n_inputs ) values(ii) = i_body_end end if c *** Set the entire domain as the default range for the dynamic model if( n_les .eq. 2 ) then if( idm_end .eq. 0 ) then idm_end = Nx+1 ii = index_param( 'idm_end', labels, n_inputs ) values(ii) = idm_end end if if( jdm_end .eq. 0 ) then jdm_end = Ny+1 ii = index_param( 'jdm_end', labels, n_inputs ) values(ii) = jdm_end end if if( kdm_end .eq. 0 ) then kdm_end = Nz+1 ii = index_param( 'kdm_end', labels, n_inputs ) values(ii) = kdm_end end if end if n_shear = 0 i = index_param( 'i_shear_dir1', labels, n_inputs ) if( found(i) ) n_shear = 1 i = index_param( 'i_shear_dir2', labels, n_inputs ) if( found(i) ) n_shear = 2 i = index_param( 'i_shear_dir3', labels, n_inputs ) if( found(i) ) n_shear = 3 c *** Set noise shaping parameters i = index_param( 'z_noise_width', labels, n_inputs ) if( found(i) ) shape_noise = .true. return end subroutine set_labels( lab, val, req, fix, Np, & n_inputs_r, n_inputs_i, n_inputs, & n_dynpar_r, n_dynpar_i, n_dynpar, & n_params, ierr ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com character( 2) ext2 character(18) lab(Np), label real val(Np) logical req(Np), fix(Np), F, done c *** Input parameters and their attributes are defined here. The logical c *** variable fix is true if the parameter should remain fixed upon restart c *** and req is true for a parameter requiring an explicit input value. c *** By default both attributes are true. val is the default value of an c *** optional parameter not listed in the input file. Real-values c *** parameters must be listed first, followed by integers. There are c *** two parameter blocks, one for inputs (static), and a second for c *** run-time (dynamic). Each block is divided into real and integer lists. ierr = 0 F= .false. req(1:Np) = .true. fix(1:Np) = .true. val(1:Np) = 0.0 c *** Static input file parameters. i = 1 c-------------------------------List reals first------------------------ lab(i)='yL2' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='z0' ; val(i)=0.0 ; i=i+1 lab(i)='z1_fade' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='dt0' ; fix(i)=F; req(i)=F; val(i)=0.01 ; i=i+1 lab(i)='cfl0' ; fix(i)=F; req(i)=F; val(i)=0.8 ; i=i+1 lab(i)='Uo' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='Vo' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='rho_o' ; i=i+1 lab(i)='To' ; i=i+1 lab(i)='grav' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='grav_s' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='vis' ; fix(i)=F; ; i=i+1 lab(i)='Pr' ; fix(i)=F; val(i)=1.0 ; i=i+1 lab(i)='gamma_gas' ; val(i)=1.4 ; i=i+1 lab(i)='R_gas' ; val(i)=287.0 ; i=i+1 lab(i)='lapse' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='lapse_pot' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='flct_u' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='flct_t' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='char_L' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='zc_damp_wind' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='zw_damp_wind' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='dt_ramp' ; req(i)=F; val(i)=1.0e+9; i=i+1 lab(i)='t_adjust_wind' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='U_wind0' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='U_wind1' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='wind_width_z' ; req(i)=F; val(i)=1.0e+9; i=i+1 lab(i)='wind_cent_z' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='lambda_x' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='Gam' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='omega' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='amplitude' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_cent_x' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_cent_y' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_cent_z' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_width_x' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_width_y' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='packet_width_z' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_cent_x' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_cent_y' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_cent_z' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_width_x' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_width_y' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='convct_width_z' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='beam_force_len' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='beam_fade_len' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='beam_cent_x' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='beam_cent_y' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='beam_cent_z' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='duct_cent_z' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='duct_width_z' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='duct_N2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='eps_filt' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='filt0' ; fix(i)=F; req(i)=F; val(i)=0.02 ; i=i+1 lab(i)='filt_inc' ; fix(i)=F; req(i)=F; val(i)=0.02 ; i=i+1 lab(i)='filt_dec' ; fix(i)=F; req(i)=F; val(i)=0.02 ; i=i+1 lab(i)='osc_max' ; fix(i)=F; req(i)=F; val(i)=0.02 ; i=i+1 lab(i)='eps_diss' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_force' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='c_smag' ; fix(i)=F; req(i)=F; val(i)=0.01 ; i=i+1 lab(i)='rot_coords' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='pos_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+9 lab(i)='w_spongex1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongex1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_spongex2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongex2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_spongey1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongey1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_spongey2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongey2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_spongez1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongez1' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_spongez2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_spongez2' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='t_initial' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='shear_vel(3)' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+3 lab(i)='shear_cent(3)' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+3 lab(i)='shear_width(3)' ; fix(i)=F; req(i)=F; val(i)=1.0+9 ; i=i+3 lab(i)='z_gauss_avg1' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='z_gauss_avg2' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_gauss_FWHM1' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='w_gauss_FWHM2' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='lat_ref' ; req(i)=F; val(i)=40.0 ; i=i+1 lab(i)='lon_ref' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='x_ref' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='y_ref' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='R_earth' ; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='omega_earth' ; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='z_noise_cent' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='z_noise_width' ; fix(i)=F; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='pert_amp' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='pert_x_len' ; fix(i)=F; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='pert_z_cent' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='pert_z_fwhm' ; fix(i)=F; req(i)=F; val(i)=1.0 ; i=i+1 lab(i)='res_max' ; fix(i)=F; req(i)=F; val(i)=1.0e+3; i=i+1 lab(i)='res_min' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 lab(i)='z_smooth_max' ; fix(i)=F; req(i)=F; val(i)=0.0 ; i=i+1 n_inputs_r=i-1 c-------------------------------Integers below this line---------------- lab(i)='Nx' ; i=i+1 lab(i)='Ny' ; i=i+1 lab(i)='Nz' ; i=i+1 lab(i)='i_decomp' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_pencil_y' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_pencil_z' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='k_interface(9)' ; req(i)=F; val(i)=0 ; i=i+9 lab(i)='i_grid' ; val(i)=0 ; i=i+1 lab(i)='i_body_beg' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_body_end' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_scalar' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_restart' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='nt_restart' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='nt_scale' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_les' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_steps' ; fix(i)=F; ; i=i+1 lab(i)='n_skip_h' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_p' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_v' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_b' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_stat' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_cfl' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='nrk_max' ; fix(i)=F; req(i)=F; val(i)=3 ; i=i+1 lab(i)='i_interp' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='i_hydrostat' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_backg_flux' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_read_b' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_read_w' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_background' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_wind' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_damp_wind' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_smooth_wind' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_smooth_pass' ; req(i)=F; val(i)=4 ; i=i+1 lab(i)='i_atmos_var' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_mass_consv' ; req(i)=F; val(i)=1 ; i=i+1 lab(i)='i_initial' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_prob' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_noise' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_diff_meanT' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_diff_meanU' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_filt' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_filt_off' ; fix(i)=F; req(i)=F; val(i)=5 ; i=i+1 lab(i)='i_diss' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_sponge' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_force' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_mol' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_skip_model' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='i_dm_avg' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='idm_beg' ; fix(i)=F; req(i)=F; val(i)=2 ; i=i+1 lab(i)='idm_end' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='jdm_beg' ; fix(i)=F; req(i)=F; val(i)=2 ; i=i+1 lab(i)='jdm_end' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='kdm_beg' ; fix(i)=F; req(i)=F; val(i)=2 ; i=i+1 lab(i)='kdm_end' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_dm_avg_x2' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_dm_avg_y2' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='n_dm_avg_z2' ; fix(i)=F; req(i)=F; val(i)=1 ; i=i+1 lab(i)='iubc_x1' ; i=i+1 lab(i)='itbc_x1' ; i=i+1 lab(i)='iubc_x2' ; i=i+1 lab(i)='itbc_x2' ; i=i+1 lab(i)='iubc_y1' ; i=i+1 lab(i)='itbc_y1' ; i=i+1 lab(i)='iubc_y2' ; i=i+1 lab(i)='itbc_y2' ; i=i+1 lab(i)='iubc_z1' ; i=i+1 lab(i)='itbc_z1' ; i=i+1 lab(i)='iubc_z2' ; fix(i)=F; ; i=i+1 lab(i)='itbc_z2' ; fix(i)=F; ; i=i+1 lab(i)='n_skip_s' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_beg_sub' ; req(i)=F; val(i)=2 ; i=i+1 lab(i)='i_end_sub' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_skip_sub' ; req(i)=F; val(i)=1 ; i=i+1 lab(i)='j_beg_sub' ; req(i)=F; val(i)=2 ; i=i+1 lab(i)='j_end_sub' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='j_skip_sub' ; req(i)=F; val(i)=1 ; i=i+1 lab(i)='k_beg_sub' ; req(i)=F; val(i)=2 ; i=i+1 lab(i)='k_end_sub' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='k_skip_sub' ; req(i)=F; val(i)=1 ; i=i+1 lab(i)='i_yz_plane(10)' ; req(i)=F; val(i)=0 ;i=i+10 lab(i)='j_xz_plane(10)' ; req(i)=F; val(i)=0 ;i=i+10 lab(i)='k_xy_plane(16)' ; req(i)=F; val(i)=0 ;i=i+16 lab(i)='i_plane_dev' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_skip_t' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_beg_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='i_end_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='j_beg_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='j_end_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='k_beg_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='k_end_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='i_dir_ts(9)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+9 lab(i)='i_shear_dir(3)' ; fix(i)=F; req(i)=F; val(i)=0 ; i=i+3 lab(i)='i_beam' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_cycles_beam' ; req(i)=F; val(i)=2 ; i=i+1 lab(i)='i_clean_restart'; fix(i)=F; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_coriolis' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='i_stretch_z_bak' ; req(i)=F; val(i)=0 ; i=i+1 lab(i)='n_filt_pass_xy' ; fix(i)=F; req(i)=F; val(i)=6 ; i=i+1 lab(i)='n_filt_pass_z' ; fix(i)=F; req(i)=F; val(i)=6 ; i=i+1 lab(i)='i_pert_wind' ; fix(i)=F; req(i)=F; val(i)=6 ; i=i+1 lab(i)='i_pert_temp' ; fix(i)=F; req(i)=F; val(i)=6 ; i=i+1 lab(i)='i_ground_ref' ; req(i)=F; val(i)=0 ; i=i+1 n_inputs = i-1 c *** Dynamic parameters for the header file c-------------------------------List reals first------------------------ lab(i) = 'time' ; i=i+1 lab(i) = 'dt' ; i=i+1 lab(i) = 't_stat' ; i=i+1 lab(i) = 'filt_wt_x(8)' ; i=i+8 lab(i) = 'filt_wt_y(8)' ; i=i+8 lab(i) = 'filt_wt_z(8)' ; i=i+8 n_dynpar_r = i-n_inputs-1 c-------------------------------Integers below this line---------------- lab(i) = 'nt' ; i=i+1 lab(i) = 'n_frame_p' ; i=i+1 lab(i) = 'n_frame_s' ; i=i+1 lab(i) = 'n_frame_t' ; i=i+1 lab(i) = 'n_hist' ; i=i+1 n_dynpar = i-n_inputs-1 n_params = n_inputs + n_dynpar n_inputs_i = n_inputs - n_inputs_r n_dynpar_i = n_dynpar - n_dynpar_r if( n_params .gt. Np ) then print *, 'ERROR: The number of input parameters exceeds ', & 'the value of L_params' ierr = 1 end if c *** Rename array style labels to a list of scalar labels. i = 1 done = .false. do while( .not. done ) ind_b = index(lab(i), '(') if( ind_b .ne. 0 ) then label = lab(i)(1:ind_b-1) ind_e = index(lab(i), ')') read(lab(i)(ind_b+1:ind_e-1),*) n_remap do j=1, n_remap write(ext2,'(i2.2)') j i_ext_start = 2 if( j .gt. 9 ) i_ext_start=1 lab(i+j-1) = trim(label) // ext2(i_ext_start:2) val(i+j-1) = val(i) req(i+j-1) = req(i) fix(i+j-1) = fix(i) end do i = i + n_remap else if( i .le. n_params ) then i = i + 1 else done = .true. end if end do return end function index_param( param, labels, n_params ) character( *) param character(18) labels(n_params), label_lower character(len(param)) param_lower call to_lower( param, param_lower, len(param) ) do i=1, n_params call to_lower( labels(i), label_lower, 18 ) if( param_lower .eq. label_lower ) then index_param = i return end if end do index_param = 0 return end subroutine read_line( i_unit, label, value, mult, done, ierr ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com c *** Error codes: c 0 - normal exit c 2 - non-numeric character in numeric field c 4 - invalic character in label c 5 - missing numeric entry character(18) label character(80) line logical done ierr = 0 done = .false. 10 continue c *** Read a line and extract the label. Comments are flagged with # c *** Either a blank line or the string 'eof' signifies the end of file read(i_unit,'(a80)',end=50) line if( line(1:1) .eq. '#' ) goto 10 do i=1, 80 if( line(i:i) .ne. ' ' ) goto 20 end do c *** If flow gets here we found a blank line, which implies end of file goto 50 20 i_lab_beg = i i_lab_end = index(line(i_lab_beg+1:80),' ')-1 + i_lab_beg len_label = i_lab_end-i_lab_beg+1 label = line(i_lab_beg:i_lab_end) c *** Check for non-alpha characters and convert to lower case do i=1, len_label ic = ichar(label(i:i)) if( ic .lt. 48 .or. & (ic .gt. 57 .and. ic .lt. 65) .or. & (ic .gt. 90 .and. ic .lt. 95) .or. & ic .eq. 96 .or. & ic .gt. 122 ) then print *, 'ERROR: Invalid character in label = ', label ierr = 4 return end if if( ic .ge. 65 .and. ic .le. 90 ) then label(i:i) = achar(ic+32) end if end do if( label(1:len_label) .eq. 'eof' ) goto 50 c *** Isolate the numeric input and convert it from character to c *** floating point value. i_num_beg = scan( line(i_lab_end+1:80), '+-.0123456789' ) + & i_lab_end if( i_num_beg .eq. i_lab_end ) then print *, 'ERROR: No numeric entry found' ierr = 5 return end if i_num_end = index(line(i_num_beg+1:80),' ')-1 + i_num_beg c *** Look for 10^2, 10^3, and 10^6 multipliers appended to c *** numerical input mult = 1 select case(line(i_num_end:i_num_end)) case('c') mult = 100 case('k') mult = 1000 case('m') mult = 1000000 end select if( mult .ne. 1 ) i_num_end = i_num_end - 1 read( line(i_num_beg:i_num_end), *, err=40 ) value if( mult .ne. 1 ) value = value*float(mult) return 40 print *, 'ERROR: Malformed numeric entry' ierr = 2 return 50 done = .true. return end subroutine to_lower( str_in, str_out, len_str ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com character(*) str_in character(len_str) str_out c str_out = lower(str_in) ! this works instead of the lines blow str_out = str_in do i=1, len(str_in) ic = ichar(str_in(i:i)) if( ic .ge. 65 .and. ic .le. 90 ) then str_out(i:i) = achar(ic+32) end if end do return end character(*) function lower( string ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com character(*) string len_string = len(string) lower(1:len_string) = string(1:len_string) do i=1, len_string ic = ichar(string(i:i)) if( ic .ge. 65 .and. ic .le. 90 ) then lower(i:i) = achar(ic+32) end if end do return end subroutine get_fsize( fname, word_size, fsize, n_rec ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com c *** Determine the size of a file in bytes (up to 1.0 Tb max size). c *** One should be able to use the inquire intrinsic to determine the c *** file size. I tried this, but could not get it to work properly c *** on all platforms where I needed it. This routine is somewhat c *** complicated, but at least it is portable. c *** Earlier versions of this routine made use of the error code for a c *** read past the end of file. Unfortunately this code is variable c *** across different compilers and/or across different platforms and c *** thus an attempt was made to discover the error code by reading c *** something like 1 TB into the file (assuming that the file was not c *** this large and thus would give an end of file error). c *** This approach does not work (on the intel compiler at least) since c *** different error codes are returned for end of file less that 2 GB c *** in size (code 36) and greater than 2 GB in size (error code 25). c *** The work around at the moment is to assume that any error code other c *** than 0 signals a read past the end of file. c *** Taking the word_size to be > 8 can work, but may lead to c *** erroneous results. Not sure why, but iostat errors do not arise c *** consistently when a read is attempted past the end of file when c *** the record size is large. It is best to have the word_size match c *** the type of data being read. This appears to work reliably in all c *** cases. implicit none character(*) fname integer iunit, ios, word_size, eof_code, n byte, allocatable, dimension(:) :: t integer(8) max_fsize, max_rec, start, end, i, fsize, n_rec c if( word_size .gt. 8 ) then c print *, 'WARNING: In get_fize word_size = ', word_size, c & ' which is greater than 8' c print *, 'This may be ok, but check to see if this is what ', c & 'was intended' c end if allocate( t(word_size) ) max_fsize = 1.0e+12 ! 1.0 Tb max_rec = max_fsize/word_size open( newunit=iunit, file=fname, status='old', access='direct', & form='unformatted', recl=word_size ) c *** Make sure that we can read the file. Exit with an error if not. read( iunit, rec=1, iostat=ios ) t if( ios .ne. 0 ) then fsize = -1 n_rec = -1 goto 99 end if c *** Attempt to discover the end of file exit code. c *** We do not use this any longer since the error code returned can depend c *** on the file size. c read( iunit, rec=max_rec, iostat=eof_code ) t c cc print *, 'eof_code = ', eof_code c c if( eof_code .eq. 0 ) then c fsize = -99 c n_rec = -99 c goto 99 c end if c *** Check to see if the file size is equal to the input value of fsize. c *** If so, we jump to the end thereby bypassing the binary search. if( fsize .gt. 0 ) then n_rec = fsize/word_size read( iunit, rec=n_rec, iostat=ios ) t if( ios .ne. 0 ) goto 10 read( iunit, rec=n_rec+1, iostat=ios ) t if( ios .ne. 0 ) goto 99 end if 10 continue c *** Binary search for file end. if we get iostat=0, we're before the c *** end of file, if we get an iostat!=0, we're past the end of file. start = 0 end = max_rec do n=1, 100 i = ( start + end )/2 read(iunit, rec=i, iostat=ios ) t c print *, n, i, ios, i*word_size if( ios .eq. 0 ) then start = i c else if( ios .eq. eof_code ) then else end = i c else c fsize = -1 c n_rec = -1 c goto 99 end if c *** The search is done if the interval is reduced to 1 if( end-start .le. 1 ) then fsize = start*word_size n_rec = start goto 99 end if end do c *** If control gets here the binary search failed fsize = -2 n_rec = -2 99 close( iunit ) deallocate( t ) return end subroutine delete_file( fname ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com implicit none integer iunit character(*) fname logical file_exists inquire(file=fname,exist=file_exists) if( file_exists ) then open(newunit=iunit,file=fname) close(iunit,status='delete') print *, 'deleted pre-existing file ', trim(fname) end if return end subroutine move_file( fname1, fname2 ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com implicit none character(*) fname1, fname2 call copy_file( fname1, fname2, 1 ) return end subroutine copy_file( fname1, fname2, i_delete ) c © 2018 NorthWest Research Associates, Inc. All Rights Reserved. c See file license.txt in the repository root for detailed terms. c Author: Thomas S. Lund, lund@cora.nwra.com implicit none character(*) :: fname1, fname2 integer :: i_delete integer, parameter :: ibs=131072 ! blocksize = 256*512 bytes integer :: i, iunit1, iunit2, word_size, nbks, & i_left_over integer(8) :: fsize, n_rec, n_bks character(1) :: data(ibs) logical :: file_exists inquire(file=fname1,exist=file_exists) if( .not. file_exists ) then print *, 'ERROR: can not copy non-existant file ', trim(fname1) stop end if word_size = 1 call get_fsize( fname1, word_size, fsize, n_rec ) n_bks = fsize/int(ibs,kind=8) i_left_over = fsize - int(n_bks*ibs,kind=8) call delete_file( fname2 ) open(newunit=iunit1,file=fname1,form='unformatted', & access='stream',action= 'read') open(newunit=iunit2,file=fname2,form='unformatted', & access='stream',action='write') do i=1, n_bks read( iunit1) data write(iunit2) data end do read( iunit1) data(1:i_left_over) write(iunit2) data(1:i_left_over) if( i_delete .eq. 0 ) then print *, 'copied file ', trim(fname1), ' to ', trim(fname2) else print *, 'moved file ', trim(fname1), ' to ', trim(fname2) end if if( i_delete .eq. 1 ) then close(iunit1,status='delete') else close(iunit1) end if close(iunit2) return end subroutine set_dims( Nx1, Ny1, Nz1, Nx, Ny, Nz, nzp, Lx, Ly, Lz, & Lzp, Nxp1, Nxp2, Nyp1, Nyp2, Nzp1, Nzp2 ) Nx = Nx1 Ny = Ny1 Nz = Nz1 nzp = Nz Lx = Nx+2 Ly = Ny+2 Lz = Nz+2 Lzp = nzp+2 Nxp1 = Nx+1 Nyp1 = Ny+1 Nzp1 = Nz+1 Nxp2 = Nx+2 Nyp2 = Ny+2 Nzp2 = Nz+2 return end