diff --git a/GFDL_tools/fv_cmip_diag.F90 b/GFDL_tools/fv_cmip_diag.F90 index 16096ce5c..fdbdae99f 100644 --- a/GFDL_tools/fv_cmip_diag.F90 +++ b/GFDL_tools/fv_cmip_diag.F90 @@ -28,7 +28,7 @@ module fv_cmip_diag_mod use fms_io_mod, only: set_domain, nullify_domain, string use time_manager_mod, only: time_type use mpp_domains_mod, only: domain2d -use diag_manager_mod, only: register_diag_field, & +use diag_manager_mod, only: register_diag_field, diag_axis_init, & send_data, get_diag_field_id, & register_static_field, & diag_field_add_attribute, & @@ -36,13 +36,13 @@ module fv_cmip_diag_mod use diag_data_mod, only: CMOR_MISSING_VALUE, null_axis_id use tracer_manager_mod, only: get_tracer_index use field_manager_mod, only: MODEL_ATMOS -use constants_mod, only: GRAV +use constants_mod, only: GRAV, RDGAS use fv_mapz_mod, only: E_Flux use fv_arrays_mod, only: fv_atmos_type use fv_diagnostics_mod, only: interpolate_vertical, & get_height_given_pressure, & - rh_calc, get_height_field + rh_calc, get_height_field, get_vorticity use atmos_cmip_diag_mod, only: register_cmip_diag_field_2d, & register_cmip_diag_field_3d, & @@ -58,7 +58,7 @@ module fv_cmip_diag_mod public :: fv_cmip_diag_init, fv_cmip_diag, fv_cmip_diag_end -integer :: sphum +integer :: sphum, nql, nqi, nqa !----------------------------------------------------------------------- !--- namelist --- @@ -71,12 +71,15 @@ module fv_cmip_diag_mod type(cmip_diag_id_type) :: ID_ta, ID_ua, ID_va, ID_hus, ID_hur, ID_wap, ID_zg, & ID_u2, ID_v2, ID_t2, ID_wap2, ID_uv, ID_ut, ID_vt, & - ID_uwap, ID_vwap, ID_twap + ID_uwap, ID_vwap, ID_twap, ID_wa, & + ID_cls, ID_clws, ID_clis + integer :: id_ps, id_orog integer :: id_ua200, id_va200, id_ua850, id_va850, & id_ta500, id_ta700, id_ta850, id_zg500, & id_zg100, id_zg10, id_zg1000, & id_hus850, id_wap500, id_ua10 +integer :: id_rv200, id_rv500, id_rv850, id_vortmean character(len=5) :: mod_name = 'atmos' @@ -107,6 +110,7 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) integer, dimension(7) :: id_plevels integer, parameter :: id_p10=1, id_p100=2, id_p200=3, id_p500=4, id_p700=5, id_p850=6, id_p1000=7 character(len=4) :: plabel +integer :: id_pl700, id_pl700_bnds, id_nv !----------------------------------------------------------------------- if (module_is_initialized) then @@ -151,6 +155,9 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) !----------------------------------------------------------------------- sphum = get_tracer_index (MODEL_ATMOS, 'sphum') + nql = get_tracer_index (MODEL_ATMOS, 'liq_wat') + nqi = get_tracer_index (MODEL_ATMOS, 'ice_wat') + nqa = get_tracer_index (MODEL_ATMOS, 'cld_amt') !----------------------------------------------------------------------- ! register cmip 3D variables (on model levels and pressure levels) @@ -163,7 +170,7 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) ID_va = register_cmip_diag_field_3d (mod_name, 'va', Time, & 'Northward Wind', 'm s-1', standard_name='northward_wind') - + ID_hus = register_cmip_diag_field_3d (mod_name, 'hus', Time, & 'Specific Humidity', '1.0', standard_name='specific_humidity') @@ -173,6 +180,9 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) ID_hur = register_cmip_diag_field_3d (mod_name, 'hur', Time, & 'Relative Humidity', '%', standard_name='relative_humidity') + ID_wa = register_cmip_diag_field_3d (mod_name, 'wa', Time, & + 'Upward Air Velocity', 'm s-1', standard_name='upward_air_velocity') + ID_zg = register_cmip_diag_field_3d (mod_name, 'zg', Time, & 'Geopotential Height', 'm', standard_name='geopotential_height', axis='half') @@ -215,6 +225,24 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) 'Air Temperature times Omega', 'K m s-1', & standard_name='product_of_omega_and_air_temperature') +!----------------------------------------------------------------------- +! stratiform cloud tracers + + if (nql > 0) then + ID_clws = register_cmip_diag_field_3d (mod_name, 'clws', Time, & + 'Mass Fraction of Stratiform Cloud Liquid Water', '1.0', & + standard_name='mass_fraction_of_stratiform_cloud_liquid_water_in_air') + endif + if (nqi > 0) then + ID_clis = register_cmip_diag_field_3d (mod_name, 'clis', Time, & + 'Mass Fraction of Stratiform Cloud Ice', '1.0', & + standard_name='mass_fraction_of_convective_cloud_ice_in_air') + endif + if (nqa > 0) then + ID_cls = register_cmip_diag_field_3d (mod_name, 'cls', Time, & + 'Percentage Cover of Stratiform Cloud', '%', & + standard_name='stratiform_cloud_area_fraction_in_atmosphere_layer') + endif !----------------------------------------------------------------------- ! 2D fields @@ -232,14 +260,14 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) id_orog = register_static_field (mod_name, 'orog', axes(1:2), & 'Surface Altitude', 'm', & standard_name='surface_altitude', & - area=area_id) + area=area_id, interp_method='conserve_order1') if (id_orog > 0) used = send_data (id_orog, Atm(n)%phis(isc:iec,jsc:jec)/GRAV, Time) #else !--- for now output this as 'zsurf' from fv_diagnostics --- ! id_orog = register_diag_field (mod_name, 'orog', axes(1:2), Time, & ! 'Surface Altitude', 'm', & ! standard_name='surface_altitude', & -! area=area_id) +! area=area_id, interp_method='conserve_order1') #endif !----------------------------------------------------------------------- @@ -259,6 +287,24 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) endif enddo + id_pl700 = register_static_field (mod_name, 'pl700', (/null_axis_id/), & + '700 hPa Average', 'Pa', standard_name='air_pressure') + if (id_pl700 > 0) then + call diag_field_add_attribute (id_pl700, 'axis', 'Z') + call diag_field_add_attribute (id_pl700, 'positive', 'down' ) + call diag_field_add_attribute (id_pl700, 'comment', 'average at levels 600,700,850 hPa' ) + ! add bounds + id_nv = diag_axis_init('nv', (/1.,2./), 'none', 'N', 'vertex number', set_name='nv') + id_pl700_bnds = register_static_field (mod_name, 'pl700_bnds', (/id_nv,null_axis_id/), & + '700 hPa boundaries', 'Pa', standard_name='air_pressure') + if (id_pl700_bnds > 0) then + call diag_field_add_attribute (id_pl700, 'bounds', 'pl700_bnds' ) + used = send_data (id_pl700_bnds, (/850.e2,600.e2/), Time) + endif + used = send_data (id_pl700, 700.e2, Time) + endif + + !---- register field on single pressure levels ---- id_ua10 = register_cmip_diag_field_2d (mod_name, 'ua10', Time, & @@ -311,6 +357,30 @@ subroutine fv_cmip_diag_init ( Atm, axes, Time ) if (id_hus850 > 0 .and. id_plevels(id_p850) > 0) & call diag_field_add_attribute (id_hus850, 'coordinates', 'p850') + !---- relative vorticity at 200, 500, 850 hPa ---- + id_rv200 = register_cmip_diag_field_2d (mod_name, 'rv200', Time, & + 'Relative Vorticity at 200 hPa', 's-1', standard_name='atmosphere_relative_vorticity') + if (id_rv200 > 0 .and. id_plevels(id_p200) > 0) & + call diag_field_add_attribute (id_rv200, 'coordinates', 'p200') + + id_rv500 = register_cmip_diag_field_2d (mod_name, 'rv500', Time, & + 'Relative Vorticity at 500 hPa', 's-1', standard_name='atmosphere_relative_vorticity') + if (id_rv500 > 0 .and. id_plevels(id_p500) > 0) & + call diag_field_add_attribute (id_rv500, 'coordinates', 'p500') + + id_rv850 = register_cmip_diag_field_2d (mod_name, 'rv850', Time, & + 'Relative Vorticity at 850 hPa', 's-1', standard_name='atmosphere_relative_vorticity') + if (id_rv850 > 0 .and. id_plevels(id_p850) > 0) & + call diag_field_add_attribute (id_rv850, 'coordinates', 'p850') + + !---- mean relative vorticity 600, 700, 850 hPa ---- + + id_vortmean = register_cmip_diag_field_2d (mod_name, 'vortmean', Time, & + 'Mean Relative Vorticity over 600-850 hPa', 's-1', & + standard_name='atmosphere_relative_vorticity') + if (id_vortmean > 0 .and. id_pl700 > 0) & + call diag_field_add_attribute (id_vortmean, 'coordinates', 'pl700') + !---- omega at 500 hPa ---- id_wap500 = register_cmip_diag_field_2d (mod_name, 'wap500', Time, & @@ -357,15 +427,18 @@ subroutine fv_cmip_diag ( Atm, zvir, Time ) integer :: isc, iec, jsc, jec, n, i, j, k, id integer :: ngc, npz logical :: used +logical :: compute_wa , compute_rh real, dimension(Atm(1)%bd%isc:Atm(1)%bd%iec, & - Atm(1)%bd%jsc:Atm(1)%bd%jec) :: pfull, dat2 + Atm(1)%bd%jsc:Atm(1)%bd%jec) :: pfull, dat2, & + rv850, rv700, rv600 + real, dimension(Atm(1)%bd%isc:Atm(1)%bd%iec, & Atm(1)%bd%jsc:Atm(1)%bd%jec,1) :: dat3 real, dimension(Atm(1)%bd%isc:Atm(1)%bd%iec, & Atm(1)%bd%jsc:Atm(1)%bd%jec, & - Atm(1)%npz) :: rhum + Atm(1)%npz) :: rhum, wa, rv real, dimension(Atm(1)%bd%isc:Atm(1)%bd%iec, & Atm(1)%bd%jsc:Atm(1)%bd%jec, & @@ -384,26 +457,45 @@ subroutine fv_cmip_diag ( Atm, zvir, Time ) call set_domain(Atm(n)%domain) + ! set flags for computing quantities + compute_rh = .false. + compute_wa = .false. + if (count(ID_hur%field_id(:)>0) > 0) compute_rh = .true. + if (count(ID_wa%field_id(:)>0) > 0) compute_wa = .true. + ! compute relative humidity at model levels (if needed) - if (count(ID_hur%field_id(:)>0) > 0) then + if (compute_rh .or. compute_wa) then do k=1,npz do j=jsc,jec do i=isc,iec pfull(i,j) = Atm(n)%delp(i,j,k)/(Atm(n)%peln(i,k+1,j)-Atm(n)%peln(i,k,j)) enddo enddo - call rh_calc (pfull, Atm(n)%pt(isc:iec,jsc:jec,k), & + ! compute relative humidity + if (compute_rh) then + call rh_calc (pfull, Atm(n)%pt(isc:iec,jsc:jec,k), & Atm(n)%q(isc:iec,jsc:jec,k,sphum), rhum(isc:iec,jsc:jec,k), do_cmip=.true.) + endif + ! compute vertical velocity + if (compute_wa) then + wa(isc:iec,jsc:jec,k) = -(Atm(n)%omga(isc:iec,jsc:jec,k)*Atm(n)%pt(isc:iec,jsc:jec,k)/ & + pfull(isc:iec,jsc:jec))*(RDGAS/GRAV) + endif enddo endif - ! height field (wz) if needed if (count(ID_zg%field_id(:)>0) > 0 .or. any((/id_zg10,id_zg100,id_zg500,id_zg1000/) > 0)) then call get_height_field(isc, iec, jsc, jec, ngc, npz, Atm(n)%flagstruct%hydrostatic, Atm(n)%delz, & wz, Atm(n)%pt, Atm(n)%q, Atm(n)%peln, zvir) endif + ! relative vorticity + if (any((/id_rv200,id_rv500,id_rv850,id_vortmean/) > 0)) then + call get_vorticity(isc, iec, jsc, jec, Atm(n)%bd%isd, Atm(n)%bd%ied, Atm(n)%bd%jsd, Atm(n)%bd%jed, npz, & + Atm(n)%u, Atm(n)%v, rv, Atm(n)%gridstruct%dx, Atm(n)%gridstruct%dy, Atm(n)%gridstruct%rarea) + endif + !---------------------------------------------------------------------- ! process 2D fields @@ -431,6 +523,10 @@ subroutine fv_cmip_diag ( Atm, zvir, Time ) if (query_cmip_diag_id(ID_hur)) & used = send_cmip_data_3d (ID_hur, rhum(isc:iec,jsc:jec,:), Time, phalf=Atm(n)%peln, opt=1) + ! vertical velocity + if (query_cmip_diag_id(ID_wa)) & + used = send_cmip_data_3d (ID_wa, wa(isc:iec,jsc:jec,:), Time, phalf=Atm(n)%peln, opt=1) + ! geopotential height if (query_cmip_diag_id(ID_zg)) & used = send_cmip_data_3d (ID_zg, wz, Time, phalf=Atm(n)%peln, opt=1, ext=.true.) @@ -478,6 +574,13 @@ subroutine fv_cmip_diag ( Atm, zvir, Time ) used = send_cmip_data_3d (ID_twap, Atm(n)%pt (isc:iec,jsc:jec,:)*Atm(n)%omga(isc:iec,jsc:jec,:), & Time, phalf=Atm(n)%peln, opt=1) +!---------------------------------------------------------------------- +! stratiform cloud tracers (only on model levels) + + if (query_cmip_diag_id(ID_cls)) used = send_cmip_data_3d (ID_cls, Atm(n)%q(isc:iec,jsc:jec,:,nqa)*100., Time) + if (query_cmip_diag_id(ID_clws)) used = send_cmip_data_3d (ID_clws, Atm(n)%q(isc:iec,jsc:jec,:,nql), Time) + if (query_cmip_diag_id(ID_clis)) used = send_cmip_data_3d (ID_clis, Atm(n)%q(isc:iec,jsc:jec,:,nqi), Time) + !---------------------------------------------------------------------- ! process 2D fields on specific pressure levels ! @@ -541,6 +644,26 @@ subroutine fv_cmip_diag ( Atm, zvir, Time ) used = send_data (id_wap500, dat2, Time) endif + if (id_rv200 > 0) then + call interpolate_vertical (isc, iec, jsc, jec, npz, 200.e2, Atm(n)%peln, rv, dat2) + used = send_data (id_rv200, dat2, Time) + endif + + if (id_rv500 > 0) then + call interpolate_vertical (isc, iec, jsc, jec, npz, 500.e2, Atm(n)%peln, rv, dat2) + used = send_data (id_rv500, dat2, Time) + endif + + if (id_rv850 > 0 .or. id_vortmean > 0 ) then + call interpolate_vertical (isc, iec, jsc, jec, npz, 850.e2, Atm(n)%peln, rv, rv850) + if (id_rv850 > 0) used = send_data (id_rv850, rv850, Time) + if (id_vortmean > 0) then + call interpolate_vertical (isc, iec, jsc, jec, npz, 600.e2, Atm(n)%peln, rv, rv600) + call interpolate_vertical (isc, iec, jsc, jec, npz, 700.e2, Atm(n)%peln, rv, rv700) + used = send_data (id_vortmean, (rv600+rv700+rv850)/3., Time) + endif + endif + if (id_zg10 > 0) then call get_height_given_pressure (isc, iec, jsc, jec, ngc, npz, wz, 1, (/id_zg10/), & (/log(10.e2)/), Atm(n)%peln, dat3) diff --git a/driver/GFDL/atmosphere.F90 b/driver/GFDL/atmosphere.F90 index 9168ddc5b..276f9c69b 100644 --- a/driver/GFDL/atmosphere.F90 +++ b/driver/GFDL/atmosphere.F90 @@ -150,7 +150,8 @@ module atmosphere_mod real, allocatable :: qtend(:,:,:,:) real :: mv = -1.e10 !miz - type(cmip_diag_id_type) :: ID_tnta, ID_tnhusa + type(cmip_diag_id_type) :: ID_tnta, ID_tnhusa, ID_tnt, ID_tnhus + integer :: nqv, nql, nqi, nqa integer :: mytile = 1 integer :: p_split = 1 @@ -360,6 +361,12 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Surf_diff, Grid_box) ID_tnhusa = register_cmip_diag_field_3d (mod_name, 'tnhusa', Time, & 'Tendency of Specific Humidity due to Advection', 's-1', & standard_name='tendency_of_specific_humidity_due_to_advection') + ID_tnt = register_cmip_diag_field_3d (mod_name, 'tnt', Time, & + 'Tendency of Air Temperature', 'K s-1', & + standard_name='tendency_of_air_temperature') + ID_tnhus = register_cmip_diag_field_3d (mod_name, 'tnhus', Time, & + 'Tendency of Specific Humidity', 's-1', & + standard_name='tendency_of_specific_humidity') !---allocate id_tracer_* allocate (id_tracerdt_dyn (num_tracers)) @@ -374,11 +381,24 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Surf_diff, Grid_box) endif enddo if (any(id_tracerdt_dyn(:)>0)) allocate(qtendyyf(isc:iec, jsc:jec,1:npz,num_tracers)) - if ( id_tdt_dyn>0 .or. query_cmip_diag_id(ID_tnta) ) allocate(ttend(isc:iec, jsc:jec, 1:npz)) + if ( id_tdt_dyn>0 .or. query_cmip_diag_id(ID_tnta) .or. query_cmip_diag_id(ID_tnt) ) & + allocate(ttend(isc:iec, jsc:jec, 1:npz)) if ( any((/ id_qdt_dyn, id_qldt_dyn, id_qidt_dyn, id_qadt_dyn /) > 0) .or. & - query_cmip_diag_id(ID_tnhusa) ) allocate(qtend(isc:iec, jsc:jec, 1:npz, 4)) + query_cmip_diag_id(ID_tnhusa) .or. query_cmip_diag_id(ID_tnhus) ) allocate(qtend(isc:iec, jsc:jec, 1:npz, 4)) !miz +! get tracer number for common moisture tracers + nqv = get_tracer_index(MODEL_ATMOS,'sphum') + nql = get_tracer_index(MODEL_ATMOS,'liq_wat') + nqi = get_tracer_index(MODEL_ATMOS,'ice_wat') + nqa = get_tracer_index(MODEL_ATMOS,'cld_amt') +! could zero out diagnostics if nXX = 0 + if (any((/nqv,nql,nqi,nqa/)==0)) call error_mesg ('atmosphere_mod', & + 'at least one moisture tracer (sphum,liq_wat,ice_wat,cld_amt) does not exist', FATAL ) + if (nqv > size(qtend,4)) id_qdt_dyn = 0 + if (nql > size(qtend,4)) id_qldt_dyn = 0 + if (nqi > size(qtend,4)) id_qidt_dyn = 0 + if (nqa > size(qtend,4)) id_qadt_dyn = 0 ! --- initialize clocks for dynamics, physics_down and physics_up id_dynam = mpp_clock_id ('FV dy-core', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) id_subgridz = mpp_clock_id ('FV subgrid_z',flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) @@ -406,15 +426,17 @@ subroutine atmosphere_dynamics ( Time, surf_diff ) #ifndef use_AM3_physics Surf_diff%ddp_dyn(:,:,:) = Atm(mytile)%delp(isc:iec, jsc:jec, :) Surf_diff%tdt_dyn(:,:,:) = Atm(mytile)%pt(isc:iec, jsc:jec, :) - Surf_diff%qdt_dyn(:,:,:) = Atm(mytile)%q (isc:iec, jsc:jec, :, 1) + & - Atm(mytile)%q (isc:iec, jsc:jec, :, 2) + & - Atm(mytile)%q (isc:iec, jsc:jec, :, 3) + Surf_diff%qdt_dyn(:,:,:) = Atm(mytile)%q (isc:iec, jsc:jec, :, nqv) + & + Atm(mytile)%q (isc:iec, jsc:jec, :, nql) + & + Atm(mytile)%q (isc:iec, jsc:jec, :, nqi) #endif !miz[M d0 - if ( id_tdt_dyn>0 .or. query_cmip_diag_id(ID_tnta) ) ttend(:, :, :) = Atm(mytile)%pt(isc:iec, jsc:jec, :) + if ( id_tdt_dyn>0 .or. query_cmip_diag_id(ID_tnta) .or. query_cmip_diag_id(ID_tnt) ) & + ttend(:, :, :) = Atm(mytile)%pt(isc:iec, jsc:jec, :) if ( any((/ id_qdt_dyn, id_qldt_dyn, id_qidt_dyn, id_qadt_dyn /) > 0) .or. & - query_cmip_diag_id(ID_tnhusa) ) qtend(:, :, :, :) = Atm(mytile)%q (isc:iec, jsc:jec, :, :) + query_cmip_diag_id(ID_tnhusa) .or. query_cmip_diag_id(ID_tnhus) ) & + qtend(:, :, :, :) = Atm(mytile)%q (isc:iec, jsc:jec, :, 1:size(qtend,4)) !miz do itrac = 1, num_tracers if (id_tracerdt_dyn (itrac) >0 ) & @@ -458,27 +480,25 @@ subroutine atmosphere_dynamics ( Time, surf_diff ) #ifndef use_AM3_physics Surf_diff%ddp_dyn(:,:,:) =(Atm(mytile)%delp(isc:iec,jsc:jec,:)-Surf_diff%ddp_dyn(:,:,:))/dt_atmos Surf_diff%tdt_dyn(:,:,:) =(Atm(mytile)%pt(isc:iec,jsc:jec,:) -Surf_diff%tdt_dyn(:,:,:))/dt_atmos - Surf_diff%qdt_dyn(:,:,:) =(Atm(mytile)%q (isc:iec,jsc:jec,:,1) + & - Atm(mytile)%q (isc:iec,jsc:jec,:,2) + & - Atm(mytile)%q (isc:iec,jsc:jec,:,3) - Surf_diff%qdt_dyn(:,:,:))/dt_atmos + Surf_diff%qdt_dyn(:,:,:) =(Atm(mytile)%q (isc:iec,jsc:jec,:,nqv) + & + Atm(mytile)%q (isc:iec,jsc:jec,:,nql) + & + Atm(mytile)%q (isc:iec,jsc:jec,:,nqi) - Surf_diff%qdt_dyn(:,:,:))/dt_atmos #endif !miz - if ( id_udt_dyn>0 ) used = send_data( id_udt_dyn, 2.0/dt_atmos*Atm(mytile)%ua(isc:iec,jsc:jec,:), Time) - if ( id_vdt_dyn>0 ) used = send_data( id_vdt_dyn, 2.0/dt_atmos*Atm(mytile)%va(isc:iec,jsc:jec,:), Time) - if ( id_tdt_dyn>0 .or. query_cmip_diag_id(ID_tnta) ) then - ttend = (Atm(mytile)%pt(isc:iec, jsc:jec, :) - ttend(:, :, : ))/dt_atmos - if (id_tdt_dyn>0) used = send_data(id_tdt_dyn, ttend(:,:,:), Time) - if (query_cmip_diag_id(ID_tnta)) used = send_cmip_data_3d (ID_tnta, ttend(:,:,:), Time) - endif + if (id_udt_dyn>0) used = send_data( id_udt_dyn, 2.0/dt_atmos*Atm(mytile)%ua(isc:iec,jsc:jec,:), Time) + if (id_vdt_dyn>0) used = send_data( id_vdt_dyn, 2.0/dt_atmos*Atm(mytile)%va(isc:iec,jsc:jec,:), Time) + if (id_tdt_dyn > 0) used = send_data( id_tdt_dyn, (Atm(mytile)%pt(isc:iec,jsc:jec,:)-ttend(:,:,:))/dt_atmos, Time) + if (query_cmip_diag_id(ID_tnta)) & + used = send_cmip_data_3d ( ID_tnta, (Atm(mytile)%pt(isc:iec,jsc:jec,:)-ttend(:,:,:))/dt_atmos, Time) + + if (id_qdt_dyn > 0) used = send_data( id_qdt_dyn , (Atm(mytile)%q(isc:iec,jsc:jec,:,nqv)-qtend(:,:,:,nqv))/dt_atmos, Time) + if (id_qldt_dyn > 0) used = send_data( id_qldt_dyn, (Atm(mytile)%q(isc:iec,jsc:jec,:,nql)-qtend(:,:,:,nql))/dt_atmos, Time) + if (id_qidt_dyn > 0) used = send_data( id_qidt_dyn, (Atm(mytile)%q(isc:iec,jsc:jec,:,nqi)-qtend(:,:,:,nqi))/dt_atmos, Time) + if (id_qadt_dyn > 0) used = send_data( id_qadt_dyn, (Atm(mytile)%q(isc:iec,jsc:jec,:,nqa)-qtend(:,:,:,nqa))/dt_atmos, Time) + if (query_cmip_diag_id(ID_tnhusa)) & + used = send_cmip_data_3d (ID_tnhusa, (Atm(mytile)%q(isc:iec,jsc:jec,:,nqv)-qtend(:,:,:,nqv))/dt_atmos, Time) + - if ( any((/ id_qdt_dyn, id_qldt_dyn, id_qidt_dyn, id_qadt_dyn /) > 0) .or. query_cmip_diag_id(ID_tnhusa) ) then - qtend = (Atm(mytile)%q (isc:iec, jsc:jec, :, :)- qtend(:, :, :, :))/dt_atmos - if (id_qdt_dyn > 0) used = send_data(id_qdt_dyn, qtend(:,:,:,1), Time) - if (id_qldt_dyn > 0) used = send_data(id_qldt_dyn, qtend(:,:,:,2), Time) - if (id_qidt_dyn > 0) used = send_data(id_qidt_dyn, qtend(:,:,:,3), Time) - if (id_qadt_dyn > 0) used = send_data(id_qadt_dyn, qtend(:,:,:,4), Time) - if (query_cmip_diag_id(ID_tnhusa)) used = send_cmip_data_3d (ID_tnhusa, qtend(:,:,:,1), Time) - endif !miz do itrac = 1, num_tracers @@ -733,7 +753,7 @@ subroutine get_bottom_mass ( t_bot, tr_bot, p_bot, z_bot, p_surf, slp ) p_surf(i,j) = Atm(mytile)%ps(i,j) t_bot(i,j) = Atm(mytile)%pt(i,j,npz) p_bot(i,j) = Atm(mytile)%delp(i,j,npz)/(Atm(mytile)%peln(i,npz+1,j)-Atm(mytile)%peln(i,npz,j)) - z_bot(i,j) = rrg*t_bot(i,j)*(1.+zvir*Atm(mytile)%q(i,j,npz,1)) * & + z_bot(i,j) = rrg*t_bot(i,j)*(1.+zvir*Atm(mytile)%q(i,j,npz,nqv)) * & (1. - Atm(mytile)%pe(i,npz,j)/p_bot(i,j)) enddo enddo @@ -819,9 +839,9 @@ subroutine get_stock_pe(index, value) do k=1,npz do i=isc,iec ! Warning: the following works only with AM2 physics: water vapor; cloud water, cloud ice. - wm(i,j) = wm(i,j) + Atm(mytile)%delp(i,j,k) * ( Atm(mytile)%q(i,j,k,1) + & - Atm(mytile)%q(i,j,k,2) + & - Atm(mytile)%q(i,j,k,3) ) + wm(i,j) = wm(i,j) + Atm(mytile)%delp(i,j,k) * ( Atm(mytile)%q(i,j,k,nqv) + & + Atm(mytile)%q(i,j,k,nql) + & + Atm(mytile)%q(i,j,k,nqi) ) enddo enddo enddo @@ -930,6 +950,12 @@ subroutine atmosphere_state_update (Time, Physics_tendency, Physics, Atm_block) call timing_off('TWOWAY_UPDATE') endif +!--- cmip6 total tendencies of temperature and specific humidity + if (query_cmip_diag_id(ID_tnt)) & + used = send_cmip_data_3d ( ID_tnt, (Atm(mytile)%pt(isc:iec,jsc:jec,:)-ttend(:,:,:))/dt_atmos, Time) + if (query_cmip_diag_id(ID_tnhus)) & + used = send_cmip_data_3d (ID_tnhus, (Atm(mytile)%q(isc:iec,jsc:jec,:,nqv)-qtend(:,:,:,nqv))/dt_atmos, Time) + #if !defined(ATMOS_NUDGE) && !defined(CLIMATE_NUDGE) && !defined(ADA_NUDGE) if ( .not.forecast_mode .and. Atm(mytile)%flagstruct%nudge .and. Atm(mytile)%flagstruct%na_init>0 ) then if(mod(seconds, 21600)==0) call adiabatic_init_drv (Time_prev, Time_next) @@ -1168,6 +1194,7 @@ subroutine atmos_radiation_driver_inputs (Time, Radiation, Atm_block) ! phase due to the way in which MPI interacts with nested OpenMP !---------------------------------------------------------------------- call compute_g_avg(Time, 'co2', Radiation, Atm_block) + call compute_g_avg(Time, 'ch4', Radiation, Atm_block) end subroutine atmos_radiation_driver_inputs diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 77e122ab5..3c549658d 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -54,9 +54,14 @@ module fv_arrays_mod id_f15, id_f25, id_f35, id_f45, id_ctp, & id_ppt, id_ts, id_tb, id_ctt, id_pmask, id_pmaskv2, & id_delp, id_delz, id_zratio, id_ws, id_iw, id_lw, & - id_pfhy, id_pfnh, & + id_pfhy, id_pfnh, & id_qn, id_qn200, id_qn500, id_qn850, id_qp, id_mdt, id_qdt, id_aam, id_amdt, & - id_acly, id_acl, id_acl2, id_dbz, id_maxdbz, id_basedbz, id_dbz4km + id_acly, id_acl, id_acl2, id_dbz, id_maxdbz, id_basedbz, id_dbz4km, & + id_uq, id_vq, id_wq, id_iuq, id_ivq, id_iwq, & ! moisture flux & vertical integral + id_ut, id_vt, id_wt, id_iut, id_ivt, id_iwt, & ! heat flux + id_uu, id_uv, id_uw, id_vv, id_vw, id_ww, & ! momentum flux + id_iuu, id_iuv, id_iuw, id_ivv, id_ivw, id_iww ! vertically integral of momentum flux + ! Selected p-level fields from 3D variables: integer :: id_vort200, id_vort500, id_w500, id_w700 @@ -68,8 +73,8 @@ module fv_arrays_mod integer:: id_u_plev, id_v_plev, id_t_plev, id_h_plev, id_q_plev, id_omg_plev ! IPCC diag - integer :: id_rh10, id_rh50, id_rh100, id_rh200, id_rh250, id_rh300, & - id_rh500, id_rh700, id_rh850, id_rh925, id_rh1000 + integer :: id_rh10, id_rh50, id_rh100, id_rh200, id_rh250, id_rh300, & + id_rh500, id_rh700, id_rh850, id_rh925, id_rh1000 integer :: id_rh1000_cmip, id_rh925_cmip, id_rh850_cmip, id_rh700_cmip, id_rh500_cmip, & id_rh300_cmip, id_rh250_cmip, id_rh100_cmip, id_rh50_cmip, id_rh10_cmip diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 3cf4cf8f0..556aaa8c0 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -2198,7 +2198,7 @@ subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, endif call prt_maxmin('ZS_FV3', Atm%phis, is, ie, js, je, 3, 1, 1./grav) - call prt_maxmin('ZS_GFS', gzc, is, ie, js, je, 0, 1, 1.) + call prt_maxmin('ZS_GFS', gzc, is, ie, js, je, 0, 1, 1./grav) call prt_maxmin('PS_Data', psc, is, ie, js, je, 0, 1, 0.01) call prt_maxmin('T_Data', ta, is, ie, js, je, 0, km, 1.) call prt_maxmin('q_Data', qa(is:ie,js:je,1:km,1), is, ie, js, je, 0, km, 1.) @@ -2234,7 +2234,7 @@ subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, ! gzc is height ! Note the following line, gz is actully Z (from Jeff's data). - gz(km+1) = gzc(i,j)*grav + gz(km+1) = gzc(i,j) do k=km,1,-1 gz(k) = gz(k+1) + rdgas*tp(i,k)*(pn0(i,k+1)-pn0(i,k)) enddo diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index f307195af..0fdd3c2c2 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -27,7 +27,8 @@ module fv_diagnostics_mod use time_manager_mod, only: time_type, get_date, get_time use mpp_domains_mod, only: domain2d, mpp_update_domains, DGRID_NE use diag_manager_mod, only: diag_axis_init, register_diag_field, & - register_static_field, send_data, diag_grid_init + register_static_field, send_data, diag_grid_init, & + diag_field_add_attribute use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_diag_type, fv_grid_bounds_type, & R_GRID !!! CLEANUP needs rem oval? @@ -78,7 +79,7 @@ module fv_diagnostics_mod public :: fv_diag_init, fv_time, fv_diag, prt_mxm, prt_maxmin, range_check!, id_divg, id_te public :: prt_mass, prt_minmax, ppme, fv_diag_init_gn, z_sum, sphum_ll_fix, eqv_pot, qcly0, gn - public :: get_height_given_pressure, interpolate_vertical, rh_calc, get_height_field + public :: get_height_given_pressure, interpolate_vertical, rh_calc, get_height_field, get_vorticity integer, parameter :: nplev = 31 integer :: levs(nplev) @@ -149,13 +150,13 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) vsrange = (/ -200., 200. /) ! surface (lowest layer) winds - vrange = (/ -330., 330. /) ! winds + vrange = (/ -300., 300. /) ! winds wrange = (/ -100., 100. /) ! vertical wind rhrange = (/ -10., 150. /) ! RH #ifdef HIWPP trange = (/ 5., 350. /) ! temperature #else - trange = (/ 100., 350. /) ! temperature + trange = (/ 100., 400. /) ! temperature #endif slprange = (/800., 1200./) ! sea-level-pressure @@ -289,9 +290,13 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) 'latitude', 'degrees_N' ) id_area = register_static_field ( trim(field), 'area', axes(1:2), & 'cell area', 'm**2' ) + if (id_area > 0) then + call diag_field_add_attribute (id_area, 'cell_methods', 'area: sum') + endif + #ifndef DYNAMICS_ZS idiag%id_zsurf = register_static_field ( trim(field), 'zsurf', axes(1:2), & - 'surface height', 'm' ) + 'surface height', 'm', interp_method='conserve_order1' ) #endif idiag%id_zs = register_static_field ( trim(field), 'zs', axes(1:2), & 'Original Mean Terrain', 'm' ) @@ -422,13 +427,13 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) #ifdef DYNAMICS_ZS idiag%id_zsurf = register_diag_field ( trim(field), 'zsurf', axes(1:2), Time, & - 'surface height', 'm') + 'surface height', 'm', interp_method='conserve_order1') #endif !------------------- ! Surface pressure !------------------- idiag%id_ps = register_diag_field ( trim(field), 'ps', axes(1:2), Time, & - 'surface pressure', 'Pa', missing_value=missing_value ) + 'surface pressure', 'Pa', missing_value=missing_value, range=(/40000.0, 110000.0/)) !------------------- ! Mountain torque @@ -484,16 +489,17 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) all(idiag%id_h(minloc(abs(levs-100)))>0) .or. all(idiag%id_h(minloc(abs(levs-200)))>0) .or. & all(idiag%id_h(minloc(abs(levs-250)))>0) .or. all(idiag%id_h(minloc(abs(levs-300)))>0) .or. & all(idiag%id_h(minloc(abs(levs-500)))>0) .or. all(idiag%id_h(minloc(abs(levs-700)))>0) .or. & - all(idiag%id_h(minloc(abs(levs-850)))>0) .or. all(idiag%id_h(minloc(abs(levs-1000)))>0) ) then - idiag%id_hght = 1 + all(idiag%id_h(minloc(abs(levs-850)))>0) .or. all(idiag%id_h(minloc(abs(levs-925)))>0) .or. & + all(idiag%id_h(minloc(abs(levs-1000)))>0) ) then + idiag%id_hght = 1 else - idiag%id_hght = 0 + idiag%id_hght = 0 endif !----------------------------- ! mean temp between 300-500 mb !----------------------------- idiag%id_tm = register_diag_field (trim(field), 'tm', axes(1:2), Time, & - 'mean 300-500 mb temp', 'K', missing_value=missing_value ) + 'mean 300-500 mb temp', 'K', missing_value=missing_value, range=(/140.0,400.0/) ) !------------------- ! Sea-level-pressure @@ -582,6 +588,72 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) idiag%id_ws = register_diag_field ( trim(field), 'ws', axes(1:2), Time, & 'Terrain W', 'm/s', missing_value=missing_value ) !-------------------- +! 3D flux terms +!-------------------- + idiag%id_uq = register_diag_field ( trim(field), 'uq', axes(1:3), Time, & + 'zonal moisture flux', 'Kg/Kg*m/sec', missing_value=missing_value ) + idiag%id_vq = register_diag_field ( trim(field), 'vq', axes(1:3), Time, & + 'meridional moisture flux', 'Kg/Kg*m/sec', missing_value=missing_value ) + + idiag%id_ut = register_diag_field ( trim(field), 'ut', axes(1:3), Time, & + 'zonal heat flux', 'K*m/sec', missing_value=missing_value ) + idiag%id_vt = register_diag_field ( trim(field), 'vt', axes(1:3), Time, & + 'meridional heat flux', 'K*m/sec', missing_value=missing_value ) + + idiag%id_uu = register_diag_field ( trim(field), 'uu', axes(1:3), Time, & + 'zonal flux of zonal wind', '(m/sec)^2', missing_value=missing_value ) + idiag%id_uv = register_diag_field ( trim(field), 'uv', axes(1:3), Time, & + 'zonal flux of meridional wind', '(m/sec)^2', missing_value=missing_value ) + idiag%id_vv = register_diag_field ( trim(field), 'vv', axes(1:3), Time, & + 'meridional flux of meridional wind', '(m/sec)^2', missing_value=missing_value ) + + if(.not.Atm(n)%flagstruct%hydrostatic) then + idiag%id_wq = register_diag_field ( trim(field), 'wq', axes(1:3), Time, & + 'vertical moisture flux', 'Kg/Kg*m/sec', missing_value=missing_value ) + idiag%id_wt = register_diag_field ( trim(field), 'wt', axes(1:3), Time, & + 'vertical heat flux', 'K*m/sec', missing_value=missing_value ) + idiag%id_uw = register_diag_field ( trim(field), 'uw', axes(1:3), Time, & + 'zonal flux of vertical wind', '(m/sec)^2', missing_value=missing_value ) + idiag%id_vw = register_diag_field ( trim(field), 'vw', axes(1:3), Time, & + 'meridional flux of vertical wind', '(m/sec)^2', missing_value=missing_value ) + idiag%id_ww = register_diag_field ( trim(field), 'ww', axes(1:3), Time, & + 'vertical flux of vertical wind', '(m/sec)^2', missing_value=missing_value ) + endif + +!-------------------- +! vertical integral of 3D flux terms +!-------------------- + idiag%id_iuq = register_diag_field ( trim(field), 'uq_vi', axes(1:2), Time, & + 'vertical integral of uq', 'Kg/Kg*m/sec*Pa', missing_value=missing_value ) + idiag%id_ivq = register_diag_field ( trim(field), 'vq_vi', axes(1:2), Time, & + 'vertical integral of vq', 'Kg/Kg*m/sec*Pa', missing_value=missing_value ) + + idiag%id_iut = register_diag_field ( trim(field), 'ut_vi', axes(1:2), Time, & + 'vertical integral of ut', 'K*m/sec*Pa', missing_value=missing_value ) + idiag%id_ivt = register_diag_field ( trim(field), 'vt_vi', axes(1:2), Time, & + 'vertical integral of vt', 'K*m/sec*Pa', missing_value=missing_value ) + + idiag%id_iuu = register_diag_field ( trim(field), 'uu_vi', axes(1:2), Time, & + 'vertical integral of uu', '(m/sec)^2*Pa', missing_value=missing_value ) + idiag%id_iuv = register_diag_field ( trim(field), 'uv_vi', axes(1:2), Time, & + 'vertical integral of uv', '(m/sec)^2*Pa', missing_value=missing_value ) + idiag%id_ivv = register_diag_field ( trim(field), 'vv_vi', axes(1:2), Time, & + 'vertical integral of vv', '(m/sec)^2*Pa', missing_value=missing_value ) + + if(.not.Atm(n)%flagstruct%hydrostatic) then + idiag%id_iwq = register_diag_field ( trim(field), 'wq_vi', axes(1:2), Time, & + 'vertical integral of wq', 'Kg/Kg*m/sec*Pa', missing_value=missing_value ) + idiag%id_iwt = register_diag_field ( trim(field), 'wt_vi', axes(1:2), Time, & + 'vertical integral of wt', 'K*m/sec*Pa', missing_value=missing_value ) + idiag%id_iuw = register_diag_field ( trim(field), 'uw_vi', axes(1:2), Time, & + 'vertical integral of uw', '(m/sec)^2*Pa', missing_value=missing_value ) + idiag%id_ivw = register_diag_field ( trim(field), 'vw_vi', axes(1:2), Time, & + 'vertical integral of vw', '(m/sec)^2*Pa', missing_value=missing_value ) + idiag%id_iww = register_diag_field ( trim(field), 'ww_vi', axes(1:2), Time, & + 'vertical integral of ww', '(m/sec)^2*Pa', missing_value=missing_value ) + endif + +!-------------------- ! 3D Condensate !-------------------- idiag%id_qn = register_diag_field ( trim(field), 'qn', axes(1:3), Time, & @@ -650,7 +722,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) ! 850-mb vorticity !-------------------------- idiag%id_vort850 = register_diag_field ( trim(field), 'vort850', axes(1:2), Time, & - '850-mb vorticity', '1/s', missing_value=missing_value ) + '850-mb vorticity', '1/s', missing_value=missing_value) if ( .not. Atm(n)%flagstruct%hydrostatic ) & idiag%id_w200 = register_diag_field ( trim(field), 'w200', axes(1:2), Time, & @@ -658,7 +730,12 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) idiag%id_vort200 = register_diag_field ( trim(field), 'vort200', axes(1:2), Time, & '200-mb vorticity', '1/s', missing_value=missing_value ) - +!-------------------------- +! 200-mb winds: +!-------------------------- + idiag%id_w200 = register_diag_field ( trim(field), 'w200', axes(1:2), Time, & + '200-mb w-wind', '1/s', missing_value=missing_value ) +! s200: wind speed for computing KE spectrum ! Cubed_2_latlon interpolation is more accurate, particularly near the poles, using ! winds speed (a scalar), rather than wind vectors or kinetic energy directly. idiag%id_s200 = register_diag_field ( trim(field), 's200', axes(1:2), Time, & @@ -698,7 +775,11 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) idiag%id_w2500m = register_diag_field ( trim(field), 'w2500m', axes(1:2), Time, & '2.5-km AGL w-wind', 'm/s', missing_value=missing_value ) endif - +!-------------------------- +! 850-mb winds: +!-------------------------- + idiag%id_w850 = register_diag_field ( trim(field), 'w850', axes(1:2), Time, & + '850-mb w-wind', '1/s', missing_value=missing_value ) ! helicity idiag%id_x850 = register_diag_field ( trim(field), 'x850', axes(1:2), Time, & '850-mb vertical comp. of helicity', 'm/s**2', missing_value=missing_value ) @@ -923,7 +1004,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) integer :: isd, ied, jsd, jed, npz, itrac integer :: ngc, nwater - real, allocatable :: a2(:,:),a3(:,:,:), wk(:,:,:), wz(:,:,:), ucoor(:,:,:), vcoor(:,:,:) + real, allocatable :: a2(:,:),a3(:,:,:), a4(:,:,:), wk(:,:,:), wz(:,:,:), ucoor(:,:,:), vcoor(:,:,:) real, allocatable :: slp(:,:), depress(:,:), ws_max(:,:), tc_count(:,:) real, allocatable :: u2(:,:), v2(:,:), x850(:,:), var1(:,:), var2(:,:), var3(:,:) real, allocatable :: dmmr(:,:,:), dvmr(:,:,:) @@ -2193,7 +2274,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif if ( idiag%id_u100m>0 .or. idiag%id_v100m>0 .or. idiag%id_w100m>0 .or. idiag%id_w5km>0 .or. idiag%id_w2500m>0 & - & .or. idiag%id_basedbz>0 .or. idiag%id_dbz4km>0) then + .or. idiag%id_basedbz.ne.0 .or. idiag%id_dbz4km.ne.0) then !! idiag%id_basedbz and idiag%id_dbz4km are INTEGERS if (.not.allocated(wz)) allocate ( wz(isc:iec,jsc:jec,npz+1) ) if ( Atm(n)%flagstruct%hydrostatic) then rgrav = 1. / grav @@ -2233,6 +2314,12 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) used=send_data(idiag%id_rain5km, a2, Time) if(prt_minmax) call prt_maxmin('rain5km', a2, isc, iec, jsc, jec, 0, 1, 1.) endif + if ( idiag%id_w200>0 ) then + call interpolate_vertical(isc, iec, jsc, jec, npz, & + 200.e2, Atm(n)%peln, Atm(n)%w(isc:iec,jsc:jec,:), a2) + used=send_data(idiag%id_w200, a2, Time) + endif + ! 250-mb if ( idiag%id_w5km>0 ) then call interpolate_z(isc, iec, jsc, jec, npz, 5.e3, wz, Atm(n)%w(isc:iec,jsc:jec,:), a2) used=send_data(idiag%id_w5km, a2, Time) @@ -2259,7 +2346,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) if(prt_minmax) call prt_maxmin('v100m', a2, isc, iec, jsc, jec, 0, 1, 1.) endif - if ( rainwat > 0 .and. (idiag%id_dbz>0 .or. idiag%id_maxdbz>0 .or. idiag%id_basedbz>0 .or. idiag%id_dbz4km>0)) then + if ( rainwat > 0 .and. (idiag%id_dbz>0 .or. idiag%id_maxdbz>0 .or. idiag%id_basedbz>0 .or. idiag%id_dbz4km.ne.0)) then if (.not. allocated(a3)) allocate(a3(isc:iec,jsc:jec,npz)) @@ -2404,6 +2491,11 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) pout, wz, Atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1) used=send_data(idiag%id_omg_plev, a3(isc:iec,jsc:jec,:), Time) endif + if ( idiag%id_x850>0 .and. idiag%id_vort850>0 ) then + x850(:,:) = x850(:,:)*a2(:,:) + used=send_data(idiag%id_x850, x850, Time) + deallocate ( x850 ) + endif if( allocated(a3) ) deallocate (a3) ! *** End cs_intp @@ -2473,7 +2565,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo else - call eqv_pot(a3, Atm(n)%pt, Atm(n)%delp, Atm(n)%delz, Atm(n)%peln, Atm(n)%pkz, (/Atm(n)%q(isd,jsd,1,sphum)/),& + call eqv_pot(a3, Atm(n)%pt, Atm(n)%delp, Atm(n)%delz, Atm(n)%peln, Atm(n)%pkz, (/Atm(n)%q(isd,jsd,1,sphum)/), & isc, iec, jsc, jec, ngc, npz, Atm(n)%flagstruct%hydrostatic, Atm(n)%flagstruct%moist_phys) endif @@ -2576,9 +2668,201 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif endif enddo - - - +!---------------------------------- +! compute 3D flux terms +!---------------------------------- + allocate ( a4(isc:iec,jsc:jec,npz) ) + + ! zonal moisture flux + if(idiag%id_uq > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%q(i,j,k,sphum) + enddo + enddo + enddo + used=send_data(idiag%id_uq, a4, Time) + if(idiag%id_iuq > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iuq, a2, Time) + endif + endif + ! meridional moisture flux + if(idiag%id_vq > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%q(i,j,k,sphum) + enddo + enddo + enddo + used=send_data(idiag%id_vq, a4, Time) + if(idiag%id_ivq > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_ivq, a2, Time) + endif + endif + + ! zonal heat flux + if(idiag%id_ut > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%pt(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_ut, a4, Time) + if(idiag%id_iut > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iut, a2, Time) + endif + endif + ! meridional heat flux + if(idiag%id_vt > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%pt(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_vt, a4, Time) + if(idiag%id_ivt > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_ivt, a2, Time) + endif + endif + + ! zonal flux of u + if(idiag%id_uu > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%ua(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_uu, a4, Time) + if(idiag%id_iuu > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iuu, a2, Time) + endif + endif + ! zonal flux of v + if(idiag%id_uv > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%va(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_uv, a4, Time) + if(idiag%id_iuv > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iuv, a2, Time) + endif + endif + ! meridional flux of v + if(idiag%id_vv > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%va(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_vv, a4, Time) + if(idiag%id_ivv > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_ivv, a2, Time) + endif + endif + + ! terms related with vertical wind ( Atm(n)%w ): + if(.not.Atm(n)%flagstruct%hydrostatic) then + ! vertical moisture flux + if(idiag%id_wq > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%w(i,j,k) * Atm(n)%q(i,j,k,sphum) + enddo + enddo + enddo + used=send_data(idiag%id_wq, a4, Time) + if(idiag%id_iwq > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iwq, a2, Time) + endif + endif + ! vertical heat flux + if(idiag%id_wt > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%w(i,j,k) * Atm(n)%pt(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_wt, a4, Time) + if(idiag%id_iwt > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iwt, a2, Time) + endif + endif + ! zonal flux of w + if(idiag%id_uw > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%w(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_uw, a4, Time) + if(idiag%id_iuw > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iuw, a2, Time) + endif + endif + ! meridional flux of w + if(idiag%id_vw > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%w(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_vw, a4, Time) + if(idiag%id_ivw > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_ivw, a2, Time) + endif + endif + ! vertical flux of w + if(idiag%id_ww > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%w(i,j,k) * Atm(n)%w(i,j,k) + enddo + enddo + enddo + used=send_data(idiag%id_ww, a4, Time) + if(idiag%id_iww > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(idiag%id_iww, a2, Time) + endif + endif + endif + + deallocate ( a4 ) + + #endif ! enddo ! end ntileMe do-loop diff --git a/tools/fv_eta.F90 b/tools/fv_eta.F90 index 673789bbc..9f2cf78d6 100644 --- a/tools/fv_eta.F90 +++ b/tools/fv_eta.F90 @@ -404,6 +404,7 @@ subroutine set_eta(km, ks, ptop, ak, bk) real a33(34),b33(34) ! miz: grid with enhanced surface-layer resolution real a47(48),b47(48) real a48(49),b48(49) + real a49(50),b49(50) real a52(53),b52(53) real a54(55),b54(55) real a56(57),b56(57) @@ -670,6 +671,45 @@ subroutine set_eta(km, ks, ptop, ak, bk) 0.95958, 0.97747, 0.99223, & 1.00000 / + data a49/ & + 1.00000, 2.69722, 5.17136, & + 8.89455, 14.24790, 22.07157, & + 33.61283, 50.48096, 74.79993, & + 109.40055, 158.00460, 225.44108, & + 317.89560, 443.19350, 611.11558, & + 833.74392, 1125.83405, 1505.20759, & + 1993.15829, 2614.86254, 3399.78420, & + 4382.06240, 5600.87014, 7100.73115, & + 8931.78242, 11149.97021, 13817.16841, & + 17001.20930, 20775.81856, 23967.33875, & + 25527.64563, 25671.22552, 24609.29622, & + 22640.51220, 20147.13482, 17477.63530, & + 14859.86462, 12414.92533, 10201.44191, & + 8241.50255, 6534.43202, 5066.178650, & + 3815.60705, 2758.60264, 1880.646310, & + 1169.33931, 618.47983, 225.000000, & + 10.00000, 0.00000 / + + data b49/ & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.01253, & + 0.04887, 0.10724, 0.18455, & + 0.27461, 0.36914, 0.46103, & + 0.54623, 0.62305, 0.69099, & + 0.75016, 0.80110, 0.84453, & + 0.88125, 0.91210, 0.93766, & + 0.95849, 0.97495, 0.98743, & + 0.99580, 1.00000 / + + ! High PBL resolution with top at 1 mb ! SJL modified May 7, 2013 to ptop ~ 100 mb data a54/100.00000, 254.83931, 729.54278, & @@ -1299,6 +1339,13 @@ subroutine set_eta(km, ks, ptop, ak, bk) bk(k) = b48(k) enddo + case (49) + ks = 28 + do k=1,km+1 + ak(k) = a49(k) + bk(k) = b49(k) + enddo + case (52) ks = 35 ! pint = 223 do k=1,km+1 @@ -1533,7 +1580,7 @@ subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate) real ep, es, alpha, beta, gama real, parameter:: akap = 2./7. !---- Tunable parameters: - real:: k_inc = 10 ! # of layers from bottom up to near const dz region + integer:: k_inc = 10 ! # of layers from bottom up to near const dz region real:: s0 = 0.8 ! lowest layer stretch factor !----------------------- real:: s_inc diff --git a/tools/fv_surf_map.F90 b/tools/fv_surf_map.F90 index 457d75ad6..294dd5c0f 100644 --- a/tools/fv_surf_map.F90 +++ b/tools/fv_surf_map.F90 @@ -477,7 +477,7 @@ subroutine FV3_zs_filter (bd, isd, ied, jsd, jed, npx, npy, npx_global, & real(kind=R_GRID), intent(in):: grid(isd:ied+1, jsd:jed+1,2) real(kind=R_GRID), intent(in):: agrid(isd:ied, jsd:jed, 2) - real, intent(IN):: sin_sg(9,isd:ied,jsd:jed) + real, intent(IN):: sin_sg(isd:ied,jsd:jed,9) real(kind=R_GRID), intent(IN):: stretch_fac logical, intent(IN) :: nested real, intent(inout):: phis(isd:ied,jsd,jed) @@ -545,7 +545,7 @@ subroutine two_delta_filter(npx, npy, q, area, dx, dy, dxa, dya, dxc, dyc, sin_s real, intent(in):: dya(bd%isd:bd%ied, bd%jsd:bd%jed) real, intent(in):: dxc(bd%isd:bd%ied+1,bd%jsd:bd%jed) real, intent(in):: dyc(bd%isd:bd%ied, bd%jsd:bd%jed+1) - real, intent(in):: sin_sg(9,bd%isd:bd%ied,bd%jsd:bd%jed) + real, intent(in):: sin_sg(bd%isd:bd%ied,bd%jsd:bd%jed,9) real, intent(in):: oro(bd%isd:bd%ied, bd%jsd:bd%jed) ! 0==water, 1==land logical, intent(in):: zero_ocean, check_slope logical, intent(in):: nested @@ -693,10 +693,10 @@ subroutine two_delta_filter(npx, npy, q, area, dx, dy, dxa, dya, dxc, dyc, sin_s do i=is,ie+1 ddx(i,j) = (q(i-1,j)-q(i,j))/dxc(i,j) if ( extm(i-1).and.extm(i) ) then - ddx(i,j) = 0.5*(sin_sg(3,i-1,j)+sin_sg(1,i,j))*dy(i,j)*ddx(i,j) + ddx(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*ddx(i,j) elseif ( abs(ddx(i,j)) > m_slope ) then fac = min(1., max(0.1,(abs(ddx(i,j))-m_slope)/m_slope ) ) - ddx(i,j) = fac*0.5*(sin_sg(3,i-1,j)+sin_sg(1,i,j))*dy(i,j)*ddx(i,j) + ddx(i,j) = fac*0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*ddx(i,j) else ddx(i,j) = 0. endif @@ -754,10 +754,10 @@ subroutine two_delta_filter(npx, npy, q, area, dx, dy, dxa, dya, dxc, dyc, sin_s do i=is,ie ddy(i,j) = (q(i,j-1)-q(i,j))/dyc(i,j) if ( ext2(i,j-1) .and. ext2(i,j) ) then - ddy(i,j) = 0.5*(sin_sg(4,i,j-1)+sin_sg(2,i,j))*dx(i,j)*ddy(i,j) + ddy(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*ddy(i,j) elseif ( abs(ddy(i,j))>m_slope ) then fac = min(1., max(0.1,(abs(ddy(i,j))-m_slope)/m_slope)) - ddy(i,j) = fac*0.5*(sin_sg(4,i,j-1)+sin_sg(2,i,j))*dx(i,j)*ddy(i,j) + ddy(i,j) = fac*0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*ddy(i,j) else ddy(i,j) = 0. endif