module bgrid_advection_mod !----------------------------------------------------------------------- ! ! performs vertical and horizontal advection ! and negative tracer borrowing/filling ! !----------------------------------------------------------------------- !--------------------------- modules ----------------------------------- !----------------------------------------------------------------------- use bgrid_horiz_mod , only: horiz_grid_type use bgrid_vert_mod , only: vert_grid_type, compute_pres_depth use bgrid_masks_mod , only: grid_mask_type use bgrid_prog_var_mod , only: prog_var_type use bgrid_polar_filter_mod, only: pfilt_control_type, polar_filter, & polar_filter_wind, TGRID use bgrid_halo_mod , only: update_halo, vel_flux_boundary, & EAST, WEST, SOUTH, NORTH, NOPOLE, & TEMP, UWND, VWND, WIND, POLEONLY use bgrid_change_grid_mod , only: change_grid, EQUAL, AREA, & UFLX_GRID, VFLX_GRID, TEMP_GRID, WIND_GRID use vert_advection_mod , only: vert_advection, vert_advection_end, & SECOND_CENTERED, FOURTH_CENTERED, & SECOND_CENTERED_WTS, FOURTH_CENTERED_WTS, & FINITE_VOLUME_LINEAR, FINITE_VOLUME_PARABOLIC, & FLUX_FORM, WEIGHTED_TENDENCY !!!use horiz_advection_mod , only: horiz_advection, & !!! FINITE_VOLUME_LINEAR_HORIZ=>FINITE_VOLUME_LINEAR, & !!! FINITE_VOLUME_PARABOLIC_HORIZ=>FINITE_VOLUME_PARABOLIC use fms_mod, only: error_mesg, FATAL, write_version_number, & file_exist, open_namelist_file, & check_nml_error, close_file, & mpp_pe, mpp_root_pe, stdlog, uppercase, & mpp_clock_id, mpp_clock_begin, mpp_clock_end, & MPP_CLOCK_SYNC, CLOCK_MODULE, CLOCK_ROUTINE, CLOCK_LOOP use mpp_mod, only: input_nml_file, mpp_max use field_manager_mod, only: MODEL_ATMOS, parse use tracer_manager_mod, only: query_method, get_tracer_names, get_number_tracers !----------------------------------------------------------------------- implicit none private public :: advection_init, advection, advection_end !----------------------------------------------------------------------- !------------ namelist: bgrid_advection_nml ------------- ! vert_advec_scheme_wind The vertical advection scheme. ! vert_advec_scheme_temp Possible values are NONE, SECOND_CENTERED, FOURTH_CENTERED, ! vert_advec_scheme_tracer FINITE_VOLUME_LINEAR, or FINITE_VOLUME_PARABOLIC. ! Using finite volume schemes for momentum has not been ! tested and may produce poor results. character(len=24) :: vert_advec_scheme_wind = 'SECOND_CENTERED' character(len=24) :: vert_advec_scheme_temp = 'SECOND_CENTERED' character(len=24) :: vert_advec_scheme_tracer = 'SECOND_CENTERED' ! horiz_advec_scheme_wind The horizontal advection scheme. ! horiz_advec_scheme_temp Possible values are NONE, SECOND_CENTERED, FOURTH_CENTERED. ! horiz_advec_scheme_tracer character(len=24) :: horiz_advec_scheme_wind = 'SECOND_CENTERED' character(len=24) :: horiz_advec_scheme_temp = 'SECOND_CENTERED' character(len=24) :: horiz_advec_scheme_tracer = 'SECOND_CENTERED' ! advec_weight_wind Weights used for modified Euler-backward time differencing ! advec_weight_temp (i.e., when the scheme is SECOND_CENTERED or FOURTH_CENTERED). ! advec_weight_tracer 0.0 = full Euler-forward (not recommended) ! 1.0 = full Euler-backward real :: advec_weight_wind = 0.7 real :: advec_weight_temp = 0.7 real :: advec_weight_tracer = 0.7 ! num_fill_pass The number of successive repetitions of the tracer borrowing scheme. ! This value applies to both the horizontal and vertical schemes. integer :: num_fill_pass = 1 ! temporary undocumented developer flags logical :: compute_vert_wind_flux = .false. character(len=16) :: vert_vel_flux = 'area_weight' character(len=16) :: horiz_vel_flux = 'equal_weight' namelist /bgrid_advection_nml/ horiz_advec_scheme_wind, vert_advec_scheme_wind, & horiz_advec_scheme_temp, vert_advec_scheme_temp, & horiz_advec_scheme_tracer, vert_advec_scheme_tracer, & advec_weight_wind, advec_weight_temp, advec_weight_tracer, & num_fill_pass ,& compute_vert_wind_flux, vert_vel_flux, horiz_vel_flux !----------------------------------------------------------------------- !----- private data ----- ! derived-type containing tracer filling parameters type trfill_control_type integer, pointer :: fill_scheme(:) =>NULL(), & npass_horiz(:) =>NULL(), & npass_vert(:) =>NULL() end type trfill_control_type ! derived-type containing advection parameters type advec_control_type type(trfill_control_type) :: Fill integer, pointer :: scheme(:,:) =>NULL() real , pointer :: weight(:) =>NULL() logical :: do_mask4_tmp, do_mask4_vel, & do_finite_volume_tmp, do_finite_volume_trs end type advec_control_type integer, parameter :: NONE = 7000, & ! may apply to advection or filling ! filling schemes EQUAL_BORROW = 8005, & BOTTOM_BORROW = 8006, & ! not implemented GLOBAL_BORROW = 8007 ! not implemented ! advection control parameters ! set via values in namelist or field table type(advec_control_type),save :: Control real, parameter :: c4 = -1./6. ! weights for 4th order centered schemes real, parameter :: c2 = 1. - 2.*c4 logical :: do_log = .true. logical :: stability_check = .false. ! perform stability check for horizontal ! finite volume (Lin-Rood) schemes ! performance timing info integer :: id_tmp, id_trs, id_vel, id_fill logical :: do_clock_init = .true. character(len=128) :: version='$Id: bgrid_advection.F90,v 19.0 2012/01/06 19:52:53 fms Exp $' character(len=128) :: tagname='$Name: tikal $' !----------------------------------------------------------------------- contains !####################################################################### subroutine advection ( Pfilt, Hgrid, Vgrid, Masks, & pfilt_opt, dt, dpde_old, dpde, & few, fns, etadot, u, v, t, Var, Var_dt ) type(pfilt_control_type), intent(in) :: Pfilt type(horiz_grid_type), intent(inout) :: Hgrid type (vert_grid_type), intent(in) :: Vgrid type (grid_mask_type), intent(in) :: Masks integer, intent(in) :: pfilt_opt real, intent(in) :: dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: & dpde_old, dpde, u, v, t real, intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: & few, fns, etadot type (prog_var_type), intent(in) :: Var type (prog_var_type), intent(inout) :: Var_dt !----------------------------------------------------------------------- ! ! Pfilt = polar filter constants ! Hgrid = horizontal grid constants ! Vgrid = vertical grid constants ! Masks = grid masking constants ! pf_opt = polar filter option flag ! dt = adjustment time step ! dpde_old = pressure thickness of model layers at the end of the ! last advective time step ! dpde = current pressure thickness of model layers ! few, fns = zonal, meridional mass fluxes (Pa-m2/s) ! (a summation since the last call to advection) ! etadot = vertical mass flux (Pa/s) (summation since the last call to advection) ! u, v, t = prognostic variables at the end of the last advective time ! step, note that tracers have not been updated since the last ! advective time step, therefore, r = Var%r + dt*Var_dt%r ! Var = prognostic variables at the end of the last dynamics time ! step, current values would be, uc = Var%u + dt*Var_dt%u ! Var_dt = prognostic variable tendencies, accumulated since the ! variable were updated in Var ! !----------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & Vgrid%nlev) :: dpdt, dpde_xy, r, uc, vc real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & Vgrid%nlev,2) :: mask4 integer :: i, j, k, n, is, ie, js, je !----------------------------------------------------------------------- ! ---- update halos for fluxes ---- call update_halo (Hgrid, TEMP, few) call update_halo (Hgrid, VWND, fns) ! ---- compute pressure tendency in each layer ---- dpdt = (dpde - dpde_old) / dt ! ---- stability diagnostic (CFL for horizontal finite volume schemes) ---- if (stability_check) call stability_diag ( Hgrid, dt, few, dpde_old, dpde ) !----------------------------------------------------------------------- !------------------ advection of mass fields --------------------------- !----------------------------------------------------------------------- !------ initialize fourth-order mask ---- if (Control%do_mask4_tmp) & call mask4_init (Hgrid, 1, Masks%sigma, Masks%Tmp%mask, mask4) ! compute wind components for horizontal finite-volume advection schemes if (Control%do_finite_volume_tmp .or. Control%do_finite_volume_trs) then call compute_advecting_wind (Hgrid, dpde_old, dpde, few, fns, uc, vc) else uc = 0. vc = 0. endif !------ temperature advection ------ call mpp_clock_begin (id_tmp) if (Masks%sigma) then call advect_mass (Hgrid, Pfilt, Control % scheme(:,0), & Control % weight(0), pfilt_opt, & dt, dpdt, dpde, few, fns, etadot, & uc, vc, Var%t, t, Var_dt%t) else call advect_mass (Hgrid, Pfilt, Control % scheme(:,0), & Control % weight(0), pfilt_opt, & dt, dpdt, dpde, few, fns, etadot, & uc, vc, Var%t, t, Var_dt%t, & mask=Masks%Tmp%mask, mask4=mask4) endif call mpp_clock_end (id_tmp) !------ tracer advection ----------- ! (need to pass all tracers to update pressure tendency part) call mpp_clock_begin (id_trs) do n = 1, Var_dt%ntrace r = Var%r(:,:,:,n) + dt*Var_dt%r(:,:,:,n) if (Masks%sigma) then call advect_mass (Hgrid, Pfilt, Control%scheme(:,n), & Control % weight(n), & pfilt_opt, dt, dpdt, dpde, & few, fns, etadot, uc, vc, & Var%r(:,:,:,n), r, Var_dt%r(:,:,:,n) ) else call advect_mass (Hgrid, Pfilt, Control%scheme(:,n), & Control % weight(n), & pfilt_opt, dt, dpdt, dpde, & few, fns, etadot, uc, vc, & Var%r(:,:,:,n), r, Var_dt%r(:,:,:,n),& mask=Masks%Tmp%mask, mask4=mask4 ) endif enddo call mpp_clock_end (id_trs) !------ tracer hole filling ------- ! (remove negative tracer values with vert/horiz borrowing) call mpp_clock_begin (id_fill) call vert_borrow ( dt, dpde, Var%r(:,:,:,1:Var_dt%ntrace), & Var_dt%r(:,:,:,1:Var_dt%ntrace), & iters = Control % Fill % npass_vert ) if (Masks%sigma) then call horiz_borrow ( Hgrid, dt, dpde, & Var%r(:,:,:,1:Var_dt%ntrace), & Var_dt%r(:,:,:,1:Var_dt%ntrace), & iters = Control % Fill % npass_horiz ) else call horiz_borrow ( Hgrid, dt, dpde, & Var%r(:,:,:,1:Var_dt%ntrace), & Var_dt%r(:,:,:,1:Var_dt%ntrace), & mask = Masks%Tmp%mask, & iters = Control % Fill % npass_horiz ) endif call mpp_clock_end (id_fill) !----------------------------------------------------------------------- !------------------ advection of momentum fields ----------------------- !----------------------------------------------------------------------- call mpp_clock_begin (id_vel) is = Hgrid % Vel % is; ie = Hgrid % Vel % ie js = Hgrid % Vel % js; je = Hgrid % Vel % je !--- mass fluxes between velocity points --- if (trim(horiz_vel_flux(1:5)) == 'equal') then call change_grid (Hgrid, UFLX_GRID, VFLX_GRID, few, few, weight=EQUAL) call change_grid (Hgrid, VFLX_GRID, UFLX_GRID, fns, fns, weight=EQUAL) else call change_grid (Hgrid, UFLX_GRID, VFLX_GRID, few, few, weight=AREA) call change_grid (Hgrid, VFLX_GRID, UFLX_GRID, fns, fns, weight=AREA) endif ! compute vertical flux for momentum from horizonal momentum fluxes if (compute_vert_wind_flux) then call compute_etadot_vel (Hgrid, Vgrid, Masks, few, fns, etadot ) endif ! no flux across pole call vel_flux_boundary (Hgrid, fns) !--- interpolate mass fields to momentum grid --- if (.not.compute_vert_wind_flux) then if (trim(vert_vel_flux(1:4)) == 'area') then call change_grid (Hgrid, TEMP_GRID, WIND_GRID, etadot, etadot, weight=AREA) else call change_grid (Hgrid, TEMP_GRID, WIND_GRID, etadot, etadot, weight=EQUAL) endif endif ! always area-weighted (the default) call change_grid (Hgrid, TEMP_GRID, WIND_GRID, dpde, dpde_xy) call change_grid (Hgrid, TEMP_GRID, WIND_GRID, dpdt, dpdt) ! advection of momentum ! determine whether step-mountain mask are needed if (Control%do_mask4_vel .and. .not.Masks%sigma) then call mask4_init (Hgrid, 2, Masks%sigma, Masks%Vel%mask, mask4) call advect_vel (Hgrid, Pfilt, Control % scheme(:,-1), & Control % weight(-1), & pfilt_opt, dt, dpdt, dpde_xy, & few, fns, etadot, Var%u, Var%v, u, v, & Var_dt%u, Var_dt%v, & mask=Masks%Vel%mask, mask4=mask4) else ! sigma case call advect_vel (Hgrid, Pfilt, Control % scheme(:,-1), & Control % weight(-1), & pfilt_opt, dt, dpdt, dpde_xy, & few, fns, etadot, Var%u, Var%v, u, v, & Var_dt%u, Var_dt%v) endif call mpp_clock_end (id_vel) !----------------------------------------------------------------------- !------ done with fluxes - zero-out ??? ------- few = 0.0 fns = 0.0 etadot = 0.0 !----------------------------------------------------------------------- end subroutine advection !####################################################################### subroutine advect_mass ( Hgrid, Pfilt, scheme, coeff, fopt, dt, & dpdt, dpn, few, fns, fud, uc, vc, & ro, r, r_dt, mask, mask4 ) type(horiz_grid_type), intent(inout) :: Hgrid type(pfilt_control_type), intent(in) :: Pfilt integer, intent(in) :: scheme(2) integer, intent(in) :: fopt real, intent(in) :: coeff, dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dpdt, dpn,& few, fns, fud, ro, r, uc, vc real, intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: r_dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:), optional :: mask real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:), optional :: mask4 !----------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(r,3)) :: rst, rst_dt, rst_dt_v, rdpdt real :: rcoef integer :: pass, npass, j, k integer :: vert_scheme, horiz_scheme integer, parameter :: FINITE_VOLUME = 1234567 !can not match other numbers? !----------------------------------------------------------------------- ! set the vertical and horizontal differencing scheme horiz_scheme = scheme(1) vert_scheme = scheme(2) ! rename horizontal finite volume schemes to correct horizontal name ! if (horiz_scheme == FINITE_VOLUME_LINEAR) horiz_scheme = FINITE_VOLUME_LINEAR_HORIZ ! if (horiz_scheme == FINITE_VOLUME_PARABOLIC) horiz_scheme = FINITE_VOLUME_PARABOLIC_HORIZ !----------------------------------------------------------------------- rst = ro + dt * r_dt ! use updated value rdpdt = r*dpdt rcoef = coeff npass = 1 if (rcoef > 1.e-4) npass = 2 ! no advection, psdt correction only if (vert_scheme == NONE) then rst_dt_v = 0.0 if (horiz_scheme == NONE) npass = 1 endif ! both vert and horiz are using finite volume if ( (vert_scheme == FINITE_VOLUME_LINEAR .or. & vert_scheme == FINITE_VOLUME_PARABOLIC) .and. & (horiz_scheme == FINITE_VOLUME_LINEAR .or. & horiz_scheme == FINITE_VOLUME_PARABOLIC) ) then !(horiz_scheme == FINITE_VOLUME_LINEAR_HORIZ .or. & ! horiz_scheme == FINITE_VOLUME_PARABOLIC_HORIZ) ) then npass = 2 rcoef = 1. endif do pass = 1, npass !---- vertical tendency ---- ! finite volume schemes on pass 1 only if (vert_scheme == SECOND_CENTERED .or. vert_scheme == SECOND_CENTERED_WTS .or. & vert_scheme == FOURTH_CENTERED .or. vert_scheme == FOURTH_CENTERED_WTS .or. & ((vert_scheme == FINITE_VOLUME_LINEAR .or. & vert_scheme == FINITE_VOLUME_PARABOLIC) .and. pass == 1)) then call vert_advection ( dt, fud, dpn, rst, rst_dt_v, & mask=mask, scheme=vert_scheme, & form=FLUX_FORM, flags=WEIGHTED_TENDENCY ) endif !---- horizontal tendency ---- if (horiz_scheme == SECOND_CENTERED) then call advect_mass_horiz_2nd (Hgrid, few, fns, rst, rst_dt) if ( fopt >= 1 ) call polar_filter (Pfilt, rst_dt, TGRID, mask) call update_halo (Hgrid, TEMP, rst_dt) else if (horiz_scheme == FOURTH_CENTERED) then call advect_mass_horiz_4th (Hgrid, few, fns, rst, rst_dt, mask4) if ( fopt >= 1 ) call polar_filter (Pfilt, rst_dt, TGRID, mask) call update_halo (Hgrid, TEMP, rst_dt) ! compute horizontal finite volume scheme on 2nd pass else if ((horiz_scheme == FINITE_VOLUME_LINEAR .or. & horiz_scheme == FINITE_VOLUME_PARABOLIC) .and. pass == 2) then ! temporary error check call error_mesg ('bgrid_advection_mod', & 'horizontal finite volume schemes & &are not implemented with this release', FATAL) !call horiz_advection ( Hgrid%Tmp%Domain, dt, Hgrid%Tmp%dx, Hgrid%Tmp%dy, & ! uc, vc, few, fns, rst, rst_dt, scheme=horiz_scheme ) ! halos updated by finite volume scheme ! need update only poles (model b.c. differs) call update_halo (Hgrid, TEMP, rst_dt, flags=POLEONLY) else rst_dt = 0.0 endif !---- combine vertical and horizontal tendencies ---- !---- adjust for pressure tendency ---- rst_dt = (rst_dt + rst_dt_v - rdpdt) / dpn ! apply step-mountain mask? if (present(mask)) rst_dt = rst_dt*mask !---- compute new value or return tendency ---- if (pass < npass) then rst = ro + dt * (r_dt + rcoef*rst_dt) else r_dt = r_dt + rst_dt endif enddo !----------------------------------------------------------------------- end subroutine advect_mass !####################################################################### subroutine advect_vel ( Hgrid, Pfilt, scheme, coeff, fopt, dt, & dpdt, dpn, few, fns, fud, & uo, vo, u, v, u_dt, v_dt, mask, mask4 ) type(horiz_grid_type), intent(inout) :: Hgrid type(pfilt_control_type), intent(in) :: Pfilt integer, intent(in) :: scheme(2) integer, intent(in) :: fopt real, intent(in) :: coeff, dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dpdt, dpn,& few, fns, fud, uo, vo, u, v real, intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: u_dt, v_dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:), optional :: mask real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:), optional :: mask4 !----------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(u,3)) :: ust_dt, ust_dt_v, & vst_dt, vst_dt_v ! store ust & vst components together to create larger halo updates real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(u,3),2) :: uvst integer, parameter :: UCOMP=1, VCOMP=2 integer :: vert_scheme, horiz_scheme integer :: pass, npass, i, j, k, is, ie, js, je, nlev !----------------------------------------------------------------------- ! set the vertical differencing scheme horiz_scheme = scheme(1) vert_scheme = scheme(2) !----------------------------------------------------------------------- is = Hgrid%Vel%is; ie = Hgrid%Vel%ie js = Hgrid%Vel%js; je = Hgrid%Vel%je nlev = size(u,3) uvst(:,:,:,UCOMP) = uo + dt * u_dt ! use updated values uvst(:,:,:,VCOMP) = vo + dt * v_dt ! tendency halos are not up-to-date ! need to update halos call update_halo ( Hgrid, WIND, uvst ) npass = 1 if (coeff > 1.e-4) npass = 2 do pass = 1, npass !---- vertical tendency ---- ! finite volume schemes on pass 1 only if (vert_scheme == SECOND_CENTERED .or. vert_scheme == SECOND_CENTERED_WTS .or. & vert_scheme == FOURTH_CENTERED .or. vert_scheme == FOURTH_CENTERED_WTS .or. & ((vert_scheme == FINITE_VOLUME_LINEAR .or. & vert_scheme == FINITE_VOLUME_PARABOLIC) .and. pass == 1)) then call vert_advection ( dt, fud, dpn, uvst(:,:,:,UCOMP), ust_dt_v, & mask=mask, scheme=vert_scheme, & form=FLUX_FORM, flags=WEIGHTED_TENDENCY ) call vert_advection ( dt, fud, dpn, uvst(:,:,:,VCOMP), vst_dt_v, & mask=mask, scheme=vert_scheme, & form=FLUX_FORM, flags=WEIGHTED_TENDENCY ) else if (vert_scheme == NONE .and. pass == 1) then ust_dt_v = 0.0 vst_dt_v = 0.0 endif !---- horizontal tendency ---- if (horiz_scheme == SECOND_CENTERED) then call advect_vel_horiz_2nd (Hgrid, few, fns, uvst(:,:,:,UCOMP), & uvst(:,:,:,VCOMP), ust_dt, vst_dt ) else if (horiz_scheme == FOURTH_CENTERED) then call advect_vel_horiz_4th (Hgrid, few, fns, & uvst(:,:,:,UCOMP), uvst(:,:,:,VCOMP), & ust_dt, vst_dt, mask4 ) endif ! polar filter horizontal tendency if (fopt >= 2) then call polar_filter_wind (Pfilt, ust_dt, vst_dt, mask) endif !---- combine vertical and horizontal tendencies ---- !---- adjust for pressure tendency ---- do k = 1, nlev do j = js, je do i = is, ie ust_dt(i,j,k) = (ust_dt(i,j,k) + ust_dt_v(i,j,k) - & u(i,j,k)*dpdt(i,j,k)) / dpn(i,j,k) vst_dt(i,j,k) = (vst_dt(i,j,k) + vst_dt_v(i,j,k) - & v(i,j,k)*dpdt(i,j,k)) / dpn(i,j,k) enddo enddo enddo ! apply step-mountain mask? if (present(mask)) then do k = 1, nlev ust_dt(is:ie,js:je,k) = ust_dt(is:ie,js:je,k)*mask(is:ie,js:je,k) vst_dt(is:ie,js:je,k) = vst_dt(is:ie,js:je,k)*mask(is:ie,js:je,k) enddo endif !---- compute new value or return tendency ---- if (pass < npass) then do k = 1, nlev do j = js, je do i = is, ie uvst(i,j,k,UCOMP) = uo(i,j,k) + dt * (u_dt(i,j,k) + coeff*ust_dt(i,j,k)) uvst(i,j,k,VCOMP) = vo(i,j,k) + dt * (v_dt(i,j,k) + coeff*vst_dt(i,j,k)) enddo enddo enddo call update_halo (Hgrid, WIND, uvst) else do k = 1, nlev do j = js, je do i = is, ie u_dt(i,j,k) = u_dt(i,j,k) + ust_dt(i,j,k) v_dt(i,j,k) = v_dt(i,j,k) + vst_dt(i,j,k) enddo enddo enddo ! NOTE: halos not updated for output tendencies ! do this in the dynamical core for efficiency endif enddo !----------------------------------------------------------------------- end subroutine advect_vel !####################################################################### subroutine advect_mass_horiz_2nd (Hgrid, few, fns, r, rdt) type(horiz_grid_type), intent(in) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, fns, r real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: rdt real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) :: frew, frns integer :: i, j, k ! second order centered differencing on B-grid temperature grid ! constant grid spacing assumed ! minimum halo size = 1 do k = 1, size(r,3) !--- horizontal fluxes --- do j = Hgrid%Tmp%js-1, Hgrid%Tmp%je do i = Hgrid%Tmp%is-1, Hgrid%Tmp%ie frew(i,j) = few(i,j,k) * (r(i+1,j,k) + r(i,j,k)) frns(i,j) = fns(i,j,k) * (r(i,j+1,k) + r(i,j,k)) enddo enddo !--- horizontal advective tendency --- do j = Hgrid%Tmp%js, Hgrid%Tmp%je do i = Hgrid%Tmp%is, Hgrid%Tmp%ie rdt(i,j,k) = -0.5*Hgrid%Tmp%rarea(j)* & ((frew(i ,j)+frns(i,j )) & -(frew(i-1,j)+frns(i,j-1))) enddo enddo enddo end subroutine advect_mass_horiz_2nd !####################################################################### subroutine advect_mass_horiz_4th (Hgrid, few, fns, r, rdt, mask4) type(horiz_grid_type), intent(inout) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, fns, r real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: rdt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:), optional :: & mask4 real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(r,3)) :: rew, rns, frew, frns real, dimension(Hgrid%ilb:Hgrid%iub) :: x2, x4 real, dimension(Hgrid%ilb:Hgrid%iub) :: y2, y4 integer :: i, j, k, is, ie, js, je ! fourth order centered differencing on B-grid temperature grid ! constant grid spacing assumed ! assumed halo size = 1 is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie js = Hgrid % Tmp % js; je = Hgrid % Tmp % je rns(:,Hgrid%jub,:) = 0.0 do k = 1, size(r,3) do j = Hgrid%jlb,Hgrid%jub do i = Hgrid%ilb,Hgrid%iub-1 rew(i,j,k) = r(i+1,j,k)+r(i,j,k) enddo enddo do j = Hgrid%jlb,Hgrid%jub-1 do i = Hgrid%ilb,Hgrid%iub rns(i,j,k) = r(i,j+1,k)+r(i,j,k) enddo enddo enddo ! assumed halosize = 1 call update_halo (Hgrid,TEMP,rew,halos=EAST) call update_halo (Hgrid,UWND,rns,halos=NORTH+NOPOLE) ! ---- horizontal fluxes ---- if (present(mask4)) then frew=0.0; frns=0.0 do k = 1, size(r,3) do j = js, je do i = Hgrid%ilb+1, ie x4(i) = mask4(i,j,k,1) x2(i) = 1.0 - 2.*x4(i) frew(i,j,k) = few(i,j,k) * ( x2(i) * rew(i,j,k) & + x4(i) * (rew(i-1,j,k) + rew(i+1,j,k)) ) enddo enddo do j = Hgrid%jlb+1, je do i = is, ie y4(i) = mask4(i,j,k,2) y2(i) = 1.0 - 2.*y4(i) frns(i,j,k) = fns(i,j,k) * ( y2(i) * rns(i,j,k) & + y4(i) * (rns(i,j-1,k) + rns(i,j+1,k)) ) enddo enddo enddo else frew=0.0; frns=0.0 do k = 1, size(r,3) do j = js, je do i = Hgrid%ilb+1, ie frew(i,j,k) = few(i,j,k) * ( c2 * rew(i,j,k) & + c4 * (rew(i-1,j,k) + rew(i+1,j,k)) ) enddo enddo do j = Hgrid%jlb+1, je do i = is, ie frns(i,j,k) = fns(i,j,k) * ( c2 * rns(i,j,k) & + c4 * (rns(i,j-1,k) + rns(i,j+1,k)) ) enddo enddo enddo endif ! assumed halosize = 1 call update_halo (Hgrid,TEMP,frew,halos=WEST) call update_halo (Hgrid,VWND,frns,halos=SOUTH) ! ---- horizontal advective tendency ---- do k = 1, size(r,3) do j = js, je do i = is, ie rdt(i,j,k) = -0.5*Hgrid%Tmp%rarea(j)* & ((frew(i ,j,k)+frns(i,j ,k)) & -(frew(i-1,j,k)+frns(i,j-1,k))) enddo enddo enddo end subroutine advect_mass_horiz_4th !####################################################################### subroutine advect_vel_horiz_2nd (Hgrid, few, fns, u, v, udt, vdt) type(horiz_grid_type), intent(in) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, fns, & u, v real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: udt, vdt real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) :: & fuew, funs, fvew, fvns integer :: i, j, k ! second order centered differencing on B-grid momentum grid ! constant grid spacing assumed ! minimum halo size = 1 do k = 1, size(u,3) !--- horizontal fluxes --- do j = Hgrid%Vel%js, Hgrid%Vel%je do i = Hgrid%Vel%is, Hgrid%Vel%ie+1 fuew(i,j) = few(i,j,k) * (u(i-1,j,k) + u(i,j,k)) fvew(i,j) = few(i,j,k) * (v(i-1,j,k) + v(i,j,k)) enddo enddo do j = Hgrid%Vel%js, Hgrid%Vel%je+1 do i = Hgrid%Vel%is, Hgrid%Vel%ie funs(i,j) = fns(i,j,k) * (u(i,j-1,k) + u(i,j,k)) fvns(i,j) = fns(i,j,k) * (v(i,j-1,k) + v(i,j,k)) enddo enddo !--- horizontal advective tendency ---- do j = Hgrid%Vel%js, Hgrid%Vel%je do i = Hgrid%Vel%is, Hgrid%Vel%ie udt(i,j,k) = -0.5*Hgrid%Vel%rarea(j)* & ((fuew(i+1,j)+funs(i,j+1)) & -(fuew(i ,j)+funs(i,j ))) vdt(i,j,k) = -0.5*Hgrid%Vel%rarea(j)* & ((fvew(i+1,j)+fvns(i,j+1)) & -(fvew(i ,j)+fvns(i,j ))) enddo enddo enddo end subroutine advect_vel_horiz_2nd !####################################################################### subroutine advect_vel_horiz_4th (Hgrid, few, fns, u, v, udt, vdt, & mask4) type(horiz_grid_type), intent(inout) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, fns, u, v real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: udt, vdt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:), optional :: & mask4 ! compress u & v components into 4th dimension ! this will create larger and more efficient halo updates real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(u,3),2) :: uvew, uvns, fuvew, fuvns integer, parameter :: UCOMP=1, VCOMP=2 real, dimension(Hgrid%ilb:Hgrid%iub) :: x2, x4 real, dimension(Hgrid%ilb:Hgrid%iub) :: y2, y4 real :: z2, z4 integer :: i, j, k, is, ie, js, je ! fourth order centered differencing on B-grid momentum grid ! constant grid spacing assumed ! assumed halo size = 1 is = Hgrid % Vel % is; ie = Hgrid % Vel % ie js = Hgrid % Vel % js; je = Hgrid % Vel % je do k = 1, size(u,3) do j = Hgrid%jlb ,Hgrid%jub do i = Hgrid%ilb+1,Hgrid%iub uvew(i,j,k,UCOMP) = u(i-1,j,k)+u(i,j,k) uvew(i,j,k,VCOMP) = v(i-1,j,k)+v(i,j,k) enddo enddo do j = Hgrid%jlb+1,Hgrid%jub do i = Hgrid%ilb ,Hgrid%iub uvns(i,j,k,UCOMP) = u(i,j-1,k)+u(i,j,k) uvns(i,j,k,VCOMP) = v(i,j-1,k)+v(i,j,k) enddo enddo enddo ! assumed halosize = 1 call update_halo (Hgrid,WIND,uvew,halos=WEST) call update_halo (Hgrid,TEMP,uvns,halos=SOUTH) !---- horizontal fluxes ---- if (present(mask4)) then do k = 1, size(u,3) do j = js, je do i = is, Hgrid%iub-1 x4(i) = mask4(i,j,k,1) x2(i) = 1.0 - 2.*x4(i) fuvew(i,j,k,UCOMP) = few(i,j,k) * ( x2(i) * uvew(i,j,k,UCOMP) & + x4(i) * (uvew(i-1,j,k,UCOMP) + uvew(i+1,j,k,UCOMP)) ) fuvew(i,j,k,VCOMP) = few(i,j,k) * ( x2(i) * uvew(i,j,k,VCOMP) & + x4(i) * (uvew(i-1,j,k,VCOMP) + uvew(i+1,j,k,VCOMP)) ) enddo enddo do j = js, Hgrid%jub-1 do i = is, ie y4(i) = mask4(i,j,k,2) y2(i) = 1.0 - 2.*y4(i) fuvns(i,j,k,UCOMP) = fns(i,j,k) * ( y2(i) * uvns(i,j,k,UCOMP) & + y4(i) * (uvns(i,j-1,k,UCOMP) + uvns(i,j+1,k,UCOMP)) ) fuvns(i,j,k,VCOMP) = fns(i,j,k) * ( y2(i) * uvns(i,j,k,VCOMP) & + y4(i) * (uvns(i,j-1,k,VCOMP) + uvns(i,j+1,k,VCOMP)) ) enddo enddo enddo else do k = 1, size(u,3) do j = js, je do i = is, Hgrid%iub-1 fuvew(i,j,k,UCOMP) = few(i,j,k) * ( c2 * uvew(i,j,k,UCOMP) & + c4 * (uvew(i-1,j,k,UCOMP) + uvew(i+1,j,k,UCOMP)) ) fuvew(i,j,k,VCOMP) = few(i,j,k) * ( c2 * uvew(i,j,k,VCOMP) & + c4 * (uvew(i-1,j,k,VCOMP) + uvew(i+1,j,k,VCOMP)) ) enddo enddo do j = js, Hgrid%jub-1 ! second order near poles z4 = c4 if (j <= Hgrid%Vel%jsg+1) z4 = 0. if (j >= Hgrid%Vel%jeg ) z4 = 0. z2 = 1.-2.*z4 do i = is, ie fuvns(i,j,k,UCOMP) = fns(i,j,k) * ( z2 * uvns(i,j,k,UCOMP) & + z4 * (uvns(i,j-1,k,UCOMP) + uvns(i,j+1,k,UCOMP)) ) fuvns(i,j,k,VCOMP) = fns(i,j,k) * ( z2 * uvns(i,j,k,VCOMP) & + z4 * (uvns(i,j-1,k,VCOMP) + uvns(i,j+1,k,VCOMP)) ) enddo enddo enddo endif ! assumed halosize = 1 call update_halo (Hgrid,WIND,fuvew,halos=EAST) call update_halo (Hgrid,TEMP,fuvns,halos=NORTH) !---- horizontal advective tendency ---- do k = 1, size(u,3) do j = js, je do i = is, ie udt(i,j,k) = -0.5*Hgrid%Vel%rarea(j)* & ((fuvew(i+1,j,k,UCOMP)+fuvns(i,j+1,k,UCOMP)) & -(fuvew(i ,j,k,UCOMP)+fuvns(i,j ,k,UCOMP))) vdt(i,j,k) = -0.5*Hgrid%Vel%rarea(j)* & ((fuvew(i+1,j,k,VCOMP)+fuvns(i,j+1,k,VCOMP)) & -(fuvew(i ,j,k,VCOMP)+fuvns(i,j ,k,VCOMP))) enddo enddo enddo end subroutine advect_vel_horiz_4th !####################################################################### ! computes advecting wind for horizontal finite-volume advection subroutine compute_advecting_wind ( Hgrid, dpo, dpn, few, fns, uc, vc ) type(horiz_grid_type), intent(inout) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dpo, dpn, few, fns real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: uc, vc real :: dp integer :: i, j, k do k = 1, size(uc,3) do j = Hgrid%jlb, Hgrid%jub do i = Hgrid%ilb, Hgrid%iub-1 dp = (dpo(i,j,k)+dpn(i,j,k)+dpo(i+1,j,k)+dpn(i+1,j,k))*0.25 uc(i,j,k) = few(i,j,k)/(dp*Hgrid%Tmp%dy) enddo enddo enddo do k = 1, size(uc,3) do j = Hgrid%jlb, Hgrid%jub-1 do i = Hgrid%ilb, Hgrid%iub dp = (dpo(i,j,k)+dpn(i,j,k)+dpo(i,j+1,k)+dpn(i,j+1,k))*0.25 vc(i,j,k) = fns(i,j,k)/(dp*Hgrid%Vel%dx(j)) enddo enddo enddo call update_halo ( Hgrid, TEMP, uc, halos=EAST ) call update_halo ( Hgrid, VWND, vc, halos=NORTH ) end subroutine compute_advecting_wind !####################################################################### !--- stability diagnostic (CFL for horizontal finite volume scheme) ---- subroutine stability_diag ( Hgrid, dt, few, dpo, dpn ) type(horiz_grid_type), intent(inout) :: Hgrid real, intent(in) :: dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, dpo, dpn real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: uc real :: cflmax, cfl, dpa integer :: i, j, k cflmax = 0. do k = 1, size(few,3) ! compute zonal wind do j = Hgrid%Tmp%js , Hgrid%Tmp%je do i = Hgrid%Tmp%is-1, Hgrid%Tmp%ie dpa = (dpo(i,j,k)+dpn(i,j,k)+dpo(i+1,j,k)+dpn(i+1,j,k))*0.25 uc(i,j) = few(i,j,k)/(dpa*Hgrid%Tmp%dy) enddo enddo do j = Hgrid%Tmp%js, Hgrid%Tmp%je do i = Hgrid%Tmp%is, Hgrid%Tmp%ie cfl = abs(uc(i,j)-uc(i-1,j))*dt/Hgrid%Tmp%dx(j) cflmax = max(cflmax,cfl) enddo enddo enddo call mpp_max (cflmax) if (mpp_pe() == mpp_root_pe()) then if (cflmax > 1.0) then print *, 'x-axis stability violated, cfl = ', cflmax endif endif end subroutine stability_diag !####################################################################### subroutine mask4_init (Hgrid, grid, sigma, mask, mask4) type(horiz_grid_type), intent(inout) :: Hgrid integer, intent(in) :: grid logical, intent(in) :: sigma real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: mask real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: mask4 real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub, & size(mask,3)) :: mew, mns integer :: i,j,k,is,ie,js,je,n ! horizontal masks select case (grid) case (1) if (sigma) then mask4(:,:,:,1:2) = c4 else ! initialize mask4 = 0.0 mns(:,Hgrid%jub,:) = 0.0 do k = 1, size(mask,3) do j = Hgrid%jlb,Hgrid%jub do i = Hgrid%ilb,Hgrid%iub-1 mew(i,j,k) = mask(i+1,j,k)*mask(i,j,k) enddo enddo do j = Hgrid%jlb,Hgrid%jub-1 do i = Hgrid%ilb,Hgrid%iub mns(i,j,k) = mask(i,j+1,k)*mask(i,j,k) enddo enddo enddo ! assumed halosize = 1 call update_halo (Hgrid,TEMP,mew,halos=EAST) call update_halo (Hgrid,UWND,mns,halos=NORTH,flags=NOPOLE) do k = 1, size(mask,3) ! x-axis do j = Hgrid%Tmp%js, Hgrid%Tmp%je do i = Hgrid%ilb+1, Hgrid%Tmp%ie !mask4(i,j,k,1) = mew(i-1,j,k)<.01 .or. mew(i,j,k)<.01 .or. mew(i+1,j,k)<.01 !mask4(i,j,k,1) = mew(i-1,j,k) * mew(i,j,k) * mew(i+1,j,k) mask4(i,j,k,1) = c4 * mew(i-1,j,k) * mew(i+1,j,k) enddo enddo ! y-axis do j = Hgrid%jlb+1, Hgrid%Tmp%je do i = Hgrid%Tmp%is, Hgrid%Tmp%ie !mask4(i,j,k,2) = mns(i,j-1,k)<.01 .or. mns(i,j,k)<.01 .or. mns(i,j-1,k)<.01 !mask4(i,j,k,2) = mns(i,j-1,k) * mns(i,j,k) * mns(i,j+1,k) mask4(i,j,k,2) = c4 * mns(i,j-1,k) * mns(i,j+1,k) enddo enddo enddo endif case (2) if (sigma) then mask4(:,:,:,1:2) = c4 mns = 1.0 else ! initialize mask4 = 0.0 do k = 1, size(mask,3) do j = Hgrid%jlb ,Hgrid%jub do i = Hgrid%ilb+1,Hgrid%iub mew(i,j,k) = mask(i-1,j,k)*mask(i,j,k) enddo enddo do j = Hgrid%jlb+1,Hgrid%jub do i = Hgrid%ilb ,Hgrid%iub mns(i,j,k) = mask(i,j-1,k)*mask(i,j,k) enddo enddo enddo ! assumed halosize = 1 call update_halo (Hgrid,UWND,mew,halos=WEST) call update_halo (Hgrid,TEMP,mns,halos=SOUTH,flags=NOPOLE) endif ! use second order in meridional direction near poles call vel_flux_boundary (Hgrid,mns) if (.not.sigma) then do k = 1, size(mask,3) do j = Hgrid%Vel%js, Hgrid%Vel%je do i = Hgrid%Vel%is, Hgrid%iub-1 mask4(i,j,k,1) = c4 * mew(i-1,j,k) * mew(i+1,j,k) enddo enddo enddo endif do k = 1, size(mask,3) do j = Hgrid%Vel%js, Hgrid%jub-1 do i = Hgrid%Vel%is, Hgrid%Vel%ie mask4(i,j,k,2) = c4 * mns(i,j-1,k) * mns(i,j+1,k) enddo enddo enddo end select end subroutine mask4_init !####################################################################### ! initialization routine subroutine advection_init ( Hgrid ) ! INPUT: Hgrid = horizontal grid constants ! RESULT: reads bgrid_advection_nml namelist ! defines advection control parameters (stored in advec_control_type) type(horiz_grid_type), intent(in) :: Hgrid real :: wt integer :: n, m, np, ntrace, unit, ierr, io, logunit character(len=32) :: advec_methods(2) = (/ 'advec_horiz', 'advec_vert ' /) character(len=128) :: scheme, params, name !----------------------------------------------------------------------- ! read namelist if (file_exist('input.nml')) then #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=bgrid_advection_nml, iostat=io) ierr = check_nml_error(io,'bgrid_advection_nml') #else unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=bgrid_advection_nml, iostat=io, end=5) ierr = check_nml_error (io, 'bgrid_advection_nml') enddo 5 call close_file (unit) #endif endif logunit = stdlog() ! write version, namelist info to log file if (do_log) then call write_version_number(version, tagname) if (mpp_pe() == mpp_root_pe()) write (logunit, nml=bgrid_advection_nml) do_log = .false. endif ! determine the number of prognostic tracers call get_number_tracers ( MODEL_ATMOS, num_prog=ntrace ) ! allocate space for the control parameters allocate ( Control%scheme (2,-1:ntrace), & Control%weight (-1:ntrace) ) allocate ( Control%Fill%fill_scheme (ntrace), & Control%Fill%npass_horiz (ntrace), & Control%Fill%npass_vert (ntrace) ) ! set namelist values of control parameters Control%scheme(1,-1) = set_advec_scheme( horiz_advec_scheme_wind ) Control%scheme(2,-1) = set_advec_scheme( vert_advec_scheme_wind ) Control%weight (-1) = advec_weight_wind Control%scheme(1, 0) = set_advec_scheme( horiz_advec_scheme_temp ) Control%scheme(2, 0) = set_advec_scheme( vert_advec_scheme_temp ) Control%weight ( 0) = advec_weight_temp do n = 1, ntrace Control%scheme(1, n) = set_advec_scheme( horiz_advec_scheme_tracer ) Control%scheme(2, n) = set_advec_scheme( vert_advec_scheme_tracer ) Control%weight ( n) = advec_weight_tracer ! parameters for tracer filling/borrowing Control%Fill%npass_horiz(n) = num_fill_pass Control%Fill%npass_vert (n) = num_fill_pass Control%Fill%fill_scheme(n) = NONE enddo ! currently only one possible filling scheme if (maxval(Control%Fill%npass_horiz) > 0 .or. maxval(Control%Fill%npass_vert) > 0) & Control%Fill%fill_scheme = EQUAL_BORROW ! process tracer table information for advection methods do n = 1, ntrace do m = 1, 2 if (query_method(trim(advec_methods(m)), MODEL_ATMOS, n, scheme, params)) then Control%scheme(m,n) = set_advec_scheme( scheme ) ! parse e-b weight if (Control%scheme(m,n) == SECOND_CENTERED .or. & Control%scheme(m,n) == FOURTH_CENTERED) then if (parse(params,'wt', wt) == 1) Control%weight(n) = wt endif endif enddo enddo ! error check for all Euler weights do n = -1, ntrace if (Control%weight(n) < 0.0 .or. Control%weight(n) > 1.0) & call error_mesg ('bgrid_advection_mod', & 'E-B weight out of range [0,1]', FATAL) enddo ! process tracer table information for filling method do n = 1, ntrace if (query_method('filling', MODEL_ATMOS, n, scheme, params)) then if (uppercase(trim(scheme)) == 'LOCAL') then Control%Fill%fill_scheme(n) = EQUAL_BORROW Control%Fill%npass_horiz(n) = 1 Control%Fill%npass_vert (n) = 1 !else if (uppercase(trim(scheme)) == 'GLOBAL') then ! Control%Fill%fill_scheme(n) = GLOBAL_BORROW else if (uppercase(trim(scheme)) == 'NONE') then Control%Fill%fill_scheme(n) = NONE Control%Fill%npass_horiz(n) = 0 Control%Fill%npass_vert (n) = 0 else call error_mesg ('bgrid_advection_mod', & 'invalid filling scheme, '//uppercase(trim(scheme)), FATAL) endif ! parse number of filling passes if (Control%Fill%fill_scheme(n) == EQUAL_BORROW) then if (parse(params,'hp', np) == 1) Control%Fill%npass_horiz(n) = np if (parse(params,'vp', np) == 1) Control%Fill%npass_vert (n) = np endif endif enddo ! print results if (mpp_pe() == mpp_root_pe()) then write (logunit,10) 'momentum', & (trim(echo_advec_scheme(Control%scheme(m,-1))),m=1,2), Control%weight(-1) write (logunit,10) 'temperature', & (trim(echo_advec_scheme(Control%scheme(m, 0))),m=1,2), Control%weight( 0) do n = 1, ntrace call get_tracer_names (MODEL_ATMOS, n, name) write (logunit,11) n, trim(name), & ( trim(echo_advec_scheme(Control%scheme(m,n))),m=1,2), & Control%weight(n), & Control%Fill%npass_horiz(n), & Control%Fill%npass_vert(n) 10 format (3x, a24, ', HORIZ=',a24, ', VERT=',a24, ', wt=',f7.3) 11 format (i3, a24, ', HORIZ=',a24, ', VERT=',a24, ', wt=',f7.3, & ', hp=',i2, ', vp=',i2) enddo endif ! check if fourth order centered horizontal scheme Control%do_mask4_vel = Control%scheme(1,-1) == FOURTH_CENTERED Control%do_mask4_tmp = .false. do n = 0, ntrace if (Control%scheme(1,n) == FOURTH_CENTERED) then Control%do_mask4_tmp = .true. exit endif enddo ! check if finite_volume horizontal scheme for temperature if (Control%scheme(1,0) == FINITE_VOLUME_LINEAR .or. & Control%scheme(1,0) == FINITE_VOLUME_PARABOLIC) then Control%do_finite_volume_tmp = .true. else Control%do_finite_volume_tmp = .false. endif ! check if finite_volume horizontal scheme for tracers Control%do_finite_volume_trs = .false. do n = 1, ntrace if (Control%scheme(1,n) == FINITE_VOLUME_LINEAR .or. & Control%scheme(1,n) == FINITE_VOLUME_PARABOLIC) then Control%do_finite_volume_trs = .true. exit endif enddo ! error checks ! many horizontal schemes are not support with this release do n = -1, ntrace if (Control%scheme(1,n) == FINITE_VOLUME_LINEAR .or. & Control%scheme(1,n) == FINITE_VOLUME_PARABOLIC .or. & Control%scheme(1,n) == SECOND_CENTERED_WTS .or. & Control%scheme(1,n) == FOURTH_CENTERED_WTS ) then call error_mesg ('bgrid_advection_mod', & 'advection scheme not supported', FATAL) endif enddo ! initialize code sections for performance timing if (do_clock_init) then id_tmp = mpp_clock_id ('BGRID: temperature advection', & flags=MPP_CLOCK_SYNC, grain=CLOCK_ROUTINE) id_trs = mpp_clock_id ('BGRID: tracer advection', & flags=MPP_CLOCK_SYNC, grain=CLOCK_ROUTINE) id_vel = mpp_clock_id ('BGRID: momentum advection', & flags=MPP_CLOCK_SYNC, grain=CLOCK_ROUTINE) id_fill = mpp_clock_id ('BGRID: tracer borrowing', & flags=MPP_CLOCK_SYNC, grain=CLOCK_ROUTINE) do_clock_init = .false. endif !----------------------------------------------------------------------- end subroutine advection_init !####################################################################### subroutine advection_end ! terminate vertical advection module call vert_advection_end ! terminate finite-volume advection module ! if (Control%do_finite_volume_tmp .or. Control%do_finite_volume_trs) then ! call horiz_advection_end ! endif ! free up memory used (although not much) deallocate ( Control%scheme, Control%weight, & Control%Fill%fill_scheme, & Control%Fill%npass_horiz, & Control%Fill%npass_vert ) end subroutine advection_end !####################################################################### ! converts character string names to integer parameters ! checks for all valid names, otherwise produces an error function set_advec_scheme ( scheme ) result ( advec_scheme ) character(len=*), intent(in) :: scheme integer :: advec_scheme if (uppercase(trim(scheme)) == 'SECOND_CENTERED') then advec_scheme = SECOND_CENTERED else if (uppercase(trim(scheme)) == 'FOURTH_CENTERED') then advec_scheme = FOURTH_CENTERED else if (uppercase(trim(scheme)) == 'FOURTH_CENTERED_WTS') then advec_scheme = FOURTH_CENTERED_WTS else if (uppercase(trim(scheme)) == 'SECOND_CENTERED_WTS') then advec_scheme = SECOND_CENTERED_WTS else if (uppercase(trim(scheme)) == 'FINITE_VOLUME_LINEAR') then advec_scheme = FINITE_VOLUME_LINEAR else if (uppercase(trim(scheme)) == 'FINITE_VOLUME_PARABOLIC') then advec_scheme = FINITE_VOLUME_PARABOLIC else if (uppercase(trim(scheme)) == 'NONE') then advec_scheme = NONE else call error_mesg ('bgrid_advection_mod', & 'invalid advection scheme, '//uppercase(trim(scheme)), FATAL) endif end function set_advec_scheme !----------------------------------------------- ! converts integer parameter values to character string names ! checks for all valid schemes, otherwise produces an error function echo_advec_scheme ( scheme_number ) result ( scheme_name ) integer, intent(in) :: scheme_number character(len=24) :: scheme_name select case (scheme_number) case(NONE) scheme_name = 'none ' case(SECOND_CENTERED) scheme_name = 'second_centered ' case(FOURTH_CENTERED) scheme_name = 'fourth_centered ' case(SECOND_CENTERED_WTS) scheme_name = 'second_centered_wts ' case(FOURTH_CENTERED_WTS) scheme_name = 'fourth_centered_wts ' case(FINITE_VOLUME_LINEAR) scheme_name = 'finite_volume_linear ' case(FINITE_VOLUME_PARABOLIC) scheme_name = 'finite_volume_parabolic ' case default call error_mesg ('bgrid_advection_mod', & 'invalid advection scheme number', FATAL) end select end function echo_advec_scheme !####################################################################### ! routines to minimize negative tracer values ! conservative horizontal and vertical borrowing !####################################################################### subroutine horiz_borrow (Hgrid, dt, dpde, var, var_dt, mask, iters) !----------------------------------------------------------------------- ! ! removes negative values by borrowing from horizontal ! neighbors. (usually for specific humidity or tke ) ! !----------------------------------------------------------------------- type(horiz_grid_type), intent(inout) :: Hgrid real, intent(in) :: dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dpde real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: var real, intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: var_dt real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:), optional :: mask integer, intent(in), dimension(:), optional :: iters !----------------------------------------------------------------------- ! ( note: 0.0 < flt <= 0.125 ) real, parameter :: flt = 0.125 real, parameter :: flt4 = flt*0.25 real, parameter :: rmin = 0.0 !----------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,size(var,3)) :: & hew, hns, radp, few, fns, rcur, rdif real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: rew, rns real, dimension(Hgrid%jlb:Hgrid%jub) :: areax, areay integer :: i, j, k, n, is, ie, js, je, nlev, ntrace, knt, mxknt(size(var,4)) !----------------------------------------------------------------------- mxknt = 1; if (present(iters)) mxknt = iters if (maxval(mxknt) == 0) return nlev = size(var,3) ntrace = size(var,4) is = Hgrid%Tmp%is; ie = Hgrid%Tmp%ie js = Hgrid%Tmp%js; je = Hgrid%Tmp%je !----------------------------------------------------------------------- !------ compute flux coeff common to all variables ------- do j = js-1, je areax(j) = Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j) areay(j) = Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j+1) enddo do k = 1, nlev do j = js-1, je do i = is-1, ie hew(i,j,k) = areax(j)*(dpde(i,j,k)+dpde(i+1,j,k)) hns(i,j,k) = areay(j)*(dpde(i,j,k)+dpde(i,j+1,k)) enddo enddo enddo ! apply step-mountain mask? if (present(mask)) then do k = 1, nlev do j = js-1, je do i = is-1, ie hew(i,j,k) = hew(i,j,k)*mask(i,j,k)*mask(i+1,j,k) hns(i,j,k) = hns(i,j,k)*mask(i,j,k)*mask(i,j+1,k) enddo enddo enddo endif do k = 1, nlev do j = js, je radp(:,j,k) = flt4/(Hgrid%Tmp%area(j)*dpde(:,j,k)) enddo enddo !----------------------------------------------------------------------- !---- tracer loop ----- !---- store current value of tracer in rcur ---- do n = 1, ntrace rcur(:,:,:) = var(:,:,:,n) + var_dt(:,:,:,n)*dt !---- iteration loop ----- few = 0.0; fns = 0.0 do knt = 1, mxknt(n) !--------------2-nd order lat/lon contributions------------------------- ! --- do borrowing where adjacent values have opposite sign --- ! but do not turn off fluxes previously turned on do k = 1, nlev do j = js-1, je do i = is-1, ie ! east-west fluxes if ((rcur(i+1,j,k) < rmin .and. rcur(i,j,k) >= rmin) .or. & (rcur(i+1,j,k) >= rmin .and. rcur(i,j,k) < rmin)) & few(i,j,k) = hew(i,j,k) rew(i,j) = (rcur(i+1,j,k)-rcur(i,j,k))*few(i,j,k) ! north-south fluxes if ((rcur(i,j+1,k) < rmin .and. rcur(i,j,k) >= rmin) .or. & (rcur(i,j+1,k) >= rmin .and. rcur(i,j,k) < rmin)) & fns(i,j,k) = hns(i,j,k) rns(i,j) = (rcur(i,j+1,k)-rcur(i,j,k))*fns(i,j,k) enddo enddo do j = js, je do i = is, ie rdif(i,j,k) = (rew(i,j)-rew(i-1,j)+ & rns(i,j)-rns(i,j-1))*radp(i,j,k) enddo enddo enddo !---- halo update (assumed that halo size = 1) ---- call update_halo (Hgrid, TEMP, rdif) !---- update current value of tracer ---- rcur(:,:,:) = rcur(:,:,:) + rdif(:,:,:) !----------------------------------------------------------------------- enddo ! iteration loop !----------------------------------------------------------------------- !---- return the tendency ----- var_dt(:,:,:,n) = (rcur(:,:,:) - var(:,:,:,n)) / dt !----------------------------------------------------------------------- enddo ! tracer loop !----------------------------------------------------------------------- end subroutine horiz_borrow !####################################################################### subroutine vert_borrow (dt, dpde, var, var_dt, iters) !----------------------------------------------------------------------- ! ! This removes negative specific humidity/mixing ratios by borrowing ! from the grid boxes immediately above and below. If not enough ! is available to fill the negative then a negative will remain. ! !----------------------------------------------------------------------- real, intent(in) :: dt, dpde(:,:,:), var(:,:,:,:) real, intent(inout) :: var_dt(:,:,:,:) integer, intent(in), optional :: iters(:) !----------------------------------------------------------------------- real, dimension(size(var,3)) :: deficit, surplus, rdpdt real :: divid, ratio_dn, ratio_up, var_dp integer :: i, j, k, n, m, nlev, num_iters(size(var,4)), num, num_var !----------------------------------------------------------------------- num_iters = 4; if (present(iters)) num_iters = iters if (maxval(num_iters) == 0) return num_var = size(var,4) if (num_var == 0) return nlev = size(var,3) !---- variable loop and iteration loop ----- do n = 1, num_var do m = 1, num_iters(n) !---- existing negatives will not be corrected ---- do j = 1, size(var,2) do i = 1, size(var,1) do k = 1, nlev var_dp = (var(i,j,k,n)+dt*var_dt(i,j,k,n))*dpde(i,j,k) deficit(k) = min(var_dp,0.0) surplus(k) = max(var_dp,0.0) rdpdt(k) = 1./(dpde(i,j,k)*dt) enddo !---- top level ---- if (deficit(1) < 0.0) then divid = max(-surplus(2)/deficit(1),1.0) ratio_dn = surplus(2)/divid var_dt(i,j,1,n) = var_dt(i,j,1,n) + ratio_dn*rdpdt(1) var_dt(i,j,2,n) = var_dt(i,j,2,n) - ratio_dn*rdpdt(2) surplus(2) = max(surplus(2)-ratio_dn,0.) endif !---- interior levels ---- do k = 2, nlev-1 if (deficit(k) < 0.0) then divid = max(-(surplus(k-1)+surplus(k+1))/deficit(k),1.0) ratio_up = surplus(k-1)/divid ratio_dn = surplus(k+1)/divid var_dt(i,j,k,n) = var_dt(i,j,k,n)+(ratio_up+ratio_dn)*rdpdt(k) var_dt(i,j,k-1,n) = var_dt(i,j,k-1,n)-ratio_up*rdpdt(k-1) var_dt(i,j,k+1,n) = var_dt(i,j,k+1,n)-ratio_dn*rdpdt(k+1) surplus(k+1) = max(surplus(k+1)-ratio_dn,0.) endif enddo !---- bottom level ---- if (deficit(nlev) < 0.0) then divid = max(-surplus(nlev-1)/deficit(nlev),1.0) ratio_up = surplus(nlev-1)/divid var_dt(i,j,nlev ,n) = var_dt(i,j,nlev ,n)+ratio_up*rdpdt(nlev) var_dt(i,j,nlev-1,n) = var_dt(i,j,nlev-1,n)-ratio_up*rdpdt(nlev-1) endif enddo enddo enddo enddo !----------------------------------------------------------------------- end subroutine vert_borrow !####################################################################### subroutine compute_etadot_vel ( Hgrid, Vgrid, Masks, few, fns, etadot ) !subroutine compute_etadot_vel ( Hgrid, Vgrid, Masks, res, few, fns, etadot ) type(horiz_grid_type), intent(in) :: Hgrid type (vert_grid_type), intent(in) :: Vgrid type (grid_mask_type), intent(in) :: Masks !real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:) :: res real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: few, fns real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: etadot real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: sdiv integer :: i, j, k, is, ie, js, je, nlev ! this code will not work with the eta coordinate ! the variable res need to be added as an argument if (.not.Masks%sigma) call error_mesg ('bgrid_advection_mod', & 'etadot cannot be recomputed for the eta coordinate', FATAL) is = Hgrid%Vel%is; ie = Hgrid%Vel%ie js = Hgrid%Vel%js; je = Hgrid%Vel%je nlev = size(few,3) etadot = 0. sdiv = 0. do k = 1, nlev ! vertical integral of divergence at velocity points ! only need etadot for velocity compute domain do j = js, je do i = is, ie sdiv(i,j) = sdiv(i,j) + ((few(i+1,j,k)+fns(i,j+1,k))- & (few(i ,j,k)+fns(i,j ,k))) & *Hgrid%Vel%rarea(j) enddo enddo if (k < nlev) etadot(:,:,k+1) = -sdiv(:,:) enddo ! finally include d(ps)/dt term do k = 2, nlev !etadot(:,:,k) = etadot(:,:,k) + sdiv(:,:)*res(:,:)*Vgrid%eta(k)) etadot(:,:,k) = etadot(:,:,k) + sdiv(:,:)*Vgrid%eta(k) if (.not.Masks%sigma) etadot(:,:,k) = etadot(:,:,k)*Masks%Vel%mask(:,:,k) enddo end subroutine compute_etadot_vel !####################################################################### end module bgrid_advection_mod