From 8a93d5a550025764456ecda303d960830db57afb Mon Sep 17 00:00:00 2001 From: Geen Date: Thu, 14 May 2020 10:11:09 +0100 Subject: [PATCH 1/3] Changes to allow components of the evaporative flux to be output as diagnostics, and to allow the wind seen by the surface heat exchange to be set to a constant value --- .../driver/solo/idealized_moist_phys.F90 | 156 ++++----------- src/coupler/surface_flux.F90 | 177 ++++++------------ 2 files changed, 98 insertions(+), 235 deletions(-) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 8115c08a2..3938e8409 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -66,12 +66,6 @@ module idealized_moist_phys_mod use rrtm_vars #endif -#ifdef SOC_NO_COMPILE - ! Socrates not included -#else -use socrates_interface_mod -use soc_constants_mod -#endif implicit none private @@ -110,7 +104,6 @@ module idealized_moist_phys_mod !s Radiation options logical :: two_stream_gray = .true. logical :: do_rrtm_radiation = .false. -logical :: do_socrates_radiation = .false. !s MiMA uses damping logical :: do_damping = .false. @@ -149,8 +142,7 @@ module idealized_moist_phys_mod land_roughness_prefactor, & gp_surface, convection_scheme, & bucket, init_bucket_depth, init_bucket_depth_land, & !RG Add bucket - max_bucket_depth_land, robert_bucket, raw_bucket, & - do_socrates_radiation + max_bucket_depth_land, robert_bucket, raw_bucket integer, parameter :: num_time_levels = 2 !RG Add bucket - number of time levels added to allow timestepping in this module @@ -179,6 +171,7 @@ module idealized_moist_phys_mod drag_m, & ! momentum drag coefficient drag_t, & ! heat drag coefficient drag_q, & ! moisture drag coefficient + rho, & ! air density at surface w_atm, & ! wind speed ustar, & ! friction velocity bstar, & ! buoyancy scale @@ -195,15 +188,7 @@ module idealized_moist_phys_mod rough, & ! roughness for vert_turb_driver albedo, & !s albedo now defined in mixed_layer_init coszen, & !s make sure this is ready for assignment in run_rrtmg - pbltop, & !s Used as an input to damping_driver, outputted from vert_turb_driver - ex_del_m, & !mp586 for 10m winds and 2m temp - ex_del_h, & !mp586 for 10m winds and 2m temp - ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, & !mp586 for 10m winds and 2m temp - u_10m, & !mp586 for 10m winds and 2m temp - v_10m, & !mp586 for 10m winds and 2m temp - q_2m, & ! Add 2m specific humidity - rh_2m ! Add 2m relative humidity + pbltop !s Used as an input to damping_driver, outputted from vert_turb_driver real, allocatable, dimension(:,:,:) :: & diff_m, & ! momentum diffusion coeff. @@ -264,20 +249,18 @@ module idealized_moist_phys_mod id_bucket_depth_conv, & ! bucket depth variation induced by convection - RG Add bucket id_bucket_depth_cond, & ! bucket depth variation induced by condensation - RG Add bucket id_bucket_depth_lh, & ! bucket depth variation induced by LH - RG Add bucket + id_w_atm, & ! wind speed - RG Add lh flux breakdown + id_drag_q, & ! moisture drag coefficient - RG Add lh flux breakdown + id_rho, & ! density at surface - RG Add lh flux breakdown + id_q_atm, & ! lowest level specific humidity - RG Add lh flux breakdown + id_q_surf, & ! surface humidity - RG Add lh flux breakdown id_rh, & ! Relative humidity id_diss_heat_ray,& ! Heat dissipated by rayleigh bottom drag if gp_surface=.True. id_z_tg, & ! Relative humidity id_cape, & - id_cin, & - id_flux_u, & ! surface flux of zonal mom. - id_flux_v, & ! surface flux of meridional mom. - id_temp_2m, & !mp586 for 10m winds and 2m temp - id_u_10m, & !mp586 for 10m winds and 2m temp - id_v_10m, & !mp586 for 10m winds and 2m temp - id_q_2m, & ! Add 2m specific humidity - id_rh_2m ! Add 2m relative humidity - -integer, allocatable, dimension(:,:) :: convflag ! indicates which qe convection subroutines are used + id_cin + +integer, allocatable, dimension(:,:) :: & convflag ! indicates which qe convection subroutines are used real, allocatable, dimension(:,:) :: rad_lat, rad_lon real, allocatable, dimension(:) :: pref, p_half_1d, ln_p_half_1d, p_full_1d,ln_p_full_1d !s pref is a reference pressure profile, which in 2006 MiMA is just the initial full pressure levels, and an extra level with the reference surface pressure. Others are only necessary to calculate pref. real, allocatable, dimension(:,:) :: capeflag !s Added for Betts Miller scheme (rather than the simplified Betts Miller scheme). @@ -445,6 +428,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(drag_m (is:ie, js:je)) allocate(drag_t (is:ie, js:je)) allocate(drag_q (is:ie, js:je)) +allocate(rho (is:ie, js:je)) allocate(w_atm (is:ie, js:je)) allocate(ustar (is:ie, js:je)) allocate(bstar (is:ie, js:je)) @@ -457,14 +441,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(dedq_atm (is:ie, js:je)) allocate(dtaudv_atm (is:ie, js:je)) allocate(dtaudu_atm (is:ie, js:je)) -allocate(ex_del_m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(ex_del_h (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(ex_del_q (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(temp_2m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(u_10m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(v_10m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp -allocate(q_2m (is:ie, js:je)) ! Add 2m specific humidity -allocate(rh_2m (is:ie, js:je)) ! Add 2m relative humidity allocate(land (is:ie, js:je)); land = .false. allocate(land_ones (is:ie, js:je)); land_ones = 0.0 allocate(avail (is:ie, js:je)); avail = .true. @@ -629,11 +605,20 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Convective Available Potential Energy','J/kg') id_cin = register_diag_field(mod_name, 'cin', & axes(1:2), Time, 'Convective Inhibition','J/kg') -id_flux_u = register_diag_field(mod_name, 'flux_u', & - axes(1:2), Time, 'Zonal momentum flux', 'Pa') -id_flux_v = register_diag_field(mod_name, 'flux_v', & - axes(1:2), Time, 'Meridional momentum flux', 'Pa') +if(.not.gp_surface) then + id_w_atm = register_diag_field(mod_name, 'wind_speed', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Lowest level wind speed','m/s') + id_drag_q = register_diag_field(mod_name, 'drag_q', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Moisture drag coefficient','none') + id_rho = register_diag_field(mod_name, 'rho', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Air density at lowest level','kg/m/m/m') + id_q_atm = register_diag_field(mod_name, 'q_atm', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Lowest level specific humidity','kg/kg') + id_q_surf = register_diag_field(mod_name, 'q_surf', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Surface specific humidity','kg/kg') +endif + if(bucket) then id_bucket_depth = register_diag_field(mod_name, 'bucket_depth', & ! RG Add bucket axes(1:2), Time, 'Depth of surface reservoir', 'm') @@ -645,24 +630,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Tendency of bucket depth induced by LH', 'm/s') endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! - -id_temp_2m = register_diag_field(mod_name, 'temp_2m', & !mp586 add 2m temp - axes(1:2), Time, 'Air temperature 2m above surface', 'K') -id_u_10m = register_diag_field(mod_name, 'u_10m', & !mp586 add 10m wind (u) - axes(1:2), Time, 'Zonal wind 10m above surface', 'm/s') -id_v_10m = register_diag_field(mod_name, 'v_10m', & !mp586 add 10m wind (v) - axes(1:2), Time, 'Meridional wind 10m above surface', 'm/s') - -!!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -id_q_2m = register_diag_field(mod_name, 'sphum_2m', & - axes(1:2), Time, 'Specific humidity 2m above surface', 'kg/kg') !Add 2m specific humidity -id_rh_2m = register_diag_field(mod_name, 'rh_2m', & - axes(1:2), Time, 'Relative humidity 2m above surface', 'percent') !Add 2m relative humidity - select case(r_conv_scheme) case(SIMPLE_BETTS_CONV) @@ -747,16 +714,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l endif #endif -#ifdef SOC_NO_COMPILE - if (do_socrates_radiation) then - call error_mesg('idealized_moist_phys','do_socrates_radiation is .true. but compiler flag -D SOC_NO_COMPILE used. Stopping.', FATAL) - endif -#else -if (do_socrates_radiation) then - call socrates_init(is, ie, js, je, num_levels, axes, Time, rad_lat, rad_lonb_2d, rad_latb_2d, Time_step_in) -endif -#endif - if(turb) then call vert_turb_driver_init (rad_lonb_2d, rad_latb_2d, ie-is+1,je-js+1, & num_levels,get_axis_id(),Time, doing_edt, doing_entrain) @@ -767,9 +724,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l id_diff_dt_vg = register_diag_field(mod_name, 'dt_vg_diffusion', & axes(1:3), Time, 'meridional wind tendency from diffusion','m/s^2') id_diff_dt_tg = register_diag_field(mod_name, 'dt_tg_diffusion', & - axes(1:3), Time, 'temperature diffusion tendency','K/s') + axes(1:3), Time, 'temperature diffusion tendency','T/s') id_diff_dt_qg = register_diag_field(mod_name, 'dt_qg_diffusion', & - axes(1:3), Time, 'moisture diffusion tendency','kg/kg/s') + axes(1:3), Time, 'moisture diffusion tendency','T/s') endif id_rh = register_diag_field ( mod_name, 'rh', & @@ -790,6 +747,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg real :: delta_t real, dimension(size(ug,1), size(ug,2), size(ug,3)) :: tg_tmp, qg_tmp, RH,tg_interp, mc, dt_ug_conv, dt_vg_conv + real, intent(in) , dimension(:,:,:), optional :: mask integer, intent(in) , dimension(:,:), optional :: kbot @@ -978,9 +936,8 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg !!$ t_surf = surface_temperature(tg(:,:,:,previous), p_full(:,:,:,current), p_half(:,:,:,current)) end if - if(.not.gp_surface) then - call surface_flux( & +call surface_flux( & tg(:,:,num_levels,previous), & grid_tracers(:,:,num_levels,previous,nsphum), & ug(:,:,num_levels,previous), & @@ -1012,6 +969,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg drag_m(:,:), & ! is intent(out) drag_t(:,:), & ! is intent(out) drag_q(:,:), & ! is intent(out) + rho(:,:), & w_atm(:,:), & ! is intent(out) ustar(:,:), & ! is intent(out) bstar(:,:), & ! is intent(out) @@ -1024,36 +982,16 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg dedq_atm(:,:), & ! is intent(out) dtaudu_atm(:,:), & ! is intent(out) dtaudv_atm(:,:), & ! is intent(out) - ex_del_m(:,:), & ! mp586 for 10m winds and 2m temp - ex_del_h(:,:), & ! mp586 for 10m winds and 2m temp - ex_del_q(:,:), & ! mp586 for 10m winds and 2m temp - temp_2m(:,:), & ! mp586 for 10m winds and 2m temp - u_10m(:,:), & ! mp586 for 10m winds and 2m temp - v_10m(:,:), & ! mp586 for 10m winds and 2m temp - q_2m(:,:), & ! Add 2m specific humidity - rh_2m(:,:), & ! Add 2m relative humidity - delta_t, & + delta_t, & land(:,:), & .not.land(:,:), & avail(:,:) ) - if(id_flux_u > 0) used = send_data(id_flux_u, flux_u, Time) - if(id_flux_v > 0) used = send_data(id_flux_v, flux_v, Time) - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! - - if(id_temp_2m > 0) used = send_data(id_temp_2m, temp_2m, Time) ! mp586 add 2m temp - if(id_u_10m > 0) used = send_data(id_u_10m, u_10m, Time) ! mp586 add 10m wind (u) - if(id_v_10m > 0) used = send_data(id_v_10m, v_10m, Time) ! mp586 add 10m wind (v) - - - !!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if(id_q_2m > 0) used = send_data(id_q_2m, q_2m, Time) ! Add 2m specific humidity - if(id_rh_2m > 0) used = send_data(id_rh_2m, rh_2m*1e2, Time) ! Add 2m relative humidity - +if(id_w_atm > 0) used = send_data(id_w_atm, w_atm, Time) ! RG Add lh flux breakdown +if(id_drag_q > 0) used = send_data(id_drag_q, drag_q, Time) ! RG Add lh flux breakdown +if(id_rho > 0) used = send_data(id_rho, rho, Time) ! RG Add lh flux breakdown +if(id_q_atm > 0) used = send_data(id_q_atm, grid_tracers(:,:,num_levels,previous,nsphum), Time) ! RG Add lh flux breakdown +if(id_q_surf > 0) used = send_data(id_q_surf, q_surf, Time) ! RG Add lh flux breakdown endif ! Now complete the radiation calculation by computing the upward and net fluxes. @@ -1080,19 +1018,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg endif #endif -#ifdef SOC_NO_COMPILE - if (do_socrates_radiation) then - call error_mesg('idealized_moist_phys','do_socrates_radiation is .true. but compiler flag -D SOC_NO_COMPILE used. Stopping.', FATAL) - endif -#else -if (do_socrates_radiation) then - ! Socrates interface - - call run_socrates(Time, Time+Time_step, rad_lat, rad_lon, tg(:,:,:,previous), grid_tracers(:,:,:,previous,nsphum), t_surf(:,:), p_full(:,:,:,current), & - p_half(:,:,:,current),z_full(:,:,:,current),z_half(:,:,:,current), albedo, dt_tg(:,:,:), net_surf_sw_down(:,:), surf_lw_down(:,:), delta_t) -endif -#endif if(gp_surface) then @@ -1101,7 +1027,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg call compute_rayleigh_bottom_drag( 1, ie-is+1, & 1, je-js+1, & Time, delta_t, & - rad_lat(:,:), dt_ug(:,:,: ), & + rad_lat(:,:), dt_ug(:,:,: ), & dt_vg(:,:,: ), & ug(:,:,:,previous), vg(:,:,:,previous), & p_half(:,:,:,previous), p_full(:,:,:,previous), & @@ -1205,8 +1131,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg if(mixed_layer_bc) then call mixed_layer( & Time, Time+Time_step, & - js, & - je, & t_surf(:,:), & ! t_surf is intent(inout) flux_t(:,:), & flux_q(:,:), & @@ -1310,12 +1234,6 @@ subroutine idealized_moist_phys_end if(mixed_layer_bc) call mixed_layer_end(t_surf, bucket_depth, bucket) if(do_damping) call damping_driver_end -#ifdef SOC_NO_COMPILE - !No need to end socrates -#else -if(do_socrates_radiation) call run_socrates_end -#endif - end subroutine idealized_moist_phys_end !================================================================================================================================= diff --git a/src/coupler/surface_flux.F90 b/src/coupler/surface_flux.F90 index 97a207f51..9bb775878 100644 --- a/src/coupler/surface_flux.F90 +++ b/src/coupler/surface_flux.F90 @@ -256,6 +256,7 @@ module surface_flux_mod logical :: use_mixing_ratio = .false. real :: gust_const = 1.0 real :: gust_min = 0.0 +real :: w_atm_const = 0.0 logical :: ncar_ocean_flux = .false. logical :: ncar_ocean_flux_orig = .false. ! for backwards compatibility logical :: raoult_sat_vap = .false. @@ -273,6 +274,7 @@ module surface_flux_mod alt_gustiness, & gust_const, & gust_min, & + w_atm_const, & old_dtaudv, & use_mixing_ratio, & ncar_ocean_flux, & @@ -343,13 +345,10 @@ subroutine surface_flux_1d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, & + cd_m, cd_t, cd_q, rho, & !RG add rho as output w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) ! ! slm Mar 28 2002 -- remove agument drag_q since it is just cd_q*wind @@ -367,12 +366,7 @@ subroutine surface_flux_1d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m ! Add 2m q and RH - - + cd_m, cd_t, cd_q, rho real, intent(inout), dimension(:) :: q_surf real, intent(inout), dimension(:) :: bucket_depth !RG Add bucket real, intent(inout), dimension(:) :: depth_change_lh_1d !RG Add bucket @@ -383,17 +377,14 @@ subroutine surface_flux_1d ( & ! ---- local constants ----------------------------------------------------- ! temperature increment and its reciprocal value for comp. of derivatives real, parameter:: del_temp=0.1, del_temp_inv=1.0/del_temp - real:: zrefm, zrefh !mp586 for 10m winds and 2m temp - ! ---- local vars ---------------------------------------------------------- real, dimension(size(t_atm(:))) :: & thv_atm, th_atm, tv_atm, thv_surf, & e_sat, e_sat1, q_sat, q_sat1, p_ratio, & t_surf0, t_surf1, u_dif, v_dif, & - rho_drag, drag_t, drag_m, drag_q, rho, & - q_atm, q_surf0, dw_atmdu, dw_atmdv, w_gust, & - e_sat_2m, q_sat_2m + rho_drag, drag_t, drag_m, drag_q, & + q_atm, q_surf0, dw_atmdu, dw_atmdv, w_gust integer :: i, nbad @@ -505,60 +496,6 @@ subroutine surface_flux_1d ( & rough_mom, rough_heat, rough_moist, w_atm, & cd_m, cd_t, cd_q, u_star, b_star, avail ) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! - - - zrefm = 10. !want winds at 10m - zrefh = 2. !want temp and q at 2m - - call mo_profile( zrefm, zrefh, z_atm, rough_mom, & - rough_heat, rough_moist, & - u_star, b_star, q_star, & - ex_del_m, ex_del_h, ex_del_q, avail ) - - -! adapted from https://github.com/mom-ocean/MOM5/blob/3702ad86f9653f4e315b98613eb824a47d89cf00/src/coupler/flux_exchange.F90#L1932 - - - ! ------- reference temp ----------- - where (avail) & - temp_2m = t_surf + (t_atm - t_surf) * ex_del_h !t_ca = canopy temperature, assuming that there is no canopy (no difference between land and ocean), t_ca = t_surf - - ! ------- reference u comp ----------- - where (avail) & - u_10m = u_atm * ex_del_m ! setting u at surface to 0. - - ! ------- reference v comp ----------- - where (avail) & - v_10m = v_atm * ex_del_m ! setting v at surface to 0. - -!!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Add 2m q and RH - - ! ------- reference specific humidity ----------- - where (avail) & - q_2m = q_surf + (q_atm - q_surf) * ex_del_q - - call escomp ( temp_2m, e_sat_2m ) - - if(use_mixing_ratio) then - ! surface mixing ratio at saturation - q_sat_2m = d622 * e_sat_2m / (p_surf - e_sat_2m) - elseif(do_simple) then - q_sat_2m = d622 * e_sat_2m / p_surf - else - ! surface specific humidity at saturation - q_sat_2m = d622 * e_sat_2m / (p_surf - d378*e_sat) - endif - - ! ------- reference relative humidity ----------- - where (avail) & - rh_2m = q_2m / q_sat_2m - - ! override with ocean fluxes from NCAR calculation if (ncar_ocean_flux .or. ncar_ocean_flux_orig) then call ncar_ocean_fluxes (w_atm, th_atm, t_surf0, q_atm, q_surf0, z_atm, & @@ -569,22 +506,50 @@ subroutine surface_flux_1d ( & ! scale momentum drag coefficient on orographic roughness cd_m = cd_m*(log(z_atm/rough_mom+1)/log(z_atm/rough_scale+1))**2 ! surface layer drag coefficients - drag_t = cd_t * w_atm - drag_q = cd_q * w_atm + !drag_t = cd_t * w_atm ! RG transfer to if loop below to allow w_atm to be fixed + !drag_q = cd_q * w_atm ! RG transfer to if loop below to allow w_atm to be fixed drag_m = cd_m * w_atm ! density rho = p_atm / (rdgas * tv_atm) - ! sensible heat flux - rho_drag = cp_air * drag_t * rho - flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) - dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) - dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) + ! sensible heat flux ! RG transfer to if loop below to allow w_atm to be fixed + !rho_drag = cp_air * drag_t * rho + !flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) + !dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) + !dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) ! evaporation - rho_drag = drag_q * rho + !rho_drag = drag_q * rho ! RG transfer to if loop below to allow w_atm to be fixed end where + + ! RG Add option to fix w_atm in the evaporation and sensible heat equations. + if (w_atm_const > 0.0) then + where (avail) + ! sensible heat flux + drag_t = cd_t * w_atm_const + rho_drag = cp_air * drag_t * rho + flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) + dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) + dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) + + drag_q = cd_q * w_atm_const + rho_drag = drag_q * rho + + end where + else + where (avail) + ! sensible heat flux + drag_t = cd_t * w_atm + rho_drag = cp_air * drag_t * rho + flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) + dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) + dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) + + drag_q = cd_q * w_atm + rho_drag = drag_q * rho + end where + end if !RG Add bucket - if bucket is on evaluate fluxes based on moisture availability. !RG Note changes to avail statements to allow bucket to be switched on or off @@ -648,11 +613,22 @@ subroutine surface_flux_1d ( & !RG end Add bucket changes + +! RG Add option to fix w_atm in the evaporation and sensible heat equations. +if (w_atm_const > 0.0) then + where (avail) + q_surf = q_atm + flux_q / (rho*cd_q*w_atm_const) ! surface specific humidity + end where +else + where (avail) + q_surf = q_atm + flux_q / (rho*cd_q*w_atm) ! surface specific humidity + end where +end if + + where (avail) q_star = flux_q / (u_star * rho) ! moisture scale - ! ask Chris and Steve K if we still want to keep this for diagnostics - q_surf = q_atm + flux_q / (rho*cd_q*w_atm) ! surface specific humidity ! upward long wave radiation flux_r = stefan*t_surf**4 ! (W/m**2) @@ -713,9 +689,6 @@ subroutine surface_flux_0d ( & w_atm_0, u_star_0, b_star_0, q_star_0, & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0, dtaudv_atm_0, & - ex_del_m_0, ex_del_h_0, ex_del_q_0, & !mp586 for 10m winds and 2m temp - temp_2m_0, u_10m_0, v_10m_0, & !mp586 for 10m winds and 2m temp - q_2m_0, rh_2m_0, & !2m q and RH dt, land_0, seawater_0, avail_0 ) ! ---- arguments ----------------------------------------------------------- @@ -730,10 +703,7 @@ subroutine surface_flux_0d ( & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0,dtaudv_atm_0, & w_atm_0, u_star_0, b_star_0, q_star_0, & - cd_m_0, cd_t_0, cd_q_0, & - ex_del_m_0, ex_del_h_0, ex_del_q_0, & !mp586 for 10m winds and 2m temp - temp_2m_0, u_10m_0, v_10m_0, & !mp586 for 10m winds and 2m temp - q_2m_0, rh_2m_0 + cd_m_0, cd_t_0, cd_q_0 real, intent(inout) :: q_surf_0 real, intent(in) :: dt @@ -750,11 +720,7 @@ subroutine surface_flux_0d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m !Add 2m q and RH - + cd_m, cd_t, cd_q, rho real, dimension(1) :: q_surf real, dimension(1) :: bucket_depth !RG Add bucket real, dimension(1) :: depth_change_lh_1d !RG Add bucket @@ -792,13 +758,10 @@ subroutine surface_flux_0d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, & + cd_m, cd_t, cd_q, rho, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) flux_t_0 = flux_t(1) @@ -822,14 +785,6 @@ subroutine surface_flux_0d ( & cd_m_0 = cd_m(1) cd_t_0 = cd_t(1) cd_q_0 = cd_q(1) - ex_del_m_0 = ex_del_m(1) !mp586 for 10m winds and 2m temp - ex_del_h_0 = ex_del_h(1) !mp586 for 10m winds and 2m temp - ex_del_q_0 = ex_del_q(1) !mp586 for 10m winds and 2m temp - temp_2m_0 = temp_2m(1) !mp586 for 10m winds and 2m temp - u_10m_0 = u_10m(1) !mp586 for 10m winds and 2m temp - v_10m_0 = v_10m(1) !mp586 for 10m winds and 2m temp - q_2m_0 = q_2m(1) !Add 2m q - rh_2m_0 = rh_2m(1) !Add 2m RH end subroutine surface_flux_0d @@ -841,13 +796,10 @@ subroutine surface_flux_2d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, & + cd_m, cd_t, cd_q, rho, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) ! ---- arguments ----------------------------------------------------------- @@ -862,11 +814,7 @@ subroutine surface_flux_2d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, & - ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp - temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp - q_2m, rh_2m !Add 2m q and RH - + cd_m, cd_t, cd_q, rho real, intent(inout), dimension(:,:) :: q_surf logical, intent(in) :: bucket !RG Add bucket real, intent(inout), dimension(:,:) :: bucket_depth ! RG Add bucket @@ -887,13 +835,10 @@ subroutine surface_flux_2d ( & u_surf(:,j), v_surf(:,j), & rough_mom(:,j), rough_heat(:,j), rough_moist(:,j), rough_scale(:,j), gust(:,j), & flux_t(:,j), flux_q(:,j), flux_r(:,j), flux_u(:,j), flux_v(:,j), & - cd_m(:,j), cd_t(:,j), cd_q(:,j), & + cd_m(:,j), cd_t(:,j), cd_q(:,j), rho(:,j), & w_atm(:,j), u_star(:,j), b_star(:,j), q_star(:,j), & dhdt_surf(:,j), dedt_surf(:,j), dedq_surf(:,j), drdt_surf(:,j), & dhdt_atm(:,j), dedq_atm(:,j), dtaudu_atm(:,j), dtaudv_atm(:,j), & - ex_del_m(:,j), ex_del_h(:,j), ex_del_q(:,j), & !mp586 for 10m winds and 2m temp - temp_2m(:,j), u_10m(:,j), v_10m(:,j), & !mp586 for 10m winds and 2m temp - q_2m(:,j), rh_2m(:,j), & dt, land(:,j), seawater(:,j), avail(:,j) ) end do end subroutine surface_flux_2d From 63dacbde760a8ea67e445593094c9769d0170789 Mon Sep 17 00:00:00 2001 From: Geen Date: Thu, 14 May 2020 10:38:09 +0100 Subject: [PATCH 2/3] Overlay edits to surface_flux.F90 on up-to-date version --- .../driver/solo/idealized_moist_phys.F90 | 156 +++++++++++++----- src/coupler/surface_flux.F90 | 132 ++++++++++++--- 2 files changed, 229 insertions(+), 59 deletions(-) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 3938e8409..8115c08a2 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -66,6 +66,12 @@ module idealized_moist_phys_mod use rrtm_vars #endif +#ifdef SOC_NO_COMPILE + ! Socrates not included +#else +use socrates_interface_mod +use soc_constants_mod +#endif implicit none private @@ -104,6 +110,7 @@ module idealized_moist_phys_mod !s Radiation options logical :: two_stream_gray = .true. logical :: do_rrtm_radiation = .false. +logical :: do_socrates_radiation = .false. !s MiMA uses damping logical :: do_damping = .false. @@ -142,7 +149,8 @@ module idealized_moist_phys_mod land_roughness_prefactor, & gp_surface, convection_scheme, & bucket, init_bucket_depth, init_bucket_depth_land, & !RG Add bucket - max_bucket_depth_land, robert_bucket, raw_bucket + max_bucket_depth_land, robert_bucket, raw_bucket, & + do_socrates_radiation integer, parameter :: num_time_levels = 2 !RG Add bucket - number of time levels added to allow timestepping in this module @@ -171,7 +179,6 @@ module idealized_moist_phys_mod drag_m, & ! momentum drag coefficient drag_t, & ! heat drag coefficient drag_q, & ! moisture drag coefficient - rho, & ! air density at surface w_atm, & ! wind speed ustar, & ! friction velocity bstar, & ! buoyancy scale @@ -188,7 +195,15 @@ module idealized_moist_phys_mod rough, & ! roughness for vert_turb_driver albedo, & !s albedo now defined in mixed_layer_init coszen, & !s make sure this is ready for assignment in run_rrtmg - pbltop !s Used as an input to damping_driver, outputted from vert_turb_driver + pbltop, & !s Used as an input to damping_driver, outputted from vert_turb_driver + ex_del_m, & !mp586 for 10m winds and 2m temp + ex_del_h, & !mp586 for 10m winds and 2m temp + ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, & !mp586 for 10m winds and 2m temp + u_10m, & !mp586 for 10m winds and 2m temp + v_10m, & !mp586 for 10m winds and 2m temp + q_2m, & ! Add 2m specific humidity + rh_2m ! Add 2m relative humidity real, allocatable, dimension(:,:,:) :: & diff_m, & ! momentum diffusion coeff. @@ -249,18 +264,20 @@ module idealized_moist_phys_mod id_bucket_depth_conv, & ! bucket depth variation induced by convection - RG Add bucket id_bucket_depth_cond, & ! bucket depth variation induced by condensation - RG Add bucket id_bucket_depth_lh, & ! bucket depth variation induced by LH - RG Add bucket - id_w_atm, & ! wind speed - RG Add lh flux breakdown - id_drag_q, & ! moisture drag coefficient - RG Add lh flux breakdown - id_rho, & ! density at surface - RG Add lh flux breakdown - id_q_atm, & ! lowest level specific humidity - RG Add lh flux breakdown - id_q_surf, & ! surface humidity - RG Add lh flux breakdown id_rh, & ! Relative humidity id_diss_heat_ray,& ! Heat dissipated by rayleigh bottom drag if gp_surface=.True. id_z_tg, & ! Relative humidity id_cape, & - id_cin - -integer, allocatable, dimension(:,:) :: & convflag ! indicates which qe convection subroutines are used + id_cin, & + id_flux_u, & ! surface flux of zonal mom. + id_flux_v, & ! surface flux of meridional mom. + id_temp_2m, & !mp586 for 10m winds and 2m temp + id_u_10m, & !mp586 for 10m winds and 2m temp + id_v_10m, & !mp586 for 10m winds and 2m temp + id_q_2m, & ! Add 2m specific humidity + id_rh_2m ! Add 2m relative humidity + +integer, allocatable, dimension(:,:) :: convflag ! indicates which qe convection subroutines are used real, allocatable, dimension(:,:) :: rad_lat, rad_lon real, allocatable, dimension(:) :: pref, p_half_1d, ln_p_half_1d, p_full_1d,ln_p_full_1d !s pref is a reference pressure profile, which in 2006 MiMA is just the initial full pressure levels, and an extra level with the reference surface pressure. Others are only necessary to calculate pref. real, allocatable, dimension(:,:) :: capeflag !s Added for Betts Miller scheme (rather than the simplified Betts Miller scheme). @@ -428,7 +445,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(drag_m (is:ie, js:je)) allocate(drag_t (is:ie, js:je)) allocate(drag_q (is:ie, js:je)) -allocate(rho (is:ie, js:je)) allocate(w_atm (is:ie, js:je)) allocate(ustar (is:ie, js:je)) allocate(bstar (is:ie, js:je)) @@ -441,6 +457,14 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(dedq_atm (is:ie, js:je)) allocate(dtaudv_atm (is:ie, js:je)) allocate(dtaudu_atm (is:ie, js:je)) +allocate(ex_del_m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(ex_del_h (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(ex_del_q (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(temp_2m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(u_10m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(v_10m (is:ie, js:je)) !mp586 added for 10m wind and 2m temp +allocate(q_2m (is:ie, js:je)) ! Add 2m specific humidity +allocate(rh_2m (is:ie, js:je)) ! Add 2m relative humidity allocate(land (is:ie, js:je)); land = .false. allocate(land_ones (is:ie, js:je)); land_ones = 0.0 allocate(avail (is:ie, js:je)); avail = .true. @@ -605,20 +629,11 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Convective Available Potential Energy','J/kg') id_cin = register_diag_field(mod_name, 'cin', & axes(1:2), Time, 'Convective Inhibition','J/kg') +id_flux_u = register_diag_field(mod_name, 'flux_u', & + axes(1:2), Time, 'Zonal momentum flux', 'Pa') +id_flux_v = register_diag_field(mod_name, 'flux_v', & + axes(1:2), Time, 'Meridional momentum flux', 'Pa') -if(.not.gp_surface) then - id_w_atm = register_diag_field(mod_name, 'wind_speed', & ! RG Add lh flux breakdown - axes(1:2), Time, 'Lowest level wind speed','m/s') - id_drag_q = register_diag_field(mod_name, 'drag_q', & ! RG Add lh flux breakdown - axes(1:2), Time, 'Moisture drag coefficient','none') - id_rho = register_diag_field(mod_name, 'rho', & ! RG Add lh flux breakdown - axes(1:2), Time, 'Air density at lowest level','kg/m/m/m') - id_q_atm = register_diag_field(mod_name, 'q_atm', & ! RG Add lh flux breakdown - axes(1:2), Time, 'Lowest level specific humidity','kg/kg') - id_q_surf = register_diag_field(mod_name, 'q_surf', & ! RG Add lh flux breakdown - axes(1:2), Time, 'Surface specific humidity','kg/kg') -endif - if(bucket) then id_bucket_depth = register_diag_field(mod_name, 'bucket_depth', & ! RG Add bucket axes(1:2), Time, 'Depth of surface reservoir', 'm') @@ -630,6 +645,24 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Tendency of bucket depth induced by LH', 'm/s') endif +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! + +id_temp_2m = register_diag_field(mod_name, 'temp_2m', & !mp586 add 2m temp + axes(1:2), Time, 'Air temperature 2m above surface', 'K') +id_u_10m = register_diag_field(mod_name, 'u_10m', & !mp586 add 10m wind (u) + axes(1:2), Time, 'Zonal wind 10m above surface', 'm/s') +id_v_10m = register_diag_field(mod_name, 'v_10m', & !mp586 add 10m wind (v) + axes(1:2), Time, 'Meridional wind 10m above surface', 'm/s') + +!!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +id_q_2m = register_diag_field(mod_name, 'sphum_2m', & + axes(1:2), Time, 'Specific humidity 2m above surface', 'kg/kg') !Add 2m specific humidity +id_rh_2m = register_diag_field(mod_name, 'rh_2m', & + axes(1:2), Time, 'Relative humidity 2m above surface', 'percent') !Add 2m relative humidity + select case(r_conv_scheme) case(SIMPLE_BETTS_CONV) @@ -714,6 +747,16 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l endif #endif +#ifdef SOC_NO_COMPILE + if (do_socrates_radiation) then + call error_mesg('idealized_moist_phys','do_socrates_radiation is .true. but compiler flag -D SOC_NO_COMPILE used. Stopping.', FATAL) + endif +#else +if (do_socrates_radiation) then + call socrates_init(is, ie, js, je, num_levels, axes, Time, rad_lat, rad_lonb_2d, rad_latb_2d, Time_step_in) +endif +#endif + if(turb) then call vert_turb_driver_init (rad_lonb_2d, rad_latb_2d, ie-is+1,je-js+1, & num_levels,get_axis_id(),Time, doing_edt, doing_entrain) @@ -724,9 +767,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l id_diff_dt_vg = register_diag_field(mod_name, 'dt_vg_diffusion', & axes(1:3), Time, 'meridional wind tendency from diffusion','m/s^2') id_diff_dt_tg = register_diag_field(mod_name, 'dt_tg_diffusion', & - axes(1:3), Time, 'temperature diffusion tendency','T/s') + axes(1:3), Time, 'temperature diffusion tendency','K/s') id_diff_dt_qg = register_diag_field(mod_name, 'dt_qg_diffusion', & - axes(1:3), Time, 'moisture diffusion tendency','T/s') + axes(1:3), Time, 'moisture diffusion tendency','kg/kg/s') endif id_rh = register_diag_field ( mod_name, 'rh', & @@ -747,7 +790,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg real :: delta_t real, dimension(size(ug,1), size(ug,2), size(ug,3)) :: tg_tmp, qg_tmp, RH,tg_interp, mc, dt_ug_conv, dt_vg_conv - real, intent(in) , dimension(:,:,:), optional :: mask integer, intent(in) , dimension(:,:), optional :: kbot @@ -936,8 +978,9 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg !!$ t_surf = surface_temperature(tg(:,:,:,previous), p_full(:,:,:,current), p_half(:,:,:,current)) end if + if(.not.gp_surface) then -call surface_flux( & + call surface_flux( & tg(:,:,num_levels,previous), & grid_tracers(:,:,num_levels,previous,nsphum), & ug(:,:,num_levels,previous), & @@ -969,7 +1012,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg drag_m(:,:), & ! is intent(out) drag_t(:,:), & ! is intent(out) drag_q(:,:), & ! is intent(out) - rho(:,:), & w_atm(:,:), & ! is intent(out) ustar(:,:), & ! is intent(out) bstar(:,:), & ! is intent(out) @@ -982,16 +1024,36 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg dedq_atm(:,:), & ! is intent(out) dtaudu_atm(:,:), & ! is intent(out) dtaudv_atm(:,:), & ! is intent(out) - delta_t, & + ex_del_m(:,:), & ! mp586 for 10m winds and 2m temp + ex_del_h(:,:), & ! mp586 for 10m winds and 2m temp + ex_del_q(:,:), & ! mp586 for 10m winds and 2m temp + temp_2m(:,:), & ! mp586 for 10m winds and 2m temp + u_10m(:,:), & ! mp586 for 10m winds and 2m temp + v_10m(:,:), & ! mp586 for 10m winds and 2m temp + q_2m(:,:), & ! Add 2m specific humidity + rh_2m(:,:), & ! Add 2m relative humidity + delta_t, & land(:,:), & .not.land(:,:), & avail(:,:) ) -if(id_w_atm > 0) used = send_data(id_w_atm, w_atm, Time) ! RG Add lh flux breakdown -if(id_drag_q > 0) used = send_data(id_drag_q, drag_q, Time) ! RG Add lh flux breakdown -if(id_rho > 0) used = send_data(id_rho, rho, Time) ! RG Add lh flux breakdown -if(id_q_atm > 0) used = send_data(id_q_atm, grid_tracers(:,:,num_levels,previous,nsphum), Time) ! RG Add lh flux breakdown -if(id_q_surf > 0) used = send_data(id_q_surf, q_surf, Time) ! RG Add lh flux breakdown + if(id_flux_u > 0) used = send_data(id_flux_u, flux_u, Time) + if(id_flux_v > 0) used = send_data(id_flux_v, flux_v, Time) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! + + if(id_temp_2m > 0) used = send_data(id_temp_2m, temp_2m, Time) ! mp586 add 2m temp + if(id_u_10m > 0) used = send_data(id_u_10m, u_10m, Time) ! mp586 add 10m wind (u) + if(id_v_10m > 0) used = send_data(id_v_10m, v_10m, Time) ! mp586 add 10m wind (v) + + + !!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + if(id_q_2m > 0) used = send_data(id_q_2m, q_2m, Time) ! Add 2m specific humidity + if(id_rh_2m > 0) used = send_data(id_rh_2m, rh_2m*1e2, Time) ! Add 2m relative humidity + endif ! Now complete the radiation calculation by computing the upward and net fluxes. @@ -1018,7 +1080,19 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg endif #endif +#ifdef SOC_NO_COMPILE + if (do_socrates_radiation) then + call error_mesg('idealized_moist_phys','do_socrates_radiation is .true. but compiler flag -D SOC_NO_COMPILE used. Stopping.', FATAL) + endif +#else +if (do_socrates_radiation) then + ! Socrates interface + + call run_socrates(Time, Time+Time_step, rad_lat, rad_lon, tg(:,:,:,previous), grid_tracers(:,:,:,previous,nsphum), t_surf(:,:), p_full(:,:,:,current), & + p_half(:,:,:,current),z_full(:,:,:,current),z_half(:,:,:,current), albedo, dt_tg(:,:,:), net_surf_sw_down(:,:), surf_lw_down(:,:), delta_t) +endif +#endif if(gp_surface) then @@ -1027,7 +1101,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg call compute_rayleigh_bottom_drag( 1, ie-is+1, & 1, je-js+1, & Time, delta_t, & - rad_lat(:,:), dt_ug(:,:,: ), & + rad_lat(:,:), dt_ug(:,:,: ), & dt_vg(:,:,: ), & ug(:,:,:,previous), vg(:,:,:,previous), & p_half(:,:,:,previous), p_full(:,:,:,previous), & @@ -1131,6 +1205,8 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg if(mixed_layer_bc) then call mixed_layer( & Time, Time+Time_step, & + js, & + je, & t_surf(:,:), & ! t_surf is intent(inout) flux_t(:,:), & flux_q(:,:), & @@ -1234,6 +1310,12 @@ subroutine idealized_moist_phys_end if(mixed_layer_bc) call mixed_layer_end(t_surf, bucket_depth, bucket) if(do_damping) call damping_driver_end +#ifdef SOC_NO_COMPILE + !No need to end socrates +#else +if(do_socrates_radiation) call run_socrates_end +#endif + end subroutine idealized_moist_phys_end !================================================================================================================================= diff --git a/src/coupler/surface_flux.F90 b/src/coupler/surface_flux.F90 index 9bb775878..29fa9bbad 100644 --- a/src/coupler/surface_flux.F90 +++ b/src/coupler/surface_flux.F90 @@ -345,10 +345,13 @@ subroutine surface_flux_1d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, rho, & !RG add rho as output + cd_m, cd_t, cd_q, rho, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) ! ! slm Mar 28 2002 -- remove agument drag_q since it is just cd_q*wind @@ -366,7 +369,12 @@ subroutine surface_flux_1d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, rho + cd_m, cd_t, cd_q, rho, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m ! Add 2m q and RH + + real, intent(inout), dimension(:) :: q_surf real, intent(inout), dimension(:) :: bucket_depth !RG Add bucket real, intent(inout), dimension(:) :: depth_change_lh_1d !RG Add bucket @@ -377,6 +385,8 @@ subroutine surface_flux_1d ( & ! ---- local constants ----------------------------------------------------- ! temperature increment and its reciprocal value for comp. of derivatives real, parameter:: del_temp=0.1, del_temp_inv=1.0/del_temp + real:: zrefm, zrefh !mp586 for 10m winds and 2m temp + ! ---- local vars ---------------------------------------------------------- real, dimension(size(t_atm(:))) :: & @@ -384,7 +394,8 @@ subroutine surface_flux_1d ( & e_sat, e_sat1, q_sat, q_sat1, p_ratio, & t_surf0, t_surf1, u_dif, v_dif, & rho_drag, drag_t, drag_m, drag_q, & - q_atm, q_surf0, dw_atmdu, dw_atmdv, w_gust + q_atm, q_surf0, dw_atmdu, dw_atmdv, w_gust, & + e_sat_2m, q_sat_2m integer :: i, nbad @@ -496,6 +507,60 @@ subroutine surface_flux_1d ( & rough_mom, rough_heat, rough_moist, w_atm, & cd_m, cd_t, cd_q, u_star, b_star, avail ) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!! added by mp586 for 10m winds and 2m temperature add mo_profile()!!!!!!!! + + + zrefm = 10. !want winds at 10m + zrefh = 2. !want temp and q at 2m + + call mo_profile( zrefm, zrefh, z_atm, rough_mom, & + rough_heat, rough_moist, & + u_star, b_star, q_star, & + ex_del_m, ex_del_h, ex_del_q, avail ) + + +! adapted from https://github.com/mom-ocean/MOM5/blob/3702ad86f9653f4e315b98613eb824a47d89cf00/src/coupler/flux_exchange.F90#L1932 + + + ! ------- reference temp ----------- + where (avail) & + temp_2m = t_surf + (t_atm - t_surf) * ex_del_h !t_ca = canopy temperature, assuming that there is no canopy (no difference between land and ocean), t_ca = t_surf + + ! ------- reference u comp ----------- + where (avail) & + u_10m = u_atm * ex_del_m ! setting u at surface to 0. + + ! ------- reference v comp ----------- + where (avail) & + v_10m = v_atm * ex_del_m ! setting v at surface to 0. + +!!!!!!!!!!!! end of mp586 additions !!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Add 2m q and RH + + ! ------- reference specific humidity ----------- + where (avail) & + q_2m = q_surf + (q_atm - q_surf) * ex_del_q + + call escomp ( temp_2m, e_sat_2m ) + + if(use_mixing_ratio) then + ! surface mixing ratio at saturation + q_sat_2m = d622 * e_sat_2m / (p_surf - e_sat_2m) + elseif(do_simple) then + q_sat_2m = d622 * e_sat_2m / p_surf + else + ! surface specific humidity at saturation + q_sat_2m = d622 * e_sat_2m / (p_surf - d378*e_sat) + endif + + ! ------- reference relative humidity ----------- + where (avail) & + rh_2m = q_2m / q_sat_2m + + ! override with ocean fluxes from NCAR calculation if (ncar_ocean_flux .or. ncar_ocean_flux_orig) then call ncar_ocean_fluxes (w_atm, th_atm, t_surf0, q_atm, q_surf0, z_atm, & @@ -506,22 +571,11 @@ subroutine surface_flux_1d ( & ! scale momentum drag coefficient on orographic roughness cd_m = cd_m*(log(z_atm/rough_mom+1)/log(z_atm/rough_scale+1))**2 ! surface layer drag coefficients - !drag_t = cd_t * w_atm ! RG transfer to if loop below to allow w_atm to be fixed - !drag_q = cd_q * w_atm ! RG transfer to if loop below to allow w_atm to be fixed drag_m = cd_m * w_atm - ! density rho = p_atm / (rdgas * tv_atm) - - ! sensible heat flux ! RG transfer to if loop below to allow w_atm to be fixed - !rho_drag = cp_air * drag_t * rho - !flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) - !dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) - !dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) - - ! evaporation - !rho_drag = drag_q * rho ! RG transfer to if loop below to allow w_atm to be fixed - end where + end where + ! RG Add option to fix w_atm in the evaporation and sensible heat equations. if (w_atm_const > 0.0) then @@ -551,6 +605,7 @@ subroutine surface_flux_1d ( & end where end if + !RG Add bucket - if bucket is on evaluate fluxes based on moisture availability. !RG Note changes to avail statements to allow bucket to be switched on or off if (bucket) then @@ -629,6 +684,8 @@ subroutine surface_flux_1d ( & where (avail) q_star = flux_q / (u_star * rho) ! moisture scale + ! ask Chris and Steve K if we still want to keep this for diagnostics + q_surf = q_atm + flux_q / (rho*cd_q*w_atm) ! surface specific humidity ! upward long wave radiation flux_r = stefan*t_surf**4 ! (W/m**2) @@ -689,6 +746,9 @@ subroutine surface_flux_0d ( & w_atm_0, u_star_0, b_star_0, q_star_0, & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0, dtaudv_atm_0, & + ex_del_m_0, ex_del_h_0, ex_del_q_0, & !mp586 for 10m winds and 2m temp + temp_2m_0, u_10m_0, v_10m_0, & !mp586 for 10m winds and 2m temp + q_2m_0, rh_2m_0, & !2m q and RH dt, land_0, seawater_0, avail_0 ) ! ---- arguments ----------------------------------------------------------- @@ -703,7 +763,10 @@ subroutine surface_flux_0d ( & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0,dtaudv_atm_0, & w_atm_0, u_star_0, b_star_0, q_star_0, & - cd_m_0, cd_t_0, cd_q_0 + cd_m_0, cd_t_0, cd_q_0, & + ex_del_m_0, ex_del_h_0, ex_del_q_0, & !mp586 for 10m winds and 2m temp + temp_2m_0, u_10m_0, v_10m_0, & !mp586 for 10m winds and 2m temp + q_2m_0, rh_2m_0 real, intent(inout) :: q_surf_0 real, intent(in) :: dt @@ -720,7 +783,11 @@ subroutine surface_flux_0d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, rho + cd_m, cd_t, cd_q, rho, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m !Add 2m q and RH + real, dimension(1) :: q_surf real, dimension(1) :: bucket_depth !RG Add bucket real, dimension(1) :: depth_change_lh_1d !RG Add bucket @@ -758,10 +825,13 @@ subroutine surface_flux_0d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, rho, & + cd_m, cd_t, cd_q, rho, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) flux_t_0 = flux_t(1) @@ -785,6 +855,14 @@ subroutine surface_flux_0d ( & cd_m_0 = cd_m(1) cd_t_0 = cd_t(1) cd_q_0 = cd_q(1) + ex_del_m_0 = ex_del_m(1) !mp586 for 10m winds and 2m temp + ex_del_h_0 = ex_del_h(1) !mp586 for 10m winds and 2m temp + ex_del_q_0 = ex_del_q(1) !mp586 for 10m winds and 2m temp + temp_2m_0 = temp_2m(1) !mp586 for 10m winds and 2m temp + u_10m_0 = u_10m(1) !mp586 for 10m winds and 2m temp + v_10m_0 = v_10m(1) !mp586 for 10m winds and 2m temp + q_2m_0 = q_2m(1) !Add 2m q + rh_2m_0 = rh_2m(1) !Add 2m RH end subroutine surface_flux_0d @@ -796,10 +874,13 @@ subroutine surface_flux_2d ( & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & flux_t, flux_q, flux_r, flux_u, flux_v, & - cd_m, cd_t, cd_q, rho, & + cd_m, cd_t, cd_q, rho, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m, & !Add 2m q and RH dt, land, seawater, avail ) ! ---- arguments ----------------------------------------------------------- @@ -814,7 +895,11 @@ subroutine surface_flux_2d ( & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & - cd_m, cd_t, cd_q, rho + cd_m, cd_t, cd_q, rho, & + ex_del_m, ex_del_h, ex_del_q, & !mp586 for 10m winds and 2m temp + temp_2m, u_10m, v_10m, & !mp586 for 10m winds and 2m temp + q_2m, rh_2m !Add 2m q and RH + real, intent(inout), dimension(:,:) :: q_surf logical, intent(in) :: bucket !RG Add bucket real, intent(inout), dimension(:,:) :: bucket_depth ! RG Add bucket @@ -835,10 +920,13 @@ subroutine surface_flux_2d ( & u_surf(:,j), v_surf(:,j), & rough_mom(:,j), rough_heat(:,j), rough_moist(:,j), rough_scale(:,j), gust(:,j), & flux_t(:,j), flux_q(:,j), flux_r(:,j), flux_u(:,j), flux_v(:,j), & - cd_m(:,j), cd_t(:,j), cd_q(:,j), rho(:,j), & + cd_m(:,j), cd_t(:,j), cd_q(:,j), rho(:,j), & w_atm(:,j), u_star(:,j), b_star(:,j), q_star(:,j), & dhdt_surf(:,j), dedt_surf(:,j), dedq_surf(:,j), drdt_surf(:,j), & dhdt_atm(:,j), dedq_atm(:,j), dtaudu_atm(:,j), dtaudv_atm(:,j), & + ex_del_m(:,j), ex_del_h(:,j), ex_del_q(:,j), & !mp586 for 10m winds and 2m temp + temp_2m(:,j), u_10m(:,j), v_10m(:,j), & !mp586 for 10m winds and 2m temp + q_2m(:,j), rh_2m(:,j), & dt, land(:,j), seawater(:,j), avail(:,j) ) end do end subroutine surface_flux_2d From c5b23bbdf1653e2ce1aaad52256b70a993e8317d Mon Sep 17 00:00:00 2001 From: Geen Date: Thu, 14 May 2020 10:46:34 +0100 Subject: [PATCH 3/3] Add surface flux diag edits to up-to-date idealised_moist_phys.F90 --- .../driver/solo/idealized_moist_phys.F90 | 31 +++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 8115c08a2..306a023ff 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -179,6 +179,7 @@ module idealized_moist_phys_mod drag_m, & ! momentum drag coefficient drag_t, & ! heat drag coefficient drag_q, & ! moisture drag coefficient + rho, & ! air density at surface w_atm, & ! wind speed ustar, & ! friction velocity bstar, & ! buoyancy scale @@ -264,6 +265,11 @@ module idealized_moist_phys_mod id_bucket_depth_conv, & ! bucket depth variation induced by convection - RG Add bucket id_bucket_depth_cond, & ! bucket depth variation induced by condensation - RG Add bucket id_bucket_depth_lh, & ! bucket depth variation induced by LH - RG Add bucket + id_w_atm, & ! wind speed - RG Add lh flux breakdown + id_drag_q, & ! moisture drag coefficient - RG Add lh flux breakdown + id_rho, & ! density at surface - RG Add lh flux breakdown + id_q_atm, & ! lowest level specific humidity - RG Add lh flux breakdown + id_q_surf, & ! surface humidity - RG Add lh flux breakdown id_rh, & ! Relative humidity id_diss_heat_ray,& ! Heat dissipated by rayleigh bottom drag if gp_surface=.True. id_z_tg, & ! Relative humidity @@ -445,6 +451,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(drag_m (is:ie, js:je)) allocate(drag_t (is:ie, js:je)) allocate(drag_q (is:ie, js:je)) +allocate(rho (is:ie, js:je)) allocate(w_atm (is:ie, js:je)) allocate(ustar (is:ie, js:je)) allocate(bstar (is:ie, js:je)) @@ -633,7 +640,20 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Zonal momentum flux', 'Pa') id_flux_v = register_diag_field(mod_name, 'flux_v', & axes(1:2), Time, 'Meridional momentum flux', 'Pa') - + +if(.not.gp_surface) then + id_w_atm = register_diag_field(mod_name, 'wind_speed', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Lowest level wind speed','m/s') + id_drag_q = register_diag_field(mod_name, 'drag_q', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Moisture drag coefficient','none') + id_rho = register_diag_field(mod_name, 'rho', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Air density at lowest level','kg/m/m/m') + id_q_atm = register_diag_field(mod_name, 'q_atm', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Lowest level specific humidity','kg/kg') + id_q_surf = register_diag_field(mod_name, 'q_surf', & ! RG Add lh flux breakdown + axes(1:2), Time, 'Surface specific humidity','kg/kg') +endif + if(bucket) then id_bucket_depth = register_diag_field(mod_name, 'bucket_depth', & ! RG Add bucket axes(1:2), Time, 'Depth of surface reservoir', 'm') @@ -1012,6 +1032,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg drag_m(:,:), & ! is intent(out) drag_t(:,:), & ! is intent(out) drag_q(:,:), & ! is intent(out) + rho(:,:), & ! is intent(out) w_atm(:,:), & ! is intent(out) ustar(:,:), & ! is intent(out) bstar(:,:), & ! is intent(out) @@ -1036,7 +1057,13 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg land(:,:), & .not.land(:,:), & avail(:,:) ) - + + if(id_w_atm > 0) used = send_data(id_w_atm, w_atm, Time) ! RG Add lh flux breakdown + if(id_drag_q > 0) used = send_data(id_drag_q, drag_q, Time) ! RG Add lh flux breakdown + if(id_rho > 0) used = send_data(id_rho, rho, Time) ! RG Add lh flux breakdown + if(id_q_atm > 0) used = send_data(id_q_atm, grid_tracers(:,:,num_levels,previous,nsphum), Time) ! RG Add lh flux breakdown + if(id_q_surf > 0) used = send_data(id_q_surf, q_surf, Time) ! RG Add lh flux breakdown + if(id_flux_u > 0) used = send_data(id_flux_u, flux_u, Time) if(id_flux_v > 0) used = send_data(id_flux_v, flux_v, Time)