subroutine write_inputs( base_name, nzmax, mx, & tempz, temp, cwz, cw, hwz, hw, qz, q, & rho, press, theta, bvfsq, & time_meas_p, y_meas_p, z_meas_p, & time_meas_s, y_meas_s, z_meas_s, & time_meas_cir_p, cir_meas_p, & time_meas_cir_s, cir_meas_s, & time_model, & y_model_p, z_model_p, cir_model_p, & y_model_s, z_model_s, cir_model_s, & y0, z0, v0, b0, gam0, zmf, & ac_speed, hw_avg, ntemp, ncw, nhw, nq, & nt_meas_p, nt_meas_s, & nt_meas_cir_p, nt_meas_cir_s, & nt_model, tld, tssd, & write_lin, write_star ) implicit none real y_meas_p(nzmax), z_meas_p(nzmax), cir_meas_p(nzmax), & y_meas_s(nzmax), z_meas_s(nzmax), cir_meas_s(nzmax), & time_meas_p( nzmax), time_meas_s( nzmax), & time_meas_cir_p(nzmax), time_meas_cir_s(nzmax), & tempz(nzmax), temp(nzmax), theta(nzmax), bvfsq(nzmax), & hwz(nzmax), hw(nzmax), cwz(nzmax), cw(nzmax), & qz(nzmax), q(nzmax), rho(nzmax), press(nzmax), & time_model(mx), & y_model_p( mx), z_model_p(mx), cir_model_p(mx), & y_model_s( mx), z_model_s(mx), cir_model_s(mx), & y0, z0, v0, b0, gam0, zmf, ac_speed, tld, tssd, & scale_time, scale_circ, scale_z, hw_avg, t_fac, t_lag, & tan3, pi integer nzmax, mx, ntemp, ncw, nhw, nq, & nt_meas_p, nt_meas_s, nt_meas_cir_p, nt_meas_cir_s, & nt_model logical write_lin, write_star, found_data, file_exists character(80) base_name, fname real z_min, z_max, z_tick, bvf, & eps, t_star, t_die, t_max, g_star, A_0, A_1, edr_star, & N2_star, N_star, edr_star_avg, N_star_avg, g_star_0, c13, & z_bot, t_scale, sum3, wt0, wt1, & cirp_tld, cirs_tld, cir_tld, yp_tld, ys_tld, & zp_tld, zs_tld, cirp_tssd, cirs_tssd, cir_tssd, & yp_tssd, ys_tssd, zp_tssd, zs_tssd, & meas_y_min, meas_y_max, meas_z_min, meas_z_max, & modl_y_min, modl_y_max, modl_z_min, modl_z_max, & meas_cir_min, meas_cir_max, meas_time_min, meas_time_max, & modl_cir_min, modl_cir_max, modl_time_min, modl_time_max, & meas_N_min, meas_N_max, meas_temp_min, meas_temp_max, & modl_N_min, modl_N_max, modl_temp_min, modl_temp_max, & meas_edr_min, meas_edr_max, meas_wind_min, meas_wind_max, & modl_edr_min, modl_edr_max, modl_wind_min, modl_wind_max real meas_time_maxo, modl_time_maxo, & meas_y_mino, modl_y_mino, meas_y_maxo, modl_y_maxo, & meas_z_mino, modl_z_mino, meas_z_maxo, modl_z_maxo, & meas_cir_mino, modl_cir_mino, meas_cir_maxo, modl_cir_maxo, & meas_N_mino, modl_N_mino, meas_N_maxo, modl_N_maxo, & meas_edr_mino, modl_edr_mino, meas_edr_maxo, modl_edr_maxo, & meas_temp_mino, modl_temp_mino, & meas_temp_maxo, modl_temp_maxo, & meas_wind_mino, modl_wind_mino, & meas_wind_maxo, modl_wind_maxo, & tldo, cir_tldo, yp_tldo, ys_tldo, zp_tldo, zs_tldo, & tssdo, cir_tssdo, yp_tssdo, ys_tssdo, zp_tssdo, zs_tssdo integer i, isamp, i_start, i_stop character(8) eps_ticks, eps_tickso pi = acos(-1.0) tan3 = tan(3.0/180.0*pi) c13 = 1.0/3.0 eps_ticks = 'autofreq' c *** Compute min and max values from the measurements meas_y_max = 0.0 meas_y_min = 1.0e+6*(nt_meas_p+nt_meas_s) meas_z_max = 0.0 meas_z_min = 1.0e+6*(nt_meas_p+nt_meas_s) meas_time_max = 0.0 meas_time_min = 1.0e+6*(nt_meas_p+nt_meas_s+ & nt_meas_cir_p+nt_meas_cir_s) meas_cir_max = 0.0 meas_cir_min = 1.0e+6*(nt_meas_cir_p+nt_meas_cir_s) do i=1, nt_meas_p meas_y_min = min( meas_y_min, y_meas_p(i) ) meas_y_max = max( meas_y_max, y_meas_p(i) ) meas_z_min = min( meas_z_min, z_meas_p(i) ) meas_z_max = max( meas_z_max, z_meas_p(i) ) meas_time_min = min( meas_time_min, time_meas_p(i) ) meas_time_max = max( meas_time_max, time_meas_p(i) ) end do do i=1, nt_meas_s meas_y_min = min( meas_y_min, y_meas_s(i) ) meas_y_max = max( meas_y_max, y_meas_s(i) ) meas_z_min = min( meas_z_min, z_meas_s(i) ) meas_z_max = max( meas_z_max, z_meas_s(i) ) meas_time_min = min( meas_time_min, time_meas_s(i) ) meas_time_max = max( meas_time_max, time_meas_s(i) ) end do do i=1, nt_meas_cir_p meas_cir_min = min( meas_cir_min, cir_meas_p(i) ) meas_cir_max = max( meas_cir_max, cir_meas_p(i) ) meas_time_min = min( meas_time_min, time_meas_cir_p(i) ) meas_time_max = max( meas_time_max, time_meas_cir_p(i) ) end do do i=1, nt_meas_cir_s meas_cir_min = min( meas_cir_min, cir_meas_s(i) ) meas_cir_max = max( meas_cir_max, cir_meas_s(i) ) meas_time_min = min( meas_time_min, time_meas_cir_s(i) ) meas_time_max = max( meas_time_max, time_meas_cir_s(i) ) end do c *** Compute min and max values from the model run modl_y_max = 0.0 modl_y_min = 1.0e+6 modl_z_max = 0.0 modl_z_min = 1.0e+6 modl_cir_max = 0.0 modl_cir_min = 1.0e+6 modl_time_min = time_model(1) modl_time_max = time_model(nt_model) do i=1, nt_model modl_y_min = min( modl_y_min, y_model_p(i), y_model_s(i) ) modl_y_max = max( modl_y_max, y_model_p(i), y_model_s(i) ) modl_z_min = min( modl_z_min, z_model_p(i), z_model_s(i) ) modl_z_max = max( modl_z_max, z_model_p(i), z_model_s(i) ) modl_cir_min = min( modl_cir_min, & abs(cir_model_p(i)) ,cir_model_s(i) ) modl_cir_max = max( modl_cir_max, & abs(cir_model_p(i)), cir_model_s(i) ) end do c *** Compute the min and max values of the environmental parameters c *** over the altitude ranges covered by the vortex measurements meas_N_max = -1.0e+6*meas_z_max meas_N_min = 1.0e+6*meas_z_max meas_edr_max = -1.0e+6*meas_z_max meas_edr_min = 1.0e+6*meas_z_max meas_temp_max = -1.0e+6*meas_z_max meas_temp_min = 1.0e+6*meas_z_max meas_wind_max = -1.0e+6*meas_z_max meas_wind_min = 1.0e+6*meas_z_max do i=1, ntemp if( tempz(i) .ge. meas_z_min .and. & tempz(i) .le. meas_z_max ) then bvf = bvfsq(i)/abs(bvfsq(i)+1.0e-12)*sqrt(abs(bvfsq(i))) meas_N_min = min( meas_N_min, bvf ) meas_N_max = max( meas_N_max, bvf ) meas_temp_min = min( meas_temp_min, temp(i) ) meas_temp_max = max( meas_temp_max, temp(i) ) end if end do found_data = .false. do i=1, nq if( qz(i) .ge. meas_z_min .and. & qz(i) .le. meas_z_max ) then meas_edr_min = min( meas_edr_min, q(i) ) meas_edr_max = max( meas_edr_max, q(i) ) found_data = .true. end if end do if( .not. found_data ) then meas_edr_min = minval( q, nq ) meas_edr_max = maxval( q, nq ) end if found_data = .false. do i=1, nhw if( hwz(i) .ge. meas_z_min .and. & hwz(i) .le. meas_z_max ) then meas_wind_min = min( meas_wind_min, hw(i) ) meas_wind_max = max( meas_wind_max, hw(i) ) found_data = .true. end if end do if( .not. found_data ) then meas_wind_min = minval( hw, nhw ) meas_wind_max = maxval( hw, nhw ) end if found_data = .false. do i=1, ncw if( cwz(i) .ge. meas_z_min .and. & cwz(i) .le. meas_z_max ) then meas_wind_min = min( meas_wind_min, cw(i) ) meas_wind_max = max( meas_wind_max, cw(i) ) found_data = .true. end if end do if( .not. found_data ) then meas_wind_min = min( meas_wind_min, minval( cw, ncw ) ) meas_wind_max = max( meas_wind_max, maxval( cw, ncw ) ) end if c *** Compute the min and max values of the environmental parameters c *** over the altitude ranges covered by the model run. modl_N_max = -10.0 modl_N_min = 10.0 modl_edr_max = -1.0 modl_edr_min = 1.0 modl_temp_max = -100.0 modl_temp_min = 100.0 modl_wind_max = -100.0 modl_wind_min = 100.0 do i=1, ntemp if( tempz(i) .ge. modl_z_min .and. & tempz(i) .le. modl_z_max ) then bvf = bvfsq(i)/abs(bvfsq(i)+1.0e-12)*sqrt(abs(bvfsq(i))) modl_N_min = min( modl_N_min, bvf ) modl_N_max = max( modl_N_max, bvf ) modl_temp_min = min( modl_temp_min, temp(i) ) modl_temp_max = max( modl_temp_max, temp(i) ) end if end do found_data = .false. do i=1, nq if( qz(i) .ge. modl_z_min .and. & qz(i) .le. modl_z_max ) then modl_edr_min = min( modl_edr_min, q(i) ) modl_edr_max = max( modl_edr_max, q(i) ) found_data = .true. end if end do if( .not. found_data ) then modl_edr_min = minval( q, nq ) modl_edr_max = maxval( q, nq ) end if found_data = .false. do i=1, nhw if( hwz(i) .ge. modl_z_min .and. & hwz(i) .le. modl_z_max ) then modl_wind_min = min( modl_wind_min, hw(i) ) modl_wind_max = max( modl_wind_max, hw(i) ) found_data = .true. end if end do if( .not. found_data ) then modl_wind_min = minval( hw, nhw ) modl_wind_max = maxval( hw, nhw ) end if found_data = .false. do i=1, ncw if( cwz(i) .ge. modl_z_min .and. & cwz(i) .le. modl_z_max ) then modl_wind_min = min( modl_wind_min, cw(i) ) modl_wind_max = max( modl_wind_max, cw(i) ) found_data = .true. end if end do if( .not. found_data ) then modl_wind_min = min( modl_wind_min, minval( cw, ncw ) ) modl_wind_max = max( modl_wind_max, maxval( cw, ncw ) ) end if c *** Set the edr* tics level for cases where the value is constant if( nq .eq. 2 .and. q(1) .ne. 0.0 ) eps_ticks = '0.001' c *** Determine the circulation, y, and z values corresponding to c *** the instability times c print *, 'tld, tssd = ', tld, tssd cirp_tld =0.0; cirs_tld =0.0; cir_tld=0.0 yp_tld =0.0; ys_tld =0.0; zp_tld =0.0; zs_tld =0.0; cirp_tssd=0.0; cirs_tssd=0.0; cir_tssd=0.0 yp_tssd=0.0; ys_tssd=0.0; zp_tssd=0.0; zs_tssd=0.0; do i=1, nt_model-1 if( tld .ge. time_model(i) .and. & tld .lt. time_model(i+1) ) then wt0 = ( time_model(i+1) - tld )/ & ( time_model(i+1) - time_model(i) ) wt1 = 1.0 - wt0 cirp_tld = abs(wt0*cir_model_p(i)+wt1*cir_model_p(i+1)) cirs_tld = abs(wt0*cir_model_s(i)+wt1*cir_model_s(i+1)) yp_tld = wt0* y_model_p(i)+wt1* y_model_p(i+1) ys_tld = wt0* y_model_s(i)+wt1* y_model_s(i+1) zp_tld = wt0* z_model_p(i)+wt1* z_model_p(i+1) zs_tld = wt0* z_model_s(i)+wt1* z_model_s(i+1) cir_tld = min( cirp_tld, cirs_tld ) end if if( tssd .ge. time_model(i) .and. & tssd .lt. time_model(i+1) ) then wt0 = ( time_model(i+1) - tssd )/ & ( time_model(i+1) - time_model(i) ) wt1 = 1.0 - wt0 cirp_tssd = abs(wt0*cir_model_p(i)+wt1*cir_model_p(i+1)) cirs_tssd = abs(wt0*cir_model_s(i)+wt1*cir_model_s(i+1)) yp_tssd = wt0* y_model_p(i)+wt1* y_model_p(i+1) ys_tssd = wt0* y_model_s(i)+wt1* y_model_s(i+1) zp_tssd = wt0* z_model_p(i)+wt1* z_model_p(i+1) zs_tssd = wt0* z_model_s(i)+wt1* z_model_s(i+1) cir_tssd = min( cirp_tssd, cirs_tssd ) end if end do c *** Check to see if a limits file exists. If so, read its contents and c *** take its values into account. fname = 'Output/'//trim(base_name)//'.limits' inquire(file=trim(fname),exist=file_exists) open(unit=9,file=trim(fname)) if( file_exists ) then do i=1, 4 read(9,*) end do read(9,8) meas_time_maxo, modl_time_maxo, & meas_y_mino, modl_y_mino, meas_y_maxo, modl_y_maxo, & meas_z_mino, modl_z_mino, meas_z_maxo, modl_z_maxo, & meas_cir_mino, modl_cir_mino, meas_cir_maxo, modl_cir_maxo, & meas_N_mino, modl_N_mino, meas_N_maxo, modl_N_maxo, & meas_edr_mino, modl_edr_mino, meas_edr_maxo, modl_edr_maxo, & meas_temp_mino, modl_temp_mino, & meas_temp_maxo, modl_temp_maxo, & meas_wind_mino, modl_wind_mino, & meas_wind_maxo, modl_wind_maxo, & eps_tickso, & tldo, cir_tldo, yp_tldo, ys_tldo, zp_tldo, zs_tldo, & tssdo, cir_tssdo, yp_tssdo, ys_tssdo, zp_tssdo, zs_tssdo 8 format(18(14x,f8.1),4(14x,e10.3),8(14x,f8.1),(14x,a8), & 18(12x,f8.1)) rewind(9) meas_time_max = max( meas_time_max, meas_time_maxo ) modl_time_max = max( modl_time_max, modl_time_maxo ) meas_y_min = min( meas_y_min, meas_y_mino ) modl_y_min = min( modl_y_min, modl_y_mino ) meas_y_max = max( meas_y_max, meas_y_maxo ) modl_y_max = max( modl_y_max, modl_y_maxo ) meas_z_min = min( meas_z_min, meas_z_mino ) modl_z_min = min( modl_z_min, modl_z_mino ) meas_z_max = max( meas_z_max, meas_z_maxo ) modl_z_max = max( modl_z_max, modl_z_maxo ) meas_cir_min = min( meas_cir_min, meas_cir_mino ) modl_cir_min = min( modl_cir_min, modl_cir_mino ) meas_cir_max = max( meas_cir_max, meas_cir_maxo ) modl_cir_max = max( modl_cir_max, modl_cir_maxo ) meas_N_min = min( meas_N_min, meas_N_mino ) modl_N_min = min( modl_N_min, modl_N_mino ) meas_N_max = max( meas_N_max, meas_N_maxo ) modl_N_max = max( modl_N_max, modl_N_maxo ) meas_edr_min = min( meas_edr_min, meas_edr_mino ) modl_edr_min = min( modl_edr_min, modl_edr_mino ) meas_edr_max = max( meas_edr_max, meas_edr_maxo ) modl_edr_max = max( modl_edr_max, modl_edr_maxo ) meas_temp_min = min( meas_temp_min, meas_temp_mino ) modl_temp_min = min( modl_temp_min, modl_temp_mino ) meas_temp_max = max( meas_temp_max, meas_temp_maxo ) modl_temp_max = max( modl_temp_max, modl_temp_maxo ) meas_wind_min = min( meas_wind_min, meas_wind_mino ) modl_wind_min = min( modl_wind_min, modl_wind_mino ) meas_wind_max = max( meas_wind_max, meas_wind_maxo ) modl_wind_max = max( modl_wind_max, modl_wind_maxo ) tld = max( tld , tldo ) cir_tld = max( cir_tld , cir_tldo ) zp_tld = max( zp_tld , zp_tldo ) zs_tld = max( zs_tld , zs_tldo ) cir_tssd = max( cir_tssd, cir_tssdo ) zp_tssd = max( zp_tssd, zp_tssdo ) zs_tssd = max( zs_tssd, zs_tssdo ) if( abs(yp_tldo ) .gt. abs(yp_tld ) ) yp_tld = yp_tldo if( abs(ys_tldo ) .gt. abs(ys_tld ) ) ys_tld = ys_tldo if( abs(yp_tssdo) .gt. abs(yp_tssd) ) yp_tssd = yp_tssdo if( abs(ys_tssdo) .gt. abs(ys_tssd) ) ys_tssd = ys_tssdo if( trim(eps_ticks ) .eq. '0.001' .or. & trim(eps_tickso) .eq. '0.001' ) eps_ticks = '0.001' end if c *** Write the min/max data. write(9,9) trim(base_name), b0, v0, z0, & meas_time_max, modl_time_max, & meas_y_min, modl_y_min, meas_y_max, modl_y_max, & meas_z_min, modl_z_min, meas_z_max, modl_z_max, & meas_cir_min, modl_cir_min, meas_cir_max, modl_cir_max, & meas_N_min, modl_N_min, meas_N_max, modl_N_max, & meas_edr_min, modl_edr_min, meas_edr_max, modl_edr_max, & meas_temp_min, modl_temp_min, & meas_temp_max, modl_temp_max, & meas_wind_min, modl_wind_min, & meas_wind_max, modl_wind_max, & eps_ticks, & tld, cir_tld, yp_tld, ys_tld, zp_tld, zs_tld, & tssd, cir_tssd, yp_tssd, ys_tssd, zp_tssd, zs_tssd 9 format('case',a30,/, & 'b0 ',f8.3,/, & 'v0 ',f8.3,/, & 'z0 ',f8.3,/, & 'meas_time_max ',f8.1,/, & 'modl_time_max ',f8.1,/, & 'meas_y_min ',f8.1,/, & 'modl_y_min ',f8.1,/, & 'meas_y_max ',f8.1,/, & 'modl_y_max ',f8.1,/, & 'meas_z_min ',f8.1,/, & 'modl_z_min ',f8.1,/, & 'meas_z_max ',f8.1,/, & 'modl_z_max ',f8.1,/, & 'meas_cir_min ',f8.1,/, & 'modl_cir_min ',f8.1,/, & 'meas_cir_max ',f8.1,/, & 'modl_cir_max ',f8.1,/, & 'meas_N_min ',f8.1,/, & 'modl_N_min ',f8.1,/, & 'meas_N_max ',f8.1,/, & 'modl_N_max ',f8.1,/, & 'meas_edr_min ',e10.3,/, & 'modl_edr_min ',e10.3,/, & 'meas_edr_max ',e10.3,/, & 'modl_edr_max ',e10.3,/, & 'meas_temp_min ',f8.1,/, & 'modl_temp_min ',f8.1,/, & 'meas_temp_max ',f8.1,/, & 'modl_temp_max ',f8.1,/, & 'meas_wind_min ',f8.1,/, & 'modl_wind_min ',f8.1,/, & 'meas_wind_max ',f8.1,/, & 'modl_wind_max ',f8.1,/, & 'epsticks ',a8,/, & 'tld ',f8.1,/, & 'cirtld ',f8.1,/, & 'yptld ',f8.1,/, & 'ystld ',f8.1,/, & 'zptld ',f8.1,/, & 'zstld ',f8.1,/, & 'tssd ',f8.1,/, & 'cirtssd ',f8.1,/, & 'yptssd ',f8.1,/, & 'ystssd ',f8.1,/, & 'zptssd ',f8.1,/, & 'zstssd ',f8.1 ) close(9) c *** Write various files for the input data scale_time = (v0/b0) scale_circ = 1.0/(2.0*pi*b0*V0) scale_z = 1.0/b0 if( write_star .and. nt_meas_cir_p+nt_meas_cir_s .gt. 0 ) then close(19) open(unit=19,file='Output/circ_star.dat',access='append') end if if( write_star .and. nt_meas_p+nt_meas_s .gt. 0 ) then close(20) open(unit=20,file='Output/z_star.dat',access='append') end if if( nt_meas_p .gt. 0 ) then open(unit=10,file='Output/'//trim(base_name)//'.meas_yz_p') write(10,110) 110 format('# time y_p z_p') t_fac = 1.0 - hw_avg/ac_speed do i=1, nt_meas_p write(10,40) time_meas_p(i), y_meas_p(i), z_meas_p(i) 40 format(1p,6e12.4) t_lag = time_meas_p(i)*t_fac if( write_star ) then write(20,40) t_lag*scale_time, & (z_meas_p(i)+hw_avg*tan3*t_lag-z0)*scale_z end if end do close(10) end if if( nt_meas_s .gt. 0 ) then open(unit=11,file='Output/'//trim(base_name)//'.meas_yz_s') write(11,111) 111 format('# time y_s z_s') do i=1, nt_meas_s write(11,40) time_meas_s(i), y_meas_s(i), z_meas_s(i) t_lag = time_meas_s(i)*t_fac if( write_star ) then write(20,40) t_lag*scale_time, & (z_meas_s(i)+hw_avg*tan3*t_lag-z0)*scale_z end if end do close(11) end if if( nt_meas_cir_p .gt. 0 ) then open(unit=12,file='Output/'//trim(base_name)//'.meas_cir_p') write(12,112) 112 format('# time cir_p') do i=1, nt_meas_cir_p write(12,40) time_meas_cir_p(i), cir_meas_p(i) if( write_star ) then write(19,40) time_meas_cir_p(i)*t_fac*scale_time, & cir_meas_p(i)*scale_circ end if end do close(12) end if if( nt_meas_cir_s .gt. 0 ) then open(unit=13,file='Output/'//trim(base_name)//'.meas_cir_s') write(13,113) 113 format('# time cir_s') do i=1, nt_meas_cir_s write(13,40) time_meas_cir_s(i), cir_meas_s(i) if( write_star ) then write(19,40) time_meas_cir_s(i)*t_fac*scale_time, & cir_meas_s(i)*scale_circ end if end do close(13) end if if( write_star .and. nt_meas_cir_p+nt_meas_cir_s .gt. 0) close(19) if( write_star .and. nt_meas_p +nt_meas_s .gt. 0) close(20) open(unit=14,file='Output/'//trim(base_name)//'.meas_cw') if( meas_z_min .eq. 0.0 ) then z_min = modl_z_min else z_min = min( meas_z_min, modl_z_min ) end if z_max = max( meas_z_max, modl_z_max ) z_tick = 20.0 z_min = ( int(z_min/z_tick) )*z_tick z_max = ( int(z_max/z_tick) + 1 )*z_tick do i=1, ncw if( cwz(i) .gt. z_max ) goto 202 end do 202 i_stop = min( i, ncw ) do i=ncw, 1, -1 if( cwz(i) .lt. z_min ) goto 203 end do 203 i_start = max( i, 1 ) write(14,114) 114 format('# z cross wind') if( ncw .eq. 2 .and. cw(1) .eq. cw(2) ) then write(14,40) z_min, cw(1) write(14,40) z_max, cw(2) else do i=i_start, i_stop write(14,40) cwz(i), cw(i) end do end if close(14) open(unit=15,file='Output/'//trim(base_name)//'.meas_hw') do i=1, nhw if( hwz(i) .gt. z_max ) goto 204 end do 204 i_stop = min( i, nhw ) do i=nhw, 1, -1 if( hwz(i) .lt. z_min ) goto 205 end do 205 i_start = max( i, 1 ) write(15,115) 115 format('# z head wind') if( nhw .le. 3 .and. hw(i_start) .eq. hw(i_stop) ) then write(15,40) z_min, hw(i_start) write(15,40) z_max, hw(i_start) else do i=i_start, i_stop write(15,40) hwz(i), hw(i) end do end if close(15) open(unit=16,file='Output/'//trim(base_name)//'.meas_edr') do i=1, nq if( qz(i) .gt. z_max ) goto 206 end do 206 i_stop = min( i, nq ) do i=nq, 1, -1 if( qz(i) .lt. z_min ) goto 207 end do 207 i_start = max( i, 1 ) write(16,116) 116 format('# z eddy dissipation rate') if( nq .eq. 2 .and. q(1) .eq. q(2) ) then write(16,40) z_min, q(1) write(16,40) z_max, q(2) else do i=i_start, i_stop write(16,40) qz(i), q(i) end do end if close(16) open(unit=17,file='Output/'//trim(base_name)//'.meas_state') do i=1, ntemp if( tempz(i) .gt. z_max ) goto 210 end do 210 i_stop = min( i, ntemp ) do i=ntemp, 1, -1 if( tempz(i) .lt. z_min ) goto 211 end do 211 i_start = max( i, 1 ) write(17,117) 117 format('# z rho p T ', & 'theta N^2') do i=i_start, i_stop write(17,40) tempz(i), rho(i), press(i), temp(i), theta(i), & bvfsq(i) end do close(17) c *** Write the linear model for circulation for reference if( write_lin ) then open(unit=20,file='Output/'//trim(base_name)//'.lin') t_scale = v0/b0 z_bot = 1.5*b0 sum3 = 0.0 isamp = 0 do i=1, nq if( z_bot .le. qz(i) .and. qz(i) .le. z0 ) then edr_star = ((q(i)*b0)**c13)/v0 sum3 = sum3 + edr_star isamp = isamp + 1 end if end do if( isamp .eq. 0 ) then edr_star_avg = ((q(1)*b0)**c13)/v0 else edr_star_avg = sum3/float(isamp) end if sum3 = 0.0 isamp = 0 do i=1, ntemp if( z_bot .le. tempz(i) .and. tempz(i) .le. z0 ) then N2_star = bvfsq(i)*(b0/v0)**2 N_star = N2_star/abs(N2_star+1.0e-12)* & sqrt(abs(N2_star)) sum3 = sum3 + N_star isamp = isamp + 1 end if end do if( isamp .eq. 0 ) then N2_star = bvfsq(1)*(b0/v0)**2 N_star_avg = N2_star/abs(N2_star+1.0e-12)* & sqrt(abs(N2_star)) else N_star_avg = sum3/float(isamp) end if g_star_0 = 1.0 A_0 = -0.102 A_1 = A_0 - 0.687*edr_star_avg + 0.048*N_star_avg c print *, 'edr_star_avg, N_star_avg, A1 = ' c print *, edr_star_avg, N_star_avg, A_1 t_die = 1 - ( g_star_0 + A_0 )/A_1 t_max = min( t_die, time_model(nt_model)*t_scale ) write(20,40) 0.0, g_star_0 t_star = 1.0 write(20,40) t_star, g_star_0 + A_0*t_star t_star = t_max write(20,40) t_star, g_star_0 + A_0 + A_1*(t_star-1.0) close(20) end if return end