! This module has the common variables and parameters that is shared across many subroutines module common_variables character(len=256) :: namelist_file="",exp_name="" ! Determine dt from CFL condition, and other conditions ! Smallest value from python program "bin/python/cloud_sys_modeling/project/determine_time_step_CFL_cond.py" real*4 :: dt integer :: ni,nj,nk integer :: ni_u, nj_v, nk_w real*4 :: Lx,Ly,H real*4 :: Kth, Kvel, Kqw, Cs, alpha_damp, alpha_damp_factor real*4 :: sim_hrs, sim_mins, sim_secs, sim_time,output_intvl real*4 :: theta_bot, theta_top, delta_theta,theta_swall real*4 :: qw_bot,qw_top,qw_swall,press_top, press_bot integer :: free_slip ,side_wall , same_press, print_slices !flags 0 or 1 values integer :: n_theta_bins,n_qv_bins, n_super_sat_bins !Variables integer :: ITT,ITTNOW,ITTMAX, output_itt real*4 :: T_top,T_bot, TKE real*4 :: qv_sat_bot, qv_sat_top, es_sat_bot,es_sat_top real*4 :: dz, dy, dx real*4,dimension(:),allocatable :: theta_binEdges, qv_binEdges, super_sat_binEdges !Variables with one dimension either x,y,z, or t (1D) real*4,dimension (:),allocatable :: x,y,z, theta0,u0,v0,qv0,theta_v0, press0, exner_pi0, qv_sat_swall0 !Variables with (3D) real*4,dimension(:,:,:),allocatable :: theta_l,theta,pi,qv,qc,qw, div_damp, super_sat real*4,dimension(:,:,:),allocatable :: u, v, w real*4,dimension(:,:,:),allocatable :: u_at_theta,v_at_theta, w_at_theta !Forcing at each grid cell except at boundary and ghost points - (4D) real*4,dimension(:,:,:,:),allocatable :: fth_l,fpi,fqw, fu, fv, fw !Constants real*4 :: g = 9.81, Cp =1005, Lv = 2.5 *1e6, Rd = 287.0, & pi_const = 4.D0*DATAN(1.D0), f15 = 0.376*1e-4 ! at Lat 15N/S !File handles integer :: stdout = 6,& thl_xy = 101, thl_xz = 102, u_xy = 103, u_xz = 104,& v_xy = 105, v_xz = 106, w_xy = 107, w_xz = 108,& pi_xy = 109, pi_xz = 110, qv_xy = 111, qv_xz = 112,& qc_xy = 113, qc_xz = 114, super_sat_xy = 115, super_sat_xz=116; integer :: perturb_file = 301, tke_file=302, iter_stats_max=303, iter_stats_avg = 304, iter_stats_std = 305 integer :: theta_l_profile = 401, theta_profile = 402, uvar_profile = 403, vvar_profile = 404,& wvar_profile = 405, qv_profile = 406, qc_profile = 407, pi_profile = 408,& super_sat_profile= 409 integer :: theta_hist = 501, qv_hist = 502, super_sat_hist= 503 contains subroutine model_setup !use common_variables character(len=256) :: arg="",args(0:100)="" !0-> program, 1->namelist, 2-6 model params, 7->eof character(len=256) :: to_lower,nl_filename="" integer :: error_flag = 0 integer :: i !Read command line arguments print *, "Cmd line arguments below overrides the namelist values:" i = 0 do call get_command_argument(i, arg) print *, i, ":", trim(arg) if (len_trim(arg) == 0) then exit else args(i) = to_lower(trim(arg)) end if !Extract namelist filename if ( index(trim(arg),"namelist=") >0 ) then nl_filename = args(i) nl_filename = nl_filename(len("namelist=")+1:) !Global variable namelist_file = trim(nl_filename) end if i = i+1 end do args(i) = "eof" !Validate cmd line arguments and assign them as model params if (i == 1 .or. len_trim(nl_filename) == 0 ) then print *, "***** Error: Please provide a namelist file as cmd line arg: namelist=path/to/namelist.txt" print *, " " print *, "Any model parameter can be provided at command line up to 10 parameters" print *, "Cmd line values will override the values in the file" print *, "At command line specify parameter name and value without spaces: name=value" print *, "The name and value are case insensitive" error_flag = 1 elseif (i > 10 ) then print *, "***** Error: Number of cmd lines parameters should be less than 10" print *, " " print *, "Any model parameter can be provided at command line up to 10 parameters" print *, "Cmd line values will override the values in the file" print *, "At command line specify parameter name and value without spaces: name=value" print *, "The name and value are case insensitive" error_flag = 1 else !Read model params from the namelist file call load_namelist_file(nl_filename) !Overrides model params from namelist file with command line values call override_w_cmdln_args(args) !print model params call print_model_params() end if if (error_flag == 1) stop end subroutine subroutine load_namelist_file(nl_filename) !use common_variables character(len=256) :: lines(0:100)="",ln="",nl_filename integer :: ios, nl_fhdl = 1000 character(len=256) :: to_lower integer :: i print *,"" print *,"Model params from namelist file: " , trim(nl_filename) open(nl_fhdl,file = nl_filename, status ="old", form = "formatted", iostat = ios) if (ios == 6 .or. ios == 2) then print *,"File: ", trim(nl_filename) , " not found." stop elseif (ios /= 0) then print *, "Error code: " , ios stop end if !Read model parameters into array i = 0 do 10 read(nl_fhdl,'(A)', iostat=ios) ln if (ios == -1 ) exit !EoF !Other read error if ( ios /= 0 ) then print *,"Read error code: " , ios exit endif ln = trim(ln) ln = to_lower(ln) if ( index(ln,"#") == 1 )then goto 10 !Comment elseif (len_trim(ln) == 0)then goto 10 !Empty line else print *, trim(ln) lines(i) = trim(ln) i = i+1 endif end do lines(i) = "eof" close(nl_fhdl) call extract_model_params(lines) end subroutine subroutine override_w_cmdln_args(args) character(len=256) :: args(0:100) call extract_model_params(args) end subroutine subroutine extract_model_params(lines) !use common_variables character(len=256) :: lines(0:100),ln="" integer :: ios, error_flag =0, incorrect_input =0 i = 0 do ln = lines(i) if ( index(ln,"eof") == 1) then exit elseif (index(ln,"exp_name") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) exp_name if (ios /=0 .or. len_trim(exp_name) <= 0) incorrect_input = 1 elseif (index(ln,"ni") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) ni if (ios /=0 .or. ni <= 0) incorrect_input = 1 elseif (index(ln,"nj") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) nj if (ios /=0 .or. nj <= 0) incorrect_input = 1 elseif (index(ln,"nk") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) nk if (ios /=0 .or. nk <= 0) incorrect_input = 1 elseif (index(ln,"lx") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Lx if (ios /=0 .or. Lx <= 0) incorrect_input = 1 elseif (index(ln,"ly") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Ly if (ios /=0 .or. Ly <= 0) incorrect_input = 1 elseif (index(ln,"h") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) H if (ios /=0 .or. H <= 0) incorrect_input = 1 elseif (index(ln,"sim_hrs") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) sim_hrs if (ios /=0 .or. sim_hrs < 0) incorrect_input = 1 elseif (index(ln,"sim_min") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) sim_mins if (ios /=0 .or. sim_mins < 0) incorrect_input = 1 elseif (index(ln,"sim_secs") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) sim_secs if (ios /=0 .or. sim_secs < 0) incorrect_input = 1 elseif (index(ln,"dt") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) dt if (ios /=0 .or. dt <= 0) incorrect_input = 1 elseif (index(ln,"kth") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Kth if (ios /=0 .or. Kth <= 0) incorrect_input = 1 elseif (index(ln,"kvel") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Kvel if (ios /=0 .or. Kvel <= 0) incorrect_input = 1 elseif (index(ln,"kqw") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Kqw if (ios /=0 .or. Kqw <= 0) incorrect_input = 1 elseif (index(ln,"cs") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Cs if (ios /=0 .or. Cs <= 0) incorrect_input = 1 elseif (index(ln,"alpha_damp_factor") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) alpha_damp_factor if (ios /=0 .or. alpha_damp_factor <= 0) incorrect_input = 1 elseif (index(ln,"press_bot") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) press_bot if (ios /=0 .or. press_bot <= 0) incorrect_input = 1 elseif (index(ln,"same_press") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) same_press if (ios /=0 .or. same_press < 0 .or. side_wall > 1) incorrect_input = 1 elseif (index(ln,"theta_bot") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Theta_bot if (ios /=0 .or. Theta_bot <= 0) incorrect_input = 1 elseif (index(ln,"theta_top") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Theta_top if (ios /=0 .or. Theta_top <= 0) incorrect_input = 1 elseif (index(ln,"qw_bot") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) qw_bot if (ios /=0 .or. qw_bot < 0 .or. qw_bot > 1) incorrect_input = 1 elseif (index(ln,"qw_top") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) qw_top if (ios /=0 .or. qw_top < 0 .or. qw_top > 1) incorrect_input = 1 elseif (index(ln,"side_wall") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) side_wall if (ios /=0 .or. side_wall < 0 .or. side_wall > 1) incorrect_input = 1 elseif (index(ln,"theta_swall") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) Theta_swall if (ios /=0 .or. Theta_swall < 0 ) incorrect_input = 1 elseif (index(ln,"qw_swall") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) qw_swall if (ios /=0 .or. qw_swall < -1 .or. qw_swall > 1) incorrect_input = 1 elseif (index(ln,"free_slip") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) free_slip if (ios /=0 .or. free_slip < 0 .or. free_slip > 1) incorrect_input = 1 elseif (index(ln,"output_intvl") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) output_intvl if (ios /=0 .or. output_intvl <= 0 ) incorrect_input = 1 elseif (index(ln,"print_slices") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) print_slices if (ios /=0 .or. print_slices < 0 .or. print_slices > 1) incorrect_input = 1 elseif (index(ln,"n_theta_bins") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) n_theta_bins if (ios /=0 .or. n_theta_bins <= 0) incorrect_input = 1 elseif (index(ln,"n_qv_bins") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) n_qv_bins if (ios /=0 .or. n_qv_bins <= 0) incorrect_input = 1 elseif (index(ln,"n_super_sat_bins") == 1) then read(ln(index(ln,"=")+1:),*,iostat = ios ) n_super_sat_bins if (ios /=0 .or. n_super_sat_bins <= 0 ) incorrect_input = 1 endif if (incorrect_input == 1) then print *,"" print *, "Error: Incorrect input - ", trim(ln) error_flag = 1 incorrect_input = 0 end if i = i+1 end do !Stop this program, if there is incorrect model param value if (error_flag == 1) stop end subroutine subroutine print_model_params !use common_variables print *,"" print '(2A)', "Experiment Name : " , trim(exp_name) print '(A,10I10)', "Number of interior grid points: ni,nj,nk : ", ni, nj, nk print '(A,10F10.2)' , "Domain size Lx,Ly,H : ", Lx, Ly, H print '(A,10F10.2)', "sim_time : ", sim_hrs, sim_mins, sim_secs print '(A,10F10.2)', "Time step dt : ", dt print '(A,10F10.2)', "Diffusion coeff : ", Kth, Kvel, Kqw print '(A,10F10.2)', "alpha damp factor : ", alpha_damp_factor print '(A,10F10.2)', "speed of sound : ", Cs print '(A,10I10)', "sidewall,free_slip,same_press,print_slices : ", side_wall, free_slip, same_press,print_slices print '(A,10F10.2)', "Theta bot, top, sidewall : ", Theta_bot, Theta_top, Theta_swall print '(A,10F10.2)', "qw bot, top, sidewall : ", qw_bot,qw_top, qw_swall print '(A,10F10.2)', "press bot : ", press_bot print '(A,10I10)', "Histogram bins: theta, qv, super_sat : ", n_theta_bins, n_qv_bins, n_super_sat_bins print *, " " end subroutine subroutine allocate_memory !Allocate memory and initialize them to zero !use common_variables !Allocate memory for position variables allocate( x(0:ni+1) ); x = 0.0 allocate( y(0:nj+1) ); y = 0.0 allocate( z(0:nk+1) ); z = 0.0 !Allocate memory for base state profiles at each vertical level allocate( theta0(0:nk+1) ); theta0 = 0.0 allocate( u0(0:nk+1) ); u0 = 0.0 allocate( v0(0:nk+1) ); v0 = 0.0 allocate( qv0(0:nk+1) ); qv0 = 0.0 allocate(theta_v0(0:nk+1) );theta_v0 = 0.0 allocate( press0(0:nk+1) ); press0 = 0.0 allocate(exner_pi0(0:nk+1));exner_pi0= 0.0 allocate(qv_sat_swall0(0:nk+1));qv_sat_swall0= 0.0 !Allocate variables defined at theta grid points allocate( theta_l(0:ni+1, 0:nj+1, 0:nk+1) ); theta_l = 0.0 allocate( theta(0:ni+1, 0:nj+1, 0:nk+1) ); theta = 0.0 allocate( pi(0:ni+1, 0:nj+1, 0:nk+1) ); pi = 0.0 allocate( qv(0:ni+1, 0:nj+1, 0:nk+1) ); qv = 0.0 allocate( qc(0:ni+1, 0:nj+1, 0:nk+1) ); qc = 0.0 allocate( qw(0:ni+1, 0:nj+1, 0:nk+1) ); qw = 0.0 allocate( div_damp(0:ni+1, 0:nj+1, 0:nk+1) ); div_damp = 0.0 allocate( super_sat(0:ni+1, 0:nj+1, 0:nk+1) ); super_sat = 0.0 allocate( u_at_theta(0:ni+1, 0:nj+1, 0:nk+1) ); u_at_theta = 0.0 allocate( v_at_theta(0:ni+1, 0:nj+1, 0:nk+1) ); v_at_theta = 0.0 allocate( w_at_theta(0:ni+1, 0:nj+1, 0:nk+1) ); w_at_theta = 0.0 !Allocate variables at staggerred grid points allocate( u(0:ni_u+1,0:nj+1, 0:nk+1) ); u = 0.0 allocate( v(0:ni+1, 0:nj_v+1,0:nk+1) ); v = 0.0 allocate( w(0:ni+1, 0:nj+1, 0:nk_w+1)); w = 0.0 !Allocate forcing variables allocate( fth_l(1:ni, 1:nj, 1:nk, 2) ); fth_l = 0.0 allocate( fpi(1:ni, 1:nj, 1:nk, 2) ); fpi = 0.0 allocate( fqw(1:ni, 1:nj, 1:nk, 2) ); fqw = 0.0 allocate( fu(1:ni_u, 1:nj, 1:nk, 2) ); fu = 0.0 allocate( fv(1:ni, 1:nj_v, 1:nk, 2) ); fv = 0.0 allocate( fw(1:ni, 1:nj, 1:nk_w, 2) ); fw = 0.0 end subroutine subroutine open_files() !use common_variables character(len=256) :: exp_dir,cmd,slices_dir,profiles_dir integer :: cmd_sts !===Create experiment directory print *, new_line("A"), "experiment name:" , trim(exp_name) call system("ls", cmd_sts) if (cmd_sts == 1)then !For window os exp_dir = "data\"//trim(exp_name)//"\" cmd = "rmdir /s /q "//trim(exp_dir) call system(trim(cmd), status = cmd_sts) !make dir cmd = "mkdir "//trim(exp_dir)//"slices\" call system(trim(cmd), status = cmd_sts) cmd = "mkdir "//trim(exp_dir)//"profiles\" call system(trim(cmd), status = cmd_sts) slices_dir = trim(exp_dir)//"slices\" profiles_dir = trim(exp_dir)//"profiles\" else !For unix based os exp_dir = "data/"//trim(exp_name)//"/" cmd = "rm -r "//trim(exp_dir) call system(trim(cmd), status = cmd_sts) !make dir cmd = "mkdir -p "//trim(exp_dir)//"slices/" call system(trim(cmd), status = cmd_sts) cmd = "mkdir -p "//trim(exp_dir)//"profiles/" call system(trim(cmd), status = cmd_sts) slices_dir = trim(exp_dir)//"slices/" profiles_dir = trim(exp_dir)//"profiles/" end if print *,"experiment output dir: " , trim(exp_dir) print *,"experiment output dir: " , trim(slices_dir) !=== end create directories if (print_slices == 1) then open(unit = thl_xy,file = trim(slices_dir)//"theta_l_xy.txt"); open(unit = thl_xz,file = trim(slices_dir)//"theta_l_xz.txt"); open(unit = u_xy, file = trim(slices_dir)//"u_xy.txt" ); open(unit = u_xz, file = trim(slices_dir)//"u_xz.txt" ); open(unit = v_xy, file = trim(slices_dir)//"v_xy.txt" ); open(unit = v_xz, file = trim(slices_dir)//"v_xz.txt" ); open(unit = w_xy, file = trim(slices_dir)//"w_xy.txt" ); open(unit = w_xz, file = trim(slices_dir)//"w_xz.txt" ); open(unit = pi_xy, file = trim(slices_dir)//"pi_xy.txt" ); open(unit = pi_xz, file = trim(slices_dir)//"pi_xz.txt" ); open(unit = qv_xy, file = trim(slices_dir)//"qv_xy.txt" ); open(unit = qv_xz, file = trim(slices_dir)//"qv_xz.txt" ); open(unit = qc_xy, file = trim(slices_dir)//"qc_xy.txt" ); open(unit = qc_xz, file = trim(slices_dir)//"qc_xz.txt" ); open(unit = super_sat_xy, file = trim(slices_dir)//"super_sat_xy.txt" ); open(unit = super_sat_xz, file = trim(slices_dir)//"super_sat_xz.txt" ); end if open(unit = theta_l_profile, file = trim(profiles_dir)//"theta_l_profile.txt" ); open(unit = theta_profile, file = trim(profiles_dir)//"theta_profile.txt" ); open(unit = uvar_profile, file = trim(profiles_dir)//"uvar_profile.txt" ); open(unit = vvar_profile, file = trim(profiles_dir)//"vvar_profile.txt" ); open(unit = wvar_profile, file = trim(profiles_dir)//"wvar_profile.txt" ); open(unit = qv_profile, file = trim(profiles_dir)//"qv_profile.txt" ); open(unit = qc_profile, file = trim(profiles_dir)//"qc_profile.txt" ); open(unit = pi_profile, file = trim(profiles_dir)//"pi_profile.txt" ); open(unit = super_sat_profile,file = trim(profiles_dir)//"super_sat_profile.txt" ); open(unit = perturb_file, file = trim(exp_dir)//"perturb_file.txt" ); open(unit = tke_file, file = trim(exp_dir)//"tke.txt" ); open(unit = iter_stats_max, file = trim(exp_dir)//"iter_stats_max.txt" ); write(iter_stats_max,*) "#Time(s) Nu TKE(m2/s2) u_max(m/s) v_max(m/s) w_max(m/s) theta_max(K) qv_max(g/kg) qc_max(g/kg) & &super_sat_max(%) pi_max " open(unit = iter_stats_avg, file = trim(exp_dir)//"iter_stats_avg.txt" ); write(iter_stats_avg,*) "#Time(s) Nu TKE(m2/s2) u_avg(m/s) v_avg(m/s) w_avg(m/s) theta_avg(K) qv_avg(g/kg) qc_avg(g/kg) & &super_sat_avg(%) pi_avg" open(unit = iter_stats_std, file = trim(exp_dir)//"iter_stats_std.txt" ); write(iter_stats_std,*) "#Time(s) Nu TKE(m2/s2) u_std(m/s) v_std(m/s) w_std(m/s) theta_std(K) qv_std(g/kg) qc_std(g/kg) & &super_sat_std(%) pi_std" !Histogram data open(unit = theta_hist, file = trim(exp_dir)//"theta_hist.txt" ); write(theta_hist,'(A)') "#Following are the bin edges for theta in (K)" write(theta_hist,'(A,200F10.2)') "# ", theta_binEdges(1:n_theta_bins+1) write(theta_hist,'(A,I5,A,I5)') "# no. of bin edges = ", n_theta_bins+1, "; no. of cols = " , n_theta_bins + 2 write(theta_hist,'(A)') "# First col: data < lowest bin edge," write(theta_hist,'(A)') "# Last col: data >= largest bin edge" write(theta_hist,'(A)') "# Other col: left bin edge <= data < right bin edge" open(unit = qv_hist, file = trim(exp_dir)//"qv_hist.txt" ); write(qv_hist,'(A)') "#Following are the bin edges for qv in (g/kg)" write(qv_hist,'(A,200F10.2)') "# ", qv_binEdges(1:n_qv_bins+1)*1000; write(qv_hist,'(A,I5,A,I5)') "# no. of bin edges = ", n_qv_bins+1, "; no. of cols = " , n_qv_bins + 2 write(qv_hist,'(A)') "# First col: data < lowest bin edge," write(qv_hist,'(A)') "# Last col: data >= largest bin edge" write(qv_hist,'(A)') "# Other col: left bin edge <= data < right bin edge" open(unit = super_sat_hist, file = trim(exp_dir)//"super_sat_hist.txt" ); write(super_sat_hist,'(A)') "#Following are the bin edges for super saturation " write(super_sat_hist,'(A,200F10.2)') "# ", super_sat_binEdges(1:n_super_sat_bins+1); write(super_sat_hist,'(A,I5,A,I5)') "# no. of bin edges = ", n_super_sat_bins+1, "; no. of cols = " , n_super_sat_bins + 2 write(super_sat_hist,'(A)') "# First col: data < lowest bin edge," write(super_sat_hist,'(A)') "# Last col: data >= largest bin edge" write(super_sat_hist,'(A)') "# Other col: left bin edge <= data < right bin edge" end subroutine subroutine close_files() !use common_variables if (print_slices == 1) then close(thl_xy); close(thl_xz); close(u_xy); close(u_xz); close(v_xy); close(v_xz); close(w_xy); close(w_xz); close(pi_xy); close(pi_xz); close(super_sat_xy); close(super_sat_xz); end if close(perturb_file); close(tke_file); close(iter_stats_max); close(iter_stats_avg); close(theta_l_profile); close(theta_profile); close(uvar_profile); close(vvar_profile); close(wvar_profile); close(qv_profile); close(qc_profile); close(pi_profile); close(super_sat_profile); close(theta_hist ); close(qv_hist ); close(super_sat_hist); end subroutine subroutine compute_TKE() !Subroutine to compute avg kinetic at interior grid points only !use common_variables TKE = 0.5 * sum(u_at_theta(1:ni,1:nj,1:nk)**2 & + v_at_theta(1:ni,1:nj,1:nk)**2 & + w_at_theta(1:ni,1:nj,1:nk)**2)/(1.0*ni*nj*nk) !write(6,'(A,E16.6)') "------------------------------ TKE:" , TKE write(tke_file,'(E16.6)') TKE end subroutine compute_TKE subroutine print_iteration_stats() !Stats on evolution !use common_variables real*4 :: max_u, max_v, max_w, max_qv, max_qc, max_pi, max_super_sat, max_theta real*4 :: Nu,Flux_conv,Flux_cond real*4 :: avg_u, avg_v, avg_w, avg_qv, avg_qc, avg_pi, avg_super_sat, avg_theta avg_u = sum(u_at_theta(1:ni,1:nj,1:nk)) / (ni*nj*nk) avg_v = sum(v_at_theta(1:ni,1:nj,1:nk)) / (ni*nj*nk) avg_w = sum(w_at_theta(1:ni,1:nj,1:nk)) / (ni*nj*nk) avg_theta = sum( theta(1:ni,1:nj,1:nk )) / (ni*nj*nk) avg_qv = sum( qv(1:ni,1:nj,1:nk )) / (ni*nj*nk) avg_qc = sum( qc(1:ni,1:nj,1:nk )) / (ni*nj*nk) avg_pi = sum( pi(1:ni,1:nj,1:nk )) / (ni*nj*nk) avg_super_sat = sum(super_sat(1:ni,1:nj,1:nk )) / (ni*nj*nk) std_u = sqrt(sum((u_at_theta(1:ni,1:nj,1:nk) - avg_u )**2.) / (ni*nj*nk)) std_v = sqrt(sum((v_at_theta(1:ni,1:nj,1:nk) - avg_v )**2.) / (ni*nj*nk)) std_w = sqrt(sum((w_at_theta(1:ni,1:nj,1:nk) - avg_w )**2.) / (ni*nj*nk)) std_theta = sqrt(sum(( theta(1:ni,1:nj,1:nk) - avg_theta )**2.) / (ni*nj*nk)) std_qv = sqrt(sum(( qv(1:ni,1:nj,1:nk) - avg_qv )**2.) / (ni*nj*nk)) std_qc = sqrt(sum(( qc(1:ni,1:nj,1:nk) - avg_qc )**2.) / (ni*nj*nk)) std_pi = sqrt(sum(( pi(1:ni,1:nj,1:nk) - avg_pi )**2.) / (ni*nj*nk)) std_super_sat = sqrt(sum(( super_sat(1:ni,1:nj,1:nk) - avg_super_sat)**2.) / (ni*nj*nk)) max_u = maxval(abs(u_at_theta(1:ni,1:nj,1:nk))) max_v = maxval(abs(v_at_theta(1:ni,1:nj,1:nk))) max_w = maxval(abs(w_at_theta(1:ni,1:nj,1:nk))) max_theta = maxval( theta(1:ni,1:nj,1:nk) ) max_qv = maxval( qv(1:ni,1:nj,1:nk) ) max_qc = maxval( qc(1:ni,1:nj,1:nk) ) max_pi = maxval(abs( pi(1:ni,1:nj,1:nk))) max_super_sat = maxval( super_sat(1:ni,1:nj,1:nk) ) !Nusselt number = Total vertical heat flux in convection/ ( conduction flux in domain w/o convection) !In convection near surface total flux = local conduction flux, and !near middle of the domain total flux to local convective flux Flux_cond = delta_theta/H theta_avg_lvl1 = sum(theta(1:ni,1:nj,1))/(ni*nj) Flux_conv = (theta_bot - theta_avg_lvl1)/ (0.5 * dz) Nu = Flux_conv/Flux_cond !Print important stats write(6,'(A,E12.4)') "TKE(m2/s2) :", TKE write(6,'(A,F10.2)') "Nu :", Nu write(6,'(A,F10.2)') "Max u(m/s) :", max_u write(6,'(A,F10.2)') "Max v(m/s) :", max_v write(6,'(A,F10.2)') "Max w(m/s) :", max_w write(6,'(A,E12.4)') "Max pi :", max_pi write(6,'(A,F10.2)') "Avg qv(g/kg) :", avg_qv *1000 write(6,'(A,F10.2)') "Avg qc(g/kg) :", avg_qc *1000 write(6,'(A,F12.4)') "Avg super_sat(%) :", avg_super_sat*100 write(iter_stats_max,'(F10.2,E12.4,8F10.2,E12.4)') ITT*dt, TKE,& Nu, max_u,max_v,max_w,max_theta,max_qv*1000, max_qc*1000,max_super_sat*100,max_pi write(iter_stats_avg,'(F10.2,E12.4,8F10.2,E12.4)') ITT*dt, TKE,& Nu, avg_u,avg_v,avg_w,avg_theta,avg_qv*1000, avg_qc*1000,avg_super_sat*100,avg_pi write(iter_stats_std,'(F10.2,E12.4,8F10.2,E12.4)') ITT*dt, TKE,& Nu, std_u,std_v,std_w,std_theta,std_qv*1000, std_qc*1000,std_super_sat*100,std_pi end subroutine subroutine create_hist_binEdges() !use common_variables real*4 :: delta, minbin,maxbin integer :: tmp allocate(theta_binEdges(1:n_theta_bins+1)) allocate(qv_binEdges(1:n_qv_bins+1)) allocate(super_sat_binEdges(1:n_super_sat_bins+1)) minbin = theta_top; maxbin = theta_bot; delta = maxbin-minbin theta_binEdges = (/( tmp, tmp = 0,n_theta_bins,1 )/) theta_binEdges = minbin + theta_binEdges * delta/n_theta_bins print *, new_line("A"), "theta bin edges: " print '(200F10.2)' , theta_binEdges minbin = qv_sat_top; maxbin = qv_sat_bot delta = maxbin-minbin qv_binEdges = (/( tmp, tmp = 0,n_qv_bins,1 )/) qv_binEdges = minbin + qv_binEdges * delta/n_qv_bins print *, new_line("A"), "qv bin edges: " print '(200F10.2)' , qv_binEdges*1000 minbin = -1; maxbin = 1 delta = maxbin-minbin super_sat_binEdges = (/( tmp, tmp = 0,n_super_sat_bins,1 )/) super_sat_binEdges = minbin + super_sat_binEdges * delta/n_super_sat_bins print *, new_line("A"), "super_sat bin edges: " print '(200F10.2)' , super_sat_binEdges end subroutine subroutine print_histogram() !use common_variables real*4,dimension(0:n_theta_bins+1 ) :: freq_theta real*4,dimension(0:n_qv_bins+1 ) :: freq_qv real*4,dimension(0:n_super_sat_bins+1) :: freq_super_sat real*4,dimension(:),allocatable :: vals allocate(vals(1:ni*nj*nk)) vals = pack( theta(1:ni,1:nj,1:nk), .true. ) call compute_histogram(vals,ni*nj*nk,theta_binEdges,n_theta_bins+1,freq_theta) call print_1Darray(theta_hist,freq_theta) vals = pack(qv(1:ni,1:nj,1:nk), .true.) call compute_histogram(vals,ni*nj*nk,qv_binEdges,n_qv_bins+1,freq_qv) call print_1Darray(qv_hist, freq_qv) vals = pack(super_sat(1:ni,1:nj,1:nk), .true.) call compute_histogram(vals,ni*nj*nk,super_sat_binEdges,n_super_sat_bins+1,freq_super_sat) call print_1Darray(super_sat_hist, freq_super_sat) end subroutine subroutine print_1Darray (fid,var) !use common_variables integer :: fid,col real*4,dimension(0:) :: var if ( maxval(var) < 9999 )then write(fid,'(200f10.2 )') ( var(col),col = lbound(var,1)+0,ubound(var,1)-0 ) else write(fid,'(200E10.2 )') ( var(col),col = lbound(var,1)+0,ubound(var,1)-0 ) end if end subroutine print_1Darray subroutine print_2Darray (fid,var ) !Subroutine to print 2D grid. !use common_variables integer :: fid real*4,dimension(0:,0:) :: var integer :: row, col do row = lbound(var,2), ubound(var,2) if ( maxval(var) < 9999 )then write(fid,'(100F10.2 )') ( var(col,row),col = lbound(var,1)+0,ubound(var,1)-0 ) else !for debugging purpose if instability grows larger write(fid,'(100E10.2 )') ( var(col,row),col = lbound(var,1)+0,ubound(var,1)-0 ) end if end do end subroutine print_2Darray subroutine print_3Darray (fid,var ) !Print 3D array as horizontal slices of 2D(x,y) from bottom to top !use common_variables integer :: fid real*4,dimension(0:,0:,0:):: var integer :: i,j,k do k = lbound(var,3), ubound(var,3), 1 do j = lbound(var,2), ubound(var,2),1 !Before normalizing the values are so small, e-6, hence the format f15.8 if ( maxval(var) < 9999 )then write(fid,'(100F10.2 )') ( var(i,j,k),i = lbound(var,1)+0,ubound(var,1)-0 ) else write(fid,'(100E10.2 )') ( var(i,j,k),i = lbound(var,1)+0,ubound(var,1)-0 ) end if end do end do !Empty Line !write(fid,*) " " end subroutine print_3Darray subroutine print_xy_slice(fid,var) !use common_variables integer :: fid real*4,dimension(0:,0:,0:):: var real*4 :: slice_pctile = 10 integer :: bot_idx, top_idx, mid_idx bot_idx = ceiling(nk * slice_pctile/100.0) top_idx = nk - bot_idx + 1 mid_idx = nint((bot_idx + top_idx-1) / 2.0) !write(*,'(A,10I4)') "Horz slice idx:", bot_idx, mid_idx, top_idx call print_2Darray(fid,var(:,:,top_idx)) call print_2Darray(fid,var(:,:,mid_idx)) call print_2Darray(fid,var(:,:,bot_idx)) end subroutine subroutine print_xz_slice(fid,var) !To choose a Vertical slice at the same distance from front and back !I determine back_idx from front_idx rather than using ceiling or nint !use common_variables integer :: fid real*4,dimension(0:,0:,0:):: var real*4 :: slice_pctile = 10 integer :: front_idx, back_idx, center_idx front_idx = ceiling(nj * slice_pctile/100.0) back_idx = nj - front_idx + 1 center_idx = nint((front_idx + back_idx-1) / 2.0) !write(*,'(A,10I4)') "Vert slice idx:", front_idx, center_idx, back_idx call print_2Darray(fid,var(:,back_idx,:)) call print_2Darray(fid,var(:,center_idx,:)) call print_2Darray(fid,var(:,front_idx,:)) end subroutine subroutine print_avg_profile(fid,var) !use common_variables integer :: fid real*4,dimension(0:,0:,0:) :: var real*4,dimension(:),allocatable:: profile allocate( profile(0:nk+1)) profile = 0.0 do K = 1, nk profile(K) = sum( var(1:ni,1:nj,K) )/ (ni * nj) end do call print_1darray(fid,profile) end subroutine subroutine print_variance_profile(fid,var) !use common_variables integer :: fid real*4,dimension(0:,0:,0:) :: var real*4,dimension(:),allocatable:: variance real*4 :: mean allocate(variance(0:nk+1)) variance = 0.0 do K = 1, nk mean = sum( var(1:ni,1:nj,K) )/ (ni * nj) variance(K) = sum( (var(1:ni,1:nj,K) - mean)**2.0 ) / (ni*nj) end do call print_1darray(fid,variance) end subroutine subroutine output_variables() !use common_variables real*4 :: pi_max real*4,dimension(0:ni+1,0:nj+1,0:nk+1):: pi_norm if (print_slices == 1) then call print_xy_slice(thl_xy,theta_l) call print_xy_slice(u_xy,u_at_theta) call print_xy_slice(v_xy,v_at_theta) call print_xy_slice(w_xy,w_at_theta) call print_xy_slice(qv_xy,qv*1000) call print_xy_slice(qc_xy,qc*1000) call print_xy_slice(super_sat_xy,super_sat*100) call print_xz_slice(thl_xz,theta_l) call print_xz_slice(u_xz,u_at_theta) call print_xz_slice(v_xz,v_at_theta) call print_xz_slice(w_xz,w_at_theta) call print_xz_slice(qv_xz,qv*1000) call print_xz_slice(qc_xz,qc*1000) call print_xz_slice(super_sat_xz,super_sat*100) end if call print_avg_profile(theta_l_profile,theta_l) call print_avg_profile(theta_profile,theta) call print_avg_profile(qv_profile,qv*1000) call print_avg_profile(qc_profile,qc*1000) call print_avg_profile(super_sat_profile,super_sat*100) call print_variance_profile(uvar_profile,u_at_theta) call print_variance_profile(vvar_profile,v_at_theta) call print_variance_profile(wvar_profile,w_at_theta) pi_max = maxval(abs(pi)) if (pi_max /= 0 ) then pi_norm = pi/pi_max else pi_norm = pi end if call print_avg_profile(pi_profile,pi_norm) if (print_slices == 1) then call print_xy_slice(pi_xy,pi_norm) call print_xz_slice(pi_xz,pi_norm) end if end subroutine subroutine compute_histogram(vals,n_vals, binEdges, n_binEdges, freq) !Adapted from !https://pages.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/hist-2.html real*4, dimension(1:), intent(IN) :: vals integer, intent(IN) :: n_vals real*4,dimension(1:),intent(IN) :: binEdges integer, intent(IN) :: n_binEdges integer :: i, j real*4, dimension(0:) :: freq freq = 0.0 do i = 1, n_vals ! for each input do j = 1, n_binEdges ! determine the bin if (vals(i) < binEdges(j)) then freq(j-1) = freq(j-1) + 1 exit end if end do if (vals(i)>= binEdges(n_binEdges)) freq(n_binEdges) = freq(n_binEdges) + 1 end do end subroutine end module