#include module lake_tile_mod use mpp_domains_mod, only : & domain2d, mpp_get_compute_domain, mpp_global_field #ifdef INTERNAL_FILE_NML use mpp_mod, only: input_nml_file #else use fms_mod, only: open_namelist_file #endif use fms_mod, only : & write_version_number, file_exist, check_nml_error, & read_data, close_file, stdlog use constants_mod, only : & pi, tfreeze, hlf use land_constants_mod, only : & NBANDS use land_io_mod, only : & init_cover_field use land_tile_selectors_mod, only : & tile_selector_type, SEL_LAKE, register_tile_selector implicit none private ! ==== public interfaces ===================================================== public :: lake_pars_type public :: lake_prog_type public :: lake_tile_type public :: new_lake_tile, delete_lake_tile public :: lake_tiles_can_be_merged, merge_lake_tiles public :: lake_is_selected public :: get_lake_tile_tag public :: lake_tile_stock_pe public :: lake_tile_heat public :: read_lake_data_namelist public :: lake_cover_cold_start public :: lake_data_radiation public :: lake_data_diffusion public :: lake_data_thermodynamics public :: max_lev public :: lake_width_inside_lake public :: large_lake_sill_width ! =====end of public interfaces ============================================== interface new_lake_tile module procedure lake_tile_ctor module procedure lake_tile_copy_ctor end interface ! ==== module constants ====================================================== character(len=*), private, parameter :: & version = '$Id: lake_tile.F90,v 20.0 2013/12/13 23:29:39 fms Exp $', & tagname = '$Name: tikal $', & module_name = 'lake_tile_mod' integer, parameter :: max_lev = 80 integer, parameter :: n_dim_lake_types = 1 ! size of lookup table real, parameter :: psi_wilt = -150. ! matric head at wilting real, parameter :: comp = 0.001 ! m^-1 ! from the modis brdf/albedo product user's guide: real :: g_iso = 1. real :: g_vol = 0.189184 real :: g_geo = -1.377622 real :: g0_iso = 1.0 real :: g1_iso = 0.0 real :: g2_iso = 0.0 real :: g0_vol = -0.007574 real :: g1_vol = -0.070987 real :: g2_vol = 0.307588 real :: g0_geo = -1.284909 real :: g1_geo = -0.166314 real :: g2_geo = 0.041840 ! ==== types ================================================================= type :: lake_pars_type real w_sat real awc_lm2 real k_sat_ref real psi_sat_ref real chb real alpha ! *** REPLACE LATER BY alpha(layer) real heat_capacity_ref real thermal_cond_ref real refl_dry_dir(NBANDS) real refl_dry_dif(NBANDS) real refl_sat_dir(NBANDS) real refl_sat_dif(NBANDS) real emis_dry real emis_sat real z0_momentum real z0_momentum_ice real depth_sill real width_sill real whole_area real connected_to_next real backwater real backwater_1 real tau_groundwater real rsa_exp ! riparian source-area exponent end type lake_pars_type type :: lake_prog_type real dz real wl real ws real T real K_z real groundwater real groundwater_T end type lake_prog_type type :: lake_tile_type integer :: tag ! kind of the lake type(lake_prog_type), pointer :: prog(:) type(lake_pars_type) :: pars real, pointer :: w_fc(:) real, pointer :: w_wilt(:) real :: Eg_part_ref real :: z0_scalar real :: z0_scalar_ice real :: geothermal_heat_flux real, pointer :: e(:),f(:) real, pointer :: heat_capacity_dry(:) end type lake_tile_type ! ==== module data =========================================================== real, public :: & cpw = 1952.0, & ! specific heat of water vapor at constant pressure clw = 4218.0, & ! specific heat of water (liquid) csw = 2106.0 ! specific heat of water (ice) !---- namelist --------------------------------------------------------------- real :: lake_width_inside_lake = 1.e5 real :: large_lake_sill_width = 200. real :: min_lake_frac = 0. real :: max_lake_rh = 1. real :: lake_rh_exp = 1. real :: dry_lake_depth_frac = 0. real :: k_over_B = 0.25 ! reset to 0 for MCM real :: k_over_B_ice = 0.25 real :: rate_fc = 0.1/86400 ! 0.1 mm/d drainage rate at FC real :: sfc_heat_factor = 1 real :: geothermal_heat_flux_constant = 0.0 ! true continental average is ~0.065 W/m2 integer, public :: num_l = 18 ! number of lake levels integer, public :: n_outlet = 0 logical :: use_lm2_awc = .false. logical, public :: lake_specific_width = .false. integer :: n_map_1st_lake_type = 10 integer, public :: outlet_face(100) integer, public :: outlet_i(100) integer, public :: outlet_j(100) real, public :: outlet_width(100) ! from analysis of modis data (ignoring temperature dependence): real :: f_iso_ice(NBANDS) = (/ 0.056, 0.131 /) real :: f_vol_ice(NBANDS) = (/ 0.017, 0.053 /) real :: f_geo_ice(NBANDS) = (/ 0.004, 0.010 /) real :: f_iso_liq(NBANDS) = (/ 0.056, 0.131 /) real :: f_vol_liq(NBANDS) = (/ 0.017, 0.053 /) real :: f_geo_liq(NBANDS) = (/ 0.004, 0.010 /) real :: refl_ice_dif(NBANDS), refl_liq_dif(NBANDS) ! ---- remainder are used only for cold start --------- logical :: round_frac_down = .false. ! when false, any lake_frac < min_lake_frac ! is set to min_lake_frac. ! when true, any lake_frac < min_lake_frac ! is set to 0. character(len=16):: lake_to_use = 'single-tile' ! 'multi-tile' for tiled soil [default] ! 'single-tile' for geographically varying soil with single type per ! model grid cell ! 'uniform' for global constant soil, e.g., to reproduce MCM ! 'from-rivers' to get the fraction of lakes from river module logical :: use_single_lake = .false. ! true for single global lake, ! e.g., to recover MCM logical :: use_mcm_albedo = .false. ! .true. for CLIMAP albedo inputs logical :: use_single_geo = .false. ! .true. for global gw res time, ! e.g., to recover MCM integer :: lake_index_constant = 1 ! index of global constant lake, ! used when use_single_lake real :: gw_res_time = 60.*86400 ! mean groundwater residence time, ! used when use_single_geo real :: rsa_exp_global = 1.5 real, dimension(n_dim_lake_types) :: & dat_w_sat =(/ 1.000 /),& dat_awc_lm2 =(/ 1.000 /),& dat_k_sat_ref =(/ 0.021 /),& dat_psi_sat_ref =(/ -.059 /),& dat_chb =(/ 3.5 /),& dat_heat_capacity_ref =(/ 8.4e7 /),& dat_thermal_cond_ref =(/ 8.4e7 /),& dat_emis_dry =(/ 0.950 /),& dat_emis_sat =(/ 0.980 /),& dat_z0_momentum =(/ 1.4e-4 /),& dat_z0_momentum_ice =(/ 1.4e-4 /),& dat_tf_depr =(/ 0.00 /) real, dimension(n_dim_lake_types, NBANDS) :: & dat_refl_dry_dif, dat_refl_dry_dir, & dat_refl_sat_dif, dat_refl_sat_dir data & ! VIS NIR dat_refl_dry_dif / 0.060, 0.060 /, & dat_refl_dry_dir / 0.060, 0.060 /, & dat_refl_sat_dir / 0.060, 0.060 /, & dat_refl_sat_dif / 0.060, 0.060 / integer, dimension(n_dim_lake_types) :: & input_cover_types =(/ 10 /) character(len=4), dimension(n_dim_lake_types) :: & tile_names =(/ 'lake' /) namelist /lake_data_nml/ lake_width_inside_lake, & large_lake_sill_width, & lake_specific_width, n_outlet, outlet_face, outlet_i, outlet_j, outlet_width, & min_lake_frac, round_frac_down, max_lake_rh, lake_rh_exp, & dry_lake_depth_frac, & lake_to_use,input_cover_types, tile_names, & k_over_B, k_over_B_ice, & rate_fc, sfc_heat_factor, geothermal_heat_flux_constant, & num_l, & use_lm2_awc, n_map_1st_lake_type, & use_single_lake, use_mcm_albedo, & use_single_geo, lake_index_constant, & gw_res_time, rsa_exp_global, & dat_w_sat, dat_awc_lm2, & dat_k_sat_ref, & dat_psi_sat_ref, dat_chb, & dat_heat_capacity_ref, dat_thermal_cond_ref, & dat_refl_dry_dir, dat_refl_sat_dir, & dat_refl_dry_dif, dat_refl_sat_dif, & dat_emis_dry, dat_emis_sat, & dat_z0_momentum, dat_z0_momentum_ice, dat_tf_depr, & f_iso_ice, f_vol_ice, f_geo_ice, f_iso_liq, f_vol_liq, f_geo_liq !---- end of namelist -------------------------------------------------------- contains ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ! ============================================================================ subroutine read_lake_data_namelist(lake_n_lev) integer, intent(out) :: lake_n_lev ! ---- local vars integer :: unit ! unit for namelist i/o integer :: io ! i/o status for the namelist integer :: ierr ! error code, returned by i/o routines integer :: i real :: z call write_version_number(version, tagname) #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=lake_data_nml, iostat=io) ierr = check_nml_error(io, 'lake_data_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file() ierr = 1; do while (ierr /= 0) read (unit, nml=lake_data_nml, iostat=io, end=10) ierr = check_nml_error (io, 'lake_data_nml') enddo 10 continue call close_file (unit) endif #endif unit=stdlog() write(unit, nml=lake_data_nml) ! initialize global module data here refl_ice_dif = g_iso*f_iso_ice + g_vol*f_vol_ice + g_geo*f_geo_ice refl_liq_dif = g_iso*f_iso_liq + g_vol*f_vol_liq + g_geo*f_geo_liq ! register selectors for tile-specific diagnostics do i=1, n_dim_lake_types call register_tile_selector(tile_names(i), long_name='',& tag = SEL_LAKE, idata1 = i ) enddo ! set up output arguments lake_n_lev = num_l end subroutine ! ============================================================================ function lake_tile_ctor(tag) result(ptr) type(lake_tile_type), pointer :: ptr ! return value integer, intent(in) :: tag ! kind of lake allocate(ptr) ptr%tag = tag ! allocate storage for tile data allocate(ptr%prog (num_l), & ptr%w_fc (num_l), & ptr%w_wilt (num_l), & ptr%heat_capacity_dry(num_l), & ptr%e (num_l), & ptr%f (num_l) ) call init_lake_data_0d(ptr) end function lake_tile_ctor ! ============================================================================ function lake_tile_copy_ctor(lake) result(ptr) type(lake_tile_type), pointer :: ptr ! return value type(lake_tile_type), intent(in) :: lake ! tile to copy allocate(ptr) ptr=lake ! copy all non-pointer data ! allocate storage for tile data allocate(ptr%prog (num_l), & ptr%w_fc (num_l), & ptr%w_wilt (num_l), & ptr%heat_capacity_dry(num_l), & ptr%e (num_l), & ptr%f (num_l) ) ! copy all allocatable data ptr%prog(:) = lake%prog(:) ptr%w_fc(:) = lake%w_fc(:) ptr%w_wilt(:) = lake%w_wilt(:) ptr%e(:) = lake%e(:) ptr%f(:) = lake%f(:) ptr%heat_capacity_dry(:) = lake%heat_capacity_dry(:) end function lake_tile_copy_ctor ! ============================================================================ subroutine delete_lake_tile(ptr) type(lake_tile_type), pointer :: ptr deallocate(ptr%prog, ptr%w_fc, ptr%w_wilt,ptr%heat_capacity_dry, ptr%e, ptr%f) deallocate(ptr) end subroutine delete_lake_tile subroutine init_lake_data_0d(lake) type(lake_tile_type), intent(inout) :: lake ! real tau_groundwater ! real rsa_exp ! riparian source-area exponent integer :: k k = lake%tag lake%pars%w_sat = dat_w_sat (k) lake%pars%awc_lm2 = dat_awc_lm2 (k) lake%pars%k_sat_ref = dat_k_sat_ref (k) lake%pars%psi_sat_ref = dat_psi_sat_ref (k) lake%pars%chb = dat_chb (k) lake%pars%alpha = 1 lake%pars%heat_capacity_ref = dat_heat_capacity_ref(k) lake%pars%thermal_cond_ref = dat_thermal_cond_ref (k) lake%pars%refl_dry_dir = dat_refl_dry_dir (k,:) lake%pars%refl_dry_dif = dat_refl_dry_dif (k,:) lake%pars%refl_sat_dir = dat_refl_sat_dir (k,:) lake%pars%refl_sat_dif = dat_refl_sat_dif (k,:) lake%pars%emis_dry = dat_emis_dry (k) lake%pars%emis_sat = dat_emis_sat (k) lake%pars%z0_momentum = dat_z0_momentum (k) lake%pars%z0_momentum_ice = dat_z0_momentum_ice (k) lake%pars%tau_groundwater = 86400.*30. lake%pars%rsa_exp = rsa_exp_global ! -------- derived constant lake parameters -------- ! w_fc (field capacity) set to w at which hydraulic conductivity equals ! a nominal drainage rate "rate_fc" ! w_wilt set to w at which psi is psi_wilt if (use_lm2_awc) then lake%w_wilt(:) = 0.15 lake%w_fc (:) = 0.15 + lake%pars%awc_lm2 else lake%w_wilt(:) = lake%pars%w_sat & *(lake%pars%psi_sat_ref/(psi_wilt*lake%pars%alpha))**(1/lake%pars%chb) lake%w_fc (:) = lake%pars%w_sat & *(rate_fc/(lake%pars%k_sat_ref*lake%pars%alpha**2))**(1/(3+2*lake%pars%chb)) endif ! below made use of phi_e from parlange via entekhabi lake%Eg_part_ref = (-4*lake%w_fc(1)**2*lake%pars%k_sat_ref*lake%pars%psi_sat_ref*lake%pars%chb & /(pi*lake%pars%w_sat)) * (lake%w_fc(1)/lake%pars%w_sat)**(2+lake%pars%chb) & *(2*pi/(3*lake%pars%chb**2*(1+3/lake%pars%chb)*(1+4/lake%pars%chb)))/2 lake%z0_scalar = lake%pars%z0_momentum * exp(-k_over_B) lake%z0_scalar_ice = lake%pars%z0_momentum_ice * exp(-k_over_B_ice) lake%geothermal_heat_flux = geothermal_heat_flux_constant end subroutine ! ============================================================================ function lake_cover_cold_start(land_mask, lonb, latb, domain) result (lake_frac) ! creates and initializes a field of fractional lake coverage logical, intent(in) :: land_mask(:,:) ! land mask real, intent(in) :: lonb(:,:), latb(:,:)! boundaries of the grid cells real, pointer :: lake_frac (:,:,:) ! output: map of lake fractional coverage type(domain2d), intent(in) :: domain allocate( lake_frac(size(land_mask,1),size(land_mask,2),n_dim_lake_types)) if (trim(lake_to_use)=='from-rivers') then lake_frac = 0.0 call read_data('INPUT/river_data.nc', 'lake_frac', lake_frac(:,:,1), & domain=domain) ! make sure 'missing values' don't get into the result where (lake_frac < 0) lake_frac = 0 where (lake_frac > 1) lake_frac = 1 else call init_cover_field(lake_to_use, 'INPUT/ground_type.nc', 'cover','frac', & lonb, latb, lake_index_constant, input_cover_types, lake_frac) endif if (round_frac_down) then where (lake_frac.gt.0. .and. lake_frac.lt.min_lake_frac) lake_frac = 0. else where (lake_frac.gt.0. .and. lake_frac.lt.min_lake_frac) lake_frac = min_lake_frac endif end function ! ============================================================================= function lake_tiles_can_be_merged(lake1,lake2) result(response) logical :: response type(lake_tile_type), intent(in) :: lake1,lake2 response = (lake1%tag==lake2%tag) end function ! ============================================================================= ! combine two lake tiles with specified weights; the results goes into the ! second one ! THIS NEEDS TO BE REVISED FOR TILE-DEPENDENT DZ subroutine merge_lake_tiles(t1,w1,t2,w2) type(lake_tile_type), intent(in) :: t1 type(lake_tile_type), intent(inout) :: t2 real, intent(in) :: w1, w2 ! relative weights ! ---- local vars real :: x1, x2 ! normalized relative weights real :: HEAT1, HEAT2 ! temporaries for heat real :: C1, C2 ! heat capacities real :: gw integer :: i WRITE (*,*) 'SORRY, BUT merge_lake_tiles NEEDS TO BE REVISED TO ALLOW FOR ', & 'HORIZONTALLY VARYING VERTICAL DISCRETIZATION' ! calculate normalized weights x1 = w1/(w1+w2) x2 = 1.0 - x1 ! combine state variables do i = 1,num_l ! calculate "dry" heat capacities: C1 = sfc_heat_factor*t1%pars%heat_capacity_ref C2 = sfc_heat_factor*t2%pars%heat_capacity_ref ! calculate heat content at this level for both source tiles HEAT1 = & 0. ! (C1*dz(i)+clw*t1%prog(i)%wl+csw*t1%prog(i)%ws) * (t1%prog(i)%T-tfreeze) HEAT2 = & 0. ! (C2*dz(i)+clw*t2%prog(i)%wl+csw*t2%prog(i)%ws) * (t2%prog(i)%T-tfreeze) ! calculate (and assign) combined water mass t2%prog(i)%wl = t1%prog(i)%wl*x1 + t2%prog(i)%wl*x2 t2%prog(i)%ws = t1%prog(i)%ws*x1 + t2%prog(i)%ws*x2 ! if dry heat capacity of combined lake is to be changed, update it here ! ... ! calculate combined temperature, based on total heat content and combined ! heat capacity ! t2%prog(i)%T = (HEAT1*x1+HEAT2*x2) / & ! (C2*dz(i)+clw*t2%prog(i)%wl+csw*t2%prog(i)%ws) + tfreeze t2%prog(i)%T = 0. ! calculate combined groundwater content gw = t1%prog(i)%groundwater*x1 + t2%prog(i)%groundwater*x2 ! calculate combined groundwater temperature if(gw/=0) then t2%prog(i)%groundwater_T = ( & t1%prog(i)%groundwater*x1*(t1%prog(i)%groundwater_T-tfreeze) + & t2%prog(i)%groundwater*x2*(t2%prog(i)%groundwater_T-tfreeze) & ) / gw + Tfreeze else t2%prog(i)%groundwater_T = & x1*t1%prog(i)%groundwater_T + x2*t2%prog(i)%groundwater_T endif t2%prog(i)%groundwater = gw enddo end subroutine ! ============================================================================= ! returns true if tile fits the specified selector function lake_is_selected(lake, sel) logical lake_is_selected type(tile_selector_type), intent(in) :: sel type(lake_tile_type), intent(in) :: lake lake_is_selected = (sel%idata1 == lake%tag) end function ! ============================================================================ ! returns tag of the tile function get_lake_tile_tag(lake) result(tag) integer :: tag type(lake_tile_type), intent(in) :: lake tag = lake%tag end function ! ============================================================================ ! compute bare-lake albedos and bare-lake emissivity subroutine lake_data_radiation ( lake, cosz, use_brdf, & lake_alb_dir, lake_alb_dif, lake_emis ) type(lake_tile_type), intent(in) :: lake real, intent(in) :: cosz logical, intent(in) :: use_brdf real, intent(out) :: lake_alb_dir(:), lake_alb_dif(:), lake_emis ! ---- local vars real :: lake_sfc_vlc, blend real :: liq_value_dir(NBANDS), ice_value_dir(NBANDS) real :: liq_value_dif(NBANDS), ice_value_dif(NBANDS) real :: zenith_angle, zsq, zcu ! ---- radiation properties lake_sfc_vlc = lake%prog(1)%wl/lake%prog(1)%dz blend = lake_sfc_vlc/lake%pars%w_sat if (use_brdf) then zenith_angle = acos(cosz) zsq = zenith_angle*zenith_angle zcu = zenith_angle*zsq liq_value_dir = f_iso_liq*(g0_iso+g1_iso*zsq+g2_iso*zcu) & + f_vol_liq*(g0_vol+g1_vol*zsq+g2_vol*zcu) & + f_geo_liq*(g0_geo+g1_geo*zsq+g2_geo*zcu) ice_value_dir = f_iso_ice*(g0_iso+g1_iso*zsq+g2_iso*zcu) & + f_vol_ice*(g0_vol+g1_vol*zsq+g2_vol*zcu) & + f_geo_ice*(g0_geo+g1_geo*zsq+g2_geo*zcu) liq_value_dif = refl_liq_dif ice_value_dif = refl_ice_dif else liq_value_dir = lake%pars%refl_sat_dir ice_value_dir = lake%pars%refl_dry_dir liq_value_dif = lake%pars%refl_sat_dif ice_value_dif = lake%pars%refl_dry_dif endif lake_alb_dir = ice_value_dir + blend*(liq_value_dir-ice_value_dir) lake_alb_dif = ice_value_dif + blend*(liq_value_dif-ice_value_dif) lake_emis = lake%pars%emis_dry + blend*(lake%pars%emis_sat-lake%pars%emis_dry ) end subroutine ! ============================================================================ ! compute bare-lake roughness subroutine lake_data_diffusion ( lake,lake_z0s, lake_z0m ) type(lake_tile_type), intent(in) :: lake real, intent(out) :: lake_z0s, lake_z0m ! ---- surface roughness if (lake%prog(1)%ws.le.0.) then lake_z0s = lake%z0_scalar lake_z0m = lake%pars%z0_momentum else lake_z0s = lake%z0_scalar_ice lake_z0m = lake%pars%z0_momentum_ice endif end subroutine ! ============================================================================ ! compute lake thermodynamic properties. subroutine lake_data_thermodynamics ( lake_pars, lake_depth, & lake_rh, heat_capacity_dry, thermal_cond) type(lake_pars_type), intent(in) :: lake_pars real, intent(in) :: lake_depth real, intent(out) :: lake_rh real, intent(out) :: heat_capacity_dry(:) real, intent(out) :: thermal_cond(:) ! ---- local vars integer l real lake_depth_frac, lake_rh_base ! ---------------------------------------------------------------------------- lake_depth_frac = (lake_depth/lake_pars%depth_sill) lake_rh_base = (lake_depth_frac-dry_lake_depth_frac)/(1.-dry_lake_depth_frac) if (lake_rh_base.gt.0.) then lake_rh = min(max_lake_rh, lake_rh_base**lake_rh_exp ) else lake_rh = 0. endif do l = 1, num_l heat_capacity_dry(l) = lake_pars%heat_capacity_ref thermal_cond(l) = lake_pars%thermal_cond_ref enddo end subroutine ! ============================================================================ subroutine lake_tile_stock_pe (lake, twd_liq, twd_sol ) type(lake_tile_type), intent(in) :: lake real, intent(out) :: twd_liq, twd_sol integer n twd_liq = 0. twd_sol = 0. do n=1, size(lake%prog) twd_liq = twd_liq + lake%prog(n)%wl + lake%prog(n)%groundwater twd_sol = twd_sol + lake%prog(n)%ws enddo end subroutine lake_tile_stock_pe ! ============================================================================ ! returns lake tile heat content, J/m2 function lake_tile_heat (lake) result(heat) ; real heat type(lake_tile_type), intent(in) :: lake integer :: i heat = 0 do i = 1, num_l heat = heat + & (lake%heat_capacity_dry(i)*lake%prog(i)%dz + clw*lake%prog(i)%wl & + csw*lake%prog(i)%ws)*(lake%prog(i)%T-tfreeze) + & clw*lake%prog(i)%groundwater*(lake%prog(i)%groundwater_T-tfreeze) - & hlf*lake%prog(i)%ws enddo end function end module lake_tile_mod