!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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_horiz_diff_mod !======================================================================= ! ! linear horizontal mixing ! with option for any order accuracy ! !======================================================================= use bgrid_horiz_mod , only: horiz_grid_type use bgrid_masks_mod , only: grid_mask_type use bgrid_prog_var_mod , only: prog_var_type, var_init use bgrid_halo_mod , only: update_halo, vel_flux_boundary, & EAST, NORTH, NOPOLE, POLEONLY, & TEMP, UWND, VWND use bgrid_change_grid_mod, only: change_grid, TEMP_GRID, WIND_GRID use fms_mod, only: error_mesg, FATAL, write_version_number, & open_namelist_file, check_nml_error, close_file, & mpp_clock_id, mpp_clock_begin, mpp_clock_end, & MPP_CLOCK_SYNC, CLOCK_MODULE, file_exist, & mpp_pe, mpp_root_pe, uppercase, stdlog use constants_mod, only: RADIUS 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 interfaces ---------- public horiz_diff, horiz_diff_init !----------------------------------------------------------------------- ! namelist bgrid_horiz_diff_nml ! damp_scheme_wind Determines how horizontal damping coefficients ! damp_scheme_temp vary with latitude. Possible values are: 0,1,2,3. ! [0 = no damping; 1 = coeff are uniform with latitude; ! 2 = coeffs vary as 1/(dx**2+dy**2); ! 3 = coeffs vary as 1/(dx**2); ! 4 = same as scheme 1 but with scheme 2 applied ! poleward of reflat] ! Notes: Schemes 2-4 provide increased damping in ! higher latitudes. Temperature and tracers are ! controlled using "damp_scheme_temp". integer :: damp_scheme_temp = 1 integer :: damp_scheme_wind = 1 ! reflat Latitude cutoff (in degrees) at which increased ! high latitude damping is applied. Equatorward of ! this latitude uniform damping (scheme=1) is applied; ! poleward of this latitude enhanced damping (scheme=2) ! is applied. This variable is only used when ! damp_scheme_wind=4 or damp_scheme_temp=4. real :: reflat = 85. ! damp_order_wind The horizontal damping order for momentum, ! damp_order_temp temperature, and default order for all ! damp_order_tracer prognostic tracers. Only even numbers are allowed. integer :: damp_order_wind = 4 ! use: 0, 2, 4,.... integer :: damp_order_temp = 4 integer :: damp_order_tracer = 4 ! damp_coeff_wind The horizontal damping coefficients for ! damp_coeff_temp momentum, temperature, and default value for ! damp_coeff_tracer all prognostic tracers. Only positive values ! are allowed. real :: damp_coeff_wind = 0.35 real :: damp_coeff_temp = 0.35 real :: damp_coeff_tracer = 0.35 ! slope_corr_wind The topography slope correction for horizontal ! slope_corr_temp damping of momentum, temperature, and default ! slope_corr_tracer for all prognostic tracers. real, dimension(4) :: slope_corr_wind = (/0.,0.,0.,0./) real, dimension(4) :: slope_corr_temp = (/0.,0.,0.,0./) real, dimension(4) :: slope_corr_tracer = (/0.,0.,0.,0./) namelist /bgrid_horiz_diff_nml/ damp_scheme_temp, damp_scheme_wind, & damp_order_wind, damp_order_temp, damp_order_tracer, & damp_coeff_wind, damp_coeff_temp, damp_coeff_tracer, & slope_corr_wind, slope_corr_temp, slope_corr_tracer, & reflat !----------------------------------------------------------------------- !--------- private data ---------- ! control parameters for horizontal damping ! set via the namelist or field table type hdiff_control_type integer, pointer :: order(:) =>NULL() real , pointer :: coeff(:) =>NULL() real , pointer :: slope(:,:) =>NULL() logical, pointer :: do_slope_adj(:) =>NULL() logical :: do_damping, do_slope_adj_temp integer :: damping_scheme_wind, damping_scheme_temp real, dimension(:), pointer :: areahx =>NULL(), & areahy =>NULL(), & areavx =>NULL(), & areavy => NULL() real, dimension(:,:), pointer :: wth => NULL(), & wtv =>NULL() end type hdiff_control_type type(hdiff_control_type), save :: Control integer :: nlev ! number of model levels character(len=128) :: version='$Id: bgrid_horiz_diff.f90,v 10.0 2003/10/27 23:31:04 arl Exp $' character(len=128) :: tagname='$Name: $' logical :: do_log = .true. ! timing data integer :: id_total logical :: do_clock_init = .true. !----------------------------------------------------------------------- contains !####################################################################### subroutine horiz_diff ( Hgrid, Masks, nplev, dt, dpde, pres, Var, Var_dt ) !----------------------------------------------------------------------- ! ! Hgrid = horizontal grid constants ! Masks = grid masking constants for eta coordinate ! nplev = number of "pure" pressure levels at the top of the model ! dt = adjustment time step ! dpde = pressure weight for model layers ! pres = pressure at full model levels ! Var = prognostic variables at the last updated time level ! Var_dt = tendency of prognostic variables since the last ! updated time level ! !----------------------------------------------------------------------- type(horiz_grid_type), intent(inout) :: Hgrid type (grid_mask_type), intent(in) :: Masks integer, intent(in) :: nplev real, intent(in) :: dt real, intent(in) :: dpde(Hgrid%ilb:,Hgrid%jlb:,:),& pres(Hgrid%ilb:,Hgrid%jlb:,:) type (prog_var_type), intent(in) :: Var type (prog_var_type), intent(inout) :: Var_dt !----------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(dpde,3)) :: hkew3, hkns3, hkew, hkns, & hcew, hcns, dat, vdat, hdac real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, & size(dpde,3),3) :: pterms real :: dt_inv, hsign, xcoeff, ycoeff integer :: i, j, k, n, is, ie, js, je, ntp !======================================================================= ! --- should damping be done ??? ---- if ( .not. Control%do_damping ) return if ( Control%damping_scheme_wind+Control%damping_scheme_temp == 0 ) return call mpp_clock_begin (id_total) !----------------------------------------------------------------------- ! --- check the horizontal dimensions of the input array --- nlev = size(dpde,3) if (size(dpde,1) /= Hgrid%isize .or. & size(dpde,2) /= Hgrid%jsize ) call error_mesg & ('bgrid_horiz_diff', 'input array has the wrong dimensions.', FATAL) ! ---- time step related values ---- dt_inv = 1./dt ! ---- initialize flux weights ---- hkew3 = 0.0; hkns3 = 0.0 !----------------------------------------------------------------------- !-------------setup temperature and tracer damping --------------------- !----------------------------------------------------------------------- ntp = count( Control%order(1:Var_dt%ntrace) > 0 ) if (Control%order(0) > 0 .or. ntp > 0) then ! compute grid indiecs is = Hgrid%Tmp%is; ie = Hgrid%Tmp%ie js = Hgrid%Tmp%js; je = Hgrid%Tmp%je do j = js, je hdac(:,j,:) = Hgrid%Tmp%rarea(j) enddo do k = nplev+1, nlev do j = js, je hdac(:,j,k) = hdac(:,j,k) / (2.0*dpde(:,j,k)) enddo enddo ! ----- mass-weighted fluxes ----- do j = js-1, je hkew3(:,j,1:nplev) = Control%areahx(j) hkns3(:,j,1:nplev) = Control%areahy(j) enddo do k = nplev+1, nlev do j = js-1, je do i = is-1, ie hkew3(i,j,k) = Control%areahx(j) * (dpde(i,j,k)+dpde(i+1,j,k)) hkns3(i,j,k) = Control%areahy(j) * (dpde(i,j,k)+dpde(i,j+1,k)) enddo enddo enddo ! ----- mask out fluxes below step-mountain ---- if (.not.Masks%sigma) then do k = 1, nlev do j = js-1, je do i = is-1, ie hkew3(i,j,k) = hkew3(i,j,k)*Masks%Tmp%mask(i,j,k)*Masks%Tmp%mask(i+1,j,k) hkns3(i,j,k) = hkns3(i,j,k)*Masks%Tmp%mask(i,j,k)*Masks%Tmp%mask(i,j+1,k) enddo enddo enddo hdac(:,js:je,:) = hdac(:,js:je,:) * Masks%Tmp%mask(:,js:je,:) endif !----------------------------------------------------------------------- !----- slope adjustment ------------------------------------------------ if ( Control % do_slope_adj_temp ) then call slope_correction_init ( Hgrid, Masks, nplev, pres, pterms ) else hcew = 0.0 hcns = 0.0 endif !----------------------------------------------------------------------- !------------------temperature damping---------------------------------- if (Control%order(0) > 0) then ! setting up flux weights - limit to 1/8 do k = 1, nlev do j = Hgrid%jlb, Hgrid%jub xcoeff = min(0.125,Control%coeff(0)*Control%wth(j,1)) ycoeff = min(0.125,Control%coeff(0)*Control%wth(j,2)) hkew(:,j,k) = hkew3(:,j,k) * xcoeff hkns(:,j,k) = hkns3(:,j,k) * ycoeff enddo enddo dat = Var%t + dt*Var_dt%t ! --- slope adjustent --- if ( Control%do_slope_adj(0) ) & call slope_correction ( Hgrid, Masks, nplev, Control%slope(:,0), & pterms, dat, hcew, hcns ) call diff_mass (Hgrid, dat, hcew, hcns, hkew, hkns, hdac, & Control%order(0)) hsign = 1.; if (mod(Control%order(0),4) == 0) hsign = -1. Var_dt%t = Var_dt%t + hsign * dt_inv * dat endif !----------------------------------------------------------------------- !------------ tracer damping (prognostic tracers only) --------------- if (ntp > 0) then do n = 1, Var_dt%ntrace ! setting up flux weights - limit to 1/8 do k = 1, nlev do j = Hgrid%jlb, Hgrid%jub xcoeff = min(0.125,Control%coeff(n)*Control%wth(j,1)) ycoeff = min(0.125,Control%coeff(n)*Control%wth(j,2)) hkew(:,j,k) = hkew3(:,j,k) * xcoeff hkns(:,j,k) = hkns3(:,j,k) * ycoeff enddo enddo if (Control%order(n) == 0) cycle dat = Var%r(:,:,:,n) + dt*Var_dt%r(:,:,:,n) !--- slope adjustent --- if ( Control%do_slope_adj(n) ) & call slope_correction ( Hgrid, Masks, nplev, Control%slope(:,n), & pterms, dat, hcew, hcns ) call diff_mass (Hgrid, dat, hcew, hcns, hkew, hkns, & hdac, Control%order(n)) hsign = 1.; if (mod(Control%order(n),4) == 0) hsign = -1. Var_dt%r(:,:,:,n) = Var_dt%r(:,:,:,n) + hsign * dt_inv * dat enddo endif endif !----------------------------------------------------------------------- !*********************************************************************** !-----------------setup momentum damping ------------------------------- if (Control%order(-1) > 0) then is = Hgrid%Vel%is; ie = Hgrid%Vel%ie js = Hgrid%Vel%js; je = Hgrid%Vel%je call change_grid (Hgrid, TEMP_GRID, WIND_GRID, & dpde(:,:,nplev+1:nlev), vdat(:,:,nplev+1:nlev)) call update_halo (Hgrid, UWND, vdat(:,:,nplev+1:nlev), & halos=EAST+NORTH, flags=NOPOLE) do j = js, je hdac(:,j,:) = Hgrid%Vel%rarea(j) enddo do k = nplev+1, nlev hdac(:,js:je,k) = hdac(:,js:je,k) / (2.*vdat(:,js:je,k)) enddo ! ----- mass-weighted fluxes ----- do j = js, je+1 hkew3(:,j,1:nplev) = Control%areavx(j) hkns3(:,j,1:nplev) = Control%areavy(j) enddo do k = nplev+1, nlev do j = js, je+1 do i = is, ie+1 hkew3(i,j,k) = Control%areavx(j)*(vdat(i,j,k)+vdat(i-1,j,k)) hkns3(i,j,k) = Control%areavy(j)*(vdat(i,j,k)+vdat(i,j-1,k)) enddo enddo enddo ! ----- mask out fluxes below step-mountain ---- if (.not.Masks%sigma) then do k = 1, nlev do j = js, je+1 do i = is, ie+1 hkew3(i,j,k) = hkew3(i,j,k)*Masks%Vel%mask(i,j,k)*Masks%Vel%mask(i-1,j,k) hkns3(i,j,k) = hkns3(i,j,k)*Masks%Vel%mask(i,j,k)*Masks%Vel%mask(i,j-1,k) enddo enddo enddo hdac(:,js:je,:) = hdac(:,js:je,:) * Masks%Tmp%mask(:,js:je,:) endif !----------------------------------------------------------------------- !----- slope adjustment setup ------ if ( Control % do_slope_adj(-1) ) then ! ---- pressure at velocity points ---- vdat(:,:,1:nplev) = pres(:,:,1:nplev) call change_grid ( Hgrid, TEMP_GRID, WIND_GRID, & pres(:,:,nplev+1:nlev), vdat(:,:,nplev+1:nlev) ) call update_halo (Hgrid, UWND, vdat(:,:,nplev+1:nlev), & halos=EAST+NORTH, flags=NOPOLE) call vel_slope_correction_init ( Hgrid, Masks, nplev, & Control%slope(:,-1), & vdat, pterms ) endif !----------------------------------------------------------------------- ! setting up flux weights - limit to 1/8 do k = 1, nlev do j = js, je+1 xcoeff = min(0.125,Control%coeff(-1)*Control%wtv(j,1)) ycoeff = min(0.125,Control%coeff(-1)*Control%wtv(j,2)) hkew(:,j,k) = hkew3(:,j,k) * xcoeff hkns(:,j,k) = hkns3(:,j,k) * ycoeff enddo enddo ! zero-out cross polar fluxes call vel_flux_boundary (Hgrid, hkns) !----------------------------------------------------------------------- !-------------------------momentum damping------------------------------ dat = Var%u + dt*Var_dt%u vdat = Var%v + dt*Var_dt%v if ( Control % do_slope_adj(-1) ) then if (Masks%sigma) then call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac, & Control%order(-1), pterms) else call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac, & Control%order(-1), pterms, Masks%Vel%kbot) endif else call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac, & Control%order(-1)) endif hsign = 1.; if (mod(Control%order(-1),4) == 0) hsign = -1. Var_dt%u = Var_dt%u + hsign * dt_inv * dat Var_dt%v = Var_dt%v + hsign * dt_inv * vdat endif call mpp_clock_end (id_total) !----------------------------------------------------------------------- end subroutine horiz_diff !####################################################################### !####################################################################### subroutine diff_mass (Hgrid, rdat, hcew, hcns, hkew, hkns, hdac, & order) !---------------------------------------------------------------------- ! ! diff_mass is a private interface that performs multiple ! 2nd order lapacians. ! !---------------------------------------------------------------------- type(horiz_grid_type), intent(inout) :: Hgrid integer, intent(in) :: order real , intent(inout), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & nlev) :: rdat real , intent(in), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & nlev) :: hcew, hcns, & hkew, hkns, hdac !---------------------------------------------------------------------- ! hcew, hcns are weighted corrections to the diffusive fluxes for ! sloping sigma surfaces !---------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) :: & rew, rns integer :: i, j, k, n, is, ie, js, je, nordr !----------------------------------------------------------------------- is = Hgrid%Tmp%is; ie = Hgrid%Tmp%ie js = Hgrid%Tmp%js; je = Hgrid%Tmp%je !----------------------------------------------------------------------- !------------loop for order of damping scheme------------------------- do n = 1, order/2 do k = 1, nlev !------------------contributions (fluxes) ------------------------------ if ( n == 1 ) then do j = js-1, je do i = is-1, ie rew(i,j) = (rdat(i+1,j,k)-rdat(i,j,k)+hcew(i,j,k))*hkew(i,j,k) rns(i,j) = (rdat(i,j+1,k)-rdat(i,j,k)+hcns(i,j,k))*hkns(i,j,k) enddo enddo else do j = js-1, je do i = is-1, ie rew(i,j) = (rdat(i+1,j,k)-rdat(i,j,k))*hkew(i,j,k) rns(i,j) = (rdat(i,j+1,k)-rdat(i,j,k))*hkns(i,j,k) enddo enddo endif !----------------------------------------------------------------------- do j = js, je do i = is, ie rdat(i,j,k)=(rew(i,j)-rew(i-1,j)+rns(i,j)-rns(i,j-1)) & *hdac(i,j,k) enddo enddo enddo !----------------------------------------------------------------------- !---- update all halo rows ? ---- ! do not update on last pass, halos will updated in the main program if ( n < order/2 ) then call update_halo (Hgrid, TEMP, rdat) endif !----------------------------------------------------------------------- enddo !----------------------------------------------------------------------- end subroutine diff_mass !####################################################################### subroutine diff_vel (Hgrid, udat, vdat, vkew, vkns, hdac, order, terms, kbot) !---------------------------------------------------------------------- ! ! diff_vel is a private interface that performs multiple ! 2nd order lapacians for the momentum components. ! !---------------------------------------------------------------------- type(horiz_grid_type), intent(inout) :: Hgrid integer, intent(in) :: order real , intent(inout), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & nlev) :: udat, vdat real , intent(in), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & nlev) :: vkew, vkns, hdac real , intent(in), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & nlev, 3), optional :: terms integer, intent(in), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub), & optional :: kbot !---------------------------------------------------------------------- real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) :: & uew, uns, vew, vns, dudp, dvdp, ucew, ucns, vcew, vcns integer :: i, j, k, n, is, ie, js, je, k1, k2 !----------------------------------------------------------------------- is = Hgrid%Vel%is; ie = Hgrid%Vel%ie js = Hgrid%Vel%js; je = Hgrid%Vel%je !------------loop for order of damping scheme--------------------------- do n = 1, order/2 do k = 1, nlev !------------------contributions (fluxes) ------------------------------ if ( n == 1 .and. present(terms) ) then k1 = max(k-1,1) if (present(kbot)) then do j = Hgrid%jlb, Hgrid%jub do i = Hgrid%ilb, Hgrid%iub k2 = min(k+1,kbot(i,j)) dudp(i,j) = (udat(i,j,k2)-udat(i,j,k1))*terms(i,j,k,1) dvdp(i,j) = (vdat(i,j,k2)-vdat(i,j,k1))*terms(i,j,k,1) enddo enddo else k2 = min(k+1,nlev) dudp(:,:) = (udat(:,:,k2)-udat(:,:,k1))*terms(:,:,k,1) dvdp(:,:) = (vdat(:,:,k2)-vdat(:,:,k1))*terms(:,:,k,1) endif do j = js, je+1 do i = is, ie+1 ! ---- slope correction terms ---- ucew(i,j) = (dudp(i,j)+dudp(i-1,j))*terms(i,j,k,2) ucns(i,j) = (dudp(i,j)+dudp(i,j-1))*terms(i,j,k,3) vcew(i,j) = (dvdp(i,j)+dvdp(i-1,j))*terms(i,j,k,2) vcns(i,j) = (dvdp(i,j)+dvdp(i,j-1))*terms(i,j,k,3) uew(i,j) = (udat(i,j,k)-udat(i-1,j ,k)+ucew(i,j))*vkew(i,j,k) uns(i,j) = (udat(i,j,k)-udat(i ,j-1,k)+ucns(i,j))*vkns(i,j,k) vew(i,j) = (vdat(i,j,k)-vdat(i-1,j ,k)+vcew(i,j))*vkew(i,j,k) vns(i,j) = (vdat(i,j,k)-vdat(i ,j-1,k)+vcns(i,j))*vkns(i,j,k) enddo enddo else do j = js, je+1 do i = is, ie+1 uew(i,j) = (udat(i,j,k)-udat(i-1,j ,k))*vkew(i,j,k) uns(i,j) = (udat(i,j,k)-udat(i ,j-1,k))*vkns(i,j,k) vew(i,j) = (vdat(i,j,k)-vdat(i-1,j ,k))*vkew(i,j,k) vns(i,j) = (vdat(i,j,k)-vdat(i ,j-1,k))*vkns(i,j,k) enddo enddo endif !----------------------------------------------------------------------- do j = js, je do i = is, ie udat(i,j,k) = (uew(i+1,j )-uew(i,j)+ & uns(i ,j+1)-uns(i,j))*hdac(i,j,k) vdat(i,j,k) = (vew(i+1,j )-vew(i,j)+ & vns(i ,j+1)-vns(i,j))*hdac(i,j,k) enddo enddo enddo !----------------------------------------------------------------------- !---- update all halo rows ? ---- ! do not update on last pass, halos will updated in the main program if ( n < order/2 ) then call update_halo (Hgrid, UWND, udat) call update_halo (Hgrid, VWND, vdat) endif !----------------------------------------------------------------------- enddo !----------------------------------------------------------------------- end subroutine diff_vel !####################################################################### subroutine horiz_diff_init ( Hgrid ) !----------------------------------------------------------------------- ! initialization of horizontal damping coefficients !----------------------------------------------------------------------- type(horiz_grid_type), intent(in) :: Hgrid ! horizontal grid constants !----------------------------------------------------------------------- real :: eps = 1.e-6 integer :: j integer :: n, nv, order, ntrace, unit, ierr, io real :: coeff, slope(4) character(len=128) :: scheme, params, tname !----------------------------------------------------------------------- ! read namelist if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=bgrid_horiz_diff_nml, iostat=io, end=5) ierr = check_nml_error (io, 'bgrid_horiz_diff_nml') enddo 5 call close_file (unit) endif if (do_log) then call write_version_number (version,tagname) if (mpp_pe() == mpp_root_pe()) write (stdlog(), nml=bgrid_horiz_diff_nml) do_log = .false. endif call get_number_tracers ( MODEL_ATMOS, num_prog=ntrace ) allocate ( Control%order (-1:ntrace), & Control%coeff (-1:ntrace), & Control%do_slope_adj (-1:ntrace), & Control%slope (4,-1:ntrace) ) ! defaults Control%order = 4 Control%coeff = .35 Control%slope = 0. ! namelist arguments Control%order (-1) = damp_order_wind Control%coeff (-1) = damp_coeff_wind Control%slope(:,-1) = slope_corr_wind Control%order (0) = damp_order_temp Control%coeff (0) = damp_coeff_temp Control%slope(:,0) = slope_corr_temp Control%order (1:ntrace) = damp_order_tracer Control%coeff (1:ntrace) = damp_coeff_tracer Control%slope(:,1:ntrace) = spread(slope_corr_tracer,2,ntrace) ! set the damping scheme (check for errors below) Control % damping_scheme_wind = damp_scheme_wind Control % damping_scheme_temp = damp_scheme_temp ! process tracer table information for horizontal damping methods do n = 1, ntrace if (query_method('diff_horiz', MODEL_ATMOS, n, scheme, params)) then if (uppercase(trim(scheme)) /= 'LINEAR' .and. uppercase(trim(scheme)) /= 'NONE') & call error_mesg ('bgrid_horiz_diff_mod', & 'invalid damping method, '//uppercase(trim(scheme)), FATAL) if (parse(params,'order', order) == 1) Control%order(n) = order if (uppercase(trim(scheme)) == 'NONE') Control%order(n) = 0 if (parse(params,'coeff', coeff) == 1) Control%coeff(n) = coeff nv = parse(params,'slope', slope) Control%slope (1:nv,n) = slope(1:nv) endif if (mpp_pe() == mpp_root_pe()) then call get_tracer_names (MODEL_ATMOS, n, tname) write (stdlog(),10) n, trim(tname), Control%order(n), Control%coeff(n), Control%slope(:,n) 10 format (i3, a24, ', Order=',i2, ', Coeff=',f10.5, ', Slope=',4f10.5) endif enddo ! error checking do n = -1, ntrace !-- order if (Control%order(n) < 0 .or. mod(Control%order(n),2) /= 0) & call error_mesg ('bgrid_horiz_diff_mod', 'invalid damping order', FATAL) !-- non-dimension, normalized coefficient if (Control%coeff(n) < 0. .or. Control%coeff(n) > 1.) call error_mesg & ('bgrid_horiz_diff_mod', 'invalid damping coeff', FATAL) !-- slope correction weights if (minval(Control%slope(:,n)) < 0. .or. maxval(Control%slope(:,n)) > 1.) & call error_mesg ('bgrid_horiz_diff_mod', 'invalid slope correction coeff', FATAL) enddo ! set flags Control%do_damping = maxval(Control%order) > 0 Control%do_slope_adj_temp = .false. do n = -1, ntrace Control%do_slope_adj(n) = maxval(Control%slope(:, n)) > 1.e-6 ! set flag for temp OR tracer slope adjustment if (n >= 0 .and. Control%do_slope_adj(n)) Control%do_slope_adj_temp = .true. enddo !----------------------------------------------------------------------- !----- pre-compute metric terms ------ allocate ( Control%areahx (Hgrid%jlb:Hgrid%jub), & Control%areahy (Hgrid%jlb:Hgrid%jub), & Control%areavx (Hgrid%jlb:Hgrid%jub), & Control%areavy (Hgrid%jlb:Hgrid%jub), & Control%wth (Hgrid%jlb:Hgrid%jub,2), & Control%wtv (Hgrid%jlb:Hgrid%jub,2) ) ! ---- areas averaged along axes ---- do j = Hgrid%Tmp%jsd, Hgrid%Tmp%jed-1 Control%areahx(j) = 0.5*(Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j)) Control%areahy(j) = 0.5*(Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j+1)) enddo do j = Hgrid%Vel%jsd+1, Hgrid%Vel%jed Control%areavx(j) = 0.5*(Hgrid%Vel%area(j)+Hgrid%Vel%area(j)) Control%areavy(j) = 0.5*(Hgrid%Vel%area(j)+Hgrid%Vel%area(j-1)) enddo ! ---- damping weight of x and y axis varies ---- ! ---- depending which damping scheme is used ---- call damp_scheme_init ( Hgrid, WIND_GRID, & Control%damping_scheme_wind, Control%wtv ) call damp_scheme_init ( Hgrid, TEMP_GRID, & Control%damping_scheme_temp, Control%wth ) ! initialize code sections for performance timing if (do_clock_init) then id_total = mpp_clock_id ('BGRID: horiz_diff (TOTAL)', & flags=MPP_CLOCK_SYNC, grain=CLOCK_MODULE) do_clock_init = .false. endif !----------------------------------------------------------------------- end subroutine horiz_diff_init !####################################################################### subroutine damp_scheme_init ( Hgrid, grid, scheme, weights ) type(horiz_grid_type), intent(in) :: Hgrid integer, intent(in) :: grid, scheme real, intent(out) :: weights(Hgrid%jlb:Hgrid%jub,2) real, dimension(Hgrid%jlb:Hgrid%jub) :: dx2 real :: dxdy2eq, dx2eq, dy2, factor2 integer :: j ! ---- damping weight of x and y axis varies ---- ! ---- depending which damping scheme is used ---- ! ! scheme 1: equal/constant ! scheme 2: function of diagonal grid distance ! scheme 3: function of x-axis grid distance ! select case (scheme) ! ---- uniform damping ---- case (1) weights = 0.125 case (2:5) dx2=0.0 if (grid == TEMP_GRID) then do j = Hgrid%jlb, Hgrid%jub-1 dx2(j) = Hgrid%Tmp%dx(j)**2 enddo dy2 = Hgrid%Tmp%dy**2 else if (grid == WIND_GRID) then do j = Hgrid%jlb+1, Hgrid%jub dx2(j) = Hgrid%Vel%dx(j)**2 enddo dy2 = Hgrid%Vel%dy**2 else ! error condition needed endif ! function of diagonal distance (poleward of reflat) if (scheme == 2 .or. scheme == 4) then factor2 = 1./cos(reflat*acos(0.)/90.)**2 dxdy2eq = RADIUS**2*(Hgrid%dlm**2+Hgrid%dph**2) weights(:,1) = 0.125*max(1.,dxdy2eq/(factor2*dx2+dy2)) if (scheme == 2) then weights(:,2) = weights(:,1) else weights(:,2) = 0.125 endif endif ! function of x-distance if (scheme == 3 .or. scheme == 5) then dx2eq = (RADIUS*Hgrid%dlm*cos(reflat*acos(0.)/90.))**2 where (dx2 /= 0.0) weights(:,1) = 0.125*max(1.,dx2eq/dx2) if (scheme == 2) then weights(:,2) = weights(:,1) else weights(:,2) = 0.125 endif endif case default call error_mesg ('bgrid_horiz_diff_mod', 'invalid damping scheme', & FATAL) end select end subroutine damp_scheme_init !####################################################################### subroutine slope_correction_init ( Hgrid, Masks, nplev, pres, terms ) type(horiz_grid_type), intent(in) :: Hgrid type (grid_mask_type), intent(in) :: Masks integer, intent(in) :: nplev real, intent(in) :: pres (Hgrid%ilb:,Hgrid%jlb:,:) real, intent(out) :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:) integer :: i, j, k, k1, k2, nlev ! initialization of pressure terms for the sigma slope correction ! these pressure terms do not change between mass variables ! may want make weight a function of variable and/or level ! USE ONE-HALF OF SPECIFIED WEIGHT AT LOWEST LEVEL nlev = size(pres,3) do k = nplev+1, nlev !--- reciprocal of vert gradient --- k1 = max(k-1,1) if (Masks%sigma) then k2 = min(k+1,nlev) terms(:,:,k,1) = 1.0/(pres(:,:,k2)-pres(:,:,k1)) else do j = Hgrid%jlb, Hgrid%jub do i = Hgrid%ilb, Hgrid%iub k2 = min(k+1,Masks%Tmp%kbot(i,j)) if (k2 > k1) then terms(i,j,k,1) = 1.0/(pres(i,j,k2)-pres(i,j,k1)) else terms(i,j,k,1) = 0.0 endif enddo enddo endif ! --- horiz gradients (flip sign) --- do j = Hgrid%jlb, Hgrid%jub-1 do i = Hgrid%ilb, Hgrid%iub-1 terms(i,j,k,2) = (pres(i,j,k)-pres(i+1,j,k)) terms(i,j,k,3) = (pres(i,j,k)-pres(i,j+1,k)) enddo enddo enddo end subroutine slope_correction_init !####################################################################### subroutine vel_slope_correction_init ( Hgrid, Masks, nplev, weights, pres, terms ) type(horiz_grid_type), intent(in) :: Hgrid type (grid_mask_type), intent(in) :: Masks integer, intent(in) :: nplev real, intent(in) :: weights(4) real, intent(in) :: pres (Hgrid%ilb:,Hgrid%jlb:,:) real, intent(out) :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:) real :: wt2 integer :: i, j, k, k1, k2, ks, nlev, isd, ied, jsd, jed real :: dp(size(pres,1),size(pres,2)) ! initialization of pressure terms for the sigma slope correction ! these pressure terms do not change between mass variables ! may want make weight a function of variable and/or level ! USE ONE-HALF OF SPECIFIED WEIGHT AT LOWEST LEVEL nlev = size(pres,3) isd = Hgrid%Vel%isd; ied = Hgrid%Vel%ied jsd = Hgrid%Vel%jsd; jed = Hgrid%Vel%jed do k = 1, nplev terms(:,:,k,:) = 0.0 enddo do k = nplev+1, nlev k1 = max(k-1,1) if (Masks%sigma) then k2 = min(k+1,nlev) dp = (pres(:,:,k2)-pres(:,:,k1)) else do j = Hgrid%jlb, Hgrid%jub do i = Hgrid%ilb, Hgrid%iub k2 = min(k+1,Masks%Vel%kbot(i,j)) dp(i,j) = (pres(i,j,k2)-pres(i,j,k1)) enddo enddo endif !--- reciprocal of vert gradient --- where (dp > 0.) terms(:,:,k,1) = 1.0/dp elsewhere terms(:,:,k,1) = 1.e30 ! these values should not be used where it counts endwhere ! --- horiz gradients (flip sign) --- ks = max(1,k-nlev+4) wt2 = 0.5*weights(ks) do j = Hgrid%jlb+1, Hgrid%jub do i = Hgrid%ilb+1, Hgrid%iub terms(i,j,k,2) = wt2*(pres(i-1,j,k)-pres(i,j,k)) terms(i,j,k,3) = wt2*(pres(i,j-1,k)-pres(i,j,k)) enddo enddo enddo end subroutine vel_slope_correction_init !####################################################################### subroutine slope_correction ( Hgrid, Masks, nplev, weights, terms, temp, cew, cns ) type(horiz_grid_type), intent(in) :: Hgrid type (grid_mask_type), intent(in) :: Masks integer, intent(in) :: nplev real, intent(in) :: weights(4) real, intent(in) :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:) real, intent(in) :: temp (Hgrid%ilb:,Hgrid%jlb:,:) real, intent(out) :: cew (Hgrid%ilb:,Hgrid%jlb:,:),& cns (Hgrid%ilb:,Hgrid%jlb:,:) integer :: i, j, k, k1, k2, ks, nlev real :: wt2 real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: dtdp ! computes weighted-corrections for the slope of sigma surfaces ! to east-west and north-south fluxes of field temp ! check for no correction do k = 1, 4 if (weights(k) > 1.e-6) go to 10 enddo cew = 0.0; cns = 0.0 return 10 nlev = size(temp,3) do k = 1, nplev cew(:,:,k) = 0.0 cns(:,:,k) = 0.0 enddo do k = nplev+1, nlev k1 = max(k-1,1) if (Masks%sigma) then k2 = min(k+1,nlev) ! one-sided derivative at surface dtdp(:,:) = (temp(:,:,k2)-temp(:,:,k1))*terms(:,:,k,1) else do j = Hgrid%jlb, Hgrid%jub do i = Hgrid%ilb, Hgrid%iub k2 = min(k+1,Masks%Tmp%kbot(i,j)) dtdp(i,j) = (temp(i,j,k2)-temp(i,j,k1))*terms(i,j,k,1) enddo enddo endif ks = max(1,k-nlev+4) wt2 = 0.5*weights(ks) do j = Hgrid%jlb, Hgrid%jub-1 do i = Hgrid%ilb, Hgrid%iub-1 ! note: terms are grouped for reproducibility with previous version cew(i,j,k) = (dtdp(i,j)+dtdp(i+1,j))*(terms(i,j,k,2)*wt2) cns(i,j,k) = (dtdp(i,j)+dtdp(i,j+1))*(terms(i,j,k,3)*wt2) enddo enddo enddo end subroutine slope_correction !####################################################################### end module bgrid_horiz_diff_mod