diff --git a/CMakeLists.txt b/CMakeLists.txt index d7d4eb52d..888a622d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -59,6 +59,7 @@ endif() list(APPEND model_srcs model/a2b_edge.F90 model/multi_gases.F90 + model/molecular_diffusion.F90 model/boundary.F90 model/dyn_core.F90 model/fv_arrays.F90 @@ -169,6 +170,7 @@ add_library(FV3::fv3 ALIAS fv3) set_property(SOURCE model/nh_utils.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${FAST}") set_property(SOURCE model/fv_mapz.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${FAST}") +set_property(SOURCE tools/fv_diagnostics.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "-O0") set_target_properties(fv3 PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include/fv3) target_include_directories(fv3 INTERFACE $ diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index 414d91d8f..da7a176bd 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -122,9 +122,7 @@ module fv_nggps_diags_mod logical :: use_wrtgridcomp_output=.false. integer :: sphum, liq_wat, ice_wat !< GFDL physics integer :: rainwat, snowwat, graupel - real :: vrange(2) = (/ -330., 330. /) !< winds - real :: wrange(2) = (/ -100., 100. /) !< vertical wind - real :: trange(2) = (/ 100., 350. /) !< temperature + real :: vrange(2), wrange(2), trange(2) real :: skrange(2) = (/ -10000000.0, 10000000.0 /) !< dissipation estimate for SKEB ! file name @@ -181,6 +179,20 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) kstt_tracer(:) = 0 kend_tracer(:) = 0 + if ( Atm(n)%flagstruct%molecular_diffusion ) then + vrange = (/ -830., 830. /) !< winds for wam + wrange = (/ -300., 300. /) !< vertical wind for wam + trange = (/ 100., 3000. /) !< temperature for wam + else + vrange = (/ -330., 330. /) !< winds + wrange = (/ -100., 100. /) !< vertical wind +#ifdef MULTI_GASES + trange = (/ 100., 3000. /) !< temperature +#else + trange = (/ 100., 350. /) !< temperature +#endif + endif + if (Atm(n)%flagstruct%write_3d_diags) then !------------------- ! A grid winds (lat-lon) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index ee22ef355..2e6173f72 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -111,7 +111,9 @@ module dyn_core_mod use fv_mp_mod, only: is_master use fv_mp_mod, only: start_group_halo_update, complete_group_halo_update use fv_mp_mod, only: group_halo_update_type - use sw_core_mod, only: c_sw, d_sw + use molecular_diffusion_mod, & + only: md_time, md_layers, md_consv_te, md_tadj_layers + use sw_core_mod, only: c_sw, d_sw, d_md use a2b_edge_mod, only: a2b_ord2, a2b_ord4 use nh_core_mod, only: Riem_Solver3, Riem_Solver_C, update_dz_c, update_dz_d, nh_bc use tp_core_mod, only: copy_corners @@ -139,7 +141,7 @@ module dyn_core_mod use test_cases_mod, only: test_case, case9_forcing1, case9_forcing2 #endif #ifdef MULTI_GASES - use multi_gases_mod, only: virqd, vicvqd + use multi_gases_mod, only: virqd, vicpqd, vicvqd, virq, vicvq #endif use fv_regional_mod, only: dump_field, exch_uv, H_STAGGER, U_STAGGER, V_STAGGER use fv_regional_mod, only: a_step, p_step, k_step, n_step @@ -199,7 +201,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real, intent(inout), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz):: v !< D grid meridional wind (m/s) real, intent(inout) :: w( bd%isd:,bd%jsd:,1:) !< vertical vel. (m/s) real, intent(inout) :: delz(bd%is:,bd%js:,1:) !< delta-height (m, negative) - real, intent(inout) :: cappa(bd%isd:,bd%jsd:,1:) !< moist kappa + real, intent(inout) :: cappa(bd%isd:bd%ied,bd%jsd:bd%jed,1:npz) !< moist kappa #ifdef MULTI_GASES real, intent(inout) :: kapad(bd%isd:bd%ied,bd%jsd:bd%jed,1:npz) !< multi_gases kappa #endif @@ -1222,6 +1224,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, ! Prevent accumulation of rounding errors at overlapped domain edges: call mpp_get_boundary(u, v, domain, ebuffery=ebuffer, & nbufferx=nbuffer, gridtype=DGRID_NE ) + call timing_off('COMM_TOTAL') !$OMP parallel do default(none) shared(is,ie,js,je,npz,u,nbuffer,v,ebuffer) do k=1,npz do i=is,ie @@ -1234,6 +1237,16 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, endif + if ( flagstruct%molecular_diffusion .and. md_time ) then +! ----------------------------------------------------- +! direct explicit molecular diffusion +! ----------------------------------------------------- + call molecular_diffusion_run(u, v, w, delp, pt, pkz, cappa, q, bd, & + gridstruct, flagstruct, domain, i_pack, npx, npy, npz, nq, dt, it, akap, zvir, cv_air) + endif +! ------------------------------------------------- + + call timing_on('COMM_TOTAL') #ifndef ROT3 if ( it/=n_split) & call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE) @@ -1243,9 +1256,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, #ifdef SW_DYNAMICS endif #endif - if ( gridstruct%nested ) then + if ( gridstruct%nested ) then neststruct%nest_timestep = neststruct%nest_timestep + 1 - endif + endif #ifdef SW_DYNAMICS #else @@ -1298,13 +1311,12 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, if (gridstruct%nested) then - #ifndef SW_DYNAMICS - if (.not. hydrostatic) then + if (.not. hydrostatic) then call nested_grid_BC_apply_intT(w, & 0, 0, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), & neststruct%w_BC, bctype=neststruct%nestbctype ) - end if + end if #endif SW_DYNAMICS call nested_grid_BC_apply_intT(u, & 0, 1, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), & @@ -1313,9 +1325,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, 1, 0, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), & neststruct%v_BC, bctype=neststruct%nestbctype ) - end if + end if - if (flagstruct%regional) then + if (flagstruct%regional) then #ifndef SW_DYNAMICS if (.not. hydrostatic) then @@ -1340,7 +1352,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, reg_bc_update_time ) call mpp_update_domains(u, v, domain, gridtype=DGRID_NE) - end if + + endif if ( do_diag_debug_dyn ) then call debug_column_dyn( pt, delp, delz, u, v, w, q, heat_source, cappa, akap, & @@ -2784,5 +2797,181 @@ subroutine gz_bc(gz,delzBC,bd,npx,npy,npz,step,split) end subroutine gz_bc + subroutine molecular_diffusion_run(u,v,w,delp,pt,pkz,cappa,q,bd, & + gridstruct,flagstruct,domain,i_pack,npx,npy,npz, & + nq,dt,it,akap,zvir,cv_air) + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(inout):: pkz(bd%is:bd%ie,bd%js:bd%je,npz) + real, intent(inout):: pt(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real, intent(inout):: delp(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real, intent(inout):: cappa(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real, intent(inout):: w(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real, intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz) + real, intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz) + real, intent(inout):: q( bd%isd:bd%ied,bd%jsd:bd%jed,npz, nq) + + integer, intent(in):: npx,npy,npz,nq,it + real, intent(in)::akap,zvir,dt,cv_air + type(fv_grid_type), intent(INOUT), target :: gridstruct + type(fv_flags_type), intent(IN), target :: flagstruct + type(domain2d), intent(inout) :: domain + type(group_halo_update_type), intent(inout) :: i_pack(*) + + real :: pkzf(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real, dimension (bd%isd:bd%ied,bd%jsd:bd%jed) :: p, t, e + integer :: i, j, k + 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 + +! ----------------------------------------------------- +! ------- update halo to prepare for diffusion ------- + pkzf = 0.0 + do k=1,npz + do j=js,je + pkzf(is:ie,j,k) = pkz(is:ie,j,k) + enddo + enddo + + call timing_on('COMM_TOTAL') + call start_group_halo_update(i_pack(1),delp, domain, complete=.false.) + call start_group_halo_update(i_pack(1), pt, domain, complete=.true.) + call start_group_halo_update(i_pack(2),pkzf, domain) + call start_group_halo_update(i_pack(7), w, domain) + call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE) + if ( nq > 0 ) then + call timing_on('COMM_TRACER') + call start_group_halo_update(i_pack(10), q, domain) + call timing_off('COMM_TRACER') + endif +#ifdef MOIST_CAPPA + call start_group_halo_update(i_pack(12), cappa, domain) +#endif + + call complete_group_halo_update(i_pack(1), domain) ! delp. pt + call complete_group_halo_update(i_pack(2), domain) ! pkzf + call complete_group_halo_update(i_pack(7), domain) ! w + call complete_group_halo_update(i_pack(8), domain) + if ( nq>0 ) then + call timing_on('COMM_TRACER') + call complete_group_halo_update(i_pack(10), domain) + call timing_off('COMM_TRACER') + endif +#ifdef MOIST_CAPPA + call complete_group_halo_update(i_pack(12), domain) +#endif + call timing_off('COMM_TOTAL') + + if( flagstruct%nord>0 .and. (.not. (flagstruct%regional))) then + i=mod(it-1,2)+1 ! alternatively to avoid bias + do k=1,npz + call copy_corners(pt(isd,jsd,k), npx, npy, i, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(pkzf(isd,jsd,k), npx, npy, i, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(cappa(isd,jsd,k), npx, npy, i, gridstruct%nested,bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + enddo + endif + call timing_on('d_md') + +!$OMP parallel do default(none) shared(npz,flagstruct,gridstruct,bd, & +!$OMP it,dt,is,ie,js,je,isd,ied,jsd,jed, & +!$OMP pt,u,v,w,q,pkz,pkzf,cappa,akap,nq, & +!$OMP zvir,cv_air,md_layers,md_consv_te) & +!$OMP private(k,i,j,p,t,e) +! ---------------- + do k=1, md_layers +! ---------------- + +! ------- prepare p and t for molecular diffusion coefficients + + do j=jsd,jed + do i=isd,ied + t(i,j) = pt(i,j,k) * pkzf(i,j,k) +#ifdef MULTI_GASES + t(i,j) = t(i,j) / virq(q(i,j,k,:)) +#else + t(i,j) = t(i,j) / (1+zvir*q(i,j,k,1)) +#endif +#ifdef MOIST_CAPPA + p(i,j) = exp( log(pkzf(i,j,k)) / cappa(i,j,k) ) +#else +#ifdef MULTI_GASES + p(i,j) = exp( log(pkzf(i,j,k)) / & + (akap*virqd(q(i,j,k,:))/vicpqd(q(i,j,k,:))) ) +#else + p(i,j) = exp( log(pkzf(i,j,k)) / akap ) +#endif +#endif + if( md_consv_te .gt. 0.0 ) & + e(i,j) = 0.5*(w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*( & + u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -& + (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))* & + gridstruct%cosa_s(i,j))) + enddo + enddo + +! compute molecular diffusion with implicit time and dimensional splits + + call d_md( p(isd,jsd), t(isd,jsd), & + u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), q, & + it, nq, k, npz, dt, & + gridstruct, flagstruct, bd) + + do j=js,je + do i=is,ie + if( md_consv_te .gt. 0.0 ) then + e(i,j) = e(i,j) - & + 0.5*(w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*( & + u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -& + (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))* & + gridstruct%cosa_s(i,j))) +#ifdef MOIST_CAPPA + t(i,j) = t(i,j) + e(i,j) / & +#ifdef MULTI_GASES + (rdgas* virq(q(i,j,k,:)) *(1./cappa(i,j,k)-1.)) +#else + (rdgas*(1+zvir*q(i,j,k,1))*(1./cappa(i,j,k)-1.)) +#endif + endif + pkz(i,j,k) = exp( log(p(i,j)) * cappa(i,j,k) ) +#else +#ifdef MULTI_GASES + t(i,j) = t(i,j) + e(i,j) / (cv_air*vicvqd(q(i,j,k,:))) + endif + pkz(i,j,k) = exp( log(p(i,j)) * & + (akap*virqd(q(i,j,k,:))/vicpqd(q(i,j,k,:))) ) +#else + t(i,j) = t(i,j) + e(i,j)/cv_air + endif + pkz(i,j,k) = exp( log(p(i,j)) * akap ) +#endif +#endif +#ifdef MULTI_GASES + t(i,j) = t(i,j) * virq(q(i,j,k,:)) +#else + t(i,j) = t(i,j) * (1+zvir*q(i,j,k,1)) +#endif + pt(i,j,k) = t(i,j) / pkz(i,j,k) + enddo + enddo + +! ------------------------------------------------- + enddo ! k loop of 2d molecular diffusion +! ------------------------------------------------- + return + end subroutine molecular_diffusion_run end module dyn_core_mod diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 8590cade1..a353e0a65 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -88,6 +88,17 @@ module fv_arrays_mod real, allocatable, dimension(:,:) :: rdxc, rdyc real, allocatable, dimension(:,:) :: rdxa, rdya +! MOLECULAR_DIFFUSION + real(kind=R_GRID), allocatable, dimension(:,:) :: area_u_64, area_v_64 + real(kind=R_GRID), allocatable, dimension(:,:) :: dx6_64, dy6_64 + real, allocatable, dimension(:,:) :: area_u, area_v + real, allocatable, dimension(:,:) :: rarea_u, rarea_v + real, allocatable, dimension(:,:) :: dx6, dy6 + real, allocatable, dimension(:,:) :: rdx6, rdy6 + real, allocatable, dimension(:,:) :: sina_6 + real, allocatable, dimension(:,:) :: delu_6, delv_6 + real, allocatable, dimension(:,:) :: delu_5, delv_5 + ! Scalars: real(kind=R_GRID), allocatable :: edge_s(:) real(kind=R_GRID), allocatable :: edge_n(:) @@ -671,6 +682,10 @@ module fv_arrays_mod !< original value before entering the physics; a value of 0.7 roughly !< causes the energy fixer to compensate for the amount of energy changed !< by the physics in GFDL HiRAM or AM3. + real :: tau_w = 0. !< Time scale (in days) for Rayleigh friction applied to vertical winds + !< This option allows the vertical and horizontal winds use different time + !< scales for Rayleigh friction. The default value is 0.0, then tau_w=tau, + !< the same time scale appiled to horizontal and vertical winds. real :: tau = 0. !< Time scale (in days) for Rayleigh friction applied to horizontal !< and vertical winds; lost kinetic energy is converted to heat, except !< on nested grids. The default value is 0.0, which disables damping. @@ -863,6 +878,8 @@ module fv_arrays_mod logical :: butterfly_effect = .false. !< Flip the least-significant-bit of the lowest level temperature !< at the center of the domain (the center of tile 1), if set to .true. !< The default value is .false. + logical :: molecular_diffusion = .false. !< Apply Whole Atmosphere Model (WAM) molecular diffusion + !< developed by Henry Juang real :: dz_min = 2 !< Minimum thickness depth to to enforce monotonicity of height to prevent blowup. !< 2 by default @@ -1609,6 +1626,26 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie allocate ( Atm%gridstruct% dya_64(isd_2d:ied_2d ,jsd_2d:jed_2d ) ) allocate ( Atm%gridstruct%rdya(isd_2d:ied_2d ,jsd_2d:jed_2d ) ) + if ( Atm%flagstruct%molecular_diffusion ) then + allocate ( Atm%gridstruct% area_u_64(isd_2d:ied_2d ,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct% area_v_64(isd_2d:ied_2d+1,jsd_2d:jed_2d ) ) + allocate ( Atm%gridstruct% dx6_64(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct% dy6_64(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct% area_u(isd_2d:ied_2d ,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct% area_v(isd_2d:ied_2d+1,jsd_2d:jed_2d ) ) + allocate ( Atm%gridstruct% dx6(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct% dy6(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct%rarea_u(isd_2d:ied_2d ,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct%rarea_v(isd_2d:ied_2d+1,jsd_2d:jed_2d ) ) + allocate ( Atm%gridstruct%rdx6(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct%rdy6(isd_2d:ied_2d+1,jsd_2d:jed_2d+1) ) + allocate ( Atm%gridstruct%sina_6(isd_2d:ied_2d,jsd_2d:jed_2d) ) + allocate ( Atm%gridstruct%delu_6(isd_2d:ied_2d,jsd_2d:jed_2d) ) + allocate ( Atm%gridstruct%delv_6(isd_2d:ied_2d,jsd_2d:jed_2d) ) + allocate ( Atm%gridstruct%delu_5(isd_2d:ied_2d,jsd_2d:jed_2d) ) + allocate ( Atm%gridstruct%delv_5(isd_2d:ied_2d,jsd_2d:jed_2d) ) + endif + allocate ( Atm%gridstruct%grid (isd_2d:ied_2d+1,jsd_2d:jed_2d+1,1:ndims_2d) ) allocate ( Atm%gridstruct%grid_64 (isd_2d:ied_2d+1,jsd_2d:jed_2d+1,1:ndims_2d) ) allocate ( Atm%gridstruct%agrid(isd_2d:ied_2d ,jsd_2d:jed_2d ,1:ndims_2d) ) @@ -1858,6 +1895,22 @@ subroutine deallocate_fv_atmos_type(Atm) deallocate ( Atm%gridstruct% dya ) deallocate ( Atm%gridstruct%rdya ) + if ( Atm%flagstruct%molecular_diffusion ) then + deallocate ( Atm%gridstruct% area_u ) + deallocate ( Atm%gridstruct% area_v ) + deallocate ( Atm%gridstruct%rarea_u ) + deallocate ( Atm%gridstruct%rarea_v ) + deallocate ( Atm%gridstruct% dx6 ) + deallocate ( Atm%gridstruct% dy6 ) + deallocate ( Atm%gridstruct%rdx6 ) + deallocate ( Atm%gridstruct%rdy6 ) + deallocate ( Atm%gridstruct%sina_6 ) + deallocate ( Atm%gridstruct%delu_6 ) + deallocate ( Atm%gridstruct%delv_6 ) + deallocate ( Atm%gridstruct%delu_5 ) + deallocate ( Atm%gridstruct%delv_5 ) + endif + deallocate ( Atm%gridstruct%grid ) deallocate ( Atm%gridstruct%agrid ) deallocate ( Atm%gridstruct%sina ) ! SIN(angle of intersection) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 7f900f23e..b58558119 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -158,6 +158,8 @@ module fv_control_mod use multi_gases_mod, only: multi_gases_init, & read_namelist_multi_gases_nml #endif + use molecular_diffusion_mod, only: molecular_diffusion_init, & + read_namelist_molecular_diffusion_nml implicit none private @@ -184,12 +186,15 @@ module fv_control_mod !------------------------------------------------------------------------------- - subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) + subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, & + nml_filename_in, skip_nml_read_in) type(fv_atmos_type), allocatable, intent(inout), target :: Atm(:) real, intent(in) :: dt_atmos integer, intent(OUT) :: this_grid logical, allocatable, intent(OUT) :: grids_on_this_pe(:) + character(len=32), optional, intent(in) :: nml_filename_in ! alternate nml + logical, optional, intent(in) :: skip_nml_read_in ! use previously loaded nml integer, intent(INOUT) :: p_split character(100) :: pe_list_name, errstring @@ -214,6 +219,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) real :: sdt integer :: unit, ens_root_pe, tile_id(1) + character(len=32) :: nml_filename = 'input.nml' + logical :: skip_nml_read = .false. + !!!!!!!!!! POINTERS FOR READING NAMELISTS !!!!!!!!!! !------------------------------------------ @@ -321,6 +329,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) real , pointer :: ke_bg real , pointer :: consv_te real , pointer :: tau + real , pointer :: tau_w real , pointer :: rf_cutoff logical , pointer :: filter_phys logical , pointer :: dwind_2d @@ -372,6 +381,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) logical , pointer :: nudge_qv real, pointer :: add_noise logical , pointer :: butterfly_effect + logical , pointer :: molecular_diffusion real, pointer :: dz_min integer, pointer :: psm_bc @@ -408,6 +418,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) call mp_assign_gid ens_root_pe = mpp_root_pe() + if (present(nml_filename_in)) nml_filename = nml_filename_in + if (present(skip_nml_read_in)) skip_nml_read = skip_nml_read_in + ! 1. read nesting namelists call read_namelist_nest_nml call read_namelist_fv_nest_nml @@ -482,9 +495,10 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) if (n > 1) then Atm(n)%nml_filename = 'input_'//trim(pe_list_name)//'.nml' else - Atm(n)%nml_filename = 'input.nml' +! Atm(n)%nml_filename = 'input.nml' + Atm(n)%nml_filename = trim(nml_filename) endif - if (.not. file_exists(Atm(n)%nml_filename)) then + if (.not. file_exists(Atm(n)%nml_filename) .and. .not. skip_nml_read) then call mpp_error(FATAL, "Could not find nested grid namelist "//Atm(n)%nml_filename) endif enddo @@ -550,7 +564,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) else Atm(this_grid)%nml_filename = '' endif - call read_input_nml(Atm(this_grid)%nml_filename) !re-reads into internal namelist + if (.not. skip_nml_read) then + call read_input_nml(Atm(this_grid)%nml_filename) !re-reads into internal namelist + endif #endif call read_namelist_fv_grid_nml call read_namelist_fv_core_nml(Atm(this_grid)) ! do options processing here too? @@ -558,6 +574,10 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) call read_namelist_multi_gases_nml(Atm(this_grid)%nml_filename, & Atm(this_grid)%flagstruct%ncnst, Atm(this_grid)%flagstruct%nwat) #endif + if ( Atm(this_grid)%flagstruct%molecular_diffusion ) then + call read_namelist_molecular_diffusion_nml(Atm(this_grid)%nml_filename, & + Atm(this_grid)%flagstruct%ncnst, Atm(this_grid)%flagstruct%nwat) + endif call read_namelist_test_case_nml(Atm(this_grid)%nml_filename) call mpp_get_current_pelist(Atm(this_grid)%pelist, commID=commID) ! for commID call mp_start(commID,halo_update_type) @@ -865,6 +885,7 @@ subroutine set_namelist_pointers(Atm) ke_bg => Atm%flagstruct%ke_bg consv_te => Atm%flagstruct%consv_te tau => Atm%flagstruct%tau + tau_w => Atm%flagstruct%tau_w rf_cutoff => Atm%flagstruct%rf_cutoff filter_phys => Atm%flagstruct%filter_phys dwind_2d => Atm%flagstruct%dwind_2d @@ -916,6 +937,7 @@ subroutine set_namelist_pointers(Atm) nudge_qv => Atm%flagstruct%nudge_qv add_noise => Atm%flagstruct%add_noise butterfly_effect => Atm%flagstruct%butterfly_effect + molecular_diffusion => Atm%flagstruct%molecular_diffusion dz_min => Atm%flagstruct%dz_min psm_bc => Atm%flagstruct%psm_bc a2b_ord => Atm%flagstruct%a2b_ord @@ -1029,13 +1051,13 @@ subroutine read_namelist_fv_core_nml(Atm) dry_mass, grid_type, do_Held_Suarez, do_reed_physics, reed_cond_only, & consv_te, fill, filter_phys, fill_dp, fill_wz, fill_gfs, consv_am, RF_fast, & range_warn, dwind_2d, inline_q, z_tracer, reproduce_sum, adiabatic, do_vort_damp, no_dycore, & - tau, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, & + tau, tau_w, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, & na_init, nudge_dz, hybrid_z, Make_NH, n_zs_filter, nord_zs_filter, full_zs_filter, reset_eta, & pnats, dnats, dnrts, a2b_ord, remap_t, p_ref, d2_bg_k1, d2_bg_k2, & 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, & - dz_min, psm_bc, nested, twowaynest, nudge_qv, & + molecular_diffusion, 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/fv_dynamics.F90 b/model/fv_dynamics.F90 index 794207930..2dc440c49 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -149,8 +149,10 @@ module fv_dynamics_mod use boundary_mod, only: nested_grid_BC_apply_intT use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_atmos_type, fv_nest_type, fv_diag_type, fv_grid_bounds_type, inline_mp_type use fv_nwp_nudge_mod, only: do_adiabatic_init + use time_manager_mod, only: get_time + #ifdef MULTI_GASES - use multi_gases_mod, only: virq, virqd, vicpqd + use multi_gases_mod, only: virq, vicpq, virqd, vicpqd #endif implicit none @@ -158,6 +160,7 @@ module fv_dynamics_mod logical :: pt_initialized = .false. logical :: bad_range = .false. real, allocatable :: rf(:) + real, allocatable :: rw(:) !< related to tau_w integer :: kmax=1 integer :: k_rf = 0 @@ -180,7 +183,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, & ak, bk, mfx, mfy, cx, cy, ze0, hybrid_z, & gridstruct, flagstruct, neststruct, idiag, bd, & - parent_grid, domain, diss_est, inline_mp, time_total) + parent_grid, domain, diss_est, inline_mp) use mpp_mod, only: FATAL, mpp_error use ccpp_static_api, only: ccpp_physics_timestep_init, & @@ -189,11 +192,13 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, use CCPP_data, only: cdata => cdata_tile use CCPP_data, only: CCPP_interstitial + use molecular_diffusion_mod, only: md_time, md_wait_sec, md_tadj_layers, & + thermosphere_adjustment + real, intent(IN) :: bdt !< Large time-step real, intent(IN) :: consv_te real, intent(IN) :: kappa, cp_air real, intent(IN) :: zvir, ptop - real, intent(IN), optional :: time_total integer, intent(IN) :: npx integer, intent(IN) :: npy @@ -270,6 +275,9 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #ifdef MULTI_GASES real, allocatable :: kapad(:,:,:) #endif + real, save :: time_offset = -1.00 + real t(bd%isd:bd%ied,bd%jsd:bd%jed) + real :: correct, aave_t1, ttsave, tdsave real:: akap, rdg, ph1, ph2, mdt, gam, amdt, u0 real:: recip_k_split,reg_bc_update_time integer:: kord_tracer(ncnst) @@ -284,6 +292,8 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, integer :: isd, ied, jsd, jed real :: dt2 integer :: ierr + real :: time_total + integer :: seconds, days ccpp_associate: associate( cappa => CCPP_interstitial%cappa, & dp1 => CCPP_interstitial%te0, & @@ -325,6 +335,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, call init_ijk_mem(isd,ied, jsd,jed, npz, kapad, kappa) #endif + if (flagstruct%molecular_diffusion ) then + call get_time(fv_time, seconds, days) + time_total = seconds + 86400. * days + if( time_offset .lt. 0.0) time_offset = time_total + endif + !We call this BEFORE converting pt to virtual potential temperature, !since we interpolate on (regular) temperature rather than theta. if (gridstruct%nested .or. ANY(neststruct%child_grids)) then @@ -525,9 +541,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, ! call Ray_fast(abs(dt), npx, npy, npz, pfull, flagstruct%tau, u, v, w, & ! dp_ref, ptop, hydrostatic, flagstruct%rf_cutoff, bd) ! else - call Rayleigh_Super(abs(bdt), npx, npy, npz, ks, pfull, phis, flagstruct%tau, u, v, w, pt, & + call Rayleigh_Super(abs(bdt), npx, npy, npz, ks, pfull, phis, flagstruct%tau, flagstruct%tau_w, u, v, w, pt, & +#ifdef MULTI_GASES + q, ncnst, & +#endif ua, va, delz, gridstruct%agrid, cp_air, rdgas, ptop, hydrostatic, & - .not. gridstruct%bounded_domain, flagstruct%rf_cutoff, gridstruct, domain, bd) + .not. gridstruct%bounded_domain, flagstruct%molecular_diffusion, consv_te, flagstruct%rf_cutoff, gridstruct, domain, bd) ! endif else call Rayleigh_Friction(abs(bdt), npx, npy, npz, ks, pfull, flagstruct%tau, u, v, w, pt, & @@ -751,8 +770,15 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, inline_mp, flagstruct%c2l_ord, bd, flagstruct%fv_debug, & flagstruct%moist_phys) + if ( flagstruct%molecular_diffusion ) then +! do thermosphere adjustment if it is turned on and at last_step. + if( md_tadj_layers .gt.0 .and. md_time .and. last_step ) then + call thermosphere_adjustment(domain,gridstruct,npz,bd,ng,pt) + endif ! md_tadj_layers>0 and md_time and last_step + endif + if ( flagstruct%fv_debug ) then - if (is_master()) write(*,'(A, I3, A1, I3)') 'finished k_split ', n_map, '/', k_split + if (is_master()) write(*,'(A, I3, A1, I3)') 'finished k_split ', n_map, '/', k_split call prt_mxm('T_dyn_a4', pt, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain) if (sphum > 0) call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) if (liq_wat > 0) call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) @@ -802,6 +828,13 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #endif enddo ! n_map loop call timing_off('FV_DYN_LOOP') + if ( flagstruct%molecular_diffusion ) then + if( .not. md_time .and. time_total - time_offset .gt. md_wait_sec ) then + md_time= .true. + if( is_master() ) write(*,*) 'Molecular diffusion is on with explicit scheme ' + endif + endif + if ( idiag%id_mdt > 0 .and. (.not.do_adiabatic_init) ) then ! Output temperature tendency due to inline moist physics: #ifdef __GFORTRAN__ @@ -961,6 +994,16 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, endif if ( flagstruct%range_warn ) then + if ( flagstruct%molecular_diffusion ) then + call range_check('UA_dyn', ua, is, ie, js, je, ng, npz,gridstruct%agrid,& + -880., 880., bad_range) + call range_check('VA_dyn', ua, is, ie, js, je, ng, npz,gridstruct%agrid,& + -880., 880., bad_range) + call range_check('TA_dyn', pt, is, ie, js, je, ng, npz,gridstruct%agrid,& + 150.,3350., bad_range) + call range_check('W_dyn', w, is, ie, js, je, ng, npz,gridstruct%agrid, & + -250., 250., bad_range) + else call range_check('UA_dyn', ua, is, ie, js, je, ng, npz, gridstruct%agrid, & -280., 280., bad_range, fv_time) call range_check('VA_dyn', ua, is, ie, js, je, ng, npz, gridstruct%agrid, & @@ -970,6 +1013,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, if ( .not. hydrostatic ) & call range_check('W_dyn', w, is, ie, js, je, ng, npz, gridstruct%agrid, & -50., 100., bad_range, fv_time) + endif endif ! Call CCPP timestep finalize @@ -1101,16 +1145,21 @@ end subroutine Rayleigh_fast - subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & + subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, tau_w, u, v, w, pt, & +#ifdef MULTI_GASES + q, ncnst, & +#endif ua, va, delz, agrid, cp, rg, ptop, hydrostatic, & - conserve, rf_cutoff, gridstruct, domain, bd) + conserve, molecular_diffusion, consv_te, rf_cutoff, gridstruct, domain, bd) real, intent(in):: dt real, intent(in):: tau !< time scale (days) + real, intent(in):: tau_w !< time scale (days) for w real, intent(in):: cp, rg, ptop, rf_cutoff real, intent(in), dimension(npz):: pm integer, intent(in):: npx, npy, npz, ks logical, intent(in):: hydrostatic - logical, intent(in):: conserve + logical, intent(in):: conserve, molecular_diffusion + real, intent(in):: consv_te type(fv_grid_bounds_type), intent(IN) :: bd real, intent(inout):: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) !< D grid zonal wind (m/s) real, intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz) !< D grid meridional wind (m/s) @@ -1119,31 +1168,48 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & real, intent(inout):: ua(bd%isd:bd%ied,bd%jsd:bd%jed,npz) ! real, intent(inout):: va(bd%isd:bd%ied,bd%jsd:bd%jed,npz) ! real, intent(inout):: delz(bd%is: ,bd%js: ,1: ) !< delta-height (m); non-hydrostatic only +#ifdef MULTI_GASES + integer, intent(in):: ncnst + real, intent(in) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,ncnst) ! +#endif real, intent(in) :: agrid(bd%isd:bd%ied, bd%jsd:bd%jed,2) real, intent(in) :: phis(bd%isd:bd%ied,bd%jsd:bd%jed) !< Surface geopotential (g*Z_surf) type(fv_grid_type), intent(IN) :: gridstruct type(domain2d), intent(INOUT) :: domain ! real, allocatable :: u2f(:,:,:) + real, allocatable :: w2f(:,:,:) !< related to tau_w real, parameter:: u0 = 60. !< scaling velocity real, parameter:: sday = 86400. - real rcv, tau0 + real rcp, rcv, tau0, tau1 integer i, j, k + logical convert 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 + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + rcp = 1. / cp rcv = 1. / (cp - rg) + convert=conserve + if ( molecular_diffusion )then + if ( consv_te>0 ) then + convert=.true. + else + convert=.false. + endif + endif + + if ( .not. RF_initialized ) then #ifdef HIWPP allocate ( u00(is:ie, js:je+1,npz) ) @@ -1164,23 +1230,30 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & #endif #ifdef SMALL_EARTH_TEST ! changed!!! tau0 = abs( tau ) + tau1 = abs( tau_w ) #else tau0 = abs( tau * sday ) + tau1 = abs( tau_w * sday ) #endif - if( is_master() ) write(6,*) 'Rayleigh_Super tau=',tau + if( tau_w .eq. 0.0 ) tau1 = tau0 + if( is_master() ) write(6,*) 'Rayleigh_Super in sec tau=',tau0,' tau_w=',tau1 allocate( rf(npz) ) + allocate( rw(npz) ) rf(:) = 0. + rw(:) = 0. do k=1, ks+1 if( is_master() ) write(6,*) k, 0.01*pm(k) enddo - if( is_master() ) write(6,*) 'Rayleigh friction E-folding time (days):' + if( is_master() ) write(6,*) 'Rayleigh_Super E-folding time (mb days):' do k=1, npz if ( pm(k) < rf_cutoff ) then - if ( tau < 0 ) then ! GSM formula + if ( tau < 0 ) then !< GSM formula rf(k) = dt/tau0*( log(rf_cutoff/pm(k)) )**2 + rw(k) = dt/tau1*( log(rf_cutoff/pm(k)) )**2 else rf(k) = dt/tau0*sin(0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop))**2 + rw(k) = dt/tau1*sin(0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop))**2 endif if( is_master() ) write(6,*) k, 0.01*pm(k), dt/(rf(k)*sday) kmax = k @@ -1194,32 +1267,40 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & call c2l_ord2(u, v, ua, va, gridstruct, npz, gridstruct%grid_type, bd, gridstruct%bounded_domain) allocate( u2f(isd:ied,jsd:jed,kmax) ) + allocate( w2f(isd:ied,jsd:jed,kmax) ) !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,hydrostatic,ua,va,agrid, & -!$OMP u2f,rf,w) +!$OMP u2f,w2f,rf,rw,w) do k=1,kmax if ( pm(k) < rf_cutoff ) then u2f(:,:,k) = 1. / (1.+rf(k)) + w2f(:,:,k) = 1. / (1.+rw(k)) else u2f(:,:,k) = 1. + w2f(:,:,k) = 1. endif enddo call timing_on('COMM_TOTAL') call mpp_update_domains(u2f, domain) + call mpp_update_domains(w2f, domain) call timing_off('COMM_TOTAL') !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,w,rf,u,v, & #ifdef HIWPP !$OMP u00,v00, & #endif -!$OMP conserve,hydrostatic,pt,ua,va,u2f,cp,rg,ptop,rcv) +#ifdef MULTI_GASES +!$OMP q, & +#endif +!$OMP convert,molecular_diffusion, consv_te, & +!$OMP hydrostatic,pt,ua,va,u2f,w2f,cp,rg,ptop,rcp,rcv) do k=1,kmax if ( pm(k) < rf_cutoff ) then #ifdef HIWPP if (.not. hydrostatic) then do j=js,je do i=is,ie - w(i,j,k) = w(i,j,k)/(1.+rf(k)) + w(i,j,k) = w(i,j,k)/(1.+rw(k)) enddo enddo endif @@ -1235,7 +1316,8 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & enddo #else ! Add heat so as to conserve TE - if ( conserve ) then + + if ( convert ) then if ( hydrostatic ) then do j=js,je do i=is,ie @@ -1245,7 +1327,12 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & else do j=js,je do i=is,ie - pt(i,j,k) = pt(i,j,k) + 0.5*(ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2)*(1.-u2f(i,j,k)**2)*rcv +#ifdef MULTI_GASES + pt(i,j,k) = pt(i,j,k) + 0.5*((ua(i,j,k)**2+va(i,j,k)**2)*(1.-u2f(i,j,k)**2)+w(i,j,k)**2*(1.-w2f(i,j,k)**2)) & + /(cp*vicpq(q(i,j,k,:)) - rg*virq(q(i,j,k,:))) +#else + pt(i,j,k) = pt(i,j,k) + 0.5*((ua(i,j,k)**2+va(i,j,k)**2)*(1.-u2f(i,j,k)**2)+w(i,j,k)**2*(1.-w2f(i,j,k)**2))*rcv +#endif enddo enddo endif @@ -1264,7 +1351,7 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & if ( .not. hydrostatic ) then do j=js,je do i=is,ie - w(i,j,k) = u2f(i,j,k)*w(i,j,k) + w(i,j,k) = w2f(i,j,k)*w(i,j,k) enddo enddo endif @@ -1273,6 +1360,7 @@ subroutine Rayleigh_Super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, & enddo deallocate ( u2f ) + deallocate ( w2f ) end subroutine Rayleigh_Super diff --git a/model/fv_grid_utils.F90 b/model/fv_grid_utils.F90 index 3b89bc209..3ca58f88b 100644 --- a/model/fv_grid_utils.F90 +++ b/model/fv_grid_utils.F90 @@ -159,6 +159,12 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, grid real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c real(kind=R_GRID), pointer, dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya +! For MOLECULAR_DIFFUSION + real(kind=R_GRID), pointer, dimension(:,:) :: area_u, area_v + real(kind=R_GRID), pointer, dimension(:,:) :: dx6, dy6 + real, pointer, dimension(:,:) :: sina_6 + real, pointer, dimension(:,:) :: delu_6, delv_6, delu_5, delv_5 + real, pointer, dimension(:,:) :: del6_u, del6_v real, pointer, dimension(:,:) :: divg_u, divg_v real, pointer, dimension(:,:) :: cosa_u, cosa_v, cosa_s @@ -194,6 +200,18 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) dya => Atm%gridstruct%dya_64 sina => Atm%gridstruct%sina_64 cosa => Atm%gridstruct%cosa_64 +! MOLECULAR_DIFFUSION + if ( Atm%flagstruct%molecular_diffusion ) then + area_u => Atm%gridstruct%area_u_64 + area_v => Atm%gridstruct%area_v_64 + dx6 => Atm%gridstruct%dx6_64 + dy6 => Atm%gridstruct%dy6_64 + sina_6 => Atm%gridstruct%sina_6 + delu_6 => Atm%gridstruct%delu_6 + delv_6 => Atm%gridstruct%delv_6 + delu_5 => Atm%gridstruct%delu_5 + delv_5 => Atm%gridstruct%delv_5 + endif divg_u => Atm%gridstruct%divg_u divg_v => Atm%gridstruct%divg_v @@ -513,6 +531,10 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) rsin2 = big_number cosa = big_number sina = big_number +! MOLECULAR_DIFFUSION + if ( Atm%flagstruct%molecular_diffusion ) then + sina_6 = big_number + endif do j=js,je+1 do i=is,ie+1 @@ -571,6 +593,14 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) rsin_v(i,j) = 1. / max(tiny_number, sina_v(i,j)**2) enddo enddo + if ( Atm%flagstruct%molecular_diffusion ) then + do j=jsd+1,jed + do i=isd+1,ied + sina_6(i,j) = 0.25*(sin_sg(i-1,j ,7)+sin_sg(i,j ,6)+ & + sin_sg(i-1,j-1,8)+sin_sg(i,j-1,9)) + enddo + enddo + endif do j=jsd,jed do i=isd,ied @@ -700,6 +730,9 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) cosa_s = 0. rsin_u = 1. rsin_v = 1. + if ( Atm%flagstruct%molecular_diffusion ) then + sina_6 = 1. + endif endif if ( grid_type < 3 ) then @@ -786,6 +819,20 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) del6_v(ie+1,j) = 0.5*(sin_sg(npx,j,1)+sin_sg(npx-1,j,3))*dy(ie+1,j)/dxc(ie+1,j) endif enddo + if ( Atm%flagstruct%molecular_diffusion ) then + do j=jsd,jed + do i=isd,ied + delu_5(i,j) = sin_sg(i,j,5)*dxa(i,j)/dya(i,j) + delv_5(i,j) = sin_sg(i,j,5)*dya(i,j)/dxa(i,j) + enddo + enddo + do j=jsd+1,jed + do i=isd+1,ied + delu_6(i,j) = sina_6(i,j)*dy6(i,j)/dx6(i,j) + delv_6(i,j) = sina_6(i,j)*dx6(i,j)/dy6(i,j) + enddo + enddo + endif ! Initialize cubed_sphere to lat-lon transformation: call init_cubed_to_latlon( Atm%gridstruct, Atm%flagstruct%hydrostatic, agrid, grid_type, c2l_order, Atm%bd ) @@ -806,6 +853,12 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) gridtype=CGRID_NE_PARAM, complete=.true.) call mpp_update_domains(del6_v, del6_u, Atm%domain, flags=SCALAR_PAIR, & gridtype=CGRID_NE_PARAM, complete=.true.) + if ( Atm%flagstruct%molecular_diffusion ) then + call mpp_update_domains(delu_5, Atm%domain, complete=.true.) + call mpp_update_domains(delv_5, Atm%domain, complete=.true.) + call mpp_update_domains(delu_6, Atm%domain, complete=.true.) + call mpp_update_domains(delv_6, Atm%domain, complete=.true.) + endif call edge_factors (Atm%gridstruct%edge_s, Atm%gridstruct%edge_n, Atm%gridstruct%edge_w, & Atm%gridstruct%edge_e, non_ortho, grid, agrid, npx, npy, Atm%bd) call efactor_a2c_v(Atm%gridstruct%edge_vect_s, Atm%gridstruct%edge_vect_n, & @@ -840,6 +893,12 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) Atm%gridstruct%dyc = Atm%gridstruct%dyc_64 Atm%gridstruct%cosa = Atm%gridstruct%cosa_64 Atm%gridstruct%sina = Atm%gridstruct%sina_64 + if ( Atm%flagstruct%molecular_diffusion ) then + Atm%gridstruct%area_u = Atm%gridstruct%area_u_64 + Atm%gridstruct%area_v = Atm%gridstruct%area_v_64 + Atm%gridstruct%dx6 = Atm%gridstruct%dx6_64 + Atm%gridstruct%dy6 = Atm%gridstruct%dy6_64 + endif !--- deallocate the higher-order gridstruct arrays !rab deallocate ( Atm%gridstruct%grid_64 ) @@ -854,6 +913,21 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) deallocate ( Atm%gridstruct%dyc_64 ) deallocate ( Atm%gridstruct%cosa_64 ) deallocate ( Atm%gridstruct%sina_64 ) + if ( Atm%flagstruct%molecular_diffusion ) then + deallocate ( Atm%gridstruct%area_u_64 ) + deallocate ( Atm%gridstruct%area_v_64 ) + deallocate ( Atm%gridstruct%dx6_64 ) + deallocate ( Atm%gridstruct%dy6_64 ) + nullify(area_u) + nullify(area_v) + nullify(dx6) + nullify(dy6) + nullify(sina_6) + nullify(delu_6) + nullify(delv_6) + nullify(delu_5) + nullify(delv_5) + endif nullify(agrid) nullify(grid) @@ -865,6 +939,7 @@ subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order) nullify(dyc) nullify(dxa) nullify(dya) +! nullify(sina) nullify(cosa) nullify(divg_u) diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index 56617eedf..d3313a444 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -529,7 +529,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & !------------ ! Compute pkz -!hmhj pk is pe**kappa(=rgas/cp_air), but pkz=plyr**kappa(=r/cp) +!< pk is pe**kappa(=rgas/cp_air), but pkz=plyr**kappa(=r/cp) !------------ if ( hydrostatic ) then do k=1,km @@ -3522,7 +3522,11 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai integer, intent(in):: is, ie, isd,ied, jsd,jed, km, nwat, j, k integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel +#ifdef MULTI_GASES + real, intent(in), dimension(isd:ied,jsd:jed,km,num_gas):: q +#else real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q +#endif real, intent(out), dimension(is:ie):: cpm, qd real, intent(in), optional:: t1(is:ie) ! diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 15bb41f0b..f2e589e24 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -3645,7 +3645,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & enddo enddo -! map shpum, o3mr, liq_wat tracers +! map shpum, o3mr/spo3, liq_wat tracers do iq=1,ncnst ! if (iq == sphum .or. iq == liq_wat .or. iq == o3mr) then ! only remap if the data is already set if (iq /= cld_amt) then ! don't remap cld_amt diff --git a/model/molecular_diffusion.F90 b/model/molecular_diffusion.F90 new file mode 100644 index 000000000..c34432634 --- /dev/null +++ b/model/molecular_diffusion.F90 @@ -0,0 +1,333 @@ +!*********************************************************************** +!* 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 Molecular diffusion coefficients of +!> viscosity (mur), conductivity (lam), and diffusivity (d12) +!> and their efolding time effectiveness +!>@author H.-M. H. Juang, NOAA/NWS/NCEP/EMC + +module molecular_diffusion_mod + +! Modules Included: +! +! +! +! +! +! +! +! +! +!
Module NameFunctions Included
constants_modrdgas, cp_air
+ + use constants_mod, only: rdgas, cp_air + use fv_mp_mod, only: is_master + use mpp_mod, only: stdlog, input_nml_file + use fms_mod, only: check_nml_error, open_namelist_file, close_file + use fv_grid_utils_mod, only: g_sum + use mpp_domains_mod, only: domain2d + use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type + + +#ifdef MULTI_GASES + use multi_gases_mod, only: ind_gas, num_gas, virqd, vicpqd, vicvq, vicvqd, virq +#endif + + + implicit none + + integer :: ind_gas_str, ind_gas_end + integer :: md_layers, md_tadj_layers + real tau_visc, tau_cond, tau_diff + real md_consv_te + real md_wait_hr + real md_wait_sec + logical md_time + real, parameter:: amo=15.9994, amo2=31.9999, amo3=47.9982 !g/mol + real, parameter:: amn2=28.013, amh2o=18.0154 !g/mol +!< muo3 and muh2o are not precise, correct later + real, parameter:: muo=3.9e-7, muo2=4.03e-7, muo3=4.03e-7 !kg/m/s + real, parameter:: mun2=3.43e-7, muh2o=3.43e-7 !kg/m/s +!< lao3 is not precise values, but o3_n is very small + real, parameter:: lao=75.9e-5, lao2=56.e-5, lao3=36.e-5 !kg/m/s + real, parameter:: lan2=56.e-5, lah2o=55.e-5 !kg/m/s + real, parameter:: cpo=1299.185, cpo2=918.0969, cpo3=820.2391 + real, parameter:: cpn2=1031.108, cph2o=1846.00 + real, parameter:: avgd=6.0221415e23 ! Avogadro constant + real, parameter:: bz=1.3806505e-23 ! Boltzmann constant J/K + real, parameter:: a12=9.69e18 ! O-O2 diffusion params + real, parameter:: s12=0.774 +! + public molecular_diffusion_init, read_namelist_molecular_diffusion_nml + public molecular_diffusion_coefs, thermosphere_adjustment + + CONTAINS +! -------------------------------------------------------- + subroutine molecular_diffusion_init(ncnst,nwat) +!-------------------------------------------- +! molecular diffusion control for each effect from namelist +! Input: +! ncnst : number of all prognostic tracers +! nwat : number of water and the end location of water +!-------------------------------------------- + integer, intent(in):: ncnst, nwat +! + md_wait_sec = md_wait_hr * 3600.0 + md_time = .false. + + if( is_master() ) then + write(*,*) ' molecular_diffusion is on' + write(*,*) ' molecular_diffusion initial wait seconds ',md_wait_sec + write(*,*) ' molecular_diffusion number of layers ',md_layers + write(*,*) ' thermosphere adjustment number of layers ',md_tadj_layers + write(*,*) ' energy conservation of MD (0:off, 1:on) ',md_consv_te + write(*,*) ' viscosity day ',tau_visc,' with effect ',tau_visc + write(*,*) ' conductivity day ',tau_cond,' with effect ',tau_cond + write(*,*) ' diffusivity day ',tau_diff,' with effect ',tau_diff + endif + +#ifdef MULTI_GASES + ind_gas_str = ind_gas + ind_gas_end = num_gas +#else + ind_gas_str = nwat + 1 + ind_gas_end = ncnst +#endif + + return + end subroutine molecular_diffusion_init + +! -------------------------------------------------------- + subroutine read_namelist_molecular_diffusion_nml(nml_filename,ncnst,nwat) + + character(*), intent(IN) :: nml_filename + integer, intent(IN) :: ncnst, nwat + integer :: ierr, f_unit, unit, ios + + namelist /molecular_diffusion_nml/ tau_visc, tau_cond, tau_diff, & + md_layers, md_tadj_layers, & + md_wait_hr, md_consv_te + + unit = stdlog() + +#ifdef INTERNAL_FILE_NML + + ! Read molecular_diffusion namelist + read (input_nml_file,molecular_diffusion_nml,iostat=ios) + ierr = check_nml_error(ios,'molecular_diffusion_nml') + +#else + ! Read molecular_diffusion namelist + f_unit = open_namelist_file(nml_filename) + rewind (f_unit) + read (f_unit,molecular_diffusion_nml,iostat=ios) + ierr = check_nml_error(ios,'molecular_diffusion_nml') + call close_file(f_unit) +#endif + write(unit, nml=molecular_diffusion_nml) + call molecular_diffusion_init(ncnst,nwat) + + return + end subroutine read_namelist_molecular_diffusion_nml + +! -------------------------------------------------------- + subroutine molecular_diffusion_coefs(dim,plyr,temp,q,mur,lam,d12) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Input: dim : length of field +! plyr: pressure at given layer (pascal) +! temp: temperature (K) +! q : tracers, needs only gas (kg/kg) +! Ouput: mur : viscosity for momemtum +! lam : conductivity for thermodynamics +! d12 : diffusivity for tracers +!-------------------------------------------- + integer, intent(in):: dim + real, intent(in):: plyr(dim), temp(dim), q(dim,*) + real, intent(out):: mur(dim), lam(dim), d12(dim) +! Local: + integer i, n, spo3, spo, spo2 + real am, fgas, a12bz, avgdbz + real qo, qo2, qo3, qn2, qh2o + real o_n, o2_n, o3_n, n2_n, h2o_n + real mu, la, t69, rho, cpx, cvx + +!constants + a12bz = a12 * bz + avgdbz= avgd * bz + spo3 = ind_gas_str + spo = spo3 + 1 + spo2 = spo + 1 + if( spo.gt.ind_gas_end ) spo=0 + if( spo2.gt.ind_gas_end ) spo2=0 + + do n=1,dim +!check + if(plyr(n).le.0.0) then + write(*,*) 'ERROR non positive value of plyr ',plyr(n) + stop + endif + if(temp(n).le.0.0) then + write(*,*) 'ERROR non positive value of temp ',temp(n) + stop + endif + + d12(n) = a12bz*temp(n)**s12 * temp(n)/plyr(n) + if( d12(n).lt.0.0 .and. is_master() ) then + write(*,*) 'ERROR negative d12 a12bz temp plyr s12 ', & + d12(n),a12bz,temp(n),plyr(n),s12 + endif + + fgas = 1.0 + do i=2,ind_gas_str-1 + fgas = fgas - max(0.0,q(n,i)) + enddo + fgas = max(min(fgas,1.0),1.0E-20) ! reasonable assured + fgas = 1./fgas + + qh2o = max(0.0,q(n,1))*fgas + qo = 0.0 + qo2 = 0.0 + qo3 = 0.0 + if(spo .ne.0) qo = max(0.0,q(n,spo ))*fgas + if(spo2.ne.0) qo2 = max(0.0,q(n,spo2))*fgas + if(spo3.ne.0) qo3 = max(0.0,q(n,spo3))*fgas + qn2 = 1.0 - qo - qo2 - qo3 - qh2o + +! reasonable values assure + qo = max(min(qo ,1.0),0.0) + qo2 = max(min(qo2,1.0),0.0) + qo3 = max(min(qo3,1.0),0.0) + qn2 = max(min(qn2,1.0),0.0) + qh2o = max(min(qh2o,1.0),0.0) + cpx = ( cpo*qo + cpo2*qo2 + cpo3*qo3 + cpn2*qn2 + cph2o*qh2o ) / & + ( qo + qo2 + qo3 + qn2 + qh2o ) + cvx = cpx - rdgas + + am = qo/amo + qo2/amo2 + qo3/amo3 + qn2/amn2 +qh2o/amh2o + am = 1.0 / am ! g/mol + o_n = qo * am / amo + o2_n = qo2 * am / amo2 + o3_n = qo3 * am / amo3 + n2_n = qn2 * am / amn2 + h2o_n = qh2o * am / amh2o + mu = o_n*muo + o2_n*muo2 + o3_n*muo3 + n2_n*mun2 + h2o_n*muh2o + la = o_n*lao + o2_n*lao2 + o3_n*lao3 + n2_n*lan2 + h2o_n*lah2o + + t69 = temp(n) ** 0.69 + rho = 1e-3 * am * plyr(n)/temp(n) / avgdbz ! km/m^3 + + mur(n) = mu * t69 / rho + lam(n) = la * t69 / rho / cvx +! lam(n) = la * t69 / rho / cpx + +!reasonable assured + mur(n) = min(mur(n),1.0e15) ! viscosity + lam(n) = min(lam(n),1.0e15) ! conductivity + d12(n) = min(d12(n),1.0e15) ! diffusivity + + enddo + + return + end subroutine molecular_diffusion_coefs + subroutine thermosphere_adjustment(domain,gridstruct,npz,bd,ng,pt) + type(fv_grid_bounds_type), intent(IN) :: bd + real , intent(INOUT) :: pt(bd%is:bd%ie ,bd%js:bd%je, npz) + type(fv_grid_type), intent(IN), target :: gridstruct + type(domain2d), intent(INOUT) :: domain + integer, intent(IN) :: ng, npz + real :: correct, aave_t1, ttsave, tdsave + integer :: k, j, i + real :: t(bd%is:bd%ie,bd%js:bd%je) + 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 + + + + + k=md_tadj_layers + do j=js,je + do i=is,ie + t(i,j) = pt(i,j,k) + enddo + enddo + aave_t1 = g_sum(domain,t,is,ie,js,je,ng,gridstruct%area_64, 1) + k=md_tadj_layers-1 + do j=js,je + do i=is,ie + t(i,j) = pt(i,j,k) + enddo + enddo + ttsave = g_sum(domain,t,is,ie,js,je,ng,gridstruct%area_64, 1) + tdsave = ttsave - aave_t1 +! if( is_master() ) write(*,*) ' k t1 ts ',k,aave_t1,ttsave +! correct_sum = 0.0 + do k=md_tadj_layers-2,1,-1 + do j=js,je + do i=is,ie + t(i,j) = pt(i,j,k) + enddo + enddo + aave_t1 = g_sum(domain,t,is,ie,js,je,ng,gridstruct%area_64, 1) +! if( is_master() ) write(*,*) ' k ts t1 ',k,ttsave,aave_t1 + if( aave_t1 .gt. ttsave ) then + tdsave = min(tdsave,aave_t1-ttsave) + ttsave = ttsave + tdsave + else + tdsave = 0.0 + endif + correct = ttsave - aave_t1 +! if( is_master() ) write(*,*) ' k t1 ts correct ',k,aave_t1,ttsave,correct + if( correct .ne. 0.0 ) then + do j=js,je + do i=is,ie + pt(i,j,k) = t(i,j) + correct + enddo + enddo + endif +! correct_sum = correct_sum + correct + enddo +! adjust for energy conserving by evenly distribute total correction +! if( correct_sum .ne. 0.0 ) then +! correct = - correct_sum / (md_tadj_layers + 10) +! if( is_master() ) write(*,*) ' sum=',correct_sum,' correct=',correct +! do k=1,md_tadj_layers+10 +! do j=js,je +! do i=is,ie +! pt(i,j,k) = t(i,j) + correct +! enddo +! enddo +! enddo +! endif + return + end subroutine thermosphere_adjustment + +end module molecular_diffusion_mod diff --git a/model/multi_gases.F90 b/model/multi_gases.F90 index b8ed58f4d..082a9c0f0 100644 --- a/model/multi_gases.F90 +++ b/model/multi_gases.F90 @@ -36,6 +36,9 @@ module multi_gases_mod ! ! +!* --------- Default gas values for wam ---------------------------------------- +!* R_n2+= 295.3892, R_h2o= 461.50, R_o3=173.2247, R_o= 519.674, R_o2=259.8370 +!* Cpn2+=1031.1083, Cph2o=1846.00, Cpo3=820.2391, Cpo=1299.185, Cpo2=918.0969 use constants_mod, only: rdgas, rvgas, cp_air use fv_mp_mod, only: is_master use mpp_mod, only: stdlog, input_nml_file @@ -59,10 +62,12 @@ module multi_gases_mod public virq public virq_max public virqd + public vicpq public vicpqd public virq_nodq public virq_qpz public vicpqd_qpz + public vicvq public vicvqd public vicvqd_qpz @@ -132,11 +137,11 @@ subroutine multi_gases_init(ngas, nwat) end subroutine multi_gases_init subroutine read_namelist_multi_gases_nml(nml_filename,ncnst,nwat) - character(*), intent(IN) :: nml_filename - integer, intent(IN) :: ncnst, nwat - integer :: ierr, f_unit, unit, ios + character(*), intent(IN) :: nml_filename + integer, intent(IN) :: ncnst, nwat + integer :: ierr, f_unit, unit, ios - namelist /multi_gases_nml/ ri,cpi + namelist /multi_gases_nml/ ri,cpi unit = stdlog() @@ -190,7 +195,7 @@ pure real function virq(q) virq = vir(0)+virq/(1.0-sum(q(sphump1:sphum+num_wat-1))) return - end function + end function virq !-------------------------------------------- ! -------------------------------------------------------- @@ -209,7 +214,7 @@ pure real function virq_nodq(q) end do return - end function + end function virq_nodq !-------------------------------------------- ! -------------------------------------------------------- @@ -231,7 +236,7 @@ pure real function virq_max(q, qmin) return - end function + end function virq_max !-------------------------------------------- ! -------------------------------------------------------- @@ -253,7 +258,7 @@ pure real function virq_qpz(q, qpz) return - end function + end function virq_qpz !-------------------------------------------- ! -------------------------------------------------------- @@ -273,7 +278,27 @@ pure real function virqd(q) virqd = vir(0)+virqd/(1.0-sum(q(sphum:sphum+num_wat-1))) return - end function + end function virqd +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicpq(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir/(1-qc) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + vicpq = vicp(sphum)*q(sphum) + do n=ind_gas,num_gas + vicpq = vicpq+vicp(n)*q(sphum+n-1) + end do + vicpq = vicp(0)+vicpq/(1.0-sum(q(sphump1:sphum+num_wat-1))) + + return + end function vicpq !-------------------------------------------- ! -------------------------------------------------------- @@ -293,7 +318,7 @@ pure real function vicpqd(q) vicpqd = vicp(0)+vicpqd/(1.0-sum(q(sphum:sphum+num_wat-1))) return - end function + end function vicpqd !-------------------------------------------- ! -------------------------------------------------------- @@ -314,7 +339,7 @@ pure real function vicpqd_qpz(q, qpz) vicpqd_qpz = vicp(0)+vicpqd_qpz/(1.0-qpz) return - end function + end function vicpqd_qpz !-------------------------------------------- ! -------------------------------------------------------- @@ -334,7 +359,27 @@ pure real function vicvqd(q) vicvqd = vicv(0)+vicvqd/(1.0-sum(q(sphum:sphum+num_wat-1))) return - end function + end function vicvqd +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicvq(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir/(1-qc) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + vicvq = vicv(sphum)*q(sphum) + do n=ind_gas,num_gas + vicvq = vicvq+vicv(n)*q(sphum+n-1) + end do + vicvq = vicv(0)+vicvq/(1.0-sum(q(sphump1:sphum+num_wat-1))) + + return + end function vicvq !-------------------------------------------- ! -------------------------------------------------------- @@ -355,7 +400,7 @@ pure real function vicvqd_qpz(q,qpz) vicvqd_qpz = vicv(0)+vicvqd_qpz/(1.0-qpz) return - end function + end function vicvqd_qpz !-------------------------------------------- end module multi_gases_mod diff --git a/model/nh_core.F90 b/model/nh_core.F90 index a197f3bf7..21c6790fd 100644 --- a/model/nh_core.F90 +++ b/model/nh_core.F90 @@ -155,7 +155,7 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, & 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 -!hmhj pk3 at interface , interface pk is using constant akap +!< pk3 at interface , interface pk is using constant akap pk3(i,j,k) = exp(akap*peln2(i,k)) enddo enddo diff --git a/model/sw_core.F90 b/model/sw_core.F90 index fc97b060d..6ed4a2e50 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -44,8 +44,11 @@ module sw_core_mod ! ! - use tp_core_mod, only: fv_tp_2d, pert_ppm, copy_corners - use fv_mp_mod, only: fill_corners, XDir, YDir + use molecular_diffusion_mod, only: molecular_diffusion_coefs, & + tau_visc, tau_cond, tau_diff + use tp_core_mod, only: fv_tp_2d, pert_ppm, copy_corners, & + deln_flux_explm, deln_flux_explm_udvd + use fv_mp_mod, only: is_master, fill_corners, XDir, YDir use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type, fv_flags_type use a2b_edge_mod, only: a2b_ord4 @@ -93,7 +96,8 @@ module sw_core_mod private - public :: c_sw, d_sw, fill_4corners, del6_vt_flux, divergence_corner, divergence_corner_nest + public :: c_sw, d_sw, d_md, fill_4corners, & + del6_vt_flux, divergence_corner, divergence_corner_nest contains @@ -1574,6 +1578,124 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & end subroutine d_sw + +! d_md :: D-Grid 2D molecular diffusion + +!>@brief The subroutine 'd_md' peforms D-grid molecular diffusion + subroutine d_md(p, t, u, v, w, q, & + it, nq, k, km, dt, & + gridstruct, flagstruct, bd) + + integer, intent(IN):: nq, k, km, it + real , intent(IN):: dt + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed ):: p, t + real, intent(INOUT), dimension(bd%isd: , bd%jsd: ):: w + real, intent(INOUT), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed+1):: u + real, intent(INOUT), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ):: v + real, intent(INOUT):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km,nq) + type(fv_grid_type), intent(IN), target :: gridstruct + type(fv_flags_type), intent(IN), target :: flagstruct +!--- + real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed,1:nq):: qtra + real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: temp, plyr + real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: visc, cond, diff + + real :: coefmax + + integer :: i, j, iq, n + integer :: ijm,idir + + integer :: is, ie, js, je + integer :: isd, ied, jsd, jed + integer :: npx, npy, nord + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + ijm = (ied-isd+1)*(jed-jsd+1) + + npx = flagstruct%npx + npy = flagstruct%npy + nord = flagstruct%nord + +! +! prepare molecular diffusion coefficients +! + do j=jsd,jed + do i=isd,ied + plyr(i,j) = p(i,j) + temp(i,j) = t(i,j) + qtra(i,j,1:nq) = q(i,j,k,1:nq) + enddo + enddo + +! fill up corner for coefficient computation with values + if( flagstruct%nord>0 .and. (.not. (flagstruct%regional))) then + idir = mod(it-1,2)+1 ! alternated by 1 and 2 to avoid bias + call copy_corners(plyr, npx, npy, idir, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(temp, npx, npy, idir, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + do iq=1,nq + call copy_corners(qtra(isd,jsd,iq), npx, npy, idir, & + gridstruct%nested, bd , & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + enddo + endif + + call molecular_diffusion_coefs(ijm,plyr(isd,jsd ),temp(isd,jsd),& + qtra(isd,jsd,1),visc(isd,jsd),& + cond(isd,jsd ),diff(isd,jsd)) + +! time scale and options + + coefmax=0.125*min(gridstruct%da_min,gridstruct%da_min_c)/abs(dt) +! if( is_master() ) & +! write(*,'(a,i5,5g13.6)') & +! ' k coefmax da_c da dx dy',k,coefmax,gridstruct%da_min_c,& +! gridstruct%da_min,gridstruct%dx(is,js),gridstruct%dy(is,js) + + do j=jsd,jed + do i=isd,ied + visc(i,j) = min(coefmax,visc(i,j))*abs(dt)*tau_visc + cond(i,j) = min(coefmax,cond(i,j))*abs(dt)*tau_cond + diff(i,j) = min(coefmax,diff(i,j))*abs(dt)*tau_diff + enddo + enddo + +! +! compute diffusion with dimensional split alternatively for implicit +! explicit diffusion has direct computation not dimensional splitting. + idir = mod(it-1,2)+1 + +! t + call deln_flux_explm(nord,is,ie,js,je,npx,npy,cond,t,gridstruct,bd) + +! q + do iq=1,nq + call deln_flux_explm(nord,is,ie,js,je,npx,npy,diff,q(isd,jsd,k,iq),gridstruct,bd) + enddo + +! u v + call deln_flux_explm_udvd(nord,is,ie,js,je,npx,npy,visc, & + u,v,gridstruct,bd) +! w + call deln_flux_explm(nord,is,ie,js,je,npx,npy,visc,w,gridstruct,bd) + + return + + end subroutine d_md + + !>@brief The subroutine 'del6_vt_flux' applies 2nd, 4th, or 6th-order damping !! to fluxes ("vorticity damping") subroutine del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd) diff --git a/model/tp_core.F90 b/model/tp_core.F90 index 7ea6cc4cd..f095d79e9 100644 --- a/model/tp_core.F90 +++ b/model/tp_core.F90 @@ -56,6 +56,8 @@ module tp_core_mod private public fv_tp_2d, pert_ppm, copy_corners + public deln_flux_explm, deln_flux_explm_udvd + public copy_corners_for_udvd real, parameter:: ppm_fac = 1.5 !< nonlinear scheme limiter: between 1 and 2 real, parameter:: r3 = 1./3. @@ -1385,4 +1387,331 @@ subroutine deln_flux(nord,is,ie,js,je, npx, npy, damp, q, fx, fy, gridstruct, bd end subroutine deln_flux +!=========================================== + subroutine deln_flux_explm(nord,is,ie,js,je,npx,npy,damp,q,gridstruct,bd) +!> Del-2 time split explicit damping for the cell-mean values (A grid) +!------------------ +!> nord = used for extra grid as side +!------------------ + type(fv_grid_bounds_type), intent(IN) :: bd + integer, intent(in):: nord !< del-n + integer, intent(in):: is,ie,js,je, npx, npy + real, intent(in):: damp(bd%isd:bd%ied, bd%jsd:bd%jed) + real, intent(inout):: q(bd%isd:bd%ied, bd%jsd:bd%jed) ! q ghosted on input + type(fv_grid_type), intent(IN), target :: gridstruct +! local: + real fx2(bd%isd:bd%ied,bd%jsd:bd%jed), fy2(bd%isd:bd%ied,bd%jsd:bd%jed) + real xv(bd%isd:bd%ied,bd%jsd:bd%jed) + real d2(bd%isd:bd%ied,bd%jsd:bd%jed) + integer i, j, m, n, nt, i1, i2, j1, j2 + +#ifdef USE_SG + real, pointer, dimension(:,:) :: dx, dy, rdxc, rdyc + real, pointer, dimension(:,:,:) :: sin_sg + + dx => gridstruct%dx + dy => gridstruct%dy + rdxc => gridstruct%rdxc + rdyc => gridstruct%rdyc + sin_sg => gridstruct%sin_sg +#endif + + i1 = is-1-nord; i2 = ie+1+nord + j1 = js-1-nord; j2 = je+1+nord + + do j=bd%jsd,bd%jed + do i=bd%isd,bd%ied + xv(i,j) = damp(i,j) + d2(i,j) = q(i,j) + enddo + enddo + +! explicit scheme +! for dadx * dy + call copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(xv, npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) +! + do j=j1,j2 + do i=i1+1,i2 +#ifdef USE_SG + fx2(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*(d2(i-1,j)-d2(i,j))*rdxc(i,j) +#else + fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i-1,j)-d2(i,j)) +#endif + fx2(i,j) = 0.5*(xv(i-1,j)+xv(i,j))*fx2(i,j) + enddo + enddo +! for dady * dx + call copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(xv, npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) +! + do i=i1,i2 + do j=j1+1,j2 +#ifdef USE_SG + fy2(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*(d2(i,j-1)-d2(i,j))*rdyc(i,j) +#else + fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j-1)-d2(i,j)) +#endif + fy2(i,j) = 0.5*(xv(i,j-1)+xv(i,j))*fy2(i,j) + enddo + enddo +! + do j=j1+1,j2-1 + do i=i1+1,i2-1 + q(i,j) = q(i,j) + & + (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j) + enddo + enddo + + return + end subroutine deln_flux_explm + + subroutine deln_flux_explm_udvd(nord,is,ie,js,je,npx,npy,damp,u,v,gridstruct,bd) +!> Del-2 time split explicit damping for the cell-edge values (D grid) +!------------------ +!> nord = used for extra grid as side +!------------------ + type(fv_grid_bounds_type), intent(IN) :: bd + integer, intent(in):: nord !< del-n + integer, intent(in):: is,ie,js,je, npx, npy + real, intent(in):: damp(bd%isd:bd%ied, bd%jsd:bd%jed) + real, intent(inout):: u(bd%isd:bd%ied , bd%jsd:bd%jed+1) + real, intent(inout):: v(bd%isd:bd%ied+1, bd%jsd:bd%jed ) + type(fv_grid_type), intent(IN), target :: gridstruct +! local: + real fx1(bd%isd:bd%ied,bd%jsd:bd%jed), fy1(bd%isd:bd%ied,bd%jsd:bd%jed) + real fx2(bd%isd:bd%ied,bd%jsd:bd%jed), fy2(bd%isd:bd%ied,bd%jsd:bd%jed) + real x5(bd%isd:bd%ied,bd%jsd:bd%jed) + real x6(bd%isd:bd%ied,bd%jsd:bd%jed) + real d1(bd%isd:bd%ied ,bd%jsd:bd%jed+1) + real d2(bd%isd:bd%ied+1,bd%jsd:bd%jed ) + integer i, j, m, n, nt, i1, i2, j1, j2 + +#ifdef USE_SG + real, pointer, dimension(:,:) :: dxa, dya, dx6, dy6 + real, pointer, dimension(:,:) :: rdxa, rdya, rdx6, rdy6 + real, pointer, dimension(:,:) :: sina_6 + real, pointer, dimension(:,:,:) :: sin_sg + + sin_sg => gridstruct%sin_sg + sina_6 => gridstruct%sina_6 + dxa => gridstruct%dxa + dya => gridstruct%dya + rdxa => gridstruct%rdxa + rdya => gridstruct%rdya + dx6 => gridstruct%dx6 + dy6 => gridstruct%dy6 + rdx6 => gridstruct%rdx6 + rdy6 => gridstruct%rdy6 +#endif + + i1 = is-1-nord; i2 = ie+1+nord + j1 = js-1-nord; j2 = je+1+nord + + do j=bd%jsd,bd%jed + do i=bd%isd,bd%ied + x5(i,j) = damp(i,j) + enddo + enddo + do j=bd%jsd,bd%jed+1 + do i=bd%isd,bd%ied + d1(i,j) = u(i,j) + enddo + enddo + do j=bd%jsd,bd%jed + do i=bd%isd,bd%ied+1 + d2(i,j) = v(i,j) + enddo + enddo + +! explicit scheme + +! + call copy_corners_for_udvd(d1,d2, npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(x5, npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + do j=bd%jsd+1,bd%jed + do i=bd%isd+1,bd%ied + x6(i,j) = 0.25*(x5(i,j)+x5(i,j-1)+x5(i-1,j)+x5(i-1,j-1)) + enddo + enddo +! +! dux + do j=j1+1,j2 + do i=i1+1,i2 +#ifdef USE_SG + fx1(i,j) = sina_6(i,j)*dy6(i,j)*rdx6(i,j)*(d1(i-1,j)-d1(i,j)) +#else + fx1(i,j) = gridstruct%delu_6(i,j)*(d1(i-1,j)-d1(i,j)) +#endif + fx1(i,j) = x6(i,j)*fx1(i,j) + enddo + enddo +! dvx + do j=j1,j2 + do i=i1,i2 +#ifdef USE_SG + fx2(i,j) = sin_sg(i,j,5)*dya(i,j)*rdxa(i,j)*(d2(i,j)-d2(i+1,j)) +#else + fx2(i,j) = gridstruct%delv_5(i,j)*(d2(i,j)-d2(i+1,j)) +#endif + fx2(i,j) = x5(i,j)*fx2(i,j) + enddo + enddo +! + call copy_corners_for_udvd(d1,d2, npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + call copy_corners(x5, npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, & + gridstruct%nw_corner, gridstruct%ne_corner) + do j=bd%jsd+1,bd%jed + do i=bd%isd+1,bd%ied + x6(i,j) = 0.25*(x5(i,j)+x5(i,j-1)+x5(i-1,j)+x5(i-1,j-1)) + enddo + enddo +! +! duy + do j=j1,j2 + do i=i1,i2 +#ifdef USE_SG + fy1(i,j) = sin_sg(i,j,5)*dxa(i,j)*rdya(i,j)*(d1(i,j)-d1(i,j+1)) +#else + fy1(i,j) = gridstruct%delu_5(i,j)*(d1(i,j)-d1(i,j+1)) +#endif + fy1(i,j) = x5(i,j)*fy1(i,j) + enddo + enddo +! dvy + do j=j1+1,j2 + do i=i1+1,i2 +#ifdef USE_SG + fy2(i,j) = sina_6(i,j)*dx6(i,j)*rdy6(i,j)*(d2(i,j-1)-d2(i,j)) +#else + fy2(i,j) = gridstruct%delv_6(i,j)*(d2(i,j-1)-d2(i,j)) +#endif + fy2(i,j) = x6(i,j)*fy2(i,j) + enddo + enddo +! ------------------ + do j=j1+1,j2 + do i=i1+1,i2-1 + u(i,j) = u(i,j) + & + (fx1(i,j)-fx1(i+1,j)+fy1(i,j-1)-fy1(i,j)) * & + gridstruct%rarea_u(i,j) + enddo + enddo + do j=j1+1,j2-1 + do i=i1+1,i2 + v(i,j) = v(i,j) + & + (fx2(i-1,j)-fx2(i,j)+fy2(i,j)-fy2(i,j+1)) * & + gridstruct%rarea_v(i,j) + enddo + enddo + + return + end subroutine deln_flux_explm_udvd + + + subroutine copy_corners_for_udvd(ud, vd, npx, npy, dir, nested, bd, & + sw_corner, se_corner, nw_corner, ne_corner) +type(fv_grid_bounds_type), intent(IN) :: bd + integer, intent(in):: npx, npy, dir + real, intent(inout):: ud(bd%isd:bd%ied ,bd%jsd:bd%jed+1) + real, intent(inout):: vd(bd%isd:bd%ied+1,bd%jsd:bd%jed ) + logical, intent(IN) :: nested, sw_corner, se_corner, nw_corner, ne_corner + integer i,j,ng + + ng = bd%ng + + if (nested) return + + if ( dir == 1 ) then +! XDir: D grid (note that npx = npy at grid corners + if ( sw_corner ) then + do j=1-ng,0 + do i=1-ng,0 + ud(i,j) = -vd(j,1-i) + vd(i,j) = ud(j,2-i) + enddo + enddo + endif + if ( se_corner ) then + do j=1-ng,0 + do i=npx,npx+ng-1 + ud(i ,j) = vd(npy-j+1,i-npx+1) + vd(i+1,j) = -ud(npy-j ,i-npx+2) + enddo + enddo + endif + if ( ne_corner ) then + do j=npy,npy+ng-1 + do i=npx,npx+ng-1 + ud(i ,j+1) = -vd(j+1,2*npx-1-i) + vd(i+1,j ) = ud(j ,2*npx-1-i) + enddo + enddo + endif + if ( nw_corner ) then + do j=npy,npy+ng-1 + do i=1-ng,0 + ud(i,j+1) = vd(npy-j,i-1+npx) + vd(i,j ) = -ud(npy-j,i-1+npx) + enddo + enddo + endif + + elseif ( dir == 2 ) then +! YDir: D grid + + if ( sw_corner ) then + do j=1-ng,0 + do i=1-ng,0 + ud(i,j) = vd(2-j,i) + vd(i,j) = -ud(1-j,i) + enddo + enddo + endif + if ( se_corner ) then + do j=1-ng,0 + do i=npx,npx+ng-1 + ud(i ,j) = -vd(npy+j-1,npx-i) + vd(i+1,j) = ud(npy+j-1,npx-i) + enddo + enddo + endif + if ( ne_corner ) then + do j=npy,npy+ng-1 + do i=npx,npx+ng-1 + ud(i ,j+1) = vd(2*npy-1-j,i ) + vd(i+1,j ) = -ud(2*npy-1-j,i+1) + enddo + enddo + endif + if ( nw_corner ) then + do j=npy,npy+ng-1 + do i=1-ng,0 + ud(i,j+1) = -vd(j+2-npy,npx-i ) + vd(i,j ) = ud(j+1-npy,npx-i+1) + enddo + enddo + endif + + endif + + end subroutine copy_corners_for_udvd + + + end module tp_core_mod diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index a4ad6f2c5..836e75ca8 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -227,7 +227,7 @@ subroutine get_external_ic( Atm, fv_domain, cold_start ) integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, sgs_tke, cld_amt integer :: liq_aero, ice_aero #ifdef MULTI_GASES - integer :: spfo, spfo2, spfo3 + integer :: spo, spo2, spo3 #else integer :: o3mr #endif @@ -323,9 +323,9 @@ subroutine get_external_ic( Atm, fv_domain, cold_start ) snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat') graupel = get_tracer_index(MODEL_ATMOS, 'graupel') #ifdef MULTI_GASES - spfo = get_tracer_index(MODEL_ATMOS, 'spfo') - spfo2 = get_tracer_index(MODEL_ATMOS, 'spfo2') - spfo3 = get_tracer_index(MODEL_ATMOS, 'spfo3') + spo = get_tracer_index(MODEL_ATMOS, 'spo') + spo2 = get_tracer_index(MODEL_ATMOS, 'spo2') + spo3 = get_tracer_index(MODEL_ATMOS, 'spo3') #else o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') #endif @@ -345,12 +345,12 @@ subroutine get_external_ic( Atm, fv_domain, cold_start ) if ( graupel > 0 ) & call prt_maxmin('graupel', Atm%q(:,:,:,graupel), is, ie, js, je, ng, Atm%npz, 1.) #ifdef MULTI_GASES - if ( spfo > 0 ) & - call prt_maxmin('SPFO', Atm%q(:,:,:,spfo), is, ie, js, je, ng, Atm%npz, 1.) - if ( spfo2 > 0 ) & - call prt_maxmin('SPFO2', Atm%q(:,:,:,spfo2), is, ie, js, je, ng, Atm%npz, 1.) - if ( spfo3 > 0 ) & - call prt_maxmin('SPFO3', Atm%q(:,:,:,spfo3), is, ie, js, je, ng, Atm%npz, 1.) + if ( spo > 0 ) & + call prt_maxmin('SPO', Atm%q(:,:,:,spo), is, ie, js, je, ng, Atm%npz, 1.) + if ( spo2 > 0 ) & + call prt_maxmin('SPO2', Atm%q(:,:,:,spo2), is, ie, js, je, ng, Atm%npz, 1.) + if ( spo3 > 0 ) & + call prt_maxmin('SPO3', Atm%q(:,:,:,spo3), is, ie, js, je, ng, Atm%npz, 1.) #else if ( o3mr > 0 ) & call prt_maxmin('O3MR', Atm%q(:,:,:,o3mr), is, ie, js, je, ng, Atm%npz, 1.) @@ -1955,7 +1955,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) integer :: isd, ied, jsd, jed integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, sgs_tke, cld_amt #ifdef MULTI_GASES - integer :: spfo, spfo2, spfo3 + integer :: spo, spo2, spo3 #else integer :: o3mr #endif @@ -1965,7 +1965,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) real(kind=R_GRID), dimension(3):: e1, e2, ex, ey real, allocatable:: ps_gfs(:,:), zh_gfs(:,:,:) #ifdef MULTI_GASES - real, allocatable:: spfo_gfs(:,:,:), spfo2_gfs(:,:,:), spfo3_gfs(:,:,:) + real, allocatable:: spo_gfs(:,:,:), spo2_gfs(:,:,:), spo3_gfs(:,:,:) #else real, allocatable:: o3mr_gfs(:,:,:) #endif @@ -2004,9 +2004,9 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat') graupel = get_tracer_index(MODEL_ATMOS, 'graupel') #ifdef MULTI_GASES - spfo = get_tracer_index(MODEL_ATMOS, 'spfo') - spfo2 = get_tracer_index(MODEL_ATMOS, 'spfo2') - spfo3 = get_tracer_index(MODEL_ATMOS, 'spfo3') + spo = get_tracer_index(MODEL_ATMOS, 'spo') + spo2 = get_tracer_index(MODEL_ATMOS, 'spo2') + spo3 = get_tracer_index(MODEL_ATMOS, 'spo3') #else o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') #endif @@ -2023,9 +2023,9 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) print *, 'graupel = ', graupel endif #ifdef MULTI_GASES - print *, ' spfo3 = ', spfo3 - print *, ' spfo = ', spfo - print *, ' spfo2 = ', spfo2 + print *, ' spo3 = ', spo3 + print *, ' spo = ', spo + print *, ' spo2 = ', spo2 #else print *, ' o3mr = ', o3mr #endif @@ -2072,9 +2072,9 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) !! Read in o3mr, ps and zh from GFS_data.tile?.nc #ifdef MULTI_GASES - allocate (spfo3_gfs(is:ie,js:je,levp_gfs)) - allocate ( spfo_gfs(is:ie,js:je,levp_gfs)) - allocate (spfo2_gfs(is:ie,js:je,levp_gfs)) + allocate ( spo_gfs(is:ie,js:je,levp_gfs)) + allocate (spo2_gfs(is:ie,js:je,levp_gfs)) + allocate (spo3_gfs(is:ie,js:je,levp_gfs)) #else allocate (o3mr_gfs(is:ie,js:je,levp_gfs)) #endif @@ -2084,13 +2084,13 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(GFS_restart, "lat", "y") call register_axis(GFS_restart, "lon", "x") - call register_axis(GFS_restart, "lev", size(o3mr_gfs,3)) call register_axis(GFS_restart, "levp", size(zh_gfs,3)) #ifdef MULTI_GASES - call register_restart_field(GFS_restart, 'spfo3', spfo3_gfs, dim_names_3d3, is_optional=.true.) - call register_restart_field(GFS_restart, 'spfo', spfo_gfs, dim_names_3d3, is_optional=.true.) - call register_restart_field(GFS_restart, 'spfo2', spfo2_gfs, dim_names_3d3, is_optional=.true.) + call register_restart_field(GFS_restart, 'spo3', spo3_gfs, dim_names_3d3, is_optional=.true.) + call register_restart_field(GFS_restart, 'spo', spo_gfs, dim_names_3d3, is_optional=.true.) + call register_restart_field(GFS_restart, 'spo2', spo2_gfs, dim_names_3d3, is_optional=.true.) #else + call register_axis(GFS_restart, "lev", size(o3mr_gfs,3)) call register_restart_field(GFS_restart, 'o3mr', o3mr_gfs, dim_names_3d3, is_optional=.true.) #endif call register_restart_field(GFS_restart, 'ps', ps_gfs, dim_names_2d) @@ -2117,18 +2117,18 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) if ( bk_gfs(1) < 1.E-9 ) ak_gfs(1) = max(1.e-9, ak_gfs(1)) #ifdef MULTI_GASES - iq = spfo3 - if(is_master()) write(*,*) 'Reading spfo3 from GFS_data.nc:' - if(is_master()) write(*,*) 'spfo3 =', iq - call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo3_gfs, zh_gfs, iq) - iq = spfo - if(is_master()) write(*,*) 'Reading spfo from GFS_data.nc:' - if(is_master()) write(*,*) 'spfo =', iq - call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo_gfs, zh_gfs, iq) - iq = spfo2 - if(is_master()) write(*,*) 'Reading spfo2 from GFS_data.nc:' - if(is_master()) write(*,*) 'spfo2 =', iq - call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo2_gfs, zh_gfs, iq) + iq = spo + if(is_master()) write(*,*) 'Reading spo from GFS_data.nc:' + if(is_master()) write(*,*) 'spo =', iq + call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spo_gfs, zh_gfs, iq) + iq = spo2 + if(is_master()) write(*,*) 'Reading spo2 from GFS_data.nc:' + if(is_master()) write(*,*) 'spo2 =', iq + call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spo2_gfs, zh_gfs, iq) + iq = spo3 + if(is_master()) write(*,*) 'Reading spo3 from GFS_data.nc:' + if(is_master()) write(*,*) 'spo3 =', iq + call remap_scalar_single(Atm, levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spo3_gfs, zh_gfs, iq) #else iq = o3mr if(is_master()) write(*,*) 'Reading o3mr from GFS_data.nc:' @@ -2139,9 +2139,9 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) deallocate (ak_gfs, bk_gfs) deallocate (ps_gfs, zh_gfs) #ifdef MULTI_GASES - deallocate (spfo3_gfs) - deallocate ( spfo_gfs) - deallocate (spfo2_gfs) + deallocate ( spo_gfs) + deallocate (spo2_gfs) + deallocate (spo3_gfs) #else deallocate (o3mr_gfs) #endif @@ -2959,7 +2959,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in) integer i,j,k,l,m, k2,iq integer sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt, sgs_tke #ifdef MULTI_GASES - integer spfo, spfo2, spfo3 + integer spo, spo2, spo3 #else integer o3mr #endif @@ -2978,9 +2978,9 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in) graupel = get_tracer_index(MODEL_ATMOS, 'graupel') cld_amt = get_tracer_index(MODEL_ATMOS, 'cld_amt') #ifdef MULTI_GASES - spfo = get_tracer_index(MODEL_ATMOS, 'spfo') - spfo2 = get_tracer_index(MODEL_ATMOS, 'spfo2') - spfo3 = get_tracer_index(MODEL_ATMOS, 'spfo3') + spo = get_tracer_index(MODEL_ATMOS, 'spo') + spo2 = get_tracer_index(MODEL_ATMOS, 'spo2') + spo3 = get_tracer_index(MODEL_ATMOS, 'spo3') #else o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') #endif @@ -2993,9 +2993,9 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in) print *, 'sphum = ', sphum print *, 'clwmr = ', liq_wat #ifdef MULTI_GASES - print *, 'spfo3 = ', spfo3 - print *, ' spfo = ', spfo - print *, 'spfo2 = ', spfo2 + print *, 'spo = ', spo + print *, 'spo2 = ', spo2 + print *, 'spo3 = ', spo3 #else print *, ' o3mr = ', o3mr #endif diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 859fb0dbe..0fd969277 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -176,7 +176,12 @@ module fv_diagnostics_mod logical :: prt_minmax =.false. logical :: m_calendar integer sphum, liq_wat, ice_wat, cld_amt ! GFDL physics - integer rainwat, snowwat, graupel, o3mr + integer rainwat, snowwat, graupel +#ifdef MULTI_GASES + integer spo, spo2, spo3 +#else + integer o3mr +#endif integer :: istep, mp_top real :: ptop real, parameter :: rad2deg = 180./pi @@ -211,7 +216,8 @@ module fv_diagnostics_mod integer :: yr_init, mo_init, dy_init, hr_init, mn_init, sec_init integer :: id_dx, id_dy - real :: vrange(2), vsrange(2), wrange(2), trange(2), slprange(2), rhrange(2), skrange(2) + real,dimension(2) :: vrange, vsrange, wrange, trange, slprange, rhrange, skrange + real,dimension(2) :: vrange_bad, trange_bad ! integer :: id_d_grid_ucomp, id_d_grid_vcomp ! D grid winds ! integer :: id_c_grid_ucomp, id_c_grid_vcomp ! C grid winds @@ -282,23 +288,39 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat') snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat') graupel = get_tracer_index (MODEL_ATMOS, 'graupel') +#ifdef MULTI_GASES + spo = get_tracer_index (MODEL_ATMOS, 'spo') + spo2 = get_tracer_index (MODEL_ATMOS, 'spo2') + spo3 = get_tracer_index (MODEL_ATMOS, 'spo3') +#else o3mr = get_tracer_index (MODEL_ATMOS, 'o3mr') +#endif cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt') ! valid range for some fields !!! This will need mods for more than 1 tile per pe !!! - vsrange = (/ -200., 200. /) ! surface (lowest layer) winds + if (Atm(1)%flagstruct%molecular_diffusion) then + vrange = (/ -850., 850. /) ! winds + wrange = (/ -300., 300. /) ! vertical wind + trange = (/ 5., 3500. /) ! temperature + vrange_bad = (/ -850., 850. /) ! winds + trange_bad = (/ 130., 3500. /) ! temperature + else vrange = (/ -330., 330. /) ! winds wrange = (/ -100., 100. /) ! vertical wind - rhrange = (/ -10., 150. /) ! RH + vrange_bad = (/ -250., 250. /) ! winds #ifdef HIWPP trange = (/ 5., 350. /) ! temperature + trange_bad = (/ 130., 350. /) ! temperature #else trange = (/ 100., 350. /) ! temperature + trange_bad = (/ 150., 350. /) ! temperature #endif + endif + rhrange = (/ -10., 150. /) ! RH slprange = (/800., 1200./) ! sea-level-pressure skrange = (/ -10000000.0, 10000000.0 /) ! dissipation estimate for SKEB @@ -1092,8 +1114,17 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) 'vertical liquid water flux', 'kg/m**2/s', missing_value=missing_value ) id_qiw = register_diag_field ( trim(field), 'qiw', axes(1:3), Time, & 'vertical ice water flux', 'kg/m**2/s', missing_value=missing_value ) +#ifdef MULTI_GASES + id_spow = register_diag_field ( trim(field), 'spow', axes(1:3), Time, & + 'vertical oxygen atom flux', 'kg/m**2/s', missing_value=missing_value ) + id_spo2w = register_diag_field ( trim(field), 'spo2w', axes(1:3), Time, & + 'vertical oxygen flux', 'kg/m**2/s', missing_value=missing_value ) + id_spo3w = register_diag_field ( trim(field), 'spo3w', axes(1:3), Time, & + 'vertical ozone flux', 'kg/m**2/s', missing_value=missing_value ) +#else id_o3w = register_diag_field ( trim(field), 'o3w', axes(1:3), Time, & 'vertical ozone flux', 'kg/m**2/s', missing_value=missing_value ) +#endif endif ! Total energy (only when moist_phys = .T.) @@ -1700,17 +1731,19 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) call range_check('DELP', Atm(n)%delp, isc, iec, jsc, jec, ngc, npz, Atm(n)%gridstruct%agrid, & 0.01*ptop, 200.E2, bad_range, Time) call range_check('UA', Atm(n)%ua, isc, iec, jsc, jec, ngc, npz, Atm(n)%gridstruct%agrid, & - -250., 250., bad_range, Time) + vrange_bad(1), vrange_bad(2), bad_range, Time) call range_check('VA', Atm(n)%va, isc, iec, jsc, jec, ngc, npz, Atm(n)%gridstruct%agrid, & - -250., 250., bad_range, Time) + vrange_bad(1),vrange_bad(2), bad_range, Time) + #ifndef SW_DYNAMICS call range_check('TA', Atm(n)%pt, isc, iec, jsc, jec, ngc, npz, Atm(n)%gridstruct%agrid, & -#ifdef HIWPP - 130., 350., bad_range, Time) !DCMIP ICs have very low temperatures +#ifdef MULTI_GASES + 130., 3500., bad_range, Time) #else - 150., 350., bad_range, Time) -#endif + trange_bad(1),trange_bad(2), bad_range, Time) #endif +#endif ! SW_DYNAMICS + call range_check('Qv', Atm(n)%q(:,:,:,sphum), isc, iec, jsc, jec, ngc, npz, Atm(n)%gridstruct%agrid, & -1.e-8, 1.e20, bad_range, Time) @@ -2949,7 +2982,14 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) if(id_ua > 0) used=send_data(id_ua, Atm(n)%ua(isc:iec,jsc:jec,:), Time) if(id_va > 0) used=send_data(id_va, Atm(n)%va(isc:iec,jsc:jec,:), Time) - if(id_hw > 0 .or. id_qvw > 0 .or. id_qlw > 0 .or. id_qiw > 0 .or. id_o3w > 0 ) then +#ifdef MULTI_GASES + if( id_hw > 0 .or. id_qvw > 0 .or. & + id_qlw > 0 .or. id_qiw > 0 .or. id_spo3w > 0 .or. & + id_spo2w > 0 .or. id_spow > 0 ) then +#else + if( id_hw > 0 .or. id_qvw > 0 .or. & + id_qlw > 0 .or. id_qiw > 0 .or. id_o3w > 0 ) then +#endif allocate( a3(isc:iec,jsc:jec,npz) ) do k=1,npz @@ -3015,6 +3055,47 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo used = send_data(id_qiw, a3, Time) endif +#ifdef MULTI_GASES + if (id_spow > 0) then + if (spo < 0) then + call mpp_error(FATAL, 'ow does not work without spo defined') + endif + do k=1,npz + do j=jsc,jec + do i=isc,iec + a3(i,j,k) = Atm(n)%q(i,j,k,spo)*wk(i,j,k) + enddo + enddo + enddo + used = send_data(id_spow, a3, Time) + endif + if (id_spo2w > 0) then + if (spo2 < 0) then + call mpp_error(FATAL, 'o2w does not work without spo2 defined') + endif + do k=1,npz + do j=jsc,jec + do i=isc,iec + a3(i,j,k) = Atm(n)%q(i,j,k,spo2)*wk(i,j,k) + enddo + enddo + enddo + used = send_data(id_spo2w, a3, Time) + endif + if (id_spo3w > 0) then + if (spo3 < 0) then + call mpp_error(FATAL, 'o3w does not work without spo3 defined') + endif + do k=1,npz + do j=jsc,jec + do i=isc,iec + a3(i,j,k) = Atm(n)%q(i,j,k,spo3)*wk(i,j,k) + enddo + enddo + enddo + used = send_data(id_spo3w, a3, Time) + endif +#else if (id_o3w > 0) then if (o3mr < 0) then call mpp_error(FATAL, 'o3w does not work without o3mr defined') @@ -3028,6 +3109,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo used = send_data(id_o3w, a3, Time) endif +#endif deallocate(a3) endif diff --git a/tools/fv_diagnostics.h b/tools/fv_diagnostics.h index 52107c790..9622577a8 100644 --- a/tools/fv_diagnostics.h +++ b/tools/fv_diagnostics.h @@ -96,6 +96,10 @@ integer :: id_t_dt_nudge, id_ps_dt_nudge, id_delp_dt_nudge, id_u_dt_nudge, id_v_dt_nudge ! EMC additions - integer :: id_diss, id_zratio, id_hw, id_qvw, id_qlw, id_qiw, id_o3w - + integer :: id_diss, id_zratio, id_hw, id_qvw, id_qlw, id_qiw +#ifdef MULTI_GASES + integer :: id_spo2w, id_spow, id_spo3w +#else + integer :: id_o3w +#endif #endif _FV_DIAG__ diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index 0ed76a033..d0d8a0a4b 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -576,6 +576,10 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, grid real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c + real(kind=R_GRID), pointer, dimension(:,:) :: area_u, area_v + real(kind=R_GRID), pointer, dimension(:,:) :: dx6, dy6 + real, pointer, dimension(:,:) :: rarea_u, rarea_v + real, pointer, dimension(:,:) :: rdx6, rdy6 real(kind=R_GRID), pointer, dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya real, pointer, dimension(:,:,:) :: e1, e2 @@ -615,6 +619,18 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, rarea => Atm%gridstruct%rarea rarea_c => Atm%gridstruct%rarea_c +! For MOLECULAR_DIFFUSION + if ( Atm%flagstruct%molecular_diffusion ) then + area_u => Atm%gridstruct%area_u_64 + area_v => Atm%gridstruct%area_v_64 + dx6 => Atm%gridstruct%dx6_64 + dy6 => Atm%gridstruct%dy6_64 + rarea_u => Atm%gridstruct%rarea_u + rarea_v => Atm%gridstruct%rarea_v + rdx6 => Atm%gridstruct%rdx6 + rdy6 => Atm%gridstruct%rdy6 + endif + sina => Atm%gridstruct%sina_64 cosa => Atm%gridstruct%cosa_64 dx => Atm%gridstruct%dx_64 @@ -956,6 +972,84 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, dyc(i,jed+1) = dyc(i,jed) end do + if ( Atm%flagstruct%molecular_diffusion ) then +! dx6, dy6 + do j=jsd,jed+1 + do i=isd+1,ied + call mid_pt_sphere(grid(i-1,j,1:2), grid(i ,j,1:2), p1) + call mid_pt_sphere(grid(i ,j,1:2), grid(i+1,j,1:2), p2) + dx6(i,j) = great_circle_dist( p2, p1, radius ) + enddo +!xxxxxx + !Are the following 2 lines appropriate for the regional domain? +!xxxxxx + dx6(isd ,j) = dx6(isd+1,j) + dx6(ied+1,j) = dx6(ied ,j) + enddo + + do j=jsd+1,jed + do i=isd,ied+1 + call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j ,1:2), p1) + call mid_pt_sphere(grid(i,j ,1:2), grid(i,j+1,1:2), p2) + dy6(i,j) = great_circle_dist( p2, p1, radius ) + enddo + enddo +!xxxxxx + !Are the following 2 lines appropriate for the regional domain? +!xxxxxx + do i=isd,ied+1 + dy6(i,jsd) = dy6(i,jsd+1) + dy6(i,jed+1) = dy6(i,jed) + enddo +! area_u, area_v + do j=jsd,jed + do i=isd+1,ied + call mid_pt_sphere(grid(i-1,j ,1:2), grid(i ,j ,1:2), p1) + call mid_pt_sphere(grid(i-1,j+1,1:2), grid(i ,j+1,1:2), p4) + call mid_pt_sphere(grid(i ,j ,1:2), grid(i+1,j ,1:2), p2) + call mid_pt_sphere(grid(i ,j+1,1:2), grid(i+1,j+1,1:2), p3) + area_v(i,j) = get_area(p1, p4, p2, p3, radius) + enddo +!xxxxxx + !Are the following 2 lines appropriate for the regional domain? +!xxxxxx + area_v(isd ,j) = area_v(isd+1,j) + area_v(ied+1,j) = area_v(ied ,j) + enddo + + do j=jsd+1,jed + do i=isd,ied + call mid_pt_sphere(grid(i ,j-1,1:2), grid(i ,j ,1:2), p1) + call mid_pt_sphere(grid(i ,j ,1:2), grid(i ,j+1,1:2), p4) + call mid_pt_sphere(grid(i+1,j-1,1:2), grid(i+1,j ,1:2), p2) + call mid_pt_sphere(grid(i+1,j ,1:2), grid(i+1,j+1,1:2), p3) + area_u(i,j) = get_area(p1, p4, p2, p3, radius) + enddo + enddo +!xxxxxx + !Are the following 2 lines appropriate for the regional domain? +!xxxxxx + do i=isd,ied + area_u(i,jsd) = area_u(i,jsd+1) + area_u(i,jed+1) = area_u(i,jed) + enddo + +! we are not really using the outmost, so no need to tune this +! so update and fill corners should be enough + + call mpp_update_domains( dx6, Atm%domain, position=CORNER, complete=.true.) + call mpp_update_domains( dy6, Atm%domain, position=CORNER, complete=.true.) + call mpp_update_domains( area_v, area_u, Atm%domain, flags=SCALAR_PAIR, & + gridtype=CGRID_NE_PARAM, complete=.true.) + + if (cubed_sphere .and. (.not. (Atm%neststruct%nested .or. Atm%flagstruct%regional))) then + call fill_corners( dx6, npx, npy, FILL=XDir, BGRID=.true.) + call fill_corners( dy6, npx, npy, FILL=XDir, BGRID=.true.) + call fill_corners( area_v, area_u, npx, npy, CGRID=.true.) + endif + + endif ! MOLECULAR_DIFFUSION + if( .not. stretched_grid ) & call sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, & @@ -1143,6 +1237,29 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, enddo enddo + if ( Atm%flagstruct%molecular_diffusion ) then + do j=jsd,jed+1 + do i=isd,ied + rarea_u(i,j) = 1.0/area_u(i,j) + enddo + enddo + do j=jsd,jed + do i=isd,ied+1 + rarea_v(i,j) = 1.0/area_v(i,j) + enddo + enddo + do j=jsd,jed+1 + do i=isd,ied+1 + rdx6(i,j) = 1.0/dx6(i,j) + enddo + enddo + do j=jsd,jed+1 + do i=isd,ied+1 + rdy6(i,j) = 1.0/dy6(i,j) + enddo + enddo + endif + 200 format(A,f9.2,A,f9.2,A,f9.2) 201 format(A,f9.2,A,f9.2,A,f9.2,A,f9.2) 202 format(A,A,i4.4,A,i4.4,A) @@ -1244,6 +1361,17 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, nullify( area_c) nullify(rarea_c) + if ( Atm%flagstruct%molecular_diffusion ) then + nullify( area_u) + nullify( area_v) + nullify(rarea_u) + nullify(rarea_v) + nullify( dx6) + nullify( dy6) + nullify(rdx6) + nullify(rdy6) + endif + nullify(sina) nullify(cosa) nullify(dx) diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index c4036e4f1..4fbee842d 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -197,7 +197,7 @@ subroutine read_da_inc(Atm, fv_domain, bd, npz_in, nq, u, v, q, delp, pt, delz, integer :: isd, ied, jsd, jed integer :: sphum, liq_wat #ifdef MULTI_GASES - integer :: spfo, spfo2, spfo3 + integer :: spo, spo2, spo3 #else integer :: o3mr #endif @@ -267,9 +267,9 @@ subroutine read_da_inc(Atm, fv_domain, bd, npz_in, nq, u, v, q, delp, pt, delz, sphum = get_tracer_index(MODEL_ATMOS, 'sphum') #ifdef MULTI_GASES - spfo3 = get_tracer_index(MODEL_ATMOS, 'spfo3') - spfo = get_tracer_index(MODEL_ATMOS, 'spfo') - spfo2 = get_tracer_index(MODEL_ATMOS, 'spfo2') + spo3 = get_tracer_index(MODEL_ATMOS, 'spo3') + spo = get_tracer_index(MODEL_ATMOS, 'spo') + spo2 = get_tracer_index(MODEL_ATMOS, 'spo2') #else o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') #endif @@ -287,9 +287,9 @@ subroutine read_da_inc(Atm, fv_domain, bd, npz_in, nq, u, v, q, delp, pt, delz, call apply_inc_on_3d_scalar('sphum_inc',q(:,:,:,sphum), is_in, js_in, ie_in, je_in) call apply_inc_on_3d_scalar('liq_wat_inc',q(:,:,:,liq_wat), is_in, js_in, ie_in, je_in) #ifdef MULTI_GASES - call apply_inc_on_3d_scalar('spfo3_inc',q(:,:,:,spfo3), is_in, js_in, ie_in, je_in) - call apply_inc_on_3d_scalar('spfo_inc',q(:,:,:,spfo), is_in, js_in, ie_in, je_in) - call apply_inc_on_3d_scalar('spfo2_inc',q(:,:,:,spfo2), is_in, js_in, ie_in, je_in) + call apply_inc_on_3d_scalar('spo_inc',q(:,:,:,spo), is_in, js_in, ie_in, je_in) + call apply_inc_on_3d_scalar('spo2_inc',q(:,:,:,spo2), is_in, js_in, ie_in, je_in) + call apply_inc_on_3d_scalar('spo3_inc',q(:,:,:,spo3), is_in, js_in, ie_in, je_in) #else call apply_inc_on_3d_scalar('o3mr_inc',q(:,:,:,o3mr), is_in, js_in, ie_in, je_in) #endif