!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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: 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 12.0 2005/04/14 15:33:49 fms Exp $' character(len=128) :: tagname='$Name: mom4p1_pubrel_dec2009_nnz $' !----------------------------------------------------------------------- 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 character(len=32) :: advec_methods(2) = (/ 'advec_horiz', 'advec_vert ' /) character(len=128) :: scheme, params, name !----------------------------------------------------------------------- ! read namelist if (file_exist('input.nml')) then 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 ! write version, namelist info to log file if (do_log) then call write_version_number (version,tagname) if (mpp_pe() == mpp_root_pe()) write (stdlog(), 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 (stdlog(),10) 'momentum', & (trim(echo_advec_scheme(Control%scheme(m,-1))),m=1,2), Control%weight(-1) write (stdlog(),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 (stdlog(),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