!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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_physics_mod !----------------------------------------------------------------------- ! ! interface for b-grid dynamics with atmospheric physics ! !----------------------------------------------------------------------- use bgrid_core_driver_mod, only: bgrid_dynam_type use bgrid_halo_mod, only: update_halo, NORTH, EAST, WEST, SOUTH, & TEMP, UWND, VWND use bgrid_prog_var_mod, only: prog_var_type use bgrid_vert_mod, only: vert_grid_type, compute_height, & compute_pres_full, compute_pres_half, & compute_pres_depth use bgrid_horiz_mod, only: horiz_grid_type, get_horiz_grid_bound, TGRID use bgrid_change_grid_mod, only: change_grid, TEMP_GRID, WIND_GRID use bgrid_integrals_mod, only: global_integral use bgrid_masks_mod, only: grid_mask_type use time_manager_mod, only: time_type, get_time, operator(-) use fms_mod, only: error_mesg, FATAL, write_version_number , mpp_pe use physics_driver_mod, only: physics_driver_init, & physics_driver_end, & physics_driver_moist_init, & physics_driver_moist_end, & physics_driver_down_time_vary, & physics_driver_up_time_vary, & physics_driver_down_endts, & physics_driver_up_endts, & physics_driver_down, & physics_driver_up, & surf_diff_type use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_tracer_index, NO_TRACER !----------------------------------------------------------------------- implicit none private public bgrid_physics_down, bgrid_physics_up, & bgrid_physics_init, bgrid_physics_end public surf_diff_type !----------------------------------------------------------------------- character(len=128) :: version = '$Id: bgrid_physics.F90,v 17.0.2.1 2009/10/22 17:03:10 wfc Exp $' character(len=128) :: tag = '$Name: mom4p1_pubrel_dec2009_nnz $' !----------------------------------------------------------------------- real, allocatable, dimension(:,:,:) :: u_dt, v_dt contains !####################################################################### subroutine bgrid_physics_down (window, dt_phys, & Time_prev, Time, Time_next, & Hgrid, Vgrid, Dynam, & Var, Var_dt, & frac_land, albedo, & albedo_vis_dir, albedo_nir_dir, & albedo_vis_dif, albedo_nir_dif, & rough_vel, t_surf, & u_star, b_star, q_star, & dtau_du, dtau_dv, tau_x, tau_y, & flux_sw, & flux_sw_dir, & flux_sw_dif, & flux_sw_down_vis_dir, & flux_sw_down_vis_dif, & flux_sw_down_total_dir, & flux_sw_down_total_dif, & flux_sw_vis, & flux_sw_vis_dir, & flux_sw_vis_dif, & flux_lw, coszen, & gust, Surf_diff ) !----------------------------------------------------------------------- ! ! Time_prev = time at the previous time level, tau-1 (time_type) ! Time = time at the current time level, tau (time_type) ! Time_next = time at the next time level, tau+1 (time_type) ! ! NOTE: for a two time level scheme (e.g., forward-backward scheme) ! Time_prev = Time. ! !----------------------------------------------------------------------- integer, intent(in) :: window(2) real, intent(in) :: dt_phys type(time_type),intent(in) :: Time_prev, Time, Time_next type (horiz_grid_type),intent(inout) :: Hgrid type (vert_grid_type),intent(in) :: Vgrid type(bgrid_dynam_type),intent(in) :: Dynam type (prog_var_type),intent(inout) :: Var type (prog_var_type),intent(inout) :: Var_dt ! Note: Var is intent inout because of diagnostic tracers real, intent(in), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: & frac_land, albedo, & albedo_vis_dir, albedo_nir_dir, & albedo_vis_dif, albedo_nir_dif, & rough_vel, t_surf, & u_star, b_star, & q_star, dtau_du, dtau_dv real, intent(inout), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: tau_x, tau_y real, intent(out), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: & flux_sw, & flux_sw_dir, & flux_sw_dif, & flux_lw, coszen, & flux_sw_down_vis_dir, & flux_sw_down_vis_dif, & flux_sw_down_total_dir, & flux_sw_down_total_dif, & flux_sw_vis, & flux_sw_vis_dir, & flux_sw_vis_dif, & gust type(surf_diff_type), intent(inout) :: Surf_diff !----------------------------------------------------------------------- integer :: j, k, n, is, ie, js, je, i1, i2, j1, j2, sphum, nt, ntp integer :: ix, jx, idim, jdim real :: dt integer :: sec, day !----------------------------------------------------------------------- real, dimension(window(1),window(2),Vgrid%nlev) :: p_full, z_full real, dimension(window(1),window(2),Vgrid%nlev+1) :: p_half, z_half real, dimension(1,1,Vgrid%nlev+1) :: phalfgrey real, dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & Vgrid%nlev) :: uh, vh real, dimension(window(1),window(2)) :: pssl_new, area real, dimension(size(Var%r, 4)) :: gavg_rrv real, dimension(1,1) :: psurf !----------------------------------------------------------------------- !---------------------------- do physics ------------------------------- idim = window(1) jdim = window(2) ! --- momentum and momentum tendency on mass grid --- ! note: tendency is saved for "physics_up" call call update_halo (Hgrid, UWND, Var_dt%u, halos=SOUTH+WEST) call update_halo (Hgrid, VWND, Var_dt%v, halos=SOUTH+WEST) if (Dynam%Masks%sigma) then call change_grid (Hgrid, WIND_GRID, TEMP_GRID, Var%u, Var%v, uh , vh ) call change_grid (Hgrid, WIND_GRID, TEMP_GRID, Var_dt%u, Var_dt%v, u_dt, v_dt) else call change_grid (Hgrid, WIND_GRID, TEMP_GRID, & Var%u, Var%v, uh , vh , mask_inp=Dynam%Masks%Vel%mask) call change_grid (Hgrid, WIND_GRID, TEMP_GRID, & Var_dt%u, Var_dt%v, u_dt, v_dt, mask_inp=Dynam%Masks%Vel%mask) endif ! --- loop through physics windows --- sphum = get_tracer_index ( MODEL_ATMOS, 'sphum' ) if (sphum <= 0) call error_mesg ('bgrid_physics_mod', & 'specific humidity tracer not found', FATAL) ntp = Var_dt%ntrace nt = Var %ntrace js = Hgrid%Tmp%js gavg_rrv = 0. call compute_g_avg(Hgrid, Vgrid, Var, Dynam%Masks, gavg_rrv, 'co2') !--------------------------------------------------------------------- ! compute the physics time step (from tau-1 to tau+1). !--------------------------------------------------------------------- call get_time (Time_next-Time_prev, sec, day) dt = real(sec+day*86400) !--------------------------------------------------------------------- ! call physics_driver_down_time_vary to do the time-dependent, spatially ! independent calculations before entering windows / threads loop. !--------------------------------------------------------------------- call physics_driver_down_time_vary (Time, Time_next, gavg_rrv, dt) do while ( js <= Hgrid%Tmp%je ) je = min ( js+jdim-1, Hgrid%Tmp%je ) jx = je-js+1 is = Hgrid%Tmp%is do while ( is <= Hgrid%Tmp%ie ) ie = min ( is+idim-1, Hgrid%Tmp%ie ) ix = ie-is+1 ! ---- pass updated surface pressure ---- do j = 1, jx pssl_new(1:ix,j) = Var %pssl(is:ie,js+j-1) + & Var_dt%pssl(is:ie,js+j-1) * dt_phys area (1:ix,j) = Hgrid%Tmp%area(js+j-1) enddo call compute_pres_full (Vgrid, pssl_new(1:ix,1:jx), & p_full(1:ix,1:jx,:)) call compute_pres_half (Vgrid, pssl_new(1:ix,1:jx), & p_half(1:ix,1:jx,:)) ! compute a reference profile of pressure based on psurf = 1000 hPa psurf = reshape ( (/ 100000. /), (/ 1, 1 /) ) call compute_pres_half (Vgrid, psurf, phalfgrey) ! ----- compute height (in meters) ------- ! ----- spec humidity assumed to be tracer #1 -------- call compute_height (Vgrid, Dynam%fisl(is:ie,js:je), & Var%t (is:ie,js:je,:), & Var%r (is:ie,js:je,:,sphum), & p_full (1:ix,1:jx,:), & p_half (1:ix,1:jx,:), & z_full (1:ix,1:jx,:), & z_half (1:ix,1:jx,:), & Dynam%Masks%Tmp%mask(is:ie,js:je,:) ) ! ---- j-axis indices in the global physics grid ---- j1 = js-Hgrid%Tmp%js+1; j2 = j1+(je-js) i1 = is-Hgrid%Tmp%is+1; i2 = i1+(ie-is) !----------------------------------------------------------------------- !-------------------------- call physics ------------------------------- !------------ (need to add leap-frog option for uh,vh) ----------------- !----------------------------------------------------------------------- if ( .not. Dynam%Masks%sigma ) then !------------ eta coordinate ------------------------------------------- call physics_driver_down (i1, i2, j1, j2 ,& Time_prev, Time, Time_next ,& Hgrid%Tmp%aph (is:ie,js:je), Hgrid%Tmp%alm (is:ie,js:je),& area ( 1:ix, 1:jx), p_half ( 1:ix, 1:jx,:) ,& p_full ( 1:ix, 1:jx,:), z_half ( 1:ix, 1:jx,:) ,& z_full ( 1:ix, 1:jx,:) ,& phalfgrey ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) , Var%t(is:ie,js:je,:) ,& Var%r(is:ie,js:je,:,sphum) , Var%r(is:ie,js:je,:, : ) ,& uh (is:ie,js:je,:) , vh (is:ie,js:je,:) ,& Var%t (is:ie,js:je,:) , Var%r (is:ie,js:je,:,sphum) ,& Var%r (is:ie,js:je,:, : ) , & frac_land(is:ie,js:je) , rough_vel(is:ie,js:je) ,& albedo (is:ie,js:je) , albedo_vis_dir (is:ie,js:je) ,& albedo_nir_dir (is:ie,js:je) , albedo_vis_dif (is:ie,js:je) ,& albedo_nir_dif (is:ie,js:je) , t_surf (is:ie,js:je) ,& u_star (is:ie,js:je) , b_star (is:ie,js:je) ,& q_star (is:ie,js:je) , dtau_du (is:ie,js:je), dtau_dv (is:ie,js:je),& tau_x (is:ie,js:je) , tau_y (is:ie,js:je) ,& u_dt (is:ie,js:je,:) , v_dt (is:ie,js:je,:) ,& Var_dt%t (is:ie,js:je,:) , Var_dt%r (is:ie,js:je,:,sphum),& Var_dt%r (is:ie,js:je,:,1:ntp) , flux_sw (is:ie,js:je) ,& flux_sw_dir (is:ie,js:je) ,& flux_sw_dif (is:ie,js:je) ,& flux_sw_down_vis_dir (is:ie,js:je) ,& flux_sw_down_vis_dif (is:ie,js:je) ,& flux_sw_down_total_dir (is:ie,js:je) ,& flux_sw_down_total_dif (is:ie,js:je) ,& flux_sw_vis (is:ie,js:je) ,& flux_sw_vis_dir (is:ie,js:je) ,& flux_sw_vis_dif (is:ie,js:je) ,& flux_lw (is:ie,js:je) ,& coszen (is:ie,js:je) ,& gust (is:ie,js:je) ,& Surf_diff ,& gavg_rrv ,& mask=Dynam%Masks%Tmp%mask(is:ie,js:je,:) ,& kbot=Dynam%Masks%Tmp%kbot(is:ie,js:je) ) else !------------- sigma coordinate ---------------------------------------- call physics_driver_down (i1, i2, j1, j2 ,& Time_prev, Time, Time_next ,& Hgrid%Tmp%aph (is:ie,js:je), Hgrid%Tmp%alm (is:ie,js:je),& area ( 1:ix, 1:jx), p_half ( 1:ix, 1:jx,:) ,& p_full ( 1:ix, 1:jx,:), z_half ( 1:ix, 1:jx,:) ,& z_full ( 1:ix, 1:jx,:) ,& phalfgrey ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) , Var%t(is:ie,js:je,:) ,& Var%r(is:ie,js:je,:,sphum) , Var%r(is:ie,js:je,:, : ) ,& uh (is:ie,js:je,:) , vh (is:ie,js:je,:) ,& Var%t (is:ie,js:je,:) , Var%r (is:ie,js:je,:,sphum) ,& Var%r (is:ie,js:je,:, : ) , & frac_land(is:ie,js:je) , rough_vel(is:ie,js:je) ,& albedo (is:ie,js:je) , albedo_vis_dir (is:ie,js:je) ,& albedo_nir_dir (is:ie,js:je) , albedo_vis_dif (is:ie,js:je) ,& albedo_nir_dif (is:ie,js:je) , t_surf (is:ie,js:je) ,& u_star (is:ie,js:je) , b_star (is:ie,js:je) ,& q_star (is:ie,js:je) , dtau_du (is:ie,js:je) ,dtau_dv (is:ie,js:je),& tau_x (is:ie,js:je) , tau_y (is:ie,js:je) ,& u_dt (is:ie,js:je,:) , v_dt (is:ie,js:je,:) ,& Var_dt%t (is:ie,js:je,:) , Var_dt%r (is:ie,js:je,:,sphum),& Var_dt%r (is:ie,js:je,:,1:ntp) , flux_sw (is:ie,js:je) ,& ! Time_prev, Time, Time_next ,& ! Hgrid%Tmp%aph (is:ie,js:je) ,& ! Hgrid%Tmp%alm (is:ie,js:je) ,& ! area ( 1:ix, 1:jx) ,& ! p_half ( 1:ix, 1:jx,:) ,& ! p_full ( 1:ix, 1:jx,:) ,& ! z_half ( 1:ix, 1:jx,:) ,& ! z_full ( 1:ix, 1:jx,:) ,& ! uh (is:ie,js:je,:) ,& ! vh (is:ie,js:je,:) ,& ! Var%t(is:ie,js:je,:) ,& ! Var%r(is:ie,js:je,:,sphum) ,& ! Var%r(is:ie,js:je,:,1:ntp) ,& ! uh (is:ie,js:je,:) ,& ! vh (is:ie,js:je,:) ,& ! Var%t (is:ie,js:je,:) ,& ! Var%r (is:ie,js:je,:,sphum) ,& ! Var%r (is:ie,js:je,:,1:ntp) ,& ! Var%r (is:ie,js:je,:,ntp+1:nt) ,& ! frac_land(is:ie,js:je) ,& ! rough_vel(is:ie,js:je) ,& ! albedo (is:ie,js:je) ,& ! albedo_vis_dir (is:ie,js:je) ,& ! albedo_nir_dir (is:ie,js:je) ,& ! albedo_vis_dif (is:ie,js:je) ,& ! albedo_nir_dif (is:ie,js:je) ,& ! t_surf (is:ie,js:je) ,& ! u_star (is:ie,js:je) ,& ! b_star (is:ie,js:je) ,& ! q_star (is:ie,js:je) ,& ! dtau_dv (is:ie,js:je) ,& ! tau_x (is:ie,js:je) ,& ! tau_y (is:ie,js:je) ,& ! u_dt (is:ie,js:je,:) ,& ! v_dt (is:ie,js:je,:) ,& ! Var_dt%t (is:ie,js:je,:) ,& ! Var_dt%r (is:ie,js:je,:,sphum) ,& ! Var_dt%r (is:ie,js:je,:,1:ntp) ,& ! flux_sw (is:ie,js:je) ,& flux_sw_dir (is:ie,js:je) ,& flux_sw_dif (is:ie,js:je) ,& flux_sw_down_vis_dir (is:ie,js:je) ,& flux_sw_down_vis_dif (is:ie,js:je) ,& flux_sw_down_total_dir (is:ie,js:je) ,& flux_sw_down_total_dif (is:ie,js:je) ,& flux_sw_vis (is:ie,js:je) ,& flux_sw_vis_dir (is:ie,js:je) ,& flux_sw_vis_dif (is:ie,js:je) ,& flux_lw (is:ie,js:je) ,& coszen (is:ie,js:je) ,& gust (is:ie,js:je) ,& Surf_diff ,& gavg_rrv ) endif is = is + idim enddo js = js + jdim enddo call physics_driver_down_endts (is-Hgrid%Tmp%is+1, js-Hgrid%Tmp%js+1) ! halo rows for tendencies do not need updating until after physics_up ! update halos for diagnostic tracers only if (ntp < nt) call update_halo (Hgrid, TEMP, Var%r(:,:,:,ntp+1:nt)) !----------------------------------------------------------------------- end subroutine bgrid_physics_down !####################################################################### subroutine bgrid_physics_up (window, dt_phys, & Time_prev, Time, Time_next, & Hgrid, Vgrid, Dynam, & Var, Var_dt, omega, & frac_land, Surf_diff, lprec, fprec, gust, & u_star, b_star, q_star ) !----------------------------------------------------------------------- ! ! Time_prev = time at the previous time level, tau-1 (time_type) ! Time = time at the current time level, tau (time_type) ! Time_next = time at the next time level, tau+1 (time_type) ! ! NOTE: for a two time level scheme (e.g., forward-backward scheme) ! Time_prev = Time. ! !----------------------------------------------------------------------- integer, intent(in) :: window(2) real, intent(in) :: dt_phys type(time_type),intent(in) :: Time_prev, Time, Time_next type (horiz_grid_type),intent(inout) :: Hgrid type (vert_grid_type),intent(in) :: Vgrid type(bgrid_dynam_type),intent(in) :: Dynam type (prog_var_type),intent(in) :: Var type (prog_var_type),intent(inout) :: Var_dt real, intent(in), dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & Vgrid%nlev) :: omega real, intent(in), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: frac_land type(surf_diff_type), intent(inout) :: Surf_diff real, intent(out), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: lprec, fprec real, intent(inout), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: gust real, intent(in), dimension(Hgrid%Tmp%is:,Hgrid%Tmp%js:) :: u_star, b_star, q_star !----------------------------------------------------------------------- integer :: j, k, n, is, ie, js, je, i1, i2, j1, j2, sphum, ntp, npz integer :: ix, jx, idim, jdim integer :: sec, day real :: dt !----------------------------------------------------------------------- real, dimension(window(1),window(2),Vgrid%nlev) :: p_full, z_full real, dimension(window(1),window(2),Vgrid%nlev+1) :: p_half, z_half real, dimension(Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, & Vgrid%nlev) :: uh, vh, uh_dt, vh_dt real, dimension(window(1),window(2)) :: pssl_new, area !--------------------------------------------------------------------- ! compute the physics time step (from tau-1 to tau+1). !--------------------------------------------------------------------- call get_time (Time_next-Time_prev, sec, day) dt = real(sec+day*86400) call physics_driver_up_time_vary (Time, dt) !----------------------------------------------------------------------- !---------------------------- do physics ------------------------------- idim = window(1) jdim = window(2) ! --- momentum and previous momentum tendency on mass grid --- if (Dynam%Masks%sigma) then call change_grid (Hgrid, WIND_GRID, TEMP_GRID, Var%u, Var%v, uh , vh ) call change_grid (Hgrid, WIND_GRID, TEMP_GRID, Var_dt%u, Var_dt%v, uh_dt, vh_dt) else call change_grid (Hgrid, WIND_GRID, TEMP_GRID, & Var%u, Var%v, uh , vh , mask_inp=Dynam%Masks%Vel%mask) call change_grid (Hgrid, WIND_GRID, TEMP_GRID, & Var_dt%u, Var_dt%v, uh_dt, vh_dt, mask_inp=Dynam%Masks%Vel%mask) endif ! --- loop through physics windows --- sphum = get_tracer_index ( MODEL_ATMOS, 'sphum' ) if (sphum <= 0) call error_mesg ('bgrid_physics_mod', & 'specific humidity tracer not found', FATAL) ntp = Var_dt%ntrace js = Hgrid%Tmp%js npz = size(p_full,3) call physics_driver_moist_init (window(1), window(2), npz, ntp) do while ( js <= Hgrid%Tmp%je ) je = min ( js+jdim-1, Hgrid%Tmp%je ) jx = je-js+1 is = Hgrid%Tmp%is do while ( is <= Hgrid%Tmp%ie ) ie = min ( is+idim-1, Hgrid%Tmp%ie ) ix = ie-is+1 ! ---- pass updated surface pressure ---- do j = 1, jx pssl_new(1:ix,j) = Var %pssl(is:ie,js+j-1) + & Var_dt%pssl(is:ie,js+j-1) * dt_phys area (1:ix,j) = Hgrid%Tmp%area(js+j-1) enddo call compute_pres_full (Vgrid, pssl_new(1:ix,1:jx), & p_full(1:ix,1:jx,:)) call compute_pres_half (Vgrid, pssl_new(1:ix,1:jx), & p_half(1:ix,1:jx,:)) call compute_height (Vgrid, Dynam%fisl(is:ie,js:je), & Var%t (is:ie,js:je,:), & Var%r (is:ie,js:je,:,sphum), & p_full (1:ix,1:jx,:), & p_half (1:ix,1:jx,:), & z_full (1:ix,1:jx,:), & z_half (1:ix,1:jx,:), & Dynam%Masks%Tmp%mask(is:ie,js:je,:) ) ! ---- j-axis indices in the global physics grid ---- j1 = js-Hgrid%Tmp%js+1; j2 = j1+(je-js) i1 = is-Hgrid%Tmp%is+1; i2 = i1+(ie-is) !----------------------------------------------------------------------- !-------------------------- call physics ------------------------------- !------------ (need to add leap-frog option for uh,vh) ----------------- !----------------------------------------------------------------------- if ( .not. Dynam%Masks%sigma ) then !------------ eta coordinate ------------------------------------------- call physics_driver_up (i1, i2, j1, j2 ,& Time_prev, Time, Time_next ,& Hgrid%Tmp%aph (is:ie,js:je) ,& Hgrid%Tmp%alm (is:ie,js:je) ,& area ( 1:ix, 1:jx) ,& p_half ( 1:ix, 1:jx,:) ,& p_full ( 1:ix, 1:jx,:) ,& z_half ( 1:ix, 1:jx,:) ,& z_full ( 1:ix, 1:jx,:) ,& omega (is:ie,js:je,:) ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) ,& Var%t(is:ie,js:je,:) ,& Var%r(is:ie,js:je,:,sphum) ,& Var%r(is:ie,js:je,:,1:ntp) ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) ,& Var%t (is:ie,js:je,:) ,& Var%r (is:ie,js:je,:,sphum) ,& Var%r (is:ie,js:je,:,1:ntp) ,& frac_land(is:ie,js:je) ,& u_star(is:ie,js:je) ,& b_star(is:ie,js:je) ,& q_star(is:ie,js:je) ,& u_dt (is:ie,js:je,:) ,& v_dt (is:ie,js:je,:) ,& Var_dt%t (is:ie,js:je,:) ,& Var_dt%r (is:ie,js:je,:,sphum) ,& Var_dt%r (is:ie,js:je,:,1:ntp) ,& Surf_diff ,& lprec (is:ie,js:je) ,& fprec (is:ie,js:je) ,& gust (is:ie,js:je) ,& mask=Dynam%Masks%Tmp%mask(is:ie,js:je,:) ,& kbot=Dynam%Masks%Tmp%kbot(is:ie,js:je) ) else !------------- sigma coordinate ---------------------------------------- call physics_driver_up (i1, i2, j1, j2 ,& Time_prev, Time, Time_next ,& Hgrid%Tmp%aph (is:ie,js:je) ,& Hgrid%Tmp%alm (is:ie,js:je) ,& area ( 1:ix, 1:jx) ,& p_half ( 1:ix, 1:jx,:) ,& p_full ( 1:ix, 1:jx,:) ,& z_half ( 1:ix, 1:jx,:) ,& z_full ( 1:ix, 1:jx,:) ,& omega (is:ie,js:je,:) ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) ,& Var%t(is:ie,js:je,:) ,& Var%r(is:ie,js:je,:,sphum) ,& Var%r(is:ie,js:je,:,1:ntp) ,& uh (is:ie,js:je,:) ,& vh (is:ie,js:je,:) ,& Var%t (is:ie,js:je,:) ,& Var%r (is:ie,js:je,:,sphum) ,& Var%r (is:ie,js:je,:,1:ntp) ,& frac_land(is:ie,js:je) ,& u_star(is:ie,js:je) ,& b_star(is:ie,js:je) ,& q_star(is:ie,js:je) ,& u_dt (is:ie,js:je,:) ,& v_dt (is:ie,js:je,:) ,& Var_dt%t (is:ie,js:je,:) ,& Var_dt%r (is:ie,js:je,:,sphum) ,& Var_dt%r (is:ie,js:je,:,1:ntp) ,& Surf_diff ,& lprec (is:ie,js:je) ,& fprec (is:ie,js:je) ,& gust (is:ie,js:je) ) endif is = is + idim enddo js = js + jdim enddo call physics_driver_moist_end call physics_driver_up_endts (is-Hgrid%Tmp%is+1, js-Hgrid%Tmp%js+1) !----------------------------------------------------------------------- ! compute momentum tendencies on mass grid for physics only ! udt(phys) = udt(current) - udt(before physics) uh(:,:,:) = u_dt(:,:,:) - uh_dt(:,:,:) vh(:,:,:) = v_dt(:,:,:) - vh_dt(:,:,:) ! update halos of momentum tendencies on mass grid ! then move momentum tendencies to momentum grid call update_halo (Hgrid, TEMP, uh, halos=NORTH+EAST) call update_halo (Hgrid, TEMP, vh, halos=NORTH+EAST) call change_grid (Hgrid, TEMP_GRID, WIND_GRID, uh, vh, uh, vh) uh(:,Hgrid%jub,:) = 0.0 ! zero out unused polar halo row vh(:,Hgrid%jub,:) = 0.0 ! no harm done when not polar row !---- update momentum tendencies ---- Var_dt%u = Var_dt%u + uh * Dynam%Masks%Vel%mask Var_dt%v = Var_dt%v + vh * Dynam%Masks%Vel%mask !---- update all halo rows ---- call update_halo (Hgrid, TEMP, Var_dt%t) call update_halo (Hgrid, TEMP, Var_dt%r) call update_halo (Hgrid, UWND, Var_dt%u) call update_halo (Hgrid, VWND, Var_dt%v) !----------------------------------------------------------------------- end subroutine bgrid_physics_up !####################################################################### subroutine bgrid_physics_init (axes, Time, Hgrid, Vgrid, Dynam, & Var, Surf_diff) !----------------------------------------------------------------------- ! ! axes = array of axis indices for diagnostics (x,y,pf,ph) ! Time = current time (time_type) ! Hgrid = horizontal grid constants ! Vgrid = vertical grid constants ! Dynam = current state of the dynamical core ! Var = prognostic variables ! !----------------------------------------------------------------------- integer, intent(in) :: axes(4) type (time_type), intent(in) :: Time type (horiz_grid_type),intent(inout) :: Hgrid type (vert_grid_type),intent(in) :: Vgrid type(bgrid_dynam_type),intent(in) :: Dynam type (prog_var_type),intent(inout) :: Var type (surf_diff_type), intent(inout) :: Surf_diff !----------------------------------------------------------------------- real, dimension(Hgrid%Tmp%is:Hgrid%Tmp%ie+1) :: lonb real, dimension(Hgrid%Tmp%js:Hgrid%Tmp%je+1) :: latb real, dimension(Hgrid%Tmp%is:Hgrid%Tmp%ie, & Hgrid%Tmp%js:Hgrid%Tmp%je,Vgrid%nlev+1) :: phalf real, dimension(Vgrid%nlev+1,2) :: pref real, dimension(2,1,Vgrid%nlev) :: pref_full real, dimension(2,1) :: pref_sl real, dimension(Hgrid%Tmp%is:Hgrid%Tmp%ie+1,Hgrid%Tmp%js:Hgrid%Tmp%je+1) :: lonb2d, latb2d !----------------------------------------------------------------------- integer :: i, j, k, n, nt, is, ie, js, je, unit !----------------------------------------------------------------------- nt = Var%ntrace !----- write version to logfile -------- call write_version_number (version,tag) !---------- get local grid box edges --------- call get_horiz_grid_bound (Hgrid, TGRID, lonb, latb) do j = Hgrid%Tmp%js,Hgrid%Tmp%je+1 lonb2d(:,j) = lonb(:) end do do i = Hgrid%Tmp%is,Hgrid%Tmp%ie+1 latb2d(i,:) = latb(:) end do !---------- reference profile ----------- pref_sl = reshape ( (/ 101325., 81060. /), (/ 2, 1 /) ) call compute_pres_full (Vgrid, pref_sl, pref_full) pref(1:Vgrid%nlev,1) = pref_full(1,1,:) pref(1:Vgrid%nlev,2) = pref_full(2,1,:) pref(Vgrid%nlev+1,1) = pref_sl(1,1) pref(Vgrid%nlev+1,2) = pref_sl(2,1) !------- pressure at model layer interfaces ----- is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie js = Hgrid % Tmp % js; je = Hgrid % Tmp % je call compute_pres_half (Vgrid, Var%pssl(is:ie,js:je), phalf) !---------- initialize physics ------- if (Dynam%Masks%sigma) then call physics_driver_init (Time, lonb2d, latb2d, axes, pref, & Var%r(is:ie,js:je,:,1:nt), & Surf_diff, phalf ) else call physics_driver_init (Time, lonb2d, latb2d, axes, pref, & Var%r(is:ie,js:je,:,1:nt), & Surf_diff, phalf, & mask=Dynam%Masks%Tmp%mask(is:ie,js:je,:), & kbot=Dynam%Masks%Vel%kbot(is:ie,js:je) ) endif ! ---- boundaries for tracers ---- call update_halo (Hgrid, TEMP, Var%r) ! ---- storage for global data ---- allocate (u_dt(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,Vgrid%nlev)) allocate (v_dt(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,Vgrid%nlev)) !----------------------------------------------------------------------- end subroutine bgrid_physics_init !####################################################################### subroutine bgrid_physics_end (Time) !----------------------------------------------------------------------- type(time_type), intent(in) :: Time !----------------------------------------------------------------------- ! NOTE: this is not the dynamics time !----------------------------------------------------------------------- call physics_driver_end (Time) !----------------------------------------------------------------------- end subroutine bgrid_physics_end !####################################################################### subroutine compute_g_avg(Hgrid, Vgrid, Var, Masks, rrv, tracer_name) type (horiz_grid_type),intent(inout) :: Hgrid type (vert_grid_type),intent(in) :: Vgrid type (prog_var_type),intent(in) :: Var type (grid_mask_type), intent(in) :: Masks real, dimension(:), intent(inout) :: rrv character(len=*), intent(in) :: tracer_name real, dimension (Hgrid%ilb:Hgrid%iub, & Hgrid%jlb:Hgrid%jub, Var%nlev) :: dpde, avg real avgps real psfc_sum, qp_sum, qp integer j, i, k, idx psfc_sum = 0. qp_sum = 0. idx = get_tracer_index(MODEL_ATMOS, trim(tracer_name)) if(idx /= NO_TRACER) then call compute_pres_depth (Vgrid, Var%pssl, dpde) avg = Var%r(:,:,:,idx) * dpde avgps = global_integral (Hgrid, 1, Var%ps, do_exact=.true.) rrv(idx) = global_integral(Hgrid, 1, avg, Masks, .true.) / avgps endif end subroutine compute_g_avg end module bgrid_physics_mod