From 306ff31371e74694e5d9f4a57584295c7122b9ac Mon Sep 17 00:00:00 2001 From: XiaqiongZhou-NOAA <48254930+XiaqiongZhou-NOAA@users.noreply.github.com> Date: Mon, 25 Jan 2021 16:09:19 -0500 Subject: [PATCH 1/2] Add an option to use zero-gradient BC to reconstruct interface u/v for height advection; change dz_min as a input (#53) * 1)Add an option to use original BC or zero-gradient BC to reconstruct interface u/v with PSM for height advection 2)Change dz_min as a namelist input * Change to the original version by adjust the upper interface of ZH when dz Atm%flagstruct%nudge_qv add_noise => Atm%flagstruct%add_noise butterfly_effect => Atm%flagstruct%butterfly_effect + dz_min => Atm%flagstruct%dz_min + psm_bc => Atm%flagstruct%psm_bc a2b_ord => Atm%flagstruct%a2b_ord c2l_ord => Atm%flagstruct%c2l_ord ndims => Atm%flagstruct%ndims @@ -1047,6 +1051,10 @@ subroutine read_namelist_fv_core_nml(Atm) !! !> \param[in] butterfly\_effect Logical: whether to flip the least-significant-bit of the lowest level temperature. False by default. !! +!> \param[in] dz\_min Real: Minimum layer thickness to enforce monotonicity of height to prevent blowup. Set to 2 by default. +!! +!> \param[in] psm\_bc Integer: Options to use origional BCs (0) or zero-gradient BCs (1) to reconstruct interface u/v for the advection of height. Set to 0 by default. +!! !> \param[in] adjust\_dry\_mass Logical: whether to adjust the global dry-air mass to the value set by dry\_mass. This is only done in an initialization step, particularly when using an initial condition from an external dataset, interpolated from another resolution (either horizontal or vertical), or when changing the topography, so that the global mass of the atmosphere matches some estimate of observed value. False by default. It is recommended to only set this to True when initializing the model. !! !> \param[in] breed\_vortex\_inline Logical: whether to bogus tropical cyclones into the model, which are specified by an external file. Options are set in fv\_nwp\_nudge\_nml. False by default. @@ -1438,7 +1446,7 @@ subroutine read_namelist_fv_core_nml(Atm) c2l_ord, dx_const, dy_const, umax, deglat, & deglon_start, deglon_stop, deglat_start, deglat_stop, & phys_hydrostatic, use_hydro_pressure, make_hybrid_z, old_divg_damp, add_noise, butterfly_effect, & - nested, twowaynest, nudge_qv, & + dz_min, psm_bc, nested, twowaynest, nudge_qv, & nestbctype, nestupdate, nsponge, s_weight, & check_negative, nudge_ic, halo_update_type, gfs_phil, agrid_vel_rst, & do_uni_zfull, adj_mass_vmr, fac_n_spl, fhouri, update_blend, regional, bc_update_interval, & diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 893a8725c..c5abae66b 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -1,2218 +1,2341 @@ - -!*********************************************************************** -!* GNU Lesser General Public License -!* -!* This file is part of the FV3 dynamical core. -!* -!* The FV3 dynamical core is free software: you can redistribute it -!* and/or modify it under the terms of the -!* GNU Lesser General Public License as published by the -!* Free Software Foundation, either version 3 of the License, or -!* (at your option) any later version. -!* -!* The FV3 dynamical core is distributed in the hope that it will be -!* useful, but WITHOUT ANYWARRANTY; 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 Lesser General Public -!* License along with the FV3 dynamical core. -!* If not, see . -!*********************************************************************** - -!>@brief The module 'nh_utils' peforms non-hydrostatic computations. -!>@author S. J. Lin, NOAA/GFDL - -module nh_utils_mod - -! Modules Included: -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -! -!
Module NameFunctions Included
constants_modrdgas, cp_air, grav
fv_arrays_modfv_grid_bounds_type, fv_grid_type
sw_core_modfill_4corners, del6_vt_flux
tp_core_modfv_tp_2d
- - use constants_mod, only: rdgas, cp_air, grav - use tp_core_mod, only: fv_tp_2d - use sw_core_mod, only: fill_4corners, del6_vt_flux - use fv_arrays_mod, only: fv_grid_bounds_type, fv_grid_type,fv_nest_BC_type_3d -#ifdef MULTI_GASES - use multi_gases_mod, only: vicpqd, vicvqd -#endif - - implicit none - private - - public update_dz_c, update_dz_d, nh_bc - public sim_solver, sim1_solver, sim3_solver - public sim3p0_solver, rim_2d - public Riem_Solver_c - - real, parameter:: dz_min = 6. - real, parameter:: r3 = 1./3. - -CONTAINS - - subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, & - npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type) -! !INPUT PARAMETERS: - type(fv_grid_bounds_type), intent(IN) :: bd - integer, intent(in):: is, ie, js, je, ng, km, npx, npy, grid_type - logical, intent(IN):: sw_corner, se_corner, ne_corner, nw_corner - real, intent(in):: dt - real, intent(in):: dp0(km) - real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: ut, vt - real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng):: area - real, intent(inout):: gz(is-ng:ie+ng,js-ng:je+ng,km+1) - real, intent(in ):: zs(is-ng:ie+ng, js-ng:je+ng) - real, intent( out):: ws(is-ng:ie+ng, js-ng:je+ng) -! Local Work array: - real:: gz2(is-ng:ie+ng,js-ng:je+ng) - real, dimension(is-1:ie+2,js-1:je+1):: xfx, fx - real, dimension(is-1:ie+1,js-1:je+2):: yfx, fy - real, parameter:: r14 = 1./14. - integer i, j, k - integer:: is1, ie1, js1, je1 - integer:: ie2, je2 - real:: rdt, top_ratio, bot_ratio, int_ratio -!-------------------------------------------------------------------- - - rdt = 1. / dt - - top_ratio = dp0(1 ) / (dp0( 1)+dp0(2 )) - bot_ratio = dp0(km) / (dp0(km-1)+dp0(km)) - - is1 = is - 1 - js1 = js - 1 - - ie1 = ie + 1 - je1 = je + 1 - - ie2 = ie + 2 - je2 = je + 2 - -!$OMP parallel do default(none) shared(js1,je1,is1,ie2,km,je2,ie1,ut,top_ratio,vt, & -!$OMP bot_ratio,dp0,js,je,ng,is,ie,gz,grid_type, & -!$OMP bd,npx,npy,sw_corner,se_corner,ne_corner, & -!$OMP nw_corner,area) & -!$OMP private(gz2, xfx, yfx, fx, fy, int_ratio) - do 6000 k=1,km+1 - - if ( k==1 ) then - do j=js1, je1 - do i=is1, ie2 - xfx(i,j) = ut(i,j,1) + (ut(i,j,1)-ut(i,j,2))*top_ratio - enddo - enddo - do j=js1, je2 - do i=is1, ie1 - yfx(i,j) = vt(i,j,1) + (vt(i,j,1)-vt(i,j,2))*top_ratio - enddo - enddo - elseif ( k==km+1 ) then -! Bottom extrapolation - do j=js1, je1 - do i=is1, ie2 - xfx(i,j) = ut(i,j,km) + (ut(i,j,km)-ut(i,j,km-1))*bot_ratio -! xfx(i,j) = r14*(3.*ut(i,j,km-2)-13.*ut(i,j,km-1)+24.*ut(i,j,km)) -! if ( xfx(i,j)*ut(i,j,km)<0. ) xfx(i,j) = 0. - enddo - enddo - do j=js1, je2 - do i=is1, ie1 - yfx(i,j) = vt(i,j,km) + (vt(i,j,km)-vt(i,j,km-1))*bot_ratio -! yfx(i,j) = r14*(3.*vt(i,j,km-2)-13.*vt(i,j,km-1)+24.*vt(i,j,km)) -! if ( yfx(i,j)*vt(i,j,km)<0. ) yfx(i,j) = 0. - enddo - enddo - else - int_ratio = 1./(dp0(k-1)+dp0(k)) - do j=js1, je1 - do i=is1, ie2 - xfx(i,j) = (dp0(k)*ut(i,j,k-1)+dp0(k-1)*ut(i,j,k))*int_ratio - enddo - enddo - do j=js1, je2 - do i=is1, ie1 - yfx(i,j) = (dp0(k)*vt(i,j,k-1)+dp0(k-1)*vt(i,j,k))*int_ratio - enddo - enddo - endif - - do j=js-ng, je+ng - do i=is-ng, ie+ng - gz2(i,j) = gz(i,j,k) - enddo - enddo - - if (grid_type < 3) call fill_4corners(gz2, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner) - do j=js1, je1 - do i=is1, ie2 - if( xfx(i,j) > 0. ) then - fx(i,j) = gz2(i-1,j) - else - fx(i,j) = gz2(i ,j) - endif - fx(i,j) = xfx(i,j)*fx(i,j) - enddo - enddo - - if (grid_type < 3) call fill_4corners(gz2, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner) - do j=js1,je2 - do i=is1,ie1 - if( yfx(i,j) > 0. ) then - fy(i,j) = gz2(i,j-1) - else - fy(i,j) = gz2(i,j) - endif - fy(i,j) = yfx(i,j)*fy(i,j) - enddo - enddo - - do j=js1, je1 - do i=is1,ie1 - gz(i,j,k) = (gz2(i,j)*area(i,j) + fx(i,j)- fx(i+1,j)+ fy(i,j)- fy(i,j+1)) & - / ( area(i,j) + xfx(i,j)-xfx(i+1,j)+yfx(i,j)-yfx(i,j+1)) - enddo - enddo -6000 continue - -! Enforce monotonicity of height to prevent blowup -!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km) - do j=js1, je1 - do k=2, km+1 - do i=is1, ie1 - gz(i,j,k) = min( gz(i,j,k), gz(i,j,k-1) - dz_min ) - enddo - enddo - do i=is1, ie1 - ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt - enddo - enddo - - end subroutine update_dz_c - - - subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea, & - dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac) - - type(fv_grid_bounds_type), intent(IN) :: bd - integer, intent(in):: is, ie, js, je, ng, km, npx, npy - integer, intent(in):: hord - real, intent(in) :: rdt - real, intent(in) :: dp0(km) - real, intent(in) :: area(is-ng:ie+ng,js-ng:je+ng) - real, intent(in) :: rarea(is-ng:ie+ng,js-ng:je+ng) - real, intent(inout):: damp(km+1) - integer, intent(inout):: ndif(km+1) - real, intent(in ) :: zs(is-ng:ie+ng,js-ng:je+ng) - real, intent(inout) :: zh(is-ng:ie+ng,js-ng:je+ng,km+1) - real, intent(inout), dimension(is:ie+1,js-ng:je+ng,km):: crx, xfx - real, intent(inout), dimension(is-ng:ie+ng,js:je+1,km):: cry, yfx - real, intent(out) :: ws(is:ie,js:je) - type(fv_grid_type), intent(IN), target :: gridstruct - real, intent(in) :: lim_fac -!----------------------------------------------------- -! Local array: - real, dimension(is: ie+1, js-ng:je+ng,km+1):: crx_adv, xfx_adv - real, dimension(is-ng:ie+ng,js: je+1,km+1 ):: cry_adv, yfx_adv - real, dimension(is:ie+1,js:je ):: fx - real, dimension(is:ie ,js:je+1):: fy - real, dimension(is-ng:ie+ng+1,js-ng:je+ng ):: fx2 - real, dimension(is-ng:ie+ng ,js-ng:je+ng+1):: fy2 - real, dimension(is-ng:ie+ng ,js-ng:je+ng ):: wk2, z2 - real:: ra_x(is:ie,js-ng:je+ng) - real:: ra_y(is-ng:ie+ng,js:je) -!-------------------------------------------------------------------- - integer i, j, k, isd, ied, jsd, jed - logical:: uniform_grid - - uniform_grid = .false. - - damp(km+1) = damp(km) - ndif(km+1) = ndif(km) - - isd = is - ng; ied = ie + ng - jsd = js - ng; jed = je + ng - -!$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, & -!$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv) - do j=jsd,jed - call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & - dp0, uniform_grid, 0) - if ( j<=je+1 .and. j>=js ) & - call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & - dp0, uniform_grid, 0) - enddo - -!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,area,xfx_adv,yfx_adv, & -!$OMP damp,zh,crx_adv,cry_adv,npx,npy,hord,gridstruct,bd, & -!$OMP ndif,rarea,lim_fac) & -!$OMP private(z2, fx2, fy2, ra_x, ra_y, fx, fy,wk2) - do k=1,km+1 - - do j=jsd,jed - do i=is,ie - ra_x(i,j) = area(i,j) + xfx_adv(i,j,k) - xfx_adv(i+1,j,k) - enddo - enddo - do j=js,je - do i=isd,ied - ra_y(i,j) = area(i,j) + yfx_adv(i,j,k) - yfx_adv(i,j+1,k) - enddo - enddo - - if ( damp(k)>1.E-5 ) then - do j=jsd,jed - do i=isd,ied - z2(i,j) = zh(i,j,k) - enddo - enddo - call fv_tp_2d(z2, crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, & - fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac) - call del6_vt_flux(ndif(k), npx, npy, damp(k), z2, wk2, fx2, fy2, gridstruct, bd) - do j=js,je - do i=is,ie - zh(i,j,k) = (z2(i,j)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & - / (ra_x(i,j)+ra_y(i,j)-area(i,j)) + (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*rarea(i,j) - enddo - enddo - else - call fv_tp_2d(zh(isd,jsd,k), crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, & - fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac) - do j=js,je - do i=is,ie - zh(i,j,k) = (zh(i,j,k)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & - / (ra_x(i,j) + ra_y(i,j) - area(i,j)) -! zh(i,j,k) = rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & -! + zh(i,j,k)*(3.-rarea(i,j)*(ra_x(i,j) + ra_y(i,j))) - enddo - enddo - endif - - enddo - -!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt) - do j=js, je - do k=2, km+1 - do i=is, ie -! Enforce monotonicity of height to prevent blowup - zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min ) - enddo - enddo - do i=is,ie - ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt - enddo - enddo - - end subroutine update_dz_d - - - subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & - akap, cappa, cp, & -#ifdef MULTI_GASES - kapad, & -#endif - ptop, hs, w3, pt, q_con, & - delp, gz, pef, ws, p_fac, a_imp, scale_m) - - integer, intent(in):: is, ie, js, je, ng, km - integer, intent(in):: ms - real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m - real, intent(in):: ws(is-ng:ie+ng,js-ng:je+ng) - real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delp - real, intent(in), dimension(is-ng:,js-ng:,1:):: q_con, cappa -#ifdef MULTI_GASES - real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: kapad -#endif - real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng) - real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: w3 -! OUTPUT PARAMETERS - real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: gz - real, intent( out), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: pef -! Local: - real, dimension(is-1:ie+1,km ):: dm, dz2, w2, pm2, gm2, cp2 - real, dimension(is-1:ie+1,km+1):: pem, pe2, peg -#ifdef MULTI_GASES - real, dimension(is-1:ie+1,km ):: kapad2 -#endif - real gama, rgrav - integer i, j, k - integer is1, ie1 - - gama = 1./(1.-akap) - rgrav = 1./grav - - is1 = is - 1 - ie1 = ie + 1 - -!$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, & -#ifdef MULTI_GASES -!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad) & -!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg, kapad2) -#else -!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) & -!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg) -#endif - do 2000 j=js-1, je+1 - - do k=1,km - do i=is1, ie1 - dm(i,k) = delp(i,j,k) - enddo - enddo - - do i=is1, ie1 - pef(i,j,1) = ptop ! full pressure at top - pem(i,1) = ptop -#ifdef USE_COND - peg(i,1) = ptop -#endif - enddo - - do k=2,km+1 - do i=is1, ie1 - pem(i,k) = pem(i,k-1) + dm(i,k-1) -#ifdef USE_COND -! Excluding contribution from condensates: - peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1)) -#endif - enddo - enddo - - do k=1,km - do i=is1, ie1 - dz2(i,k) = gz(i,j,k+1) - gz(i,j,k) -#ifdef USE_COND - pm2(i,k) = (peg(i,k+1)-peg(i,k))/log(peg(i,k+1)/peg(i,k)) - -#ifdef MOIST_CAPPA - cp2(i,k) = cappa(i,j,k) - gm2(i,k) = 1. / (1.-cp2(i,k)) -#endif - -#else - pm2(i,k) = dm(i,k)/log(pem(i,k+1)/pem(i,k)) -#endif -#ifdef MULTI_GASES - kapad2(i,k) = kapad(i,j,k) -#endif - dm(i,k) = dm(i,k) * rgrav - w2(i,k) = w3(i,j,k) - enddo - enddo - - - if ( a_imp < -0.01 ) then - call SIM3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, scale_m) - elseif ( a_imp <= 0.5 ) then - call RIM_2D(ms, dt, is1, ie1, km, rdgas, gama, gm2, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, & - dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.) - else - call SIM1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, & - dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac) - endif - - - do k=2,km+1 - do i=is1, ie1 - pef(i,j,k) = pe2(i,k) + pem(i,k) ! add hydrostatic full-component - enddo - enddo - -! Compute Height * grav (for p-gradient computation) - do i=is1, ie1 - gz(i,j,km+1) = hs(i,j) - enddo - - do k=km,1,-1 - do i=is1, ie1 - gz(i,j,k) = gz(i,j,k+1) - dz2(i,k)*grav - enddo - enddo - -2000 continue - - end subroutine Riem_Solver_c - - -!>GFDL - This routine will not give absoulte reproducibility when compiled with -fast-transcendentals. -!! GFDL - It is now inside of nh_core.F90 and being compiled without -fast-transcendentals. - subroutine Riem_Solver3test(ms, dt, is, ie, js, je, km, ng, & - isd, ied, jsd, jed, akap, cappa, cp, & -#ifdef MULTI_GASES - kapad, & -#endif - ptop, zs, q_con, w, delz, pt, & - delp, zh, pe, ppe, pk3, pk, peln, & - ws, scale_m, p_fac, a_imp, & - use_logp, last_call, fp_out) -!-------------------------------------------- -! !OUTPUT PARAMETERS -! Ouput: gz: grav*height at edges -! pe: full hydrostatic pressure -! ppe: non-hydrostatic pressure perturbation -!-------------------------------------------- - integer, intent(in):: ms, is, ie, js, je, km, ng - integer, intent(in):: isd, ied, jsd, jed - real, intent(in):: dt ! the BIG horizontal Lagrangian time step - real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m - real, intent(in):: zs(isd:ied,jsd:jed) - logical, intent(in):: last_call, use_logp, fp_out - real, intent(in):: ws(is:ie,js:je) - real, intent(in), dimension(isd:,jsd:,1:):: q_con, cappa - real, intent(in), dimension(isd:ied,jsd:jed,km):: delp, pt -#ifdef MULTI_GASES - real, intent(in), dimension(isd:ied,jsd:jed,km):: kapad -#endif - real, intent(inout), dimension(isd:ied,jsd:jed,km+1):: zh - real, intent(inout), dimension(isd:ied,jsd:jed,km):: w - real, intent(inout):: pe(is-1:ie+1,km+1,js-1:je+1) - real, intent(out):: peln(is:ie,km+1,js:je) ! ln(pe) - real, intent(out), dimension(isd:ied,jsd:jed,km+1):: ppe - real, intent(out):: delz(is:ie,js:je,km) - real, intent(out):: pk(is:ie,js:je,km+1) - real, intent(out):: pk3(isd:ied,jsd:jed,km+1) -! Local: - real, dimension(is:ie,km):: dm, dz2, pm2, w2, gm2, cp2 - real, dimension(is:ie,km+1)::pem, pe2, peln2, peg, pelng -#ifdef MULTI_GASES - real, dimension(is:ie,km):: kapad2 -#endif - real gama, rgrav, ptk, peln1 - integer i, j, k - - gama = 1./(1.-akap) - rgrav = 1./grav - peln1 = log(ptop) - ptk = exp(akap*peln1) - -!$OMP parallel do default(none) shared(is,ie,js,je,km,delp,ptop,peln1,pk3,ptk,akap,rgrav,zh,pt, & -!$OMP w,a_imp,dt,gama,ws,p_fac,scale_m,ms,delz,last_call, & -#ifdef MULTI_GASES -!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad ) & -!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2,kapad2) -#else -!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con ) & -!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2) -#endif - do 2000 j=js, je - - do k=1,km - do i=is, ie - dm(i,k) = delp(i,j,k) -#ifdef MOIST_CAPPA - cp2(i,k) = cappa(i,j,k) -#endif -#ifdef MULTI_GASES - kapad2(i,k) = kapad(i,j,k) -#endif - enddo - enddo - - do i=is,ie - pem(i,1) = ptop - peln2(i,1) = peln1 - pk3(i,j,1) = ptk -#ifdef USE_COND - peg(i,1) = ptop - pelng(i,1) = peln1 -#endif - enddo - do k=2,km+1 - do i=is,ie - pem(i,k) = pem(i,k-1) + dm(i,k-1) - peln2(i,k) = log(pem(i,k)) -#ifdef USE_COND -! Excluding contribution from condensates: -! peln used during remap; pk3 used only for p_grad - peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1)) - pelng(i,k) = log(peg(i,k)) -#endif - pk3(i,j,k) = exp(akap*peln2(i,k)) - enddo - enddo - - do k=1,km - do i=is, ie -#ifdef USE_COND - pm2(i,k) = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k)) - -#ifdef MOIST_CAPPA - gm2(i,k) = 1. / (1.-cp2(i,k)) -#endif - -#else - pm2(i,k) = dm(i,k)/(peln2(i,k+1)-peln2(i,k)) -#endif - dm(i,k) = dm(i,k) * rgrav - dz2(i,k) = zh(i,j,k+1) - zh(i,j,k) - w2(i,k) = w(i,j,k) - enddo - enddo - - - if ( a_imp < -0.999 ) then - call SIM3p0_solver(dt, is, ie, km, rdgas, gama, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, scale_m ) - elseif ( a_imp < -0.5 ) then - call SIM3_solver(dt, is, ie, km, rdgas, gama, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), abs(a_imp), p_fac, scale_m) - elseif ( a_imp <= 0.5 ) then - call RIM_2D(ms, dt, is, ie, km, rdgas, gama, gm2, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, & - dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.) - elseif ( a_imp > 0.999 ) then - call SIM1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac) - else - call SIM_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), & - a_imp, p_fac, scale_m) - endif - - - do k=1, km - do i=is, ie - w(i,j,k) = w2(i,k) - delz(i,j,k) = dz2(i,k) - enddo - enddo - - if ( last_call ) then - do k=1,km+1 - do i=is,ie - peln(i,k,j) = peln2(i,k) - pk(i,j,k) = pk3(i,j,k) - pe(i,k,j) = pem(i,k) - enddo - enddo - endif - - if( fp_out ) then - do k=1,km+1 - do i=is, ie - ppe(i,j,k) = pe2(i,k) + pem(i,k) - enddo - enddo - else - do k=1,km+1 - do i=is, ie - ppe(i,j,k) = pe2(i,k) - enddo - enddo - endif - - if ( use_logp ) then - do k=2,km+1 - do i=is, ie - pk3(i,j,k) = peln2(i,k) - enddo - enddo - endif - - do i=is, ie - zh(i,j,km+1) = zs(i,j) - enddo - do k=km,1,-1 - do i=is, ie - zh(i,j,k) = zh(i,j,k+1) - dz2(i,k) - enddo - enddo - -2000 continue - - end subroutine Riem_Solver3test - - subroutine imp_diff_w(j, is, ie, js, je, ng, km, cd, delz, ws, w, w3) - integer, intent(in) :: j, is, ie, js, je, km, ng - real, intent(in) :: cd - real, intent(in) :: delz(is:ie, km) !< delta-height (m) - real, intent(in) :: w(is:ie, km) !< vertical vel. (m/s) - real, intent(in) :: ws(is:ie) - real, intent(out) :: w3(is-ng:ie+ng,js-ng:je+ng,km) -! Local: - real, dimension(is:ie,km):: c, gam, dz, wt - real:: bet(is:ie) - real:: a - integer:: i, k - - do k=2,km - do i=is,ie - dz(i,k) = 0.5*(delz(i,k-1)+delz(i,k)) - enddo - enddo - do k=1,km-1 - do i=is,ie - c(i,k) = -cd/(dz(i,k+1)*delz(i,k)) - enddo - enddo - -! model top: - do i=is,ie - bet(i) = 1. - c(i,1) ! bet(i) = b - wt(i,1) = w(i,1) / bet(i) - enddo - -! Interior: - do k=2,km-1 - do i=is,ie - gam(i,k) = c(i,k-1)/bet(i) - a = cd/(dz(i,k)*delz(i,k)) - bet(i) = (1.+a-c(i,k)) + a*gam(i,k) - wt(i,k) = (w(i,k) + a*wt(i,k-1)) / bet(i) - enddo - enddo - -! Bottom: - do i=is,ie - gam(i,km) = c(i,km-1) / bet(i) - a = cd/(dz(i,km)*delz(i,km)) - wt(i,km) = (w(i,km) + 2.*ws(i)*cd/delz(i,km)**2 & - + a*wt(i,km-1))/(1. + a + (cd+cd)/delz(i,km)**2 + a*gam(i,km)) - enddo - - do k=km-1,1,-1 - do i=is,ie - wt(i,k) = wt(i,k) - gam(i,k+1)*wt(i,k+1) - enddo - enddo - - do k=1,km - do i=is,ie - w3(i,j,k) = wt(i,k) - enddo - enddo - - end subroutine imp_diff_w - - - subroutine RIM_2D(ms, bdt, is, ie, km, rgas, gama, gm2, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm2, pm2, w2, dz2, pt2, ws, c_core ) - - integer, intent(in):: ms, is, ie, km - real, intent(in):: bdt, gama, rgas - real, intent(in), dimension(is:ie,km):: dm2, pm2, gm2 - logical, intent(in):: c_core - real, intent(in ) :: pt2(is:ie,km) - real, intent(in ) :: ws(is:ie) -! IN/OUT: - real, intent(inout):: dz2(is:ie,km) - real, intent(inout):: w2(is:ie,km) - real, intent(out ):: pe2(is:ie,km+1) -#ifdef MULTI_GASES - real, intent(inout), dimension(is:ie,km):: kapad2 -#endif -! Local: - real:: ws2(is:ie) - real, dimension(km+1):: m_bot, m_top, r_bot, r_top, pe1, pbar, wbar - real, dimension(km):: r_hi, r_lo, dz, wm, dm, dts - real, dimension(km):: pf1, wc, cm , pp, pt1 - real:: dt, rdt, grg, z_frac, ptmp1, rden, pf, time_left - real:: m_surf -#ifdef MULTI_GASES - real gamax -#endif - integer:: i, k, n, ke, kt1, ktop - integer:: ks0, ks1 - - grg = gama * rgas - rdt = 1. / bdt - dt = bdt / real(ms) - - pbar(:) = 0. - wbar(:) = 0. - - do i=is,ie - ws2(i) = 2.*ws(i) - enddo - - do 6000 i=is,ie - - do k=1,km - dz(k) = dz2(i,k) - dm(k) = dm2(i,k) - wm(k) = w2(i,k)*dm(k) - pt1(k) = pt2(i,k) - enddo - - pe1(:) = 0. - wbar(km+1) = ws(i) - - ks0 = 1 - if ( ms > 1 .and. ms < 8 ) then -! Continuity of (pbar, wbar) is maintained - do k=1, km - rden = -rgas*dm(k)/dz(k) -#ifdef MOIST_CAPPA - pf1(k) = exp( gm2(i,k)*log(rden*pt1(k)) ) -! dts(k) = -dz(k)/sqrt(gm2(i,k)*rgas*pf1(k)/rden) - dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - pf1(k) = exp( gamax*log(rden*pt1(k)) ) -#else - pf1(k) = exp( gama*log(rden*pt1(k)) ) -#endif - dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden) -#endif - if ( bdt > dts(k) ) then - ks0 = k-1 - goto 222 - endif - enddo - ks0 = km -222 if ( ks0==1 ) goto 244 - - do k=1, ks0 - cm(k) = dm(k) / dts(k) - wc(k) = wm(k) / dts(k) - pp(k) = pf1(k) - pm2(i,k) - enddo - - wbar(1) = (wc(1)+pp(1)) / cm(1) - do k=2, ks0 - wbar(k) = (wc(k-1)+wc(k) + pp(k)-pp(k-1)) / (cm(k-1)+cm(k)) - pbar(k) = bdt*(cm(k-1)*wbar(k) - wc(k-1) + pp(k-1)) - pe1(k) = pbar(k) - enddo - - if ( ks0 == km ) then - pbar(km+1) = bdt*(cm(km)*wbar(km+1) - wc(km) + pp(km)) - if ( c_core ) then - do k=1,km - dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) - enddo - else - do k=1,km - dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) - w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k) - enddo - endif - pe2(i,1) = 0. - do k=2,km+1 - pe2(i,k) = pbar(k)*rdt - enddo - goto 6000 ! next i - else - if ( c_core ) then - do k=1, ks0-1 - dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) - enddo - else - do k=1, ks0-1 - dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) - w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k) - enddo - endif - pbar(ks0) = pbar(ks0) / real(ms) - endif - endif -244 ks1 = ks0 - - do 5000 n=1, ms - - do k=ks1, km - rden = -rgas*dm(k)/dz(k) -#ifdef MOIST_CAPPA - pf = exp( gm2(i,k)*log(rden*pt1(k)) ) -! dts(k) = -dz(k) / sqrt( gm2(i,k)*rgas*pf/rden ) - dts(k) = -dz(k) / sqrt( grg*pf/rden ) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - pf = exp( gamax*log(rden*pt1(k)) ) -#else - pf = exp( gama*log(rden*pt1(k)) ) -#endif - dts(k) = -dz(k) / sqrt( grg*pf/rden ) -#endif - ptmp1 = dts(k)*(pf - pm2(i,k)) - r_lo(k) = wm(k) + ptmp1 - r_hi(k) = wm(k) - ptmp1 - enddo - - ktop = ks1 - do k=ks1, km - if( dt > dts(k) ) then - ktop = k-1 - goto 333 - endif - enddo - ktop = km -333 continue - - if ( ktop >= ks1 ) then - do k=ks1, ktop - z_frac = dt/dts(k) - r_bot(k ) = z_frac*r_lo(k) - r_top(k+1) = z_frac*r_hi(k) - m_bot(k ) = z_frac* dm(k) - m_top(k+1) = m_bot(k) - enddo - if ( ktop == km ) goto 666 - endif - - do k=ktop+2, km+1 - m_top(k) = 0. - r_top(k) = 0. - enddo - - kt1 = max(1, ktop) - do 444 ke=km+1, ktop+2, -1 - time_left = dt - do k=ke-1, kt1, -1 - if ( time_left > dts(k) ) then - time_left = time_left - dts(k) - m_top(ke) = m_top(ke) + dm(k) - r_top(ke) = r_top(ke) + r_hi(k) - else - z_frac = time_left/dts(k) - m_top(ke) = m_top(ke) + z_frac*dm(k) - r_top(ke) = r_top(ke) + z_frac*r_hi(k) - go to 444 ! next level - endif - enddo -444 continue - - do k=ktop+1, km - m_bot(k) = 0. - r_bot(k) = 0. - enddo - - do 4000 ke=ktop+1, km - time_left = dt - do k=ke, km - if ( time_left > dts(k) ) then - time_left = time_left - dts(k) - m_bot(ke) = m_bot(ke) + dm(k) - r_bot(ke) = r_bot(ke) + r_lo(k) - else - z_frac = time_left/dts(k) - m_bot(ke) = m_bot(ke) + z_frac* dm(k) - r_bot(ke) = r_bot(ke) + z_frac*r_lo(k) - go to 4000 ! next interface - endif - enddo - m_surf = m_bot(ke) - do k=km, kt1, -1 - if ( time_left > dts(k) ) then - time_left = time_left - dts(k) - m_bot(ke) = m_bot(ke) + dm(k) - r_bot(ke) = r_bot(ke) - r_hi(k) - else - z_frac = time_left/dts(k) - m_bot(ke) = m_bot(ke) + z_frac* dm(k) - r_bot(ke) = r_bot(ke) - z_frac*r_hi(k) + (m_bot(ke)-m_surf)*ws2(i) - go to 4000 ! next interface - endif - enddo -4000 continue - -666 if ( ks1==1 ) wbar(1) = r_bot(1) / m_bot(1) - do k=ks1+1, km - wbar(k) = (r_bot(k)+r_top(k)) / (m_top(k)+m_bot(k)) - enddo -! pbar here is actually dt*pbar - do k=ks1+1, km+1 - pbar(k) = m_top(k)*wbar(k) - r_top(k) - pe1(k) = pe1(k) + pbar(k) - enddo - - if ( n==ms ) then - if ( c_core ) then - do k=ks1, km - dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k)) - enddo - else - do k=ks1, km - dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k)) - w2(i,k) = ( wm(k) + pbar(k+1) - pbar(k) ) / dm(k) - enddo - endif - else - do k=ks1, km - dz(k) = dz(k) + dt*(wbar(k+1)-wbar(k)) - wm(k) = wm(k) + pbar(k+1) - pbar(k) - enddo - endif - -5000 continue - pe2(i,1) = 0. - do k=2,km+1 - pe2(i,k) = pe1(k)*rdt - enddo - -6000 continue ! end i-loop - - end subroutine RIM_2D - - subroutine SIM3_solver(dt, is, ie, km, rgas, gama, kappa, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m) - integer, intent(in):: is, ie, km - real, intent(in):: dt, rgas, gama, kappa, alpha, p_fac, scale_m - real, intent(in ), dimension(is:ie,km):: dm, pt2 - real, intent(in ):: ws(is:ie) - real, intent(in ), dimension(is:ie,km+1):: pem - real, intent(out):: pe2(is:ie,km+1) - real, intent(inout), dimension(is:ie,km):: dz2, w2 -#ifdef MULTI_GASES - real, intent(inout), dimension(is:ie,km):: kapad2 -#endif -! Local - real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam - real, dimension(is:ie,km+1):: pp - real, dimension(is:ie):: p1, wk1, bet - real beta, t2, t1g, rdt, ra, capa1, r2g, r6g -#ifdef MULTI_GASES - real gamax, capa1x, t1gx -#endif - integer i, k - - beta = 1. - alpha - ra = 1. / alpha - t2 = beta / alpha - t1g = gama * 2.*(alpha*dt)**2 - rdt = 1. / dt - capa1 = kappa - 1. - r2g = grav / 2. - r6g = grav / 6. - - - do k=1,km - do i=is, ie - w1(i,k) = w2(i,k) -! Full pressure at center -#ifdef MULTI_GASES - gamax = 1. / (1.-kapad2(i,k)) - aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) -#else - aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) -#endif - enddo - enddo - - do k=1,km-1 - do i=is, ie - g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction - bb(i,k) = 2.*(1.+g_rat(i,k)) - dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1)) - enddo - enddo - -! pe2 is full p at edges - do i=is, ie -! Top: - bet(i) = bb(i,1) - pe2(i,1) = pem(i,1) - pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i) -! Bottom: - bb(i,km) = 2. - dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km) - enddo - - do k=2,km - do i=is, ie - gam(i,k) = g_rat(i,k-1) / bet(i) - bet(i) = bb(i,k) - gam(i,k) - pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i) - enddo - enddo - - do k=km, 2, -1 - do i=is, ie - pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1) - enddo - enddo -! done reconstruction of full: - -! pp is pert. p at edges - do k=1, km+1 - do i=is, ie - pp(i,k) = pe2(i,k) - pem(i,k) - enddo - enddo - - do k=2, km - do i=is, ie -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - t1gx = gamax*2.*(alpha*dt)**2 - aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) -#else - aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) -#endif - wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k)) - aa(i,k) = aa(i,k) - scale_m*dm(i,1) - enddo - enddo - do i=is, ie - bet(i) = dm(i,1) - aa(i,2) - w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2) + wk(i,2)) / bet(i) - enddo - do k=2,km-1 - do i=is, ie - gam(i,k) = aa(i,k) / bet(i) - bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) - w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) & - - aa(i,k)*w2(i,k-1)) / bet(i) - enddo - enddo - do i=is, ie -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,km)) - t1gx = gamax*2.*(alpha*dt)**2 - wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) -#else - wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) -#endif - gam(i,km) = aa(i,km) / bet(i) - bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) - w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk(i,km) + & - wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i) - enddo - do k=km-1, 1, -1 - do i=is, ie - w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) - enddo - enddo - -! pe2 is updated perturbation p at edges - do i=is, ie - pe2(i,1) = 0. - enddo - do k=1,km - do i=is, ie - pe2(i,k+1) = pe2(i,k) + ( dm(i,k)*(w2(i,k)-w1(i,k))*rdt & - - beta*(pp(i,k+1)-pp(i,k)) )*ra - enddo - enddo - -! Full non-hydro pressure at edges: - do i=is, ie - pe2(i,1) = pem(i,1) - enddo - do k=2,km+1 - do i=is, ie - pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k)) - enddo - enddo - - do i=is, ie -! Recover cell-averaged pressure - p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km) -#ifdef MULTI_GASES - capa1x = kapad2(i,km) - 1. - dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) ) -#else - dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) ) -#endif - enddo - - do k=km-1, 1, -1 - do i=is, ie - p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i) -#ifdef MULTI_GASES - capa1x = kapad2(i,k) - 1. - dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) ) -#else - dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) ) -#endif - enddo - enddo - - do k=1,km+1 - do i=is, ie - pe2(i,k) = pe2(i,k) - pem(i,k) - pe2(i,k) = pe2(i,k) + beta*(pp(i,k) - pe2(i,k)) - enddo - enddo - - end subroutine SIM3_solver - - subroutine SIM3p0_solver(dt, is, ie, km, rgas, gama, kappa, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm, & - pem, w2, dz2, pt2, ws, p_fac, scale_m) -! Sa SIM3, but for beta==0 - integer, intent(in):: is, ie, km - real, intent(in):: dt, rgas, gama, kappa, p_fac, scale_m - real, intent(in ), dimension(is:ie,km):: dm, pt2 - real, intent(in ):: ws(is:ie) - real, intent(in ):: pem(is:ie,km+1) - real, intent(out):: pe2(is:ie,km+1) - real, intent(inout), dimension(is:ie,km):: dz2, w2 -#ifdef MULTI_GASES - real, intent(inout), dimension(is:ie,km):: kapad2 -#endif -! Local - real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam - real, dimension(is:ie,km+1):: pp - real, dimension(is:ie):: p1, wk1, bet - real t1g, rdt, capa1, r2g, r6g -#ifdef MULTI_GASES - real gamax, capa1x, t1gx -#endif - integer i, k - - t1g = 2.*gama*dt**2 - rdt = 1. / dt - capa1 = kappa - 1. - r2g = grav / 2. - r6g = grav / 6. - - do k=1,km - do i=is, ie - w1(i,k) = w2(i,k) -! Full pressure at center -#ifdef MULTI_GASES - gamax = 1. / ( 1. - kapad2(i,k) ) - aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) -#else - aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) -#endif - enddo - enddo - - do k=1,km-1 - do i=is, ie - g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction - bb(i,k) = 2.*(1.+g_rat(i,k)) - dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1)) - enddo - enddo - -! pe2 is full p at edges - do i=is, ie -! Top: - bet(i) = bb(i,1) - pe2(i,1) = pem(i,1) - pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i) -! Bottom: - bb(i,km) = 2. - dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km) - enddo - - do k=2,km - do i=is, ie - gam(i,k) = g_rat(i,k-1) / bet(i) - bet(i) = bb(i,k) - gam(i,k) - pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i) - enddo - enddo - - do k=km, 2, -1 - do i=is, ie - pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1) - enddo - enddo -! done reconstruction of full: - -! pp is pert. p at edges - do k=1, km+1 - do i=is, ie - pp(i,k) = pe2(i,k) - pem(i,k) - enddo - enddo - - do k=2, km - do i=is, ie -#ifdef MULTI_GASES - gamax = 1. / (1.-kapad2(i,k)) - t1gx = 2.*gamax*dt**2 - aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1) -#else - aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1) -#endif - enddo - enddo - do i=is, ie - bet(i) = dm(i,1) - aa(i,2) - w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2)) / bet(i) - enddo - do k=2,km-1 - do i=is, ie - gam(i,k) = aa(i,k) / bet(i) - bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) - w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1))/bet(i) - enddo - enddo - do i=is, ie -#ifdef MULTI_GASES - gamax = 1. / (1.-kapad2(i,km)) - t1gx = 2.*gamax*dt**2 - wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) -#else - wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) -#endif - gam(i,km) = aa(i,km) / bet(i) - bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) - w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk1(i)*ws(i) - & - aa(i,km)*w2(i,km-1)) / bet(i) - enddo - do k=km-1, 1, -1 - do i=is, ie - w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) - enddo - enddo - -! pe2 is updated perturbation p at edges - do i=is, ie - pe2(i,1) = 0. - enddo - do k=1,km - do i=is, ie - pe2(i,k+1) = pe2(i,k) + dm(i,k)*(w2(i,k)-w1(i,k))*rdt - enddo - enddo - -! Full non-hydro pressure at edges: - do i=is, ie - pe2(i,1) = pem(i,1) - enddo - do k=2,km+1 - do i=is, ie - pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k)) - enddo - enddo - - do i=is, ie -! Recover cell-averaged pressure - p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km) -#ifdef MULTI_GASES - capa1x = kapad2(i,km) - 1. - dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) ) -#else - dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) ) -#endif - enddo - - do k=km-1, 1, -1 - do i=is, ie - p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3-g_rat(i,k)*p1(i) -#ifdef MULTI_GASES - capa1x = kapad2(i,k) - 1. - dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) ) -#else - dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) ) -#endif - enddo - enddo - - do k=1,km+1 - do i=is, ie - pe2(i,k) = pe2(i,k) - pem(i,k) - enddo - enddo - - end subroutine SIM3p0_solver - - - subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe, dm2, & - pm2, pem, w2, dz2, pt2, ws, p_fac) - integer, intent(in):: is, ie, km - real, intent(in):: dt, rgas, gama, kappa, p_fac - real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2 - real, intent(in ):: ws(is:ie) - real, intent(in ), dimension(is:ie,km+1):: pem - real, intent(out):: pe(is:ie,km+1) - real, intent(inout), dimension(is:ie,km):: dz2, w2 -#ifdef MULTI_GASES - real, intent(inout), dimension(is:ie,km):: kapad2 -#endif -! Local - real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam - real, dimension(is:ie,km+1):: pp - real, dimension(is:ie):: p1, bet - real t1g, rdt, capa1 -#ifdef MULTI_GASES - real gamax, capa1x, t1gx -#endif - integer i, k - -#ifdef MOIST_CAPPA - t1g = 2.*dt*dt -#else - t1g = gama * 2.*dt*dt -#endif - rdt = 1. / dt - capa1 = kappa - 1. - - do k=1,km - do i=is, ie -#ifdef MOIST_CAPPA - pe(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#else -#ifdef MULTI_GASES - gamax = 1. / ( 1. - kapad2(i,k) ) - pe(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#else - pe(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#endif -#endif - w1(i,k) = w2(i,k) - enddo - enddo - - do k=1,km-1 - do i=is, ie - g_rat(i,k) = dm2(i,k)/dm2(i,k+1) - bb(i,k) = 2.*(1.+g_rat(i,k)) - dd(i,k) = 3.*(pe(i,k) + g_rat(i,k)*pe(i,k+1)) - enddo - enddo - - do i=is, ie - bet(i) = bb(i,1) - pp(i,1) = 0. - pp(i,2) = dd(i,1) / bet(i) - bb(i,km) = 2. - dd(i,km) = 3.*pe(i,km) - enddo - - do k=2,km - do i=is, ie - gam(i,k) = g_rat(i,k-1) / bet(i) - bet(i) = bb(i,k) - gam(i,k) - pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i) - enddo - enddo - - do k=km, 2, -1 - do i=is, ie - pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1) - enddo - enddo - -! Start the w-solver - do k=2, km - do i=is, ie -#ifdef MOIST_CAPPA - aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - t1gx = gamax * 2.*dt*dt - aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) -#else - aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) -#endif -#endif - enddo - enddo - do i=is, ie - bet(i) = dm2(i,1) - aa(i,2) - w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2)) / bet(i) - enddo - do k=2,km-1 - do i=is, ie - gam(i,k) = aa(i,k) / bet(i) - bet(i) = dm2(i,k) - (aa(i,k) + aa(i,k+1) + aa(i,k)*gam(i,k)) - w2(i,k) = (dm2(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1)) / bet(i) - enddo - enddo - do i=is, ie -#ifdef MOIST_CAPPA - p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,km)) - t1gx = gamax * 2.*dt*dt - p1(i) = t1gx/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) -#else - p1(i) = t1g/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) -#endif -#endif - gam(i,km) = aa(i,km) / bet(i) - bet(i) = dm2(i,km) - (aa(i,km)+p1(i) + aa(i,km)*gam(i,km)) - w2(i,km) = (dm2(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-p1(i)*ws(i)-aa(i,km)*w2(i,km-1))/bet(i) - enddo - do k=km-1, 1, -1 - do i=is, ie - w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) - enddo - enddo - - do i=is, ie - pe(i,1) = 0. - enddo - do k=1,km - do i=is, ie - pe(i,k+1) = pe(i,k) + dm2(i,k)*(w2(i,k)-w1(i,k))*rdt - enddo - enddo - - do i=is, ie - p1(i) = ( pe(i,km) + 2.*pe(i,km+1) )*r3 -#ifdef MOIST_CAPPA - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#else -#ifdef MULTI_GASES - capa1x = kapad2(i,km)-1. - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#else - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#endif -#endif - enddo - - do k=km-1, 1, -1 - do i=is, ie - p1(i) = (pe(i,k) + bb(i,k)*pe(i,k+1) + g_rat(i,k)*pe(i,k+2))*r3 - g_rat(i,k)*p1(i) -#ifdef MOIST_CAPPA - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) - -#else -#ifdef MULTI_GASES - capa1x = kapad2(i,k)-1. - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) -#else - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) -#endif -#endif - enddo - enddo - - end subroutine SIM1_solver - - subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, & -#ifdef MULTI_GASES - kapad2, & -#endif - pe2, dm2, & - pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m) - integer, intent(in):: is, ie, km - real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m - real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2 - real, intent(in ):: ws(is:ie) - real, intent(in ), dimension(is:ie,km+1):: pem - real, intent(out):: pe2(is:ie,km+1) - real, intent(inout), dimension(is:ie,km):: dz2, w2 -#ifdef MULTI_GASES - real, intent(inout), dimension(is:ie,km):: kapad2 -#endif -! Local - real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam - real, dimension(is:ie,km+1):: pp - real, dimension(is:ie):: p1, wk1, bet - real beta, t2, t1g, rdt, ra, capa1 -#ifdef MULTI_GASES - real gamax, capa1x, t1gx -#endif - integer i, k - - beta = 1. - alpha - ra = 1. / alpha - t2 = beta / alpha -#ifdef MOIST_CAPPA - t1g = 2.*(alpha*dt)**2 -#else - t1g = 2.*gama*(alpha*dt)**2 -#endif - rdt = 1. / dt - capa1 = kappa - 1. - - do k=1,km - do i=is, ie - w1(i,k) = w2(i,k) -! P_g perturbation -#ifdef MOIST_CAPPA - pe2(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - pe2(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#else - pe2(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) -#endif -#endif - enddo - enddo - - do k=1,km-1 - do i=is, ie - g_rat(i,k) = dm2(i,k)/dm2(i,k+1) - bb(i,k) = 2.*(1.+g_rat(i,k)) - dd(i,k) = 3.*(pe2(i,k) + g_rat(i,k)*pe2(i,k+1)) - enddo - enddo - - do i=is, ie - bet(i) = bb(i,1) - pp(i,1) = 0. - pp(i,2) = dd(i,1) / bet(i) - bb(i,km) = 2. - dd(i,km) = 3.*pe2(i,km) - enddo - - do k=2,km - do i=is, ie - gam(i,k) = g_rat(i,k-1) / bet(i) - bet(i) = bb(i,k) - gam(i,k) - pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i) - enddo - enddo - - do k=km, 2, -1 - do i=is, ie - pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1) - enddo - enddo - - do k=1, km+1 - do i=is, ie -! pe2 is Full p - pe2(i,k) = pem(i,k) + pp(i,k) - enddo - enddo - - do k=2, km - do i=is, ie -#ifdef MOIST_CAPPA - aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,k)) - t1gx = 2.*gamax*(alpha*dt)**2 - aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) -#else - aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) -#endif -#endif - wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k)) - aa(i,k) = aa(i,k) - scale_m*dm2(i,1) - enddo - enddo -! Top: - do i=is, ie - bet(i) = dm2(i,1) - aa(i,2) - w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2) + wk(i,2)) / bet(i) - enddo -! Interior: - do k=2,km-1 - do i=is, ie - gam(i,k) = aa(i,k) / bet(i) - bet(i) = dm2(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) - w2(i,k) = (dm2(i,k)*w1(i,k) + dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) & - - aa(i,k)*w2(i,k-1)) / bet(i) - enddo - enddo -! Bottom: k=km - do i=is, ie -#ifdef MOIST_CAPPA - wk1(i) = t1g*gm2(i,km)/dz2(i,km)*pe2(i,km+1) -#else -#ifdef MULTI_GASES - gamax = 1./(1.-kapad2(i,km)) - t1gx = 2.*gamax*(alpha*dt)**2 - wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) -#else - wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) -#endif -#endif - gam(i,km) = aa(i,km) / bet(i) - bet(i) = dm2(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) - w2(i,km) = (dm2(i,km)*w1(i,km) + dt*(pp(i,km+1)-pp(i,km)) - wk(i,km) + & - wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i) - enddo - do k=km-1, 1, -1 - do i=is, ie - w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) - enddo - enddo - - do i=is, ie - pe2(i,1) = 0. - enddo - do k=1,km - do i=is, ie - pe2(i,k+1) = pe2(i,k) + ( dm2(i,k)*(w2(i,k)-w1(i,k))*rdt & - - beta*(pp(i,k+1)-pp(i,k)) ) * ra - enddo - enddo - - do i=is, ie - p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 -#ifdef MOIST_CAPPA - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#else -#ifdef MULTI_GASES - capa1x = kapad2(i,km)-1. - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#else - dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) -#endif -#endif - enddo - - do k=km-1, 1, -1 - do i=is, ie - p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i) -! delz = -dm*R*T_m / p_gas -#ifdef MOIST_CAPPA - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) -#else -#ifdef MULTI_GASES - capa1x = kapad2(i,k)-1. - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) -#else - dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) -#endif -#endif - enddo - enddo - - do k=1, km+1 - do i=is, ie - pe2(i,k) = pe2(i,k) + beta*(pp(i,k)-pe2(i,k)) - enddo - enddo - - end subroutine SIM_solver - - - subroutine edge_scalar(q1, qe, i1, i2, km, id) -! Optimized for wind profile reconstruction: - integer, intent(in):: i1, i2, km - integer, intent(in):: id ! 0: pp 1: wind - real, intent(in ), dimension(i1:i2,km):: q1 - real, intent(out), dimension(i1:i2,km+1):: qe -!----------------------------------------------------------------------- - real, parameter:: r2o3 = 2./3. - real, parameter:: r4o3 = 4./3. - real gak(km) - real bet - integer i, k - -!------------------------------------------------ -! Optimized coding for uniform grid: SJL Apr 2007 -!------------------------------------------------ - - if ( id==1 ) then - do i=i1,i2 - qe(i,1) = r4o3*q1(i,1) + r2o3*q1(i,2) - enddo - else - do i=i1,i2 - qe(i,1) = 1.E30 - enddo - endif - - gak(1) = 7./3. - do k=2,km - gak(k) = 1. / (4. - gak(k-1)) - do i=i1,i2 - qe(i,k) = (3.*(q1(i,k-1) + q1(i,k)) - qe(i,k-1)) * gak(k) - enddo - enddo - - bet = 1. / (1.5 - 3.5*gak(km)) - do i=i1,i2 - qe(i,km+1) = (4.*q1(i,km) + q1(i,km-1) - 3.5*qe(i,km)) * bet - enddo - - do k=km,1,-1 - do i=i1,i2 - qe(i,k) = qe(i,k) - gak(k)*qe(i,k+1) - enddo - enddo - - end subroutine edge_scalar - - - - subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) -! Optimized for wind profile reconstruction: - integer, intent(in):: i1, i2, j1, j2 - integer, intent(in):: j, km - integer, intent(in):: limiter - logical, intent(in):: uniform_grid - real, intent(in):: dp0(km) - real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2 - real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e -!----------------------------------------------------------------------- - real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values - real gak(km) - real bet, r2o3, r4o3 - real g0, gk, xt1, xt2, a_bot - integer i, k - - if ( uniform_grid ) then -!------------------------------------------------ -! Optimized coding for uniform grid: SJL Apr 2007 -!------------------------------------------------ - r2o3 = 2./3. - r4o3 = 4./3. - do i=i1,i2 - qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2) - qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2) - enddo - - gak(1) = 7./3. - do k=2,km - gak(k) = 1. / (4. - gak(k-1)) - do i=i1,i2 - qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k) - qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k) - enddo - enddo - - bet = 1. / (1.5 - 3.5*gak(km)) - do i=i1,i2 - qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet - qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet - enddo - - do k=km,1,-1 - do i=i1,i2 - qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1) - qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1) - enddo - enddo - else -! Assuming grid varying in vertical only - g0 = dp0(2) / dp0(1) - xt1 = 2.*g0*(g0+1. ) - bet = g0*(g0+0.5) - do i=i1,i2 - qe1(i,1) = ( xt1*q1(i,j,1) + q1(i,j,2) ) / bet - qe2(i,1) = ( xt1*q2(i,j,1) + q2(i,j,2) ) / bet - gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet - enddo - - do k=2,km - gk = dp0(k-1) / dp0(k) - do i=i1,i2 - bet = 2. + 2.*gk - gam(i,k-1) - qe1(i,k) = ( 3.*(q1(i,j,k-1)+gk*q1(i,j,k)) - qe1(i,k-1) ) / bet - qe2(i,k) = ( 3.*(q2(i,j,k-1)+gk*q2(i,j,k)) - qe2(i,k-1) ) / bet - gam(i,k) = gk / bet - enddo - enddo - - a_bot = 1. + gk*(gk+1.5) - xt1 = 2.*gk*(gk+1.) - do i=i1,i2 - xt2 = gk*(gk+0.5) - a_bot*gam(i,km) - qe1(i,km+1) = ( xt1*q1(i,j,km) + q1(i,j,km-1) - a_bot*qe1(i,km) ) / xt2 - qe2(i,km+1) = ( xt1*q2(i,j,km) + q2(i,j,km-1) - a_bot*qe2(i,km) ) / xt2 - enddo - - do k=km,1,-1 - do i=i1,i2 - qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1) - qe2(i,k) = qe2(i,k) - gam(i,k)*qe2(i,k+1) - enddo - enddo - endif - -!------------------ -! Apply constraints -!------------------ - if ( limiter/=0 ) then ! limit the top & bottom winds - do i=i1,i2 -! Top - if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. - if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0. -! Surface: - if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. - if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0. - enddo - endif - - do k=1,km+1 - do i=i1,i2 - q1e(i,j,k) = qe1(i,k) - q2e(i,j,k) = qe2(i,k) - enddo - enddo - - end subroutine edge_profile - -!TODO LMH 25may18: do not need delz defined on full compute domain; pass appropriate BCs instead - subroutine nh_bc(ptop, grav, kappa, cp, delp, delzBC, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - npx, npy, npz, bounded_domain, pkc_pertn, computepk3, fullhalo, bd) - - !INPUT: delp, delz (BC), pt - !OUTPUT: gz, pkc, pk3 (optional) - integer, intent(IN) :: npx, npy, npz - logical, intent(IN) :: pkc_pertn, computepk3, fullhalo, bounded_domain - real, intent(IN) :: ptop, kappa, cp, grav, BC_step, BC_split - type(fv_grid_bounds_type), intent(IN) :: bd - real, intent(IN) :: phis(bd%isd:bd%ied,bd%jsd:bd%jed) - real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: pt, delp - type(fv_nest_BC_type_3d), intent(IN) :: delzBC -#ifdef MULTI_GASES - real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,*):: q -#endif -#ifdef USE_COND - real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: q_con -#ifdef MOIST_CAPPA - real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: cappa -#endif -#endif - real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1):: gz, pkc, pk3 - - integer :: i,j,k - real :: gama !'gamma' - real :: ptk, rgrav, rkap, peln1, rdg - - integer :: istart, iend - - integer :: is, ie, js, je - integer :: isd, ied, jsd, jed - - is = bd%is - ie = bd%ie - js = bd%js - je = bd%je - isd = bd%isd - ied = bd%ied - jsd = bd%jsd - jed = bd%jed - - if (.not. bounded_domain) return - - if (is == 1) then - - call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%west_t0, delzBC%west_t1, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - pkc_pertn, computepk3, isd, ied, isd, 0, isd, 0, jsd, jed, jsd, jed, npz) - - endif - - if (ie == npx-1) then - - call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%east_t0, delzBC%east_t1, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - pkc_pertn, computepk3, isd, ied, npx, ied, npx, ied, jsd, jed, jsd, jed, npz) - - endif - - if (is == 1) then - istart = is - else - istart = isd - end if - if (ie == npx-1) then - iend = ie - else - iend = ied - end if - - if (js == 1) then - - call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%south_t0, delzBC%south_t1, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, jsd, 0, npz) - - end if - - if (je == npy-1) then - - call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%north_t0, delzBC%north_t1, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, npy, jed, npz) - endif - -end subroutine nh_bc - -subroutine nh_BC_k(ptop, grav, kappa, cp, delp, delzBC_t0, delzBC_t1, pt, phis, & -#ifdef MULTI_GASES - q , & -#endif -#ifdef USE_COND - q_con, & -#ifdef MOIST_CAPPA - cappa, & -#endif -#endif - pkc, gz, pk3, & - BC_step, BC_split, & - pkc_pertn, computepk3, isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz) - - integer, intent(IN) :: isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz - real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: delzBC_t0, delzBC_t1 - real, intent(IN) :: BC_step, BC_split - - logical, intent(IN) :: pkc_pertn, computepk3 - real, intent(IN) :: ptop, kappa, cp, grav - real, intent(IN) :: phis(isd:ied,jsd:jed) - real, intent(IN), dimension(isd:ied,jsd:jed,npz):: pt, delp -#ifdef MULTI_GASES - real, intent(IN), dimension(isd:ied,jsd:jed,npz,*):: q -#endif -#ifdef USE_COND - real, intent(IN), dimension(isd:ied,jsd:jed,npz):: q_con -#ifdef MOIST_CAPPA - real, intent(INOUT), dimension(isd:ied,jsd:jed,npz):: cappa -#endif -#endif - real, intent(INOUT), dimension(isd:ied,jsd:jed,npz+1):: gz, pkc, pk3 - - integer :: i,j,k - real :: gama !'gamma' - real :: ptk, rgrav, rkap, peln1, rdg, denom - - real, dimension(istart:iend, npz+1, jstart:jend ) :: pe, peln -#ifdef USE_COND - real, dimension(istart:iend, npz+1 ) :: peg, pelng -#endif - real, dimension(istart:iend, npz) :: gam, bb, dd, pkz - real, dimension(istart:iend, npz-1) :: g_rat - real, dimension(istart:iend) :: bet - real :: pm, delz_int - - - real :: pealn, pebln, rpkz -#ifdef MULTI_GASES - real gamax -#endif - rgrav = 1./grav - gama = 1./(1.-kappa) - ptk = ptop ** kappa - rkap = 1./kappa - peln1 = log(ptop) - rdg = - rdgas * rgrav - denom = 1./BC_split - - do j=jstart,jend - - !GZ - do i=istart,iend - gz(i,j,npz+1) = phis(i,j) - enddo - do k=npz,1,-1 - do i=istart,iend - delz_int = (delzBC_t0(i,j,k)*(BC_split-BC_step) + BC_step*delzBC_t1(i,j,k))*denom - gz(i,j,k) = gz(i,j,k+1) - delz_int*grav - enddo - enddo - - !Hydrostatic interface pressure - do i=istart,iend - pe(i,1,j) = ptop - peln(i,1,j) = peln1 -#ifdef USE_COND - peg(i,1) = ptop - pelng(i,1) = peln1 -#endif - enddo - do k=2,npz+1 - do i=istart,iend - pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1) - peln(i,k,j) = log(pe(i,k,j)) -#ifdef USE_COND - peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1)) - pelng(i,k) = log(peg(i,k)) -#endif - enddo - enddo - - !Perturbation nonhydro layer-mean pressure (NOT to the kappa) - do k=1,npz - do i=istart,iend - delz_int = (delzBC_t0(i,j,k)*(BC_split-BC_step) + BC_step*delzBC_t1(i,j,k))*denom - - !Full p -#ifdef MOIST_CAPPA - pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz_int*pt(i,j,k))) -#else -#ifdef MULTI_GASES - gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:))) - pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k))) -#else - pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k))) -#endif -#endif - !hydro -#ifdef USE_COND - pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k)) -#else - pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j)) -#endif - !Remove hydro cell-mean pressure - pkz(i,k) = pkz(i,k) - pm - enddo - enddo - - !pressure solver - do k=1,npz-1 - do i=istart,iend - g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1) - bb(i,k) = 2.*(1. + g_rat(i,k)) - dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1)) - enddo - enddo - - do i=istart,iend - bet(i) = bb(i,1) - pkc(i,j,1) = 0. - pkc(i,j,2) = dd(i,1)/bet(i) - bb(i,npz) = 2. - dd(i,npz) = 3.*pkz(i,npz) - enddo - do k=2,npz - do i=istart,iend - gam(i,k) = g_rat(i,k-1)/bet(i) - bet(i) = bb(i,k) - gam(i,k) - pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i) - enddo - enddo - do k=npz,2,-1 - do i=istart,iend - pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1) -#ifdef NHNEST_DEBUG - if (abs(pkc(i,j,k)) > 1.e5) then - print*, mpp_pe(), i,j,k, 'PKC: ', pkc(i,j,k) - endif -#endif - enddo - enddo - - - enddo - - do j=jstart,jend - - if (.not. pkc_pertn) then - do k=npz+1,1,-1 - do i=istart,iend - pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j) - enddo - enddo - endif - - !pk3 if necessary; doesn't require condenstate loading calculation - if (computepk3) then - do i=istart,iend - pk3(i,j,1) = ptk - enddo - do k=2,npz+1 - do i=istart,iend - pk3(i,j,k) = exp(kappa*log(pe(i,k,j))) - enddo - enddo - endif - - enddo - -end subroutine nh_BC_k - - -end module nh_utils_mod + +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the FV3 dynamical core. +!* +!* The FV3 dynamical core is free software: you can redistribute it +!* and/or modify it under the terms of the +!* GNU Lesser General Public License as published by the +!* Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* The FV3 dynamical core is distributed in the hope that it will be +!* useful, but WITHOUT ANYWARRANTY; 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 Lesser General Public +!* License along with the FV3 dynamical core. +!* If not, see . +!*********************************************************************** + +!>@brief The module 'nh_utils' peforms non-hydrostatic computations. +!>@author S. J. Lin, NOAA/GFDL + +module nh_utils_mod + +! Modules Included: +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +!
Module NameFunctions Included
constants_modrdgas, cp_air, grav
fv_arrays_modfv_grid_bounds_type, fv_grid_type
sw_core_modfill_4corners, del6_vt_flux
tp_core_modfv_tp_2d
+ + use constants_mod, only: rdgas, cp_air, grav + use tp_core_mod, only: fv_tp_2d + use sw_core_mod, only: fill_4corners, del6_vt_flux + use fv_arrays_mod, only: fv_grid_bounds_type, fv_grid_type,fv_nest_BC_type_3d +#ifdef MULTI_GASES + use multi_gases_mod, only: vicpqd, vicvqd +#endif + + implicit none + private + + public update_dz_c, update_dz_d, nh_bc + public sim_solver, sim1_solver, sim3_solver + public sim3p0_solver, rim_2d + public Riem_Solver_c + + real, parameter:: r3 = 1./3. + +CONTAINS + + subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, & + npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type, dz_min) +! !INPUT PARAMETERS: + type(fv_grid_bounds_type), intent(IN) :: bd + integer, intent(in):: is, ie, js, je, ng, km, npx, npy, grid_type + logical, intent(IN):: sw_corner, se_corner, ne_corner, nw_corner + real, intent(in):: dt, dz_min + real, intent(in):: dp0(km) + real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: ut, vt + real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng):: area + real, intent(inout):: gz(is-ng:ie+ng,js-ng:je+ng,km+1) + real, intent(in ):: zs(is-ng:ie+ng, js-ng:je+ng) + real, intent( out):: ws(is-ng:ie+ng, js-ng:je+ng) +! Local Work array: + real:: gz2(is-ng:ie+ng,js-ng:je+ng) + real, dimension(is-1:ie+2,js-1:je+1):: xfx, fx + real, dimension(is-1:ie+1,js-1:je+2):: yfx, fy + real, parameter:: r14 = 1./14. + integer i, j, k + integer:: is1, ie1, js1, je1 + integer:: ie2, je2 + real:: rdt, top_ratio, bot_ratio, int_ratio +!-------------------------------------------------------------------- + + rdt = 1. / dt + + top_ratio = dp0(1 ) / (dp0( 1)+dp0(2 )) + bot_ratio = dp0(km) / (dp0(km-1)+dp0(km)) + + is1 = is - 1 + js1 = js - 1 + + ie1 = ie + 1 + je1 = je + 1 + + ie2 = ie + 2 + je2 = je + 2 + +!$OMP parallel do default(none) shared(js1,je1,is1,ie2,km,je2,ie1,ut,top_ratio,vt, & +!$OMP bot_ratio,dp0,js,je,ng,is,ie,gz,grid_type, & +!$OMP bd,npx,npy,sw_corner,se_corner,ne_corner, & +!$OMP nw_corner,area) & +!$OMP private(gz2, xfx, yfx, fx, fy, int_ratio) + do 6000 k=1,km+1 + + if ( k==1 ) then + do j=js1, je1 + do i=is1, ie2 + xfx(i,j) = ut(i,j,1) + (ut(i,j,1)-ut(i,j,2))*top_ratio + enddo + enddo + do j=js1, je2 + do i=is1, ie1 + yfx(i,j) = vt(i,j,1) + (vt(i,j,1)-vt(i,j,2))*top_ratio + enddo + enddo + elseif ( k==km+1 ) then +! Bottom extrapolation + do j=js1, je1 + do i=is1, ie2 + xfx(i,j) = ut(i,j,km) + (ut(i,j,km)-ut(i,j,km-1))*bot_ratio +! xfx(i,j) = r14*(3.*ut(i,j,km-2)-13.*ut(i,j,km-1)+24.*ut(i,j,km)) +! if ( xfx(i,j)*ut(i,j,km)<0. ) xfx(i,j) = 0. + enddo + enddo + do j=js1, je2 + do i=is1, ie1 + yfx(i,j) = vt(i,j,km) + (vt(i,j,km)-vt(i,j,km-1))*bot_ratio +! yfx(i,j) = r14*(3.*vt(i,j,km-2)-13.*vt(i,j,km-1)+24.*vt(i,j,km)) +! if ( yfx(i,j)*vt(i,j,km)<0. ) yfx(i,j) = 0. + enddo + enddo + else + int_ratio = 1./(dp0(k-1)+dp0(k)) + do j=js1, je1 + do i=is1, ie2 + xfx(i,j) = (dp0(k)*ut(i,j,k-1)+dp0(k-1)*ut(i,j,k))*int_ratio + enddo + enddo + do j=js1, je2 + do i=is1, ie1 + yfx(i,j) = (dp0(k)*vt(i,j,k-1)+dp0(k-1)*vt(i,j,k))*int_ratio + enddo + enddo + endif + + do j=js-ng, je+ng + do i=is-ng, ie+ng + gz2(i,j) = gz(i,j,k) + enddo + enddo + + if (grid_type < 3) call fill_4corners(gz2, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner) + do j=js1, je1 + do i=is1, ie2 + if( xfx(i,j) > 0. ) then + fx(i,j) = gz2(i-1,j) + else + fx(i,j) = gz2(i ,j) + endif + fx(i,j) = xfx(i,j)*fx(i,j) + enddo + enddo + + if (grid_type < 3) call fill_4corners(gz2, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner) + do j=js1,je2 + do i=is1,ie1 + if( yfx(i,j) > 0. ) then + fy(i,j) = gz2(i,j-1) + else + fy(i,j) = gz2(i,j) + endif + fy(i,j) = yfx(i,j)*fy(i,j) + enddo + enddo + + do j=js1, je1 + do i=is1,ie1 + gz(i,j,k) = (gz2(i,j)*area(i,j) + fx(i,j)- fx(i+1,j)+ fy(i,j)- fy(i,j+1)) & + / ( area(i,j) + xfx(i,j)-xfx(i+1,j)+yfx(i,j)-yfx(i,j+1)) + enddo + enddo +6000 continue + +! Enforce monotonicity of height to prevent blowup +!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km,dz_min) + do j=js1, je1 + do k=km, 1, -1 + do i=is1, ie1 + gz(i,j,k) = max( gz(i,j,k), gz(i,j,k+1) + dz_min ) + enddo + enddo + do i=is1, ie1 + ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt + enddo + enddo + + end subroutine update_dz_c + + + subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea, & + dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac, dz_min, psm_bc) + + type(fv_grid_bounds_type), intent(IN) :: bd + integer, intent(in):: is, ie, js, je, ng, km, npx, npy + integer, intent(in):: hord, psm_bc + real, intent(in) :: rdt, dz_min + real, intent(in) :: dp0(km) + real, intent(in) :: area(is-ng:ie+ng,js-ng:je+ng) + real, intent(in) :: rarea(is-ng:ie+ng,js-ng:je+ng) + real, intent(inout):: damp(km+1) + integer, intent(inout):: ndif(km+1) + real, intent(in ) :: zs(is-ng:ie+ng,js-ng:je+ng) + real, intent(inout) :: zh(is-ng:ie+ng,js-ng:je+ng,km+1) + real, intent(inout), dimension(is:ie+1,js-ng:je+ng,km):: crx, xfx + real, intent(inout), dimension(is-ng:ie+ng,js:je+1,km):: cry, yfx + real, intent(out) :: ws(is:ie,js:je) + type(fv_grid_type), intent(IN), target :: gridstruct + real, intent(in) :: lim_fac +!----------------------------------------------------- +! Local array: + real, dimension(is: ie+1, js-ng:je+ng,km+1):: crx_adv, xfx_adv + real, dimension(is-ng:ie+ng,js: je+1,km+1 ):: cry_adv, yfx_adv + real, dimension(is:ie+1,js:je ):: fx + real, dimension(is:ie ,js:je+1):: fy + real, dimension(is-ng:ie+ng+1,js-ng:je+ng ):: fx2 + real, dimension(is-ng:ie+ng ,js-ng:je+ng+1):: fy2 + real, dimension(is-ng:ie+ng ,js-ng:je+ng ):: wk2, z2 + real:: ra_x(is:ie,js-ng:je+ng) + real:: ra_y(is-ng:ie+ng,js:je) +!-------------------------------------------------------------------- + integer i, j, k, isd, ied, jsd, jed + logical:: uniform_grid + + uniform_grid = .false. + + damp(km+1) = damp(km) + ndif(km+1) = ndif(km) + + isd = is - ng; ied = ie + ng + jsd = js - ng; jed = je + ng + + if (psm_bc == 0 )then +!$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, & +!$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv) + do j=jsd,jed + call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & + dp0, uniform_grid, 0) + if ( j<=je+1 .and. j>=js ) & + call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & + dp0, uniform_grid, 0) + enddo + else +!$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, & +!$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv) + do j=jsd,jed + call edge_profile_0grad(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & + dp0, uniform_grid, 0) + if ( j<=je+1 .and. j>=js ) & + call edge_profile_0grad(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & + dp0, uniform_grid, 0) + enddo + endif + +!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,area,xfx_adv,yfx_adv, & +!$OMP damp,zh,crx_adv,cry_adv,npx,npy,hord,gridstruct,bd, & +!$OMP ndif,rarea,lim_fac) & +!$OMP private(z2, fx2, fy2, ra_x, ra_y, fx, fy,wk2) + do k=1,km+1 + + do j=jsd,jed + do i=is,ie + ra_x(i,j) = area(i,j) + xfx_adv(i,j,k) - xfx_adv(i+1,j,k) + enddo + enddo + do j=js,je + do i=isd,ied + ra_y(i,j) = area(i,j) + yfx_adv(i,j,k) - yfx_adv(i,j+1,k) + enddo + enddo + + if ( damp(k)>1.E-5 ) then + do j=jsd,jed + do i=isd,ied + z2(i,j) = zh(i,j,k) + enddo + enddo + call fv_tp_2d(z2, crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, & + fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac) + call del6_vt_flux(ndif(k), npx, npy, damp(k), z2, wk2, fx2, fy2, gridstruct, bd) + do j=js,je + do i=is,ie + zh(i,j,k) = (z2(i,j)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & + / (ra_x(i,j)+ra_y(i,j)-area(i,j)) + (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*rarea(i,j) + enddo + enddo + else + call fv_tp_2d(zh(isd,jsd,k), crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, & + fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac) + do j=js,je + do i=is,ie + zh(i,j,k) = (zh(i,j,k)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & + / (ra_x(i,j) + ra_y(i,j) - area(i,j)) +! zh(i,j,k) = rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) & +! + zh(i,j,k)*(3.-rarea(i,j)*(ra_x(i,j) + ra_y(i,j))) + enddo + enddo + endif + + enddo + +!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt,dz_min) + do j=js, je + do k=km, 1, -1 + do i=is, ie +! Enforce monotonicity of height to prevent blowup + zh(i,j,k) = max( zh(i,j,k), zh(i,j,k+1) + dz_min ) + enddo + enddo + do i=is,ie + ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt + enddo + enddo + + end subroutine update_dz_d + + + subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & + akap, cappa, cp, & +#ifdef MULTI_GASES + kapad, & +#endif + ptop, hs, w3, pt, q_con, & + delp, gz, pef, ws, p_fac, a_imp, scale_m) + + integer, intent(in):: is, ie, js, je, ng, km + integer, intent(in):: ms + real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m + real, intent(in):: ws(is-ng:ie+ng,js-ng:je+ng) + real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delp + real, intent(in), dimension(is-ng:,js-ng:,1:):: q_con, cappa +#ifdef MULTI_GASES + real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: kapad +#endif + real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng) + real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: w3 +! OUTPUT PARAMETERS + real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: gz + real, intent( out), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: pef +! Local: + real, dimension(is-1:ie+1,km ):: dm, dz2, w2, pm2, gm2, cp2 + real, dimension(is-1:ie+1,km+1):: pem, pe2, peg +#ifdef MULTI_GASES + real, dimension(is-1:ie+1,km ):: kapad2 +#endif + real gama, rgrav + integer i, j, k + integer is1, ie1 + + gama = 1./(1.-akap) + rgrav = 1./grav + + is1 = is - 1 + ie1 = ie + 1 + +!$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, & +#ifdef MULTI_GASES +!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad) & +!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg, kapad2) +#else +!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) & +!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg) +#endif + do 2000 j=js-1, je+1 + + do k=1,km + do i=is1, ie1 + dm(i,k) = delp(i,j,k) + enddo + enddo + + do i=is1, ie1 + pef(i,j,1) = ptop ! full pressure at top + pem(i,1) = ptop +#ifdef USE_COND + peg(i,1) = ptop +#endif + enddo + + do k=2,km+1 + do i=is1, ie1 + pem(i,k) = pem(i,k-1) + dm(i,k-1) +#ifdef USE_COND +! Excluding contribution from condensates: + peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1)) +#endif + enddo + enddo + + do k=1,km + do i=is1, ie1 + dz2(i,k) = gz(i,j,k+1) - gz(i,j,k) +#ifdef USE_COND + pm2(i,k) = (peg(i,k+1)-peg(i,k))/log(peg(i,k+1)/peg(i,k)) + +#ifdef MOIST_CAPPA + cp2(i,k) = cappa(i,j,k) + gm2(i,k) = 1. / (1.-cp2(i,k)) +#endif + +#else + pm2(i,k) = dm(i,k)/log(pem(i,k+1)/pem(i,k)) +#endif +#ifdef MULTI_GASES + kapad2(i,k) = kapad(i,j,k) +#endif + dm(i,k) = dm(i,k) * rgrav + w2(i,k) = w3(i,j,k) + enddo + enddo + + + if ( a_imp < -0.01 ) then + call SIM3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, scale_m) + elseif ( a_imp <= 0.5 ) then + call RIM_2D(ms, dt, is1, ie1, km, rdgas, gama, gm2, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, & + dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.) + else + call SIM1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, & + dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac) + endif + + + do k=2,km+1 + do i=is1, ie1 + pef(i,j,k) = pe2(i,k) + pem(i,k) ! add hydrostatic full-component + enddo + enddo + +! Compute Height * grav (for p-gradient computation) + do i=is1, ie1 + gz(i,j,km+1) = hs(i,j) + enddo + + do k=km,1,-1 + do i=is1, ie1 + gz(i,j,k) = gz(i,j,k+1) - dz2(i,k)*grav + enddo + enddo + +2000 continue + + end subroutine Riem_Solver_c + + +!>GFDL - This routine will not give absoulte reproducibility when compiled with -fast-transcendentals. +!! GFDL - It is now inside of nh_core.F90 and being compiled without -fast-transcendentals. + subroutine Riem_Solver3test(ms, dt, is, ie, js, je, km, ng, & + isd, ied, jsd, jed, akap, cappa, cp, & +#ifdef MULTI_GASES + kapad, & +#endif + ptop, zs, q_con, w, delz, pt, & + delp, zh, pe, ppe, pk3, pk, peln, & + ws, scale_m, p_fac, a_imp, & + use_logp, last_call, fp_out) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: gz: grav*height at edges +! pe: full hydrostatic pressure +! ppe: non-hydrostatic pressure perturbation +!-------------------------------------------- + integer, intent(in):: ms, is, ie, js, je, km, ng + integer, intent(in):: isd, ied, jsd, jed + real, intent(in):: dt ! the BIG horizontal Lagrangian time step + real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m + real, intent(in):: zs(isd:ied,jsd:jed) + logical, intent(in):: last_call, use_logp, fp_out + real, intent(in):: ws(is:ie,js:je) + real, intent(in), dimension(isd:,jsd:,1:):: q_con, cappa + real, intent(in), dimension(isd:ied,jsd:jed,km):: delp, pt +#ifdef MULTI_GASES + real, intent(in), dimension(isd:ied,jsd:jed,km):: kapad +#endif + real, intent(inout), dimension(isd:ied,jsd:jed,km+1):: zh + real, intent(inout), dimension(isd:ied,jsd:jed,km):: w + real, intent(inout):: pe(is-1:ie+1,km+1,js-1:je+1) + real, intent(out):: peln(is:ie,km+1,js:je) ! ln(pe) + real, intent(out), dimension(isd:ied,jsd:jed,km+1):: ppe + real, intent(out):: delz(is:ie,js:je,km) + real, intent(out):: pk(is:ie,js:je,km+1) + real, intent(out):: pk3(isd:ied,jsd:jed,km+1) +! Local: + real, dimension(is:ie,km):: dm, dz2, pm2, w2, gm2, cp2 + real, dimension(is:ie,km+1)::pem, pe2, peln2, peg, pelng +#ifdef MULTI_GASES + real, dimension(is:ie,km):: kapad2 +#endif + real gama, rgrav, ptk, peln1 + integer i, j, k + + gama = 1./(1.-akap) + rgrav = 1./grav + peln1 = log(ptop) + ptk = exp(akap*peln1) + +!$OMP parallel do default(none) shared(is,ie,js,je,km,delp,ptop,peln1,pk3,ptk,akap,rgrav,zh,pt, & +!$OMP w,a_imp,dt,gama,ws,p_fac,scale_m,ms,delz,last_call, & +#ifdef MULTI_GASES +!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad ) & +!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2,kapad2) +#else +!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con ) & +!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2) +#endif + do 2000 j=js, je + + do k=1,km + do i=is, ie + dm(i,k) = delp(i,j,k) +#ifdef MOIST_CAPPA + cp2(i,k) = cappa(i,j,k) +#endif +#ifdef MULTI_GASES + kapad2(i,k) = kapad(i,j,k) +#endif + enddo + enddo + + do i=is,ie + pem(i,1) = ptop + peln2(i,1) = peln1 + pk3(i,j,1) = ptk +#ifdef USE_COND + peg(i,1) = ptop + pelng(i,1) = peln1 +#endif + enddo + do k=2,km+1 + do i=is,ie + pem(i,k) = pem(i,k-1) + dm(i,k-1) + peln2(i,k) = log(pem(i,k)) +#ifdef USE_COND +! Excluding contribution from condensates: +! peln used during remap; pk3 used only for p_grad + peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1)) + pelng(i,k) = log(peg(i,k)) +#endif + pk3(i,j,k) = exp(akap*peln2(i,k)) + enddo + enddo + + do k=1,km + do i=is, ie +#ifdef USE_COND + pm2(i,k) = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k)) + +#ifdef MOIST_CAPPA + gm2(i,k) = 1. / (1.-cp2(i,k)) +#endif + +#else + pm2(i,k) = dm(i,k)/(peln2(i,k+1)-peln2(i,k)) +#endif + dm(i,k) = dm(i,k) * rgrav + dz2(i,k) = zh(i,j,k+1) - zh(i,j,k) + w2(i,k) = w(i,j,k) + enddo + enddo + + + if ( a_imp < -0.999 ) then + call SIM3p0_solver(dt, is, ie, km, rdgas, gama, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, scale_m ) + elseif ( a_imp < -0.5 ) then + call SIM3_solver(dt, is, ie, km, rdgas, gama, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), abs(a_imp), p_fac, scale_m) + elseif ( a_imp <= 0.5 ) then + call RIM_2D(ms, dt, is, ie, km, rdgas, gama, gm2, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, & + dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.) + elseif ( a_imp > 0.999 ) then + call SIM1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac) + else + call SIM_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), & + a_imp, p_fac, scale_m) + endif + + + do k=1, km + do i=is, ie + w(i,j,k) = w2(i,k) + delz(i,j,k) = dz2(i,k) + enddo + enddo + + if ( last_call ) then + do k=1,km+1 + do i=is,ie + peln(i,k,j) = peln2(i,k) + pk(i,j,k) = pk3(i,j,k) + pe(i,k,j) = pem(i,k) + enddo + enddo + endif + + if( fp_out ) then + do k=1,km+1 + do i=is, ie + ppe(i,j,k) = pe2(i,k) + pem(i,k) + enddo + enddo + else + do k=1,km+1 + do i=is, ie + ppe(i,j,k) = pe2(i,k) + enddo + enddo + endif + + if ( use_logp ) then + do k=2,km+1 + do i=is, ie + pk3(i,j,k) = peln2(i,k) + enddo + enddo + endif + + do i=is, ie + zh(i,j,km+1) = zs(i,j) + enddo + do k=km,1,-1 + do i=is, ie + zh(i,j,k) = zh(i,j,k+1) - dz2(i,k) + enddo + enddo + +2000 continue + + end subroutine Riem_Solver3test + + subroutine imp_diff_w(j, is, ie, js, je, ng, km, cd, delz, ws, w, w3) + integer, intent(in) :: j, is, ie, js, je, km, ng + real, intent(in) :: cd + real, intent(in) :: delz(is:ie, km) !< delta-height (m) + real, intent(in) :: w(is:ie, km) !< vertical vel. (m/s) + real, intent(in) :: ws(is:ie) + real, intent(out) :: w3(is-ng:ie+ng,js-ng:je+ng,km) +! Local: + real, dimension(is:ie,km):: c, gam, dz, wt + real:: bet(is:ie) + real:: a + integer:: i, k + + do k=2,km + do i=is,ie + dz(i,k) = 0.5*(delz(i,k-1)+delz(i,k)) + enddo + enddo + do k=1,km-1 + do i=is,ie + c(i,k) = -cd/(dz(i,k+1)*delz(i,k)) + enddo + enddo + +! model top: + do i=is,ie + bet(i) = 1. - c(i,1) ! bet(i) = b + wt(i,1) = w(i,1) / bet(i) + enddo + +! Interior: + do k=2,km-1 + do i=is,ie + gam(i,k) = c(i,k-1)/bet(i) + a = cd/(dz(i,k)*delz(i,k)) + bet(i) = (1.+a-c(i,k)) + a*gam(i,k) + wt(i,k) = (w(i,k) + a*wt(i,k-1)) / bet(i) + enddo + enddo + +! Bottom: + do i=is,ie + gam(i,km) = c(i,km-1) / bet(i) + a = cd/(dz(i,km)*delz(i,km)) + wt(i,km) = (w(i,km) + 2.*ws(i)*cd/delz(i,km)**2 & + + a*wt(i,km-1))/(1. + a + (cd+cd)/delz(i,km)**2 + a*gam(i,km)) + enddo + + do k=km-1,1,-1 + do i=is,ie + wt(i,k) = wt(i,k) - gam(i,k+1)*wt(i,k+1) + enddo + enddo + + do k=1,km + do i=is,ie + w3(i,j,k) = wt(i,k) + enddo + enddo + + end subroutine imp_diff_w + + + subroutine RIM_2D(ms, bdt, is, ie, km, rgas, gama, gm2, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm2, pm2, w2, dz2, pt2, ws, c_core ) + + integer, intent(in):: ms, is, ie, km + real, intent(in):: bdt, gama, rgas + real, intent(in), dimension(is:ie,km):: dm2, pm2, gm2 + logical, intent(in):: c_core + real, intent(in ) :: pt2(is:ie,km) + real, intent(in ) :: ws(is:ie) +! IN/OUT: + real, intent(inout):: dz2(is:ie,km) + real, intent(inout):: w2(is:ie,km) + real, intent(out ):: pe2(is:ie,km+1) +#ifdef MULTI_GASES + real, intent(inout), dimension(is:ie,km):: kapad2 +#endif +! Local: + real:: ws2(is:ie) + real, dimension(km+1):: m_bot, m_top, r_bot, r_top, pe1, pbar, wbar + real, dimension(km):: r_hi, r_lo, dz, wm, dm, dts + real, dimension(km):: pf1, wc, cm , pp, pt1 + real:: dt, rdt, grg, z_frac, ptmp1, rden, pf, time_left + real:: m_surf +#ifdef MULTI_GASES + real gamax +#endif + integer:: i, k, n, ke, kt1, ktop + integer:: ks0, ks1 + + grg = gama * rgas + rdt = 1. / bdt + dt = bdt / real(ms) + + pbar(:) = 0. + wbar(:) = 0. + + do i=is,ie + ws2(i) = 2.*ws(i) + enddo + + do 6000 i=is,ie + + do k=1,km + dz(k) = dz2(i,k) + dm(k) = dm2(i,k) + wm(k) = w2(i,k)*dm(k) + pt1(k) = pt2(i,k) + enddo + + pe1(:) = 0. + wbar(km+1) = ws(i) + + ks0 = 1 + if ( ms > 1 .and. ms < 8 ) then +! Continuity of (pbar, wbar) is maintained + do k=1, km + rden = -rgas*dm(k)/dz(k) +#ifdef MOIST_CAPPA + pf1(k) = exp( gm2(i,k)*log(rden*pt1(k)) ) +! dts(k) = -dz(k)/sqrt(gm2(i,k)*rgas*pf1(k)/rden) + dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + pf1(k) = exp( gamax*log(rden*pt1(k)) ) +#else + pf1(k) = exp( gama*log(rden*pt1(k)) ) +#endif + dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden) +#endif + if ( bdt > dts(k) ) then + ks0 = k-1 + goto 222 + endif + enddo + ks0 = km +222 if ( ks0==1 ) goto 244 + + do k=1, ks0 + cm(k) = dm(k) / dts(k) + wc(k) = wm(k) / dts(k) + pp(k) = pf1(k) - pm2(i,k) + enddo + + wbar(1) = (wc(1)+pp(1)) / cm(1) + do k=2, ks0 + wbar(k) = (wc(k-1)+wc(k) + pp(k)-pp(k-1)) / (cm(k-1)+cm(k)) + pbar(k) = bdt*(cm(k-1)*wbar(k) - wc(k-1) + pp(k-1)) + pe1(k) = pbar(k) + enddo + + if ( ks0 == km ) then + pbar(km+1) = bdt*(cm(km)*wbar(km+1) - wc(km) + pp(km)) + if ( c_core ) then + do k=1,km + dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) + enddo + else + do k=1,km + dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) + w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k) + enddo + endif + pe2(i,1) = 0. + do k=2,km+1 + pe2(i,k) = pbar(k)*rdt + enddo + goto 6000 ! next i + else + if ( c_core ) then + do k=1, ks0-1 + dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) + enddo + else + do k=1, ks0-1 + dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k)) + w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k) + enddo + endif + pbar(ks0) = pbar(ks0) / real(ms) + endif + endif +244 ks1 = ks0 + + do 5000 n=1, ms + + do k=ks1, km + rden = -rgas*dm(k)/dz(k) +#ifdef MOIST_CAPPA + pf = exp( gm2(i,k)*log(rden*pt1(k)) ) +! dts(k) = -dz(k) / sqrt( gm2(i,k)*rgas*pf/rden ) + dts(k) = -dz(k) / sqrt( grg*pf/rden ) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + pf = exp( gamax*log(rden*pt1(k)) ) +#else + pf = exp( gama*log(rden*pt1(k)) ) +#endif + dts(k) = -dz(k) / sqrt( grg*pf/rden ) +#endif + ptmp1 = dts(k)*(pf - pm2(i,k)) + r_lo(k) = wm(k) + ptmp1 + r_hi(k) = wm(k) - ptmp1 + enddo + + ktop = ks1 + do k=ks1, km + if( dt > dts(k) ) then + ktop = k-1 + goto 333 + endif + enddo + ktop = km +333 continue + + if ( ktop >= ks1 ) then + do k=ks1, ktop + z_frac = dt/dts(k) + r_bot(k ) = z_frac*r_lo(k) + r_top(k+1) = z_frac*r_hi(k) + m_bot(k ) = z_frac* dm(k) + m_top(k+1) = m_bot(k) + enddo + if ( ktop == km ) goto 666 + endif + + do k=ktop+2, km+1 + m_top(k) = 0. + r_top(k) = 0. + enddo + + kt1 = max(1, ktop) + do 444 ke=km+1, ktop+2, -1 + time_left = dt + do k=ke-1, kt1, -1 + if ( time_left > dts(k) ) then + time_left = time_left - dts(k) + m_top(ke) = m_top(ke) + dm(k) + r_top(ke) = r_top(ke) + r_hi(k) + else + z_frac = time_left/dts(k) + m_top(ke) = m_top(ke) + z_frac*dm(k) + r_top(ke) = r_top(ke) + z_frac*r_hi(k) + go to 444 ! next level + endif + enddo +444 continue + + do k=ktop+1, km + m_bot(k) = 0. + r_bot(k) = 0. + enddo + + do 4000 ke=ktop+1, km + time_left = dt + do k=ke, km + if ( time_left > dts(k) ) then + time_left = time_left - dts(k) + m_bot(ke) = m_bot(ke) + dm(k) + r_bot(ke) = r_bot(ke) + r_lo(k) + else + z_frac = time_left/dts(k) + m_bot(ke) = m_bot(ke) + z_frac* dm(k) + r_bot(ke) = r_bot(ke) + z_frac*r_lo(k) + go to 4000 ! next interface + endif + enddo + m_surf = m_bot(ke) + do k=km, kt1, -1 + if ( time_left > dts(k) ) then + time_left = time_left - dts(k) + m_bot(ke) = m_bot(ke) + dm(k) + r_bot(ke) = r_bot(ke) - r_hi(k) + else + z_frac = time_left/dts(k) + m_bot(ke) = m_bot(ke) + z_frac* dm(k) + r_bot(ke) = r_bot(ke) - z_frac*r_hi(k) + (m_bot(ke)-m_surf)*ws2(i) + go to 4000 ! next interface + endif + enddo +4000 continue + +666 if ( ks1==1 ) wbar(1) = r_bot(1) / m_bot(1) + do k=ks1+1, km + wbar(k) = (r_bot(k)+r_top(k)) / (m_top(k)+m_bot(k)) + enddo +! pbar here is actually dt*pbar + do k=ks1+1, km+1 + pbar(k) = m_top(k)*wbar(k) - r_top(k) + pe1(k) = pe1(k) + pbar(k) + enddo + + if ( n==ms ) then + if ( c_core ) then + do k=ks1, km + dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k)) + enddo + else + do k=ks1, km + dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k)) + w2(i,k) = ( wm(k) + pbar(k+1) - pbar(k) ) / dm(k) + enddo + endif + else + do k=ks1, km + dz(k) = dz(k) + dt*(wbar(k+1)-wbar(k)) + wm(k) = wm(k) + pbar(k+1) - pbar(k) + enddo + endif + +5000 continue + pe2(i,1) = 0. + do k=2,km+1 + pe2(i,k) = pe1(k)*rdt + enddo + +6000 continue ! end i-loop + + end subroutine RIM_2D + + subroutine SIM3_solver(dt, is, ie, km, rgas, gama, kappa, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m) + integer, intent(in):: is, ie, km + real, intent(in):: dt, rgas, gama, kappa, alpha, p_fac, scale_m + real, intent(in ), dimension(is:ie,km):: dm, pt2 + real, intent(in ):: ws(is:ie) + real, intent(in ), dimension(is:ie,km+1):: pem + real, intent(out):: pe2(is:ie,km+1) + real, intent(inout), dimension(is:ie,km):: dz2, w2 +#ifdef MULTI_GASES + real, intent(inout), dimension(is:ie,km):: kapad2 +#endif +! Local + real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam + real, dimension(is:ie,km+1):: pp + real, dimension(is:ie):: p1, wk1, bet + real beta, t2, t1g, rdt, ra, capa1, r2g, r6g +#ifdef MULTI_GASES + real gamax, capa1x, t1gx +#endif + integer i, k + + beta = 1. - alpha + ra = 1. / alpha + t2 = beta / alpha + t1g = gama * 2.*(alpha*dt)**2 + rdt = 1. / dt + capa1 = kappa - 1. + r2g = grav / 2. + r6g = grav / 6. + + + do k=1,km + do i=is, ie + w1(i,k) = w2(i,k) +! Full pressure at center +#ifdef MULTI_GASES + gamax = 1. / (1.-kapad2(i,k)) + aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) +#else + aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) +#endif + enddo + enddo + + do k=1,km-1 + do i=is, ie + g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction + bb(i,k) = 2.*(1.+g_rat(i,k)) + dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1)) + enddo + enddo + +! pe2 is full p at edges + do i=is, ie +! Top: + bet(i) = bb(i,1) + pe2(i,1) = pem(i,1) + pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i) +! Bottom: + bb(i,km) = 2. + dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km) + enddo + + do k=2,km + do i=is, ie + gam(i,k) = g_rat(i,k-1) / bet(i) + bet(i) = bb(i,k) - gam(i,k) + pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i) + enddo + enddo + + do k=km, 2, -1 + do i=is, ie + pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1) + enddo + enddo +! done reconstruction of full: + +! pp is pert. p at edges + do k=1, km+1 + do i=is, ie + pp(i,k) = pe2(i,k) - pem(i,k) + enddo + enddo + + do k=2, km + do i=is, ie +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + t1gx = gamax*2.*(alpha*dt)**2 + aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) +#else + aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) +#endif + wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k)) + aa(i,k) = aa(i,k) - scale_m*dm(i,1) + enddo + enddo + do i=is, ie + bet(i) = dm(i,1) - aa(i,2) + w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2) + wk(i,2)) / bet(i) + enddo + do k=2,km-1 + do i=is, ie + gam(i,k) = aa(i,k) / bet(i) + bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) + w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) & + - aa(i,k)*w2(i,k-1)) / bet(i) + enddo + enddo + do i=is, ie +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,km)) + t1gx = gamax*2.*(alpha*dt)**2 + wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) +#else + wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) +#endif + gam(i,km) = aa(i,km) / bet(i) + bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) + w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk(i,km) + & + wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i) + enddo + do k=km-1, 1, -1 + do i=is, ie + w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) + enddo + enddo + +! pe2 is updated perturbation p at edges + do i=is, ie + pe2(i,1) = 0. + enddo + do k=1,km + do i=is, ie + pe2(i,k+1) = pe2(i,k) + ( dm(i,k)*(w2(i,k)-w1(i,k))*rdt & + - beta*(pp(i,k+1)-pp(i,k)) )*ra + enddo + enddo + +! Full non-hydro pressure at edges: + do i=is, ie + pe2(i,1) = pem(i,1) + enddo + do k=2,km+1 + do i=is, ie + pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k)) + enddo + enddo + + do i=is, ie +! Recover cell-averaged pressure + p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km) +#ifdef MULTI_GASES + capa1x = kapad2(i,km) - 1. + dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) ) +#else + dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) ) +#endif + enddo + + do k=km-1, 1, -1 + do i=is, ie + p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i) +#ifdef MULTI_GASES + capa1x = kapad2(i,k) - 1. + dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) ) +#else + dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) ) +#endif + enddo + enddo + + do k=1,km+1 + do i=is, ie + pe2(i,k) = pe2(i,k) - pem(i,k) + pe2(i,k) = pe2(i,k) + beta*(pp(i,k) - pe2(i,k)) + enddo + enddo + + end subroutine SIM3_solver + + subroutine SIM3p0_solver(dt, is, ie, km, rgas, gama, kappa, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm, & + pem, w2, dz2, pt2, ws, p_fac, scale_m) +! Sa SIM3, but for beta==0 + integer, intent(in):: is, ie, km + real, intent(in):: dt, rgas, gama, kappa, p_fac, scale_m + real, intent(in ), dimension(is:ie,km):: dm, pt2 + real, intent(in ):: ws(is:ie) + real, intent(in ):: pem(is:ie,km+1) + real, intent(out):: pe2(is:ie,km+1) + real, intent(inout), dimension(is:ie,km):: dz2, w2 +#ifdef MULTI_GASES + real, intent(inout), dimension(is:ie,km):: kapad2 +#endif +! Local + real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam + real, dimension(is:ie,km+1):: pp + real, dimension(is:ie):: p1, wk1, bet + real t1g, rdt, capa1, r2g, r6g +#ifdef MULTI_GASES + real gamax, capa1x, t1gx +#endif + integer i, k + + t1g = 2.*gama*dt**2 + rdt = 1. / dt + capa1 = kappa - 1. + r2g = grav / 2. + r6g = grav / 6. + + do k=1,km + do i=is, ie + w1(i,k) = w2(i,k) +! Full pressure at center +#ifdef MULTI_GASES + gamax = 1. / ( 1. - kapad2(i,k) ) + aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) +#else + aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k))) +#endif + enddo + enddo + + do k=1,km-1 + do i=is, ie + g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction + bb(i,k) = 2.*(1.+g_rat(i,k)) + dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1)) + enddo + enddo + +! pe2 is full p at edges + do i=is, ie +! Top: + bet(i) = bb(i,1) + pe2(i,1) = pem(i,1) + pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i) +! Bottom: + bb(i,km) = 2. + dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km) + enddo + + do k=2,km + do i=is, ie + gam(i,k) = g_rat(i,k-1) / bet(i) + bet(i) = bb(i,k) - gam(i,k) + pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i) + enddo + enddo + + do k=km, 2, -1 + do i=is, ie + pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1) + enddo + enddo +! done reconstruction of full: + +! pp is pert. p at edges + do k=1, km+1 + do i=is, ie + pp(i,k) = pe2(i,k) - pem(i,k) + enddo + enddo + + do k=2, km + do i=is, ie +#ifdef MULTI_GASES + gamax = 1. / (1.-kapad2(i,k)) + t1gx = 2.*gamax*dt**2 + aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1) +#else + aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1) +#endif + enddo + enddo + do i=is, ie + bet(i) = dm(i,1) - aa(i,2) + w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2)) / bet(i) + enddo + do k=2,km-1 + do i=is, ie + gam(i,k) = aa(i,k) / bet(i) + bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) + w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1))/bet(i) + enddo + enddo + do i=is, ie +#ifdef MULTI_GASES + gamax = 1. / (1.-kapad2(i,km)) + t1gx = 2.*gamax*dt**2 + wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) +#else + wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) +#endif + gam(i,km) = aa(i,km) / bet(i) + bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) + w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk1(i)*ws(i) - & + aa(i,km)*w2(i,km-1)) / bet(i) + enddo + do k=km-1, 1, -1 + do i=is, ie + w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) + enddo + enddo + +! pe2 is updated perturbation p at edges + do i=is, ie + pe2(i,1) = 0. + enddo + do k=1,km + do i=is, ie + pe2(i,k+1) = pe2(i,k) + dm(i,k)*(w2(i,k)-w1(i,k))*rdt + enddo + enddo + +! Full non-hydro pressure at edges: + do i=is, ie + pe2(i,1) = pem(i,1) + enddo + do k=2,km+1 + do i=is, ie + pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k)) + enddo + enddo + + do i=is, ie +! Recover cell-averaged pressure + p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km) +#ifdef MULTI_GASES + capa1x = kapad2(i,km) - 1. + dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) ) +#else + dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) ) +#endif + enddo + + do k=km-1, 1, -1 + do i=is, ie + p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3-g_rat(i,k)*p1(i) +#ifdef MULTI_GASES + capa1x = kapad2(i,k) - 1. + dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) ) +#else + dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) ) +#endif + enddo + enddo + + do k=1,km+1 + do i=is, ie + pe2(i,k) = pe2(i,k) - pem(i,k) + enddo + enddo + + end subroutine SIM3p0_solver + + + subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe, dm2, & + pm2, pem, w2, dz2, pt2, ws, p_fac) + integer, intent(in):: is, ie, km + real, intent(in):: dt, rgas, gama, kappa, p_fac + real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2 + real, intent(in ):: ws(is:ie) + real, intent(in ), dimension(is:ie,km+1):: pem + real, intent(out):: pe(is:ie,km+1) + real, intent(inout), dimension(is:ie,km):: dz2, w2 +#ifdef MULTI_GASES + real, intent(inout), dimension(is:ie,km):: kapad2 +#endif +! Local + real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam + real, dimension(is:ie,km+1):: pp + real, dimension(is:ie):: p1, bet + real t1g, rdt, capa1 +#ifdef MULTI_GASES + real gamax, capa1x, t1gx +#endif + integer i, k + +#ifdef MOIST_CAPPA + t1g = 2.*dt*dt +#else + t1g = gama * 2.*dt*dt +#endif + rdt = 1. / dt + capa1 = kappa - 1. + + do k=1,km + do i=is, ie +#ifdef MOIST_CAPPA + pe(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#else +#ifdef MULTI_GASES + gamax = 1. / ( 1. - kapad2(i,k) ) + pe(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#else + pe(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#endif +#endif + w1(i,k) = w2(i,k) + enddo + enddo + + do k=1,km-1 + do i=is, ie + g_rat(i,k) = dm2(i,k)/dm2(i,k+1) + bb(i,k) = 2.*(1.+g_rat(i,k)) + dd(i,k) = 3.*(pe(i,k) + g_rat(i,k)*pe(i,k+1)) + enddo + enddo + + do i=is, ie + bet(i) = bb(i,1) + pp(i,1) = 0. + pp(i,2) = dd(i,1) / bet(i) + bb(i,km) = 2. + dd(i,km) = 3.*pe(i,km) + enddo + + do k=2,km + do i=is, ie + gam(i,k) = g_rat(i,k-1) / bet(i) + bet(i) = bb(i,k) - gam(i,k) + pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i) + enddo + enddo + + do k=km, 2, -1 + do i=is, ie + pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1) + enddo + enddo + +! Start the w-solver + do k=2, km + do i=is, ie +#ifdef MOIST_CAPPA + aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + t1gx = gamax * 2.*dt*dt + aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) +#else + aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k)) +#endif +#endif + enddo + enddo + do i=is, ie + bet(i) = dm2(i,1) - aa(i,2) + w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2)) / bet(i) + enddo + do k=2,km-1 + do i=is, ie + gam(i,k) = aa(i,k) / bet(i) + bet(i) = dm2(i,k) - (aa(i,k) + aa(i,k+1) + aa(i,k)*gam(i,k)) + w2(i,k) = (dm2(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1)) / bet(i) + enddo + enddo + do i=is, ie +#ifdef MOIST_CAPPA + p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,km)) + t1gx = gamax * 2.*dt*dt + p1(i) = t1gx/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) +#else + p1(i) = t1g/dz2(i,km)*(pem(i,km+1)+pp(i,km+1)) +#endif +#endif + gam(i,km) = aa(i,km) / bet(i) + bet(i) = dm2(i,km) - (aa(i,km)+p1(i) + aa(i,km)*gam(i,km)) + w2(i,km) = (dm2(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-p1(i)*ws(i)-aa(i,km)*w2(i,km-1))/bet(i) + enddo + do k=km-1, 1, -1 + do i=is, ie + w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) + enddo + enddo + + do i=is, ie + pe(i,1) = 0. + enddo + do k=1,km + do i=is, ie + pe(i,k+1) = pe(i,k) + dm2(i,k)*(w2(i,k)-w1(i,k))*rdt + enddo + enddo + + do i=is, ie + p1(i) = ( pe(i,km) + 2.*pe(i,km+1) )*r3 +#ifdef MOIST_CAPPA + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#else +#ifdef MULTI_GASES + capa1x = kapad2(i,km)-1. + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#else + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#endif +#endif + enddo + + do k=km-1, 1, -1 + do i=is, ie + p1(i) = (pe(i,k) + bb(i,k)*pe(i,k+1) + g_rat(i,k)*pe(i,k+2))*r3 - g_rat(i,k)*p1(i) +#ifdef MOIST_CAPPA + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) + +#else +#ifdef MULTI_GASES + capa1x = kapad2(i,k)-1. + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) +#else + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) +#endif +#endif + enddo + enddo + + end subroutine SIM1_solver + + subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, & +#ifdef MULTI_GASES + kapad2, & +#endif + pe2, dm2, & + pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m) + integer, intent(in):: is, ie, km + real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m + real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2 + real, intent(in ):: ws(is:ie) + real, intent(in ), dimension(is:ie,km+1):: pem + real, intent(out):: pe2(is:ie,km+1) + real, intent(inout), dimension(is:ie,km):: dz2, w2 +#ifdef MULTI_GASES + real, intent(inout), dimension(is:ie,km):: kapad2 +#endif +! Local + real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam + real, dimension(is:ie,km+1):: pp + real, dimension(is:ie):: p1, wk1, bet + real beta, t2, t1g, rdt, ra, capa1 +#ifdef MULTI_GASES + real gamax, capa1x, t1gx +#endif + integer i, k + + beta = 1. - alpha + ra = 1. / alpha + t2 = beta / alpha +#ifdef MOIST_CAPPA + t1g = 2.*(alpha*dt)**2 +#else + t1g = 2.*gama*(alpha*dt)**2 +#endif + rdt = 1. / dt + capa1 = kappa - 1. + + do k=1,km + do i=is, ie + w1(i,k) = w2(i,k) +! P_g perturbation +#ifdef MOIST_CAPPA + pe2(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + pe2(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#else + pe2(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k) +#endif +#endif + enddo + enddo + + do k=1,km-1 + do i=is, ie + g_rat(i,k) = dm2(i,k)/dm2(i,k+1) + bb(i,k) = 2.*(1.+g_rat(i,k)) + dd(i,k) = 3.*(pe2(i,k) + g_rat(i,k)*pe2(i,k+1)) + enddo + enddo + + do i=is, ie + bet(i) = bb(i,1) + pp(i,1) = 0. + pp(i,2) = dd(i,1) / bet(i) + bb(i,km) = 2. + dd(i,km) = 3.*pe2(i,km) + enddo + + do k=2,km + do i=is, ie + gam(i,k) = g_rat(i,k-1) / bet(i) + bet(i) = bb(i,k) - gam(i,k) + pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i) + enddo + enddo + + do k=km, 2, -1 + do i=is, ie + pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1) + enddo + enddo + + do k=1, km+1 + do i=is, ie +! pe2 is Full p + pe2(i,k) = pem(i,k) + pp(i,k) + enddo + enddo + + do k=2, km + do i=is, ie +#ifdef MOIST_CAPPA + aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,k)) + t1gx = 2.*gamax*(alpha*dt)**2 + aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) +#else + aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) +#endif +#endif + wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k)) + aa(i,k) = aa(i,k) - scale_m*dm2(i,1) + enddo + enddo +! Top: + do i=is, ie + bet(i) = dm2(i,1) - aa(i,2) + w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2) + wk(i,2)) / bet(i) + enddo +! Interior: + do k=2,km-1 + do i=is, ie + gam(i,k) = aa(i,k) / bet(i) + bet(i) = dm2(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k)) + w2(i,k) = (dm2(i,k)*w1(i,k) + dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) & + - aa(i,k)*w2(i,k-1)) / bet(i) + enddo + enddo +! Bottom: k=km + do i=is, ie +#ifdef MOIST_CAPPA + wk1(i) = t1g*gm2(i,km)/dz2(i,km)*pe2(i,km+1) +#else +#ifdef MULTI_GASES + gamax = 1./(1.-kapad2(i,km)) + t1gx = 2.*gamax*(alpha*dt)**2 + wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1) +#else + wk1(i) = t1g/dz2(i,km)*pe2(i,km+1) +#endif +#endif + gam(i,km) = aa(i,km) / bet(i) + bet(i) = dm2(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km)) + w2(i,km) = (dm2(i,km)*w1(i,km) + dt*(pp(i,km+1)-pp(i,km)) - wk(i,km) + & + wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i) + enddo + do k=km-1, 1, -1 + do i=is, ie + w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1) + enddo + enddo + + do i=is, ie + pe2(i,1) = 0. + enddo + do k=1,km + do i=is, ie + pe2(i,k+1) = pe2(i,k) + ( dm2(i,k)*(w2(i,k)-w1(i,k))*rdt & + - beta*(pp(i,k+1)-pp(i,k)) ) * ra + enddo + enddo + + do i=is, ie + p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 +#ifdef MOIST_CAPPA + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#else +#ifdef MULTI_GASES + capa1x = kapad2(i,km)-1. + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#else + dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km)))) +#endif +#endif + enddo + + do k=km-1, 1, -1 + do i=is, ie + p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i) +! delz = -dm*R*T_m / p_gas +#ifdef MOIST_CAPPA + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) +#else +#ifdef MULTI_GASES + capa1x = kapad2(i,k)-1. + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) +#else + dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k)))) +#endif +#endif + enddo + enddo + + do k=1, km+1 + do i=is, ie + pe2(i,k) = pe2(i,k) + beta*(pp(i,k)-pe2(i,k)) + enddo + enddo + + end subroutine SIM_solver + + + subroutine edge_scalar(q1, qe, i1, i2, km, id) +! Optimized for wind profile reconstruction: + integer, intent(in):: i1, i2, km + integer, intent(in):: id ! 0: pp 1: wind + real, intent(in ), dimension(i1:i2,km):: q1 + real, intent(out), dimension(i1:i2,km+1):: qe +!----------------------------------------------------------------------- + real, parameter:: r2o3 = 2./3. + real, parameter:: r4o3 = 4./3. + real gak(km) + real bet + integer i, k + +!------------------------------------------------ +! Optimized coding for uniform grid: SJL Apr 2007 +!------------------------------------------------ + + if ( id==1 ) then + do i=i1,i2 + qe(i,1) = r4o3*q1(i,1) + r2o3*q1(i,2) + enddo + else + do i=i1,i2 + qe(i,1) = 1.E30 + enddo + endif + + gak(1) = 7./3. + do k=2,km + gak(k) = 1. / (4. - gak(k-1)) + do i=i1,i2 + qe(i,k) = (3.*(q1(i,k-1) + q1(i,k)) - qe(i,k-1)) * gak(k) + enddo + enddo + + bet = 1. / (1.5 - 3.5*gak(km)) + do i=i1,i2 + qe(i,km+1) = (4.*q1(i,km) + q1(i,km-1) - 3.5*qe(i,km)) * bet + enddo + + do k=km,1,-1 + do i=i1,i2 + qe(i,k) = qe(i,k) - gak(k)*qe(i,k+1) + enddo + enddo + + end subroutine edge_scalar + + + + subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) +! Optimized for wind profile reconstruction: + integer, intent(in):: i1, i2, j1, j2 + integer, intent(in):: j, km + integer, intent(in):: limiter + logical, intent(in):: uniform_grid + real, intent(in):: dp0(km) + real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2 + real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e +!----------------------------------------------------------------------- + real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values + real gak(km) + real bet, r2o3, r4o3 + real g0, gk, xt1, xt2, a_bot + integer i, k + + if ( uniform_grid ) then +!------------------------------------------------ +! Optimized coding for uniform grid: SJL Apr 2007 +!------------------------------------------------ + r2o3 = 2./3. + r4o3 = 4./3. + do i=i1,i2 + qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2) + qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2) + enddo + + gak(1) = 7./3. + do k=2,km + gak(k) = 1. / (4. - gak(k-1)) + do i=i1,i2 + qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k) + qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k) + enddo + enddo + + bet = 1. / (1.5 - 3.5*gak(km)) + do i=i1,i2 + qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet + qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1) + enddo + enddo + else +! Assuming grid varying in vertical only + g0 = dp0(2) / dp0(1) + xt1 = 2.*g0*(g0+1. ) + bet = g0*(g0+0.5) + do i=i1,i2 + qe1(i,1) = ( xt1*q1(i,j,1) + q1(i,j,2) ) / bet + qe2(i,1) = ( xt1*q2(i,j,1) + q2(i,j,2) ) / bet + gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet + enddo + + do k=2,km + gk = dp0(k-1) / dp0(k) + do i=i1,i2 + bet = 2. + 2.*gk - gam(i,k-1) + qe1(i,k) = ( 3.*(q1(i,j,k-1)+gk*q1(i,j,k)) - qe1(i,k-1) ) / bet + qe2(i,k) = ( 3.*(q2(i,j,k-1)+gk*q2(i,j,k)) - qe2(i,k-1) ) / bet + gam(i,k) = gk / bet + enddo + enddo + + a_bot = 1. + gk*(gk+1.5) + xt1 = 2.*gk*(gk+1.) + do i=i1,i2 + xt2 = gk*(gk+0.5) - a_bot*gam(i,km) + qe1(i,km+1) = ( xt1*q1(i,j,km) + q1(i,j,km-1) - a_bot*qe1(i,km) ) / xt2 + qe2(i,km+1) = ( xt1*q2(i,j,km) + q2(i,j,km-1) - a_bot*qe2(i,km) ) / xt2 + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gam(i,k)*qe2(i,k+1) + enddo + enddo + endif + +!------------------ +! Apply constraints +!------------------ + if ( limiter/=0 ) then ! limit the top & bottom winds + do i=i1,i2 +! Top + if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. + if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0. +! Surface: + if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. + if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0. + enddo + endif + + do k=1,km+1 + do i=i1,i2 + q1e(i,j,k) = qe1(i,k) + q2e(i,j,k) = qe2(i,k) + enddo + enddo + + end subroutine edge_profile + + subroutine edge_profile_0grad(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) +! Optimized for wind profile reconstruction: +! Added this option by Henry Juang and Xiaqiong Zhou 1/21/2021 + integer, intent(in):: i1, i2, j1, j2 + integer, intent(in):: j, km + integer, intent(in):: limiter + logical, intent(in):: uniform_grid + real, intent(in):: dp0(km) + real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2 + real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e +!----------------------------------------------------------------------- + real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values + real gak(km) + real bet, r2o3, r4o3 + real g0, gk, xt1, xt2, a_bot + integer i, k + + if ( uniform_grid ) then +!------------------------------------------------ +! Optimized coding for uniform grid: SJL Apr 2007 +!------------------------------------------------ + r2o3 = 2./3. + r4o3 = 4./3. + do i=i1,i2 + qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2) + qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2) + enddo + + gak(1) = 7./3. + do k=2,km + gak(k) = 1. / (4. - gak(k-1)) + do i=i1,i2 + qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k) + qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k) + enddo + enddo + + bet = 1. / (1.5 - 3.5*gak(km)) + do i=i1,i2 + qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet + qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1) + enddo + enddo + else +! Assuming grid varying in vertical only + g0 = dp0(1) / dp0(2) + bet = 1.5 + 2.*g0 + do i=i1,i2 + qe1(i,2) = 3.*( 0.5*q1(i,j,1) + g0*q1(i,j,2) ) / bet + qe2(i,2) = 3.*( 0.5*q2(i,j,1) + g0*q2(i,j,2) ) / bet + gam(i,1) = g0/bet + enddo + + +! for k=2,km + do k=2,km-1 + gk = dp0(k) / dp0(k+1) + do i=i1,i2 + bet = 2. + 2.*gk - gam(i,k-1) + qe1(i,k+1) = ( 3.*(q1(i,j,k)+gk*q1(i,j,k+1)) - qe1(i,k) ) / bet + qe2(i,k+1) = ( 3.*(q2(i,j,k)+gk*q2(i,j,k+1)) - qe2(i,k) ) / bet + gam(i,k) = gk / bet + enddo + enddo + +!km+1 + do i=i1,i2 + bet = 2.- gam(i,km-1) + qe1(i,km+1) = ( 3.*q1(i,j,km) - qe1(i,km) ) / bet + qe2(i,km+1) = ( 3.*q2(i,j,km) - qe2(i,km) ) / bet + enddo + + do i=i1,i2 + do k=km,2,-1 + qe1(i,k) = qe1(i,k) - gam(i,k-1)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gam(i,k-1)*qe2(i,k+1) + enddo + qe1(i,1)=1.5*q1(i,j,1)-0.5*qe1(i,2) + qe2(i,1)=1.5*q2(i,j,1)-0.5*qe2(i,2) + enddo + + endif + +!------------------ +! Apply constraints +!------------------ + if ( limiter/=0 ) then ! limit the top & bottom winds + do i=i1,i2 +! Top + if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. + if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0. +! Surface: + if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. + if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0. + enddo + endif + + do k=1,km+1 + do i=i1,i2 + q1e(i,j,k) = qe1(i,k) + q2e(i,j,k) = qe2(i,k) + enddo + enddo + + end subroutine edge_profile_0grad + +!TODO LMH 25may18: do not need delz defined on full compute domain; pass appropriate BCs instead + subroutine nh_bc(ptop, grav, kappa, cp, delp, delzBC, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + npx, npy, npz, bounded_domain, pkc_pertn, computepk3, fullhalo, bd) + + !INPUT: delp, delz (BC), pt + !OUTPUT: gz, pkc, pk3 (optional) + integer, intent(IN) :: npx, npy, npz + logical, intent(IN) :: pkc_pertn, computepk3, fullhalo, bounded_domain + real, intent(IN) :: ptop, kappa, cp, grav, BC_step, BC_split + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(IN) :: phis(bd%isd:bd%ied,bd%jsd:bd%jed) + real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: pt, delp + type(fv_nest_BC_type_3d), intent(IN) :: delzBC +#ifdef MULTI_GASES + real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,*):: q +#endif +#ifdef USE_COND + real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: q_con +#ifdef MOIST_CAPPA + real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: cappa +#endif +#endif + real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1):: gz, pkc, pk3 + + integer :: i,j,k + real :: gama !'gamma' + real :: ptk, rgrav, rkap, peln1, rdg + + integer :: istart, iend + + integer :: is, ie, js, je + integer :: isd, ied, jsd, jed + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + + if (.not. bounded_domain) return + + if (is == 1) then + + call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%west_t0, delzBC%west_t1, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + pkc_pertn, computepk3, isd, ied, isd, 0, isd, 0, jsd, jed, jsd, jed, npz) + + endif + + if (ie == npx-1) then + + call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%east_t0, delzBC%east_t1, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + pkc_pertn, computepk3, isd, ied, npx, ied, npx, ied, jsd, jed, jsd, jed, npz) + + endif + + if (is == 1) then + istart = is + else + istart = isd + end if + if (ie == npx-1) then + iend = ie + else + iend = ied + end if + + if (js == 1) then + + call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%south_t0, delzBC%south_t1, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, jsd, 0, npz) + + end if + + if (je == npy-1) then + + call nh_BC_k(ptop, grav, kappa, cp, delp, delzBC%north_t0, delzBC%north_t1, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, npy, jed, npz) + endif + +end subroutine nh_bc + +subroutine nh_BC_k(ptop, grav, kappa, cp, delp, delzBC_t0, delzBC_t1, pt, phis, & +#ifdef MULTI_GASES + q , & +#endif +#ifdef USE_COND + q_con, & +#ifdef MOIST_CAPPA + cappa, & +#endif +#endif + pkc, gz, pk3, & + BC_step, BC_split, & + pkc_pertn, computepk3, isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz) + + integer, intent(IN) :: isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz + real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: delzBC_t0, delzBC_t1 + real, intent(IN) :: BC_step, BC_split + + logical, intent(IN) :: pkc_pertn, computepk3 + real, intent(IN) :: ptop, kappa, cp, grav + real, intent(IN) :: phis(isd:ied,jsd:jed) + real, intent(IN), dimension(isd:ied,jsd:jed,npz):: pt, delp +#ifdef MULTI_GASES + real, intent(IN), dimension(isd:ied,jsd:jed,npz,*):: q +#endif +#ifdef USE_COND + real, intent(IN), dimension(isd:ied,jsd:jed,npz):: q_con +#ifdef MOIST_CAPPA + real, intent(INOUT), dimension(isd:ied,jsd:jed,npz):: cappa +#endif +#endif + real, intent(INOUT), dimension(isd:ied,jsd:jed,npz+1):: gz, pkc, pk3 + + integer :: i,j,k + real :: gama !'gamma' + real :: ptk, rgrav, rkap, peln1, rdg, denom + + real, dimension(istart:iend, npz+1, jstart:jend ) :: pe, peln +#ifdef USE_COND + real, dimension(istart:iend, npz+1 ) :: peg, pelng +#endif + real, dimension(istart:iend, npz) :: gam, bb, dd, pkz + real, dimension(istart:iend, npz-1) :: g_rat + real, dimension(istart:iend) :: bet + real :: pm, delz_int + + + real :: pealn, pebln, rpkz +#ifdef MULTI_GASES + real gamax +#endif + rgrav = 1./grav + gama = 1./(1.-kappa) + ptk = ptop ** kappa + rkap = 1./kappa + peln1 = log(ptop) + rdg = - rdgas * rgrav + denom = 1./BC_split + + do j=jstart,jend + + !GZ + do i=istart,iend + gz(i,j,npz+1) = phis(i,j) + enddo + do k=npz,1,-1 + do i=istart,iend + delz_int = (delzBC_t0(i,j,k)*(BC_split-BC_step) + BC_step*delzBC_t1(i,j,k))*denom + gz(i,j,k) = gz(i,j,k+1) - delz_int*grav + enddo + enddo + + !Hydrostatic interface pressure + do i=istart,iend + pe(i,1,j) = ptop + peln(i,1,j) = peln1 +#ifdef USE_COND + peg(i,1) = ptop + pelng(i,1) = peln1 +#endif + enddo + do k=2,npz+1 + do i=istart,iend + pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1) + peln(i,k,j) = log(pe(i,k,j)) +#ifdef USE_COND + peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1)) + pelng(i,k) = log(peg(i,k)) +#endif + enddo + enddo + + !Perturbation nonhydro layer-mean pressure (NOT to the kappa) + do k=1,npz + do i=istart,iend + delz_int = (delzBC_t0(i,j,k)*(BC_split-BC_step) + BC_step*delzBC_t1(i,j,k))*denom + + !Full p +#ifdef MOIST_CAPPA + pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz_int*pt(i,j,k))) +#else +#ifdef MULTI_GASES + gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:))) + pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k))) +#else + pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k))) +#endif +#endif + !hydro +#ifdef USE_COND + pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k)) +#else + pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j)) +#endif + !Remove hydro cell-mean pressure + pkz(i,k) = pkz(i,k) - pm + enddo + enddo + + !pressure solver + do k=1,npz-1 + do i=istart,iend + g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1) + bb(i,k) = 2.*(1. + g_rat(i,k)) + dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1)) + enddo + enddo + + do i=istart,iend + bet(i) = bb(i,1) + pkc(i,j,1) = 0. + pkc(i,j,2) = dd(i,1)/bet(i) + bb(i,npz) = 2. + dd(i,npz) = 3.*pkz(i,npz) + enddo + do k=2,npz + do i=istart,iend + gam(i,k) = g_rat(i,k-1)/bet(i) + bet(i) = bb(i,k) - gam(i,k) + pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i) + enddo + enddo + do k=npz,2,-1 + do i=istart,iend + pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1) +#ifdef NHNEST_DEBUG + if (abs(pkc(i,j,k)) > 1.e5) then + print*, mpp_pe(), i,j,k, 'PKC: ', pkc(i,j,k) + endif +#endif + enddo + enddo + + + enddo + + do j=jstart,jend + + if (.not. pkc_pertn) then + do k=npz+1,1,-1 + do i=istart,iend + pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j) + enddo + enddo + endif + + !pk3 if necessary; doesn't require condenstate loading calculation + if (computepk3) then + do i=istart,iend + pk3(i,j,1) = ptk + enddo + do k=2,npz+1 + do i=istart,iend + pk3(i,j,k) = exp(kappa*log(pe(i,k,j))) + enddo + enddo + endif + + enddo + +end subroutine nh_BC_k + + +end module nh_utils_mod From 4a27a041372ad3fd641a3ddf19937784f452a11c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 26 Mar 2021 17:30:28 -0600 Subject: [PATCH 2/2] Workaround to read the name of the CCPP suite from its current location (section atmos_model_nml) in the input namelist --- driver/fvGFS/atmosphere.F90 | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 3145998c9..4e53d04b4 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -293,6 +293,24 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) integer :: nthreads, ierr integer :: nlunit = 9999 character (len = 64) :: fn_nml = 'input.nml' + ! DH* 20210326 + ! This is a temporary workaround until the implementation + ! of the generic interface to call GFDL or CCPP physics is + ! completed. We need the name of the CCPP suite here in order + ! to run the adiabatic init with fast physics turned on. All + ! other vaiables are ignored (set to the same default value + ! as in fv3atm's atmos_model.F90. + integer :: io + integer :: blocksize = 1 + logical :: chksum_debug = .false. + logical :: dycore_only = .false. + logical :: debug = .false. + logical :: sync = .false. + integer, parameter :: maxhr = 4096 + real, dimension(maxhr) :: fdiag = 0. + real :: fhmax=384.0, fhmaxhf=120.0, fhout=3.0, fhouthf=1.0,avg_max_length=3600. + namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, avg_max_length + ! *DH 20210326 !For regional a_step = 0 @@ -422,6 +440,25 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics) + ! DH* 20210326 + ! First, read atmos_model_nml namelist section - this is a workaround to avoid + ! unnecessary additional changes to the input namelists, in anticipation of the + ! implementation of a generic interface for GFDL and CCPP fast physics soon +#ifdef INTERNAL_FILE_NML + read(input_nml_file, nml=atmos_model_nml, iostat=io) + ierr = check_nml_error(io, 'atmos_model_nml') +#else + unit = open_namelist_file ( ) + ierr=1 + do while (ierr /= 0) + read (unit, nml=atmos_model_nml, iostat=io, end=10) + ierr = check_nml_error(io,'atmos_model_nml') + enddo +10 call close_file (unit) +#endif + !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'" + ! *DH 20210326 + ! For fast physics running over the entire domain, block ! and thread number are not used; set to safe values cdata%blk_no = 1