! Mani Rajagopal ! This program is a modification of qcom_3d_w_swall project. ! This program is configurable through a file. ! Options available through file ! ni,nj,nk, Lx,Ly, H, ! Kth,Kvel,Kqw, alpha_damp, Cs ! time step ! Top and Bottom boundary conditions, ! Theta_top, Theta_bot, moisture: frac of saturated, velocity: free or no slip ! Sidewall = Theta, moisture flux, free slip or no slip ! Theta: insulated or fixed, moisture: insulated or fixed, velocity: free or no slip !output ! xy_slice, xz_slice program main ! pgf90 -c print.f ! pgf90 -o qcom QCOM.f90 print.o ! ./qcom use common_variables real*4 :: A, B !Adam Bashworth scheme integer :: N1,N2 real*4 :: start, finish integer :: i call cpu_time(start) !Load model configuration and params call model_setup() print *, "********** Experiment Name: ", trim(exp_name) , "*************" !Other model params ni_u = ni-side_wall; nj_v = nj-side_wall; nk_w = nk-1; sim_time = sim_hrs + sim_mins/60. + sim_secs/3600. ITTMAX = nint(sim_time * 3600 * (1.0/dt) ) output_itt = nint( output_intvl * (1/dt) ) delta_theta = theta_bot - theta_top dz = H/nk dx = Lx/ni dy = Ly/nj alpha_damp = alpha_damp_factor * Cs**2.0 * dt call allocate_memory() !Distance of theta grid points in x,y,z direction x = (/ (I, I = 0,ni+1) /) y = (/ (I, I = 0,nj+1) /) z = (/ (I, I = 0,nk+1) /) x = (x * dx) - dx/2 y = (y * dy) - dy/2 z = (z * dz) - dz/2 print '(A,F10.2)', "dx:", dx print '(A,F10.2)', "dy:", dy print '(A,F10.2)', "dz:", dz print '(A,F10.2)', "dt:", dt print '(A,F10.2)', "alpha_damp:", alpha_damp print '(A,I10)', "Max. no. of time steps :", ITTMAX print '(A,I10)', "Output every nth time_step :", output_itt print *, new_line("A"), "y position of theta grid points:" call print_1Darray(stdout,y) print *, new_line("A"), "z position of theta grid points:" call print_1Darray(stdout,z) !===Initialize with sounding, hydrostatic pressure, and boundary conditions CALL INIT !=== create histogram bins call create_hist_binEdges() !=== Open output files !Open file for writing data call open_files() !Verify initialization print *, new_line("A"), "Theta Initialization profile:" call print_1Darray(stdout,theta0) print *, new_line("A"), "Theta_v0 Initialization profile:" call print_1Darray(stdout,theta_v0) print *, new_line("A"), "Humidity Initialization profile:" call print_1Darray(stdout,qv0) print *, new_line("A"), "u wind Initialization profile:" call print_1Darray(stdout,u0) print *, new_line("A"), "v wind Initialization profile:" call print_1Darray(stdout,v0) !Verify initialization print *, new_line("A"), "Theta_l profile from 1st column:" call print_1Darray(stdout,theta_l(0,0,:)) print *, new_line("A"), "Theta profile from 1st column :" call print_1Darray(stdout,theta(0,0,:)) print *, new_line("A"), "total water profile from 1st column :" call print_1Darray(stdout,qw(0,0,:)*1000) print *, new_line("A"), "u wind profile from 1st column :" call print_1Darray(stdout,u(0,0,:)) print *, new_line("A"), "v wind profile from 1st column :" call print_1Darray(stdout,v(0,0,:)) !Output variables after Initialization call output_variables() call print_histogram() call print_iteration_stats() write(6,'(A,I5,A,F6.2,A)') "****************************** Iteration:" , ITT, " Time: " , ITT*dt, " secs" ! USE FORWARD SCHEME TO do first step ITT = 1 ! itt is time step index A = 1. B = 0. N1 = MOD ( ITT , 2 ) + 1 N2 = MOD ( ITT - 1, 2 ) + 1 CALL STEP ( N1, N2, A, B ) ! do first time step call compute_TKE() call print_iteration_stats() !ADAMS - BASHFORTH TWO - LEVEL SCHEME for time integration A = 3. / 2. B = - 1. / 2. ITTNOW = ITT + 1 DO ITT = ITTNOW, ITTMAX N1 = MOD ( ITT , 2 ) + 1 N2 = MOD ( ITT - 1, 2 ) + 1 CALL STEP ( N1, N2, A, B ) ! do subsequent time steps call compute_TKE() !Print variables every "output_itt" iteration if ( MOD (ITT,output_itt) == 0) then write(6,'(A,I7,A, F10.2, A)') "****************************** Iteration:" , ITT, " Time: " , ITT*dt, " secs" call output_variables() call print_histogram() call print_iteration_stats() !Flush stdout buffer flush(6) end if end do call output_variables() call close_files() call cpu_time(finish) print '(" Total CPU Time = ",f20.2," seconds.")',finish-start contains SUBROUTINE INIT !initialize all variables !use common_variables real*4 :: Amp,exner_pi_bot, T_swall, es_sat_swall real*4,dimension(0:ni+1,0:nj+1) :: perturb integer :: pidx !Boundary condition for top and bottom surfaces real*4 :: ES integer :: i,j,k !Initialize profiles theta0 = theta_bot - delta_theta/H * z qv0 = 0.0 theta_v0 = theta0 * (1 + 0.61 * qv0) u0 = 0.0 v0 = 0.0 !Test case - Make all variables as zero to figure out issue with equations !theta0 = 0.0 !qv0 = 0.0 !u0 = 0.0 !v0 = 0.0 !theta_v0 = 0.0 !End test case !Initialize the domain with profiles do j = 1, nj do i = 1, ni !Base state theta that includes vapor to determine buoyancy !Theta_l is the conserved variable in the moist c theta_l(i,j,:) = theta0 !Theta is a diagnostic variable theta(i,j,:) = theta0 !Moisture qv(i,j,:) = qv0 !Total water, qw (vapor + condensate ) is same as moisture profile since condensate = 0 qw(i,j,:) = qv0 end do end do !Create random perturbations close to surface at two levels pidx = 3 Amp = 0.5 call random_number(perturb) perturb = perturb * Amp !Perturbation near bottom theta_l(:,:, pidx) = theta_l(:,:,pidx) + perturb !level 1 theta_l(:,:, pidx-1) = theta_l(:,:,pidx-1) + perturb !level 2 !Perturbation near top theta_l(:,:, nk+1-pidx) = theta_l(:,:,nk+1-pidx) + perturb !level 1 theta_l(:,:, nk+1-pidx-1) = theta_l(:,:,nk+1-pidx-1) + perturb !level 2 !print *, "Perturbation:" call print_2Darray(perturb_file,perturb) !Set pressure to be same at all levels if (same_press == 1) then press0 = press_bot press_top = press_bot else exner_pi_bot = (press_bot/100000.)**(Rd/Cp) exner_pi0(0) = exner_pi_bot + ( g * (dz/2.0) ) /( cp * theta_v0(0) ) press0(0) = exner_pi0(0)**(Cp/Rd) * 100000 ! 1000 hpa standard pressure do k = 1, nk +1 exner_pi0(k) = exner_pi0(k-1) - (g * dz) / ( cp * 0.5*(theta_v0(k)+ theta_v0(k-1))) press0(k) = exner_pi0(k)**(Cp/Rd) * 100000 ! 1000 hpa end do endif print '(A30,f10.2 )', new_line("A")//"Pressure at the bottom surface:" , press_bot print '(A30,f10.2 )', "Pressure at the top surface :" , press_top print *, new_line("A"), "Profile of the base pressure in hpa" write(6,'(100f10.2)' ) ( press0/100 ) !Find saturation mixing ratio for top and bottom surfaces T_bot = theta_bot * (press_bot/100000.0)**(Rd/Cp) T_top = theta_top * (press_top/100000.0)**(Rd/Cp) es_sat_bot = es(T_bot) es_sat_top = es(T_top) qv_sat_bot = 0.622 * es_sat_bot / (press_bot - es_sat_bot) qv_sat_top = 0.622 * es_sat_top / (press_top - es_sat_top) !print *, "Theta(K) at bottom and top:" !write(stdout, '(2f8.2)') theta_bot,theta_top !print *, "Saturation es(pa) at bottom and top:" !write(stdout, '(2f10.2)') es_sat_bot,es_sat_top print *, new_line("A"), "Saturation qv (g/kg) at bottom and top:" write(stdout, '(2f8.2)') qv_sat_bot*1000,qv_sat_top*1000 !compute saturation mixing ratio at the side wall do k = 1, nk T_swall = theta_swall * (press0(k)/100000.0)**(Rd/Cp) es_sat_swall = es(T_swall) qv_sat_swall0(k) = 0.622 * es_sat_swall / (press0(k) - es_sat_swall) end do print *, new_line("A"), "Sidewall saturation mixing ratio(g/kg) profile" call print_1Darray(6,qv_sat_swall0*1000) call bound call diagnose END SUBROUTINE INIT SUBROUTINE STEP ( N1, N2, A, B ) integer :: N1,N2 real*4 :: A,B CALL RCALC ( N2 ) ! calculate forcing terms from variables at current time CALL AB ( N1, N2, A, B ) ! update variables using a time scheme CALL DIAGNOSE !determine values of diagnostic variables qv,qc, theta CALL BOUND END SUBROUTINE STEP SUBROUTINE RCALC ( N2 ) !CALCULATES FORCING TERMS on the right side of the equation !use common_variables integer :: N2,i,j,k real*4 :: fu_term1, fu_term2, fu_term3,fu_term4,fu_term5,fu_term6, fu_term7,fu_term8,fu_term9 real*4 :: fv_term1, fv_term2, fv_term3,fv_term4,fv_term5,fv_term6, fv_term7,fv_term8,fv_term9 real*4 :: fw_term1, fw_term2, fw_term3,fw_term4,fw_term5,fw_term6, fw_term7,fw_term8,fw_term9,fw_term10 real*4 :: fth_term1, fth_term2, fth_term3, fth_term4,fth_term5,fth_term6 real*4 :: fpi_term1, fpi_term2, fpi_term3, fpi_term4 real*4 :: fqw_term1, fqw_term2, fqw_term3, fqw_term4, fqw_term5, fqw_term6 real*4 :: fu_divdamp,fv_divdamp,fw_divdamp !set div_damp to zero beyond the domain boundary div_damp = 0.0 !Pressure perturbation equation as exner function pi do k=1,nk do j=1,nj do i=1,ni fpi_term1 = -1 * Cs**2.0 / (Cp * theta0(K)**2.0) !Horizontal gradient of theta_v0 is zero, Hence mean is same. fpi_term2 = (theta_v0(K) * u(I,J,K) - theta_v0(K) * u(I-1,J ,K))/dx fpi_term3 = (theta_v0(K) * v(I,J,K) - theta_v0(K) * v(I ,J-1,K))/dy fpi_term4 = ( 0.5 * (theta_v0(K+1) + theta_v0(K )) * w(I,J,K ) & - 0.5 * (theta_v0(K ) + theta_v0(K-1)) * w(I,J,K-1))/dz fpi(I,J,K,N2) = fpi_term1 * ( fpi_term2 + fpi_term3 + fpi_term4) !compute divergence damping term div_damp(I,J,K) = fpi_term2 + fpi_term3 + fpi_term4 end do end do end do !Periodic boundary is not required for div damping !because the loop is one iteration less for u,v,w in x,y,z direction respectively !compared to 3d domain without side walls !u momentum equation do k= 1, nk do j = 1, nj do i=1, ni_u !Pressure gradient term fu_term1 = -1.0 * Cp * theta_v0(K) * (pi(I+1,J,K) - pi(I,J,K)) / dx !Diffusion terms fu_term2 = Kvel * (u(I+1,J ,K) + u(I-1,J ,K) - 2. * u(I,J,K))/ dx**2.0 fu_term3 = Kvel * (u(I ,J+1,K) + u(I ,J-1,K) - 2. * u(I,J,K))/ dy**2.0 fu_term4 = Kvel * (u(I ,J ,K+1) + u(I ,J ,K-1) - 2. * u(I,J,K))/ dz**2.0 !Advection terms fu_term5 = -0.5*(u(I+1,J,K) + u(I-1,J,K)) * (u(I+1,J,K)- u(I-1,J,K))/(2.0*dx) fu_term6 = -0.5*(v(I+1,J ,K) + v(I,J ,K )) * (u(I,J+1,K ) - u(I,J ,K ))/ dy fu_term7 = -0.5*(v(I+1,J-1,K) + v(I,J-1,K )) * (u(I,J ,K ) - u(I,J-1,K ))/ dy fu_term8 = -0.5*(w(I+1,J ,K) + w(I,J ,K )) * (u(I,J ,K+1) - u(I,J ,K ))/ dz fu_term9 = -0.5*(w(I+1,J ,K-1) + w(I,J ,K-1)) * (u(I,J ,K ) - u(I,J ,K-1))/ dz !Divergence damping fu_divdamp = alpha_damp * (div_damp(I+1,J,K)-div_damp(I,J,K))/dx fu(I,J,K,N2)= fu_term1 + fu_term2 + fu_term3 + fu_term4 + fu_term5 & + 0.5*( fu_term6 + fu_term7) + 0.5*( fu_term8 + fu_term9) + fu_divdamp end do end do end do !v momentum equation do k= 1, nk do j = 1, nj_v do i=1, ni !Pressure gradient term fv_term1 = -1.0 * Cp * theta_v0(K) * (pi(I,J+1,K) - pi(I,J,K))/dy !Diffusion terms fv_term2 = Kvel * (v(I+1,J ,K) + v(I-1,J ,K ) - 2. * v(I,J,K))/ dx**2.0 fv_term3 = Kvel * (v(I ,J+1,K) + v(I ,J-1,K ) - 2. * v(I,J,K))/ dy**2.0 fv_term4 = Kvel * (v(I ,J ,K+1) + v(I ,J ,K-1) - 2. * v(I,J,K))/ dz**2.0 !Advection terms fv_term5 = -0.5*(u(I ,J+1,K ) + u(I ,J,K )) * (v(I+1,J,K ) - v(I ,J,K ))/ dx fv_term6 = -0.5*(u(I-1,J+1,K ) + u(I-1,J,K )) * (v(I ,J,K ) - v(I-1,J,K ))/ dx fv_term7 = -0.5*(v(I,J+1,K) + v(I,J-1,K+1))* ( v(I,J+1,K)- v(I,J-1,K))/(2.0*dy) fv_term8 = -0.5*(w(I ,J+1,K ) + w(I ,J,K )) * (v(I ,J,K+1) - v(I ,J,K ))/ dz fv_term9 = -0.5*(w(I ,J+1,K-1) + w(I ,J,K-1)) * (v(I ,J,K ) - v(I ,J,K-1))/ dz !Divergence damping fv_divdamp = alpha_damp * (div_damp(I,J+1,K)-div_damp(I,J,K))/dy fv(I,J,K,N2)= fv_term1 + fv_term2 + fv_term3 + fv_term4 + 0.5*( fv_term5 + fv_term6) & + fv_term7 + 0.5*( fv_term8 + fv_term9) + fv_divdamp end do end do end do !w momentum equation do k=1, nk_w do j=1, nj do i=1, ni !Pressure gradient term fw_term1 = -Cp * 0.5 * ( theta_v0(K+1) + theta_v0(K) ) * ( pi(I,J,K+1) - pi(I,J,K))/dz !Buoyancy terms fw_term2 = g * ( (theta(I,J,K+1) + theta(I,J,K))/ (theta0(K+1)+theta0(K)) - 1) & + g * 0.61 * ( 0.5*( qv(I,J,K+1)+ qv(I,J,K) ) - 0.5*( qv0(K+1)+qv0(K) ) )& - g * 0.5*( qc(I,J,K+1) + qc(I,J,K) ) !Diffusion terms fw_term3 = Kvel * ( w(I+1,J ,K) + w(I-1,J ,K) - 2. * w(I,J,K) ) / dx**2.0 fw_term4 = Kvel * ( w(I ,J+1,K) + w(I ,J-1,K) - 2. * w(I,J,K) ) / dy**2.0 fw_term5 = Kvel * ( w(I ,J ,K+1) + w(I ,J ,K-1) - 2. * w(I,J,K) ) / dz**2.0 !Advection terms fw_term6 = -0.5*( u(I ,J ,K+1) + u(I ,J ,K)) * (w(I+1,J ,K) - w(I ,J ,K))/dx fw_term7 = -0.5*( u(I-1,J ,K+1) + u(I-1,J ,K)) * (w(I ,J ,K) - w(I-1,J ,K))/dx fw_term8 = -0.5*( v(I ,J ,K+1) + v(I ,J ,K)) * (w(I ,J+1,K) - w(I ,J ,K))/dy fw_term9 = -0.5*( v(I ,J-1,K+1) + v(I ,J-1,K)) * (w(I ,J ,K) - w(I ,J-1,K))/dy fw_term10 = -0.5*(w(I,J,K+1) + w(I,J,K-1)) * (w(I,J,K+1) - w(I,J,K-1))/(2.0*dz) !Divergence damping fw_divdamp = alpha_damp * (div_damp(I,J,K+1)-div_damp(I,J,K))/dz fw(I,J,K,N2)= fw_term1 + fw_term2 + fw_term3 + fw_term4 + fw_term5 + 0.5*(fw_term6 + fw_term7) & + 0.5*(fw_term8 + fw_term9)+ fw_term10 + fw_divdamp end do end do end do !Theta_l equation do k=1,nk do j=1, nj do i=1,ni !Advection terms fth_term1 = -0.5*(u(I ,J ,K ) * (theta_l(I+1,J ,K ) - theta_l(I ,J ,K ))/ dx & + u(I-1,J ,K ) * (theta_l(I ,J ,K ) - theta_l(I-1,J ,K ))/ dx ) fth_term2 = -0.5*(v(I ,J ,K ) * (theta_l(I ,J+1,K ) - theta_l(I ,J ,K ))/ dy & + v(I ,J-1,K ) * (theta_l(I ,J ,K ) - theta_l(I ,J-1,K ))/ dy ) fth_term3 = -0.5*(w(I ,J ,K ) * (theta_l(I ,J ,K+1) - theta_l(I ,J ,K ))/ dz & + w(I ,J ,K-1) * (theta_l(I ,J ,K ) - theta_l(I ,J ,K-1))/ dz ) !Diffusion terms fth_term4 = Kth * ( theta_l(I+1,J ,K ) + theta_l(I-1,J ,K ) - 2. * theta_l(I,J,K)) / (dx**2.0) fth_term5 = Kth * ( theta_l(I ,J+1,K ) + theta_l(I ,J-1,K ) - 2. * theta_l(I,J,K)) / (dy**2.0) fth_term6 = Kth * ( theta_l(I ,J ,K+1) + theta_l(I ,J ,K-1) - 2. * theta_l(I,J,K)) / (dz**2.0) fth_l(I,J,K,N2) = fth_term1 + fth_term2 + fth_term3 + fth_term4 + fth_term5 + fth_term6 end do end do end do !Total water (vapor + condensate) equation do k=1,nk do j=1,nj do i=1,ni !Advection terms fqw_term1 = -0.5 * (u(I ,J ,K ) * (qw(I+1,J ,K ) - qw(I ,J ,K ))/dx & + u(I-1,J ,K ) * (qw(I ,J ,K ) - qw(I-1,J ,K ))/dx ) fqw_term2 = -0.5 * (v(I ,J ,K ) * (qw(I ,J+1,K ) - qw(I ,J ,K ))/dy & + v(I ,J-1,K ) * (qw(I ,J ,K ) - qw(I ,J-1,K ))/dy ) fqw_term3 = -0.5 * (w(I ,J ,K ) * (qw(I ,J ,K+1) - qw(I ,J ,K ))/dz & + w(I ,J ,K-1) * (qw(I ,J ,K ) - qw(I ,J ,K-1))/dz ) !Diffusion terms fqw_term4 = Kqw * (qw(I+1,J ,K ) + qw(I-1,J-1,K ) - 2.0 * qw(I,J,K))/(dx**2.0) fqw_term5 = Kqw * (qw(I ,J+1,K ) + qw(I ,J-1,K ) - 2.0 * qw(I,J,K))/(dy**2.0) fqw_term6 = Kqw * (qw(I ,J ,K+1) + qw(I ,J ,K-1) - 2.0 * qw(I,J,K))/(dz**2.0) fqw(I,J,K,N2) = fqw_term1 + fqw_term2 + fqw_term3 + fqw_term4 + fqw_term5 + fqw_term6 end do end do end do END SUBROUTINE RCALC SUBROUTINE AB ( N1, N2, A, B ) ! THE FOLLOWING LOOP UPDATES V USING EITHER THE FORWARD OR THE ! ADAMS-BASHFORTH ! SCHEME DEPENDING ON THE VALUES OF A, B. ! SUBSCRIPT N2 OF FV ALWAYS REFERS TO THE MOST RECENTLY CALCULATED ! VALUES FOR FV. !use common_variables integer :: N1,N2,i,j,k real*4 :: A,B do k = 1, nk do j = 1, nj do i = 1, ni theta_l(I,J,K) = theta_l(I,J,K) + dt * ( A * fth_l(I,J,K,N2) + B * fth_l(I,J,K,N1) ) pi(I,J,K) = pi(I,J,K) + dt * ( A * fpi(I,J,K,N2) + B * fpi(I,J,K,N1) ) qw(I,J,K) = qw(I,J,K) + dt * ( A * fqw(I,J,K,N2) + B * fqw(I,J,K,N1) ) end do end do end do do k = 1, nk do j =1,nj do i = 1,ni_u u(I,J,K) = u(I,J,K) + dt * ( A * fu(I,J,K,N2) + B * fu(I,J,K,N1) ) end do end do end do do k = 1, nk do j =1,nj_v do i = 1,ni v(I,J,K) = v(I,J,K) + dt * ( A * fv(I,J,K,N2) + B * fv(I,J,K,N1) ) end do end do end do do k = 1, nk_w do j =1,nj do i = 1,ni w(I,J,K) = w(I,J,K) + dt * ( A * fw(I,J,K,N2) + B * fw(I,J,K,N1) ) end do end do end do END SUBROUTINE AB subroutine bound !use common_variables call bc_top_n_bot if (side_wall == 1)then call bc_with_swall else call bc_wo_swall end if end subroutine subroutine bc_top_n_bot !use common_variables !BC for u,v if (free_slip == 1)then !At the top and bottom u(:,:,0 ) = u(:,:,1 ) u(:,:,nk+1) = u(:,:,nk) v(:,:,0 ) = v(:,:,1) v(:,:,nk+1) = v(:,:,nk) else !No Slip !At the top and bottom surface u(:,:,0 ) = -u(:,:,1 ) u(:,:,nk+1) = -u(:,:,nk) v(:,:,0 ) = -v(:,:,1) v(:,:,nk+1) = -v(:,:,nk) end if !BC for w : Rigid boundary w(:,:,0) = 0.0 w(:,:,nk_w+1) = 0.0 !BC for theta_l, !Conduction boundary at top and bottom theta_l(:,:,0) = 2. * theta_bot - theta_l(:,:,1) theta_l(:,:,nk+1) = 2. * theta_top - theta_l(:,:,nk) !BC for qw if (qw_bot > 0 .and. qw_bot < 1.01) then !Moisture as fraction (between 0,1) of saturated value qw(:,:,0 ) = 2. * qw_top * qv_sat_bot - qw(:,:,1) else !Dry convection with no moisture flux qw(:,:,0) = qw(:,:,1) end if if (qw_top > 0 .and. qw_top < 1.01) then !Moisture as fraction (between 0,1) of saturated value qw(:,:,nk+1) = 2. * qw_top * qv_sat_top - qw(:,:,nk) else qw(:,:,nk+1) = qw(:,:,nk) !Dry convection with no moisture flux end if !No boundary conditions for pi end subroutine subroutine bc_wo_swall !BC for sides : Periodic boundary !left to right !wo swall ni_u = ni, nj_v = nj u(0, :,:) = u(ni, :,:) u(ni+1,:,:) = u(1,:,:) v(0, :,:) = v(ni, :,:) v(ni+1,:,:) = v(1,:,:) w(0, :,:) = w(ni, :,:) w(ni+1,:,:) = w(1,:,:) theta_l(0, :,:) = theta_l(ni, :,:) theta_l(ni+1,:,:) = theta_l(1,:,:) qw(0, :,:) = qw(ni, :,:) qw(ni+1,:,:) = qw(1,:,:) pi(0, :,:) = pi(ni, :,:) pi(ni+1,:,:) = pi(1,:,:) !At front and back u(:,0, :) = u(:,nj, :) u(:,nj+1,:) = u(:,1, :) v(:,0, :) = v(:,nj, :) v(:,nj+1,:) = v(:,1, :) w(:,0, :) = w(:,nj, :) w(:,nj+1,:) = w(:,1, :) theta_l(:,0, :) = theta_l(:,nj, :) theta_l(:,nj+1,:) = theta_l(:,1, :) qw(:,0, :) = qw(:,nj, :) qw(:,nj+1,:) = qw(:,1, :) pi(:,0, :) = pi(:,nj, :) pi(:,nj+1,:) = pi(:,1, :) end subroutine subroutine bc_with_swall !BC for theta_l if (theta_swall > 0) then !Left and right side walls theta_l(0, :,:) = 2. * theta_swall - theta_l(1, :,:) theta_l(ni+1,:,:) = 2. * theta_swall - theta_l(ni,:,:) !Front and back walls theta_l(:,0, :) = 2. * theta_swall - theta_l(:,1, :) theta_l(:,nj+1,:) = 2. * theta_swall - theta_l(:,nj,:) else !Insulated side wall : No flux theta_l(0, :,:) = theta_l(1, :,:) theta_l(ni+1,:,:) = theta_l(ni,:,:) !Front and back walls theta_l(:,0, :) = theta_l(:,1, :) theta_l(:,nj+1,:) = theta_l(:,nj,:) endif !BC for velocities at sidewalls if (free_slip == 1) then !Free Slip at the boundary !At left and right v(0, :,:) = v(1, :,:) v(ni+1,:,:) = v(ni,:,:) w(0, :,:) = w(1, :,:) w(ni+1,:,:) = w(ni,:,:) !At front and back u(:,0, :) = u(:,1, :) u(:,nj+1,:) = u(:,nj,:) w(:,0, :) = w(:,1, :) w(:,nj+1,:) = w(:,nj,:) else !No Slip !At left and right side walls v(0, :,:) = -v(1, :,:) v(ni+1,:,:) = -v(ni,:,:) w(0, :,:) = -w(1, :,:) w(ni+1,:,:) = -w(ni,:,:) !At front and back side walls u(:,0, :) = -u(:,1, :) u(:,nj+1,:) = -u(:,nj,:) w(:,0, :) = -w(:,1, :) w(:,nj+1,:) = -w(:,nj,:) end if !Rigid boundary u(0,:,:) = 0.0 u(ni_u+1,:,:) = 0.0 v(:,0,:) = 0.0 v(:,nj_v+1,:) = 0.0 !BC for qw if (qw_swall == -1 ) then !Acts as moisture sink if gradient towards wall is negative else no flux !Left and right side wall do k = 1, nk do j=1, nj if ( qw(1,j,k) > qv_sat_swall0(k) ) then qw(0,j,k) = 2* qv_sat_swall0(k) - qw(1,j,k) else qw(0,j,k) = qw(1,j,k)!no flux endif if ( qw(ni, j,k) > qv_sat_swall0(k) ) then qw(ni+1,j,k) = 2* qv_sat_swall0(k) - qw(ni,j,k) else qw(ni+1,j,k) = qw(ni,j,k)!no flux endif end do end do !front and back side wall do k = 1, nk do i=1, ni if ( qw(i,1,k) > qv_sat_swall0(k) ) then qw(i,0,k) = 2* qv_sat_swall0(k) - qw(i,1,k) else qw(i,0,k) = qw(i,1,k)!no flux endif if ( qw(i, nj,k) > qv_sat_swall0(k) ) then qw(i,nj+1,k) = 2* qv_sat_swall0(k) - qw(i,nj,k) else qw(i,nj+1,k) = qw(i,nj,k)!no flux endif end do end do else !Insulated boundary: no flux qw(0, :,:) = qw(1, :,:) qw(ni+1,:,:) = qw(ni,:,:) !front and back side walls qw(:,0, :) = qw(:,1, :) qw(:,nj+1,:) = qw(:,nj,:) endif !No BC for pi END SUBROUTINE subroutine diagnose() real*4 :: temp,es_sat,qv_sat integer :: i,j,k !theta,qv,qc are diagnosed variables !Since there is no saturation adjustment, they are same theta = theta_l qv = qw !Find u,v,w at theta grid locations u_at_theta = 0.0 v_at_theta = 0.0 w_at_theta = 0.0 !Without sidewalls u has values from 0:ni+1 else with side walls from 0:ni; !ignore ni+1 for walls u_at_theta(1:ni,1:nj,1:nk) = 0.5 * ( u(0:ni-1,1:nj, 1:nk) + u(1:ni, 1:nj,1:nk) ) v_at_theta(1:ni,1:nj,1:nk) = 0.5 * ( v(1:ni ,0:nj-1,1:nk) + v(1:ni, 1:nj,1:nk) ) w_at_theta(1:ni,1:nj,1:nk) = 0.5 * ( w(1:ni ,1:nj, 0:nk-1) + w(1:ni, 1:nj,1:nk) ) !Compute supersaturation do k = 0, nk+1 do j = 0,nj+1 do i =0, ni+1 temp = theta(i,j,k) * (press0(K)/100000.0)**(Rd/Cp) es_sat = es(temp) qv_sat = 0.622 * es_sat / (press0(K) - es_sat) super_sat(i,j,k) = qv(i,j,k)/qv_sat - 1 end do end do end do end subroutine end program main function ES ( T ) !Function from Prof. Steve Krueger: taken from adjust.f !LOWE'S FORMULA FOR SATURATION VAPOR PRESSURE ( PA ). !T IS IN DEGREES KELVIN. real*4 :: ES,T real*4 :: TC,X real*4, dimension(1:7) :: A(7) = (/6.107800, & 4.436519E-01,& 1.428946E-02,& 2.650648E-04,& 3.031240E-06,& 2.034081E-08,& 6.136821E-11/) TC = T - 273.16 IF ( TC .LT. - 50. ) TC = - 50. X = A(7) do J = 1,6 X = X * TC + A(7-J) end do ES = X * 100. RETURN END function to_lower(strIn) result(strOut) ! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012) ! Original author: Clive Page character(len=256), intent(in) :: strIn character(len=256) :: strOut integer :: i,j do i = 1, len(strIn) j = iachar(strIn(i:i)) if (j>= iachar("A") .and. j<=iachar("Z") ) then strOut(i:i) = achar(iachar(strIn(i:i))+32) else strOut(i:i) = strIn(i:i) end if end do end function to_lower