!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module land_properties_mod ! ! Christopher Milly ! ! ! Elena Shevliakova ! ! ! Sergey Malyshev ! ! ! Contains land properties namelist variables and procedures relating ! to the properties of the land. ! ! ! Initialization and calculation of the land property data. The input ! cover field and implied glacier field are obtained and the glacier mask ! and cover type are redefined, if necessary. The properties that depend ! on cover and ground types are assigned. The land properties diagnostics ! are initialized and the static fields are sent to the diagnostic manager ! for output. ! ! Updates the rapidly changing parameters. Computes the albedo of ! the hypothetical no-snow and deep-snow surfaces and uses snow mass to ! blend snow-free and deep-snow albedo values. Regrids integer index data ! to any output grid from a uniformly spaced grid using a maximum-area ! voting scheme. Includes calculation of total area of each land cover ! type within the territory and the determination of the type occupying ! the biggest area within the territory. ! use mpp_domains_mod, only: domain2d, mpp_get_compute_domain, mpp_get_global_domain use time_manager_mod, only: time_type, set_date, print_date, get_date, & operator(-), operator(>) use mpp_io_mod, only: mpp_open, mpp_RDONLY use fms_mod, only: error_mesg, file_exist, open_namelist_file, & check_nml_error, stdlog, write_version_number, & mpp_pe, mpp_root_pe, close_file, read_data, & FATAL, NOTE, open_restart_file, set_domain, & field_size, stdout use constants_mod, only: tfreeze, pi use horiz_interp_mod, only: horiz_interp_type, horiz_interp_new, horiz_interp, & horiz_interp_del use numerics_mod, only: is_latlon, expand_cell use climap_albedo_mod, only: get_climap_albedo, get_climap_glacier, & get_climap_albedo_mcm, get_climap_glacier_mcm use diag_manager_mod, only : register_static_field, send_data, get_base_time use topography_mod, only: get_topog_stdev implicit none private ! ==== public interface ===================================================== public :: land_properties_init,land_properties_end public :: update_land_properties_slow public :: update_land_properties_fast public :: regrid_discrete_field ! ==== end public interface ================================================= ! ! n_dim_ground_types size of ground (soil) parameter lookup tables ! n_dim_cover_types size of cover (vegetation) parameter lookup tables ! n_map_ground_types number of ground types in input map file ! n_map_cover_types number of cover types in input map file ! ! Size of ground (soil) parameter lookup tables ! ! ! Size of cover (vegetation) parameter lookup tables ! ! ! Number of ground types in input map file ! ! ! Number of cover types in input map file ! integer, parameter :: n_dim_ground_types = 11 ! = 10+ssl integer, parameter :: n_dim_cover_types = 14 ! = 11+tun+ssl_veg+ssl_ice integer, parameter :: n_map_ground_types = 10 integer, parameter :: n_map_cover_types = 11 !---- namelist - with default values ---------------------------------- ! ! do_all_mcm run Manabe Climate Model land surface [logical] ! Setting this global control option to TRUE causes ! specification all the following (regardless of default or input ! settings, which are then ignored): ! veg_to_use = 'cons_ssl', ! soil_to_use = 'cons_ssl', ! use_climap = .true. (this forced by veg_to_use='cons_ssl'), ! use_climap_mcm = .true., ! do_mcm_masking = .true., ! use_single_geo = .true., ! geo_res_time = res_time_ssl, ! factor_root = 1., ! factor_rough = 1., ! z_root_min = 0., ! max_snow = max_snow_ssl ! veg_to_use choice of method for defining vegetation [character] ! 'variable' - use input map to index vegetation parameter vectors ! 'constant' - use veg_index_constant to index veg parameter vectors ! 'cons_tun' - use tuned global constant vegetation ! 'cons_ssl' - use Manabe Climate Model-like vegetation (forces use_climap to be true) ! soil_to_use choice of method for defining soil [character] ! 'variable' - use input map to index soil parameter vectors ! 'constant' - use soil_index_constant to index soil parameter vectors ! 'cons_ssl' - use Manabe Climate Model-like soil ! use_glaciers false to remove glaciers from land [logical] ! use_climap true to override default albedo by CLIMAP [logical] ! and to invoke use of CLIMAP to define glacier ! locations (if use_glaciers) ! use_desert_albedo_map true to override default snow-free albedo of desert only by ! albedo map (SRB) ! do_mcm_masking true to use Manabe Climate Model snow-albedo fn. [logical] ! use_single_basin true to avoid using basin maps [logical] ! use_single_geo true for global constant gw residence time [logical] ! i_dest0 if use_single_geo is set to .true., all the river ! discharge is put in a single grid cell. i_dest0 is the ! longitude index of this grid cell. ! j_dest0 if use_single_geo is set to .true., all the river ! discharge is put in a single grid cell. j_dest0 is the ! latitude index of this grid cell. ! veg_index_constant veg index used when veg_to_use='constant' - ! soil_index_constant soil index used when soil_to_use='constant' - ! soil_index_ice_substitute ground type to be substituted for ! ice ground type when such type is over-ruled due to ! absence of glacier cover type. if this is set to ! zero, then the ground type is temporarily marked as ! ocean, and, if land is present, ground type will then ! be assigned based on neighbor cells through regrid_discrete_field. ! geo_res_time time constant when use_single_geo is true s ! t_range T range over which snow/glacier albedo varies K ! factor_root global factor for critical_root_density - ! factor_stomata global factor for veg_rs_min_vec - ! factor_rough global factor for rough_momentum_vec - ! z_root_min lower bound for root-zone depth m ! max_snow value of snow above which 'snow runoff' occurs kg/m**3 ! sfc_heat_factor "fudge" factor for heat capacity and thermal ! conductivity in surface soil layers - ! num_sfc_layers number of surface layers for sfc_heat_factor - ! dynamic_cover_type set to true if cover type forcing data varies with ! time [logical] ! read_old_ascii_cover read the original ASCII file of cover type ! [logical] ! cover_dataset_init_year initial year in the cover_type dataset [integer] ! cover_dataset_entry beginning time for reading the cover_type ! dataset (yr, mo, dy, hr, mn, sc) [integer] ! ! Run Manabe Climate Model land surface. Setting this global control ! option to TRUE causes specification all the following (regardless of ! default or input settings, which are then ignored): ! veg_to_use = 'cons_ssl' ! soil_to_use = 'cons_ssl' ! use_climap = .true. (this forced by veg_to_use='cons_ssl') ! use_climap_mcm = .true. ! do_mcm_masking = .true. ! use_single_geo = .true. ! geo_res_time = res_time_ssl ! factor_root = 1. ! factor_rough = 1. ! z_root_min = 0. ! max_snow = max_snow_ssl ! ! ! Choice of method for defining vegetation. ! 'variable' - use input map to index vegetation parameter vectors ! 'constant' - use veg_index_constant to index veg parameter vectors ! 'cons_tun' - use tuned global constant vegetation ! 'cons_ssl' - use Manabe Climate Model-like vegetation (forces use_climap to be true) ! ! ! Choice of method for defining soil. ! 'variable' - use input map to index soil parameter vectors ! 'constant' - use soil_index_constant to index soil parameter vectors ! 'cons_ssl' - use Manabe Climate Model-like soil ! ! ! False to remove glaciers from land ! ! ! True to override default albedo by CLIMAP and to invoke use of CLIMAP ! to define glacier locations (if use_glaciers) ! ! ! true to override default snow-free albedo of desert only by ! albedo map (SRB) ! ! ! Run Manabe Climate Model land surface ! ! ! True to use Manabe Climate Model snow-albedo function ! ! ! True to avoid using basin maps ! ! ! True for global constant groundwater residence time ! ! ! If use_single_geo is set to .true., all the river discharge is put in a ! single grid cell. i_dest0 is the longitude index of this grid cell. ! ! ! If use_single_geo is set to .true., all the river discharge is put in a ! single grid cell. j_dest0 is the latitude index of this grid cell. ! ! ! Veg index used when veg_to_use='constant' ! ! ! Soil index used when soil_to_use='constant' ! ! ! Ground type to be substituted for ice ground type when such type is ! over-ruled due to absence of glacier cover type. if this is set to zero, ! then the ground type is temporarily marked as ocean, and, if land is ! present, ground type will then be assigned based on neighbor cells ! through regrid_discrete_field. ! ! ! Time constant when use_single_geo is true ! ! ! Temperature range over which snow/glacier albedo varies ! ! ! Global factor for critical_root_density ! ! ! Global factor for veg_rs_min_vec ! ! ! Global factor for rough_momentum_vec ! ! ! lower bound for root-zone depth ! ! ! Value of snow above which 'snow runoff' occurs ! ! ! "fudge" factor for heat capacity and thermal ! conductivity in surface soil layers ! ! ! number of surface layers for sfc_heat_factor ! ! ! Set to true if cover type forcing data varies with time. ! ! ! Set to true if reading the original ASCII static cover type forcing data ! ! ! The initial year in the cover_type dataset. ! ! ! Beginning time for reading the cover_type dataset (yr,mo,dy,hr,mn,sc). ! logical :: do_all_mcm = .false. character*8 :: veg_to_use = 'variable' character*8 :: soil_to_use = 'variable' logical :: use_glaciers = .true. logical :: use_climap = .false. logical :: use_desert_albedo_map = .false. logical :: use_climap_mcm = .false. logical :: do_mcm_masking = .false. logical :: use_single_basin = .false. logical :: use_single_geo = .false. integer :: i_dest0 = 1 integer :: j_dest0 = 1 integer :: veg_index_constant = 3 ! mixed broadleaf/needleleaf integer :: soil_index_constant = 2 ! medium texture mineral soil integer :: soil_index_ice_substitute = 1 logical :: reconcile_lakes = .false. ! if true, lake distribution ! in cover and ground type fields is reconciled based on the cover type integer :: soil_index_lake_substitute = 1 ! substitute for the ground lake type, ! used only if reconcile_lakes is true real :: geo_res_time = 60.*86400. ! two months real :: t_range = 10.0 real :: factor_root = 1.0 real :: factor_stomata = 1.0 real :: factor_rough = 1.0 real :: z_root_min = 0.01 real :: max_snow = 1000. ! (1 m water equiv) real :: sfc_heat_factor = 1. integer :: num_sfc_layers = 0 logical :: dynamic_cover_type = .true. ! true if cover type forcing data ! varies with time. Must also set ! cover_dataset_entry. logical :: read_old_ascii_cover = .false. ! true if reading original ! static ACII cover type forcing. integer, dimension(6) :: & cover_dataset_entry = (/ 1860, 1, 1, 0, 0, 0 /) ! beginning time for reading the cover_type data. Must be set when ! using read_old_ascii_cover = false. integer :: cover_dataset_init_year = 1860 ! initial year for the cover_type ! dataset. ! ! Minimum bulk stomatal resistance. For default values, refer to the table ! in the Public Code section below. ! ! ! Depth scale of root distribution. For default values, refer to the table ! in the Public Code section below. ! ! ! Root biomass areal density. For default values, refer to the table ! in the Public Code section below. ! ! ! Roughness length for momentum. For default values, refer to the table ! in the Public Code section below. ! ! ! ln (z_0_momentum / z_0_scalar). For default values, refer to the table ! in the Public Code section below. ! ! ! Snow amount that half hides surface. For default values, refer to the ! table in the Public Code section below. ! ! ! Snow-free albedo at freezing point. For default values, refer to ! the table in the Public Code section below. ! ! ! Snow-free albedo at t_range below freezing. For default values, refer to ! the table in the Public Code section below. ! ! ! Snow albedo at freezing point. For default values, refer to ! the table in the Public Code section below. ! ! ! Snow albedo at t_range below freezing. For default values, refer to ! the table in the Public Code section below. ! ! ! Available water capacity. For default values, refer to the table in ! the Public Code section below. ! ! ! Volumetric heat capacity. For default values, refer to the table in the ! Public Code section below. ! ! ! Thermal diffusivity. For default values, refer to the table in the ! Public Code section below. ! ! ! If true, the topographic momentum drag scaling scheme is used ! ! ! Maximum of topographic "roughness length" used for momentum drag scaling ! ! ! Scaling factor to convert topography variance to topographic ! "roughness length" ! ! ! Source of the sub-grid topography variance data for topographic momentum drag scaling. ! 'computed' means that the variance is calculated based on high-resolution ! topography data. 'input' means that the data will be provided in specified file ! (NetCDF of IEEE binary) ! ! ! Name of the file to be used as an input for sub-grid topography variance data. ! The file can be either NetCDF (in this case variable name can also be specified), or ! IEEE. ! ! ! Name of the NetCDF variable to be used as a topography variance field. Ignored if ! the file specified in topo_rough_file is not NetCDF file. ! ! ! ! the following namelist vectors contain properties indexed by cover (veg) type: ! 0 ocean ! 1 (BE) broadleaf evergreen trees ! 2 (BD) broadleaf deciduous trees ! 3 (BN) broadleaf/needleleaf trees ! 4 (NE) needleleaf evergreen trees ! 5 (ND) needleleaf deciduous trees ! 6 (G) grassland ! 7 (D) desert ! 8 (T) tundra ! 9 (A) agriculture ! 10 (I) ice ! 11 (L) lake ! 12 (TV) tuned global vegetation ! 13 (SV) Manabe Climate Model-like vegetation ! 14 (SI) Manabe Climate Model-like ice/glacier ! veg_rs_min_vec minimum bulk stomatal resistance s/m ! veg_zeta_vec depth scale of root distribution m ! veg_root_mass_vec root biomass areal density kg/m**2 ! rough_momentum_vec roughness length for momentum m ! k_over_B_vec ln (z_0_momentum / z_0_scalar) - ! crit_snowmass_vec snow amount that half hides surface kg/m**3 ! min_nosnow_alb_vec snow-free albedo at freezing point - ! max_nosnow_alb_vec snow-free albedo at t_range below freezing - ! min_snow_alb_vec snow albedo at freezing point - ! max_snow_alb_vec snow albedo at t_range below freezing - real, dimension(n_dim_cover_types) :: & ! BE BD BN NE ND G D T A I L TV SV SI veg_rs_min_vec =(/ 43.6, 131., 87.1, 69.7, 218., 56.6, .01, 170., 56.6, .01, .01, 67., .01, .01/),& veg_zeta_vec =(/ .26, .29, .35, .17, .17, .26, .35, .11, .25, 0.0, 1.0, .35, .35, 0.0/),& veg_root_mass_vec =(/ 4.9, 4.2, 4.3, 2.9, 2.9, 1.4, .762, 1.2, .15, 0.0, 1.0, .362, .762, 0.0/),& rough_momentum_vec=(/ 2.65, .90, 1.2, .90, .80, .07, .01, .07, .40, .01, 1.4e-4, 1.0, .045, .045/),& k_over_B_vec =(/ 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 0.25, 2., 0., 0. /),& crit_snowmass_vec =(/ 60., 10., 25., 40., 40., 5., 5., 5., 5., 5., 5., 100., 5., 5. /),& min_nosnow_alb_vec=(/0.149, 0.130, 0.132, 0.126, 0.143, 0.182, 0.333, 0.139, 0.160, 0.650, 0.06, 0.12, 999., 0.55/),& max_nosnow_alb_vec=(/0.149, 0.130, 0.132, 0.126, 0.143, 0.182, 0.333, 0.139, 0.160, 0.800, 0.06, 0.12, 999., 0.65/),& min_snow_alb_vec =(/0.650, 0.650, 0.650, 0.650, 0.650, 0.650, 0.650, 0.650, 0.650, 0.650, 0.06, 0.450, 0.450, 0.650/),& max_snow_alb_vec =(/0.800, 0.800, 0.800, 0.800, 0.800, 0.800, 0.800, 0.800, 0.800, 0.800, 0.06, 0.600, 0.600, 0.800/) ! BE BD BN NE ND G D T A I L TV SV SI ! the following vectors contain properties indexed by ground (soil) type: ! 0 ocean ! 1 (C) coarse soil ! 2 (M) medium soil ! 3 (F) fine soil ! 4 (CM) coarse/medium mix ! 5 (CF) coarse/fine mix ! 6 (MF) medium/fine mix ! 7 (CMF) coarse/medium/fine mix ! 8 (P) organic soil (peat) ! 9 (I) ice ! 10 (L) lake ! 11 (MCM) Manabe Climate Model ! soil_awc_vec available water capacity kg/m**3 ! soil_therm_cap_vec volumetric heat capacity J/(K m**3) ! soil_therm_dif_vec thermal diffusivity m**2/s real, dimension(n_dim_ground_types) :: & ! C M F CM CF MF CMF P I L MCM soil_awc_vec = & (/ 63., 132., 109., 98., 86., 120., 101., 445., 1000., 1000., 150./),& soil_therm_cap_vec= & (/ 1.8e6, 2.0e6, 2.6e6, 1.9e6, 2.2e6, 2.3e6, 2.1e6, 3.0e6, 1.6e6, 8.4e7, 1.0/),& soil_therm_dif_vec=& (/8.3e-7,4.0e-7,5.2e-7,6.2e-7,6.8e-7,4.6e-7,5.8e-7,1.3e-7,1.1e-6, 1.0, 2.0e-7/) ! C M F CM CF MF CMF P I L MCM ! logical :: use_topo_rough = .false. real :: max_topo_rough = 100 ! m real :: topo_rough_factor = 1.0 character(len=16) :: topo_rough_source = 'computed' character(len=256):: topo_rough_file = 'INPUT/mg_drag.data.nc' character(len=128):: topo_rough_var = 'ghprime' namelist /land_properties_nml/ do_all_mcm, veg_to_use, & soil_to_use, use_glaciers, & use_climap, use_desert_albedo_map, & do_mcm_masking, & use_single_basin, use_single_geo, & i_dest0, j_dest0, & veg_index_constant, soil_index_constant,& soil_index_ice_substitute, & geo_res_time, t_range, & factor_root, factor_stomata, & factor_rough, z_root_min, & max_snow, & veg_rs_min_vec, veg_zeta_vec, & veg_root_mass_vec, & rough_momentum_vec, k_over_B_vec, & crit_snowmass_vec, & min_nosnow_alb_vec, max_nosnow_alb_vec, & min_snow_alb_vec, max_snow_alb_vec, & soil_awc_vec, & soil_therm_cap_vec, soil_therm_dif_vec, & use_topo_rough, max_topo_rough, topo_rough_factor, & sfc_heat_factor, num_sfc_layers, & dynamic_cover_type, read_old_ascii_cover, & cover_dataset_init_year, cover_dataset_entry, & reconcile_lakes, soil_index_lake_substitute ! ---- private data ---------------------------------------------------------- logical :: module_is_initialized =.FALSE. character(len=128) :: version = '$Id: land_properties.F90,v 15.0 2007/08/14 04:00:11 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !---- other module variables and named constants ! nlon, nlat numbers of longitudes, latitudes in land grid ! npart number of partitions of each land grid cell ! ground_type lon/lat/tile array of ground type index ! albedo_min_no_snow lon/lat/tile array of min snow-free land albedo ! albedo_max_no_snow lon/lat/tile array of max snow-free land albedo ! critical_root_density root density at the bottom of the root zone (kg/m**3) ! max_snow_ssl Manabe Climate Model value of max_snow (kg/m**3) ! res_time_ssl Manabe Climate Model value of groundwater res. time (s) ! veg_index_desert location of desert parameters in veg vecs ! veg_index_tundra location of tundra parameters in veg vecs ! veg_index_ice location of standard ice parameters in veg vecs ! veg_index_tuned location of tuned parameters in veg vecs ! veg_index_ssl_veg location of ssl veg parameters in veg vecs ! veg_index_ssl_ice location of ssl ice parameters in veg vecs ! soil_index_coarse location of coarse soil parameters in soil vecs ! soil_index_ice location of ice parameters in soil vecs ! soil_index_ssl location of ssl soil parameters in soil vecs ! topo_stdev standard deviation of topography ! model_init_time model base time - used to calculate cover type offset ! cover_entry cover type dataset entry time integer :: nlon, nlat, npart integer, pointer, dimension(:,:,:) :: ground_type =>NULL() real, pointer, dimension(:,:,:) :: albedo_min_no_snow =>NULL() real, pointer, dimension(:,:,:) :: albedo_max_no_snow =>NULL() real :: critical_root_density = 0.125 real :: critical_glacier_albedo = 0.79 real :: max_snow_ssl = 200. real :: res_time_ssl = 60.*86400. integer :: veg_index_desert = 7 integer :: veg_index_tundra = 8 integer :: veg_index_ice = 10 integer :: veg_index_lake = 11 integer :: veg_index_tuned = 12 integer :: veg_index_ssl_veg = 13 integer :: veg_index_ssl_ice = 14 integer :: soil_index_coarse = 1 integer :: soil_index_ice = 9 integer :: soil_index_lake = 10 integer :: soil_index_ssl = 11 real, allocatable :: topo_stdev(:,:) type(time_type) :: model_init_time, cover_entry type(horiz_interp_type) :: Interp_gw ! ---- name of the diagnostic "module" --------------------------------------- character(len=*), parameter :: diag_mod_name = "soil" ! ---- diagnostic fields ids ------------------------------------------------- integer :: id_alb_max, id_alb_min integer :: id_ground_type integer :: id_topo_stdev include 'netcdf.inc' contains !####################################################################### ! ! ! Initialize land property data. ! ! ! Reads and re-grids each land input data field and allocates arrays. ! The input cover field and implied glacier field are obtained, if either ! will be needed. Then, the ice cover type is reset to tundra. It ! may be changed back to ice according to actual glacier specification ! later. The glacier mask and cover type are redefined, if necessary. ! ! The appropriate cover type to glacier cells is assigned, if any. The ! ground (soil) type is defined and the groundwater residence time is ! calculated. The properties that depend on cover and ground types are ! assigned. If requested, the snow-free albedo is changed to the CLIMAP ! array. ! ! The land properties diagnostics are initialized and the static fields ! are sent to the diagnostic manager for output. ! ! ! subroutine land_properties_init ( & lonb, latb, land, time, domain, & glacier, & lake, & rough_momentum, & rough_heat, & rough_scale, & ! topographic roughness for drag scaling soil_therm_con, & soil_therm_cap, & veg_rs_min, & max_water, & tau_groundwater, & max_snow_out, & cover_type, & id_lon, id_lat ) real, intent(in) :: lonb(:,:), latb(:,:) ! corners of the cells, ! radian logical, intent(in) :: land(:,:,:) ! land mask type(time_type), intent(in) :: time ! current time type(domain2d), intent(in) :: domain ! our domain logical, intent(out), dimension(:,:,:) :: glacier, lake ! glacier, lake masks real , intent(out), dimension(:,:,:) :: & rough_momentum, & ! roughness length for momentum rough_heat, & rough_scale, & veg_rs_min, & max_water, & tau_groundwater ! groundwater residence time real , intent(out), dimension(:,:,:,:) :: & soil_therm_con, & soil_therm_cap real , intent(out) :: max_snow_out integer, intent(out), dimension(:,:,:) :: cover_type integer, intent(in) :: id_lon, id_lat ! ids of diagnostic axes ! ! ---- local vars ---------------------------------------------------------- integer, allocatable :: integer_input_1(:,:) real, allocatable :: real_input (:,:) real, allocatable :: tmp(:,:,:) ! temp array for read data for cover_type integer, allocatable :: cover_type_input (:,:,:) ! cover_type input data integer, allocatable :: nlon_input (:) ! input data lon integer, allocatable :: nlat_input (:) ! input data lat real, allocatable :: nblon_input (:) ! input data western boundary real, dimension(size(land,1),size(land,2)):: stdev ! standard deviation of ! orography logical got_stdev character*80 input_format integer :: unit, io, ierr, i_cover_type, i_ground_type integer :: ii, jj, num_lon_input, num_lat_input, k real :: wb_degrees ! western boundary of the input grid cells real :: sb, wb, dx, dy, z_root integer :: is,ie, js,je ! boundaries of our domain integer :: gnlon, gnlat integer :: siz(4) ! size of field in field_size call integer :: model_yr, model_mo, model_dy, model_hr, model_mn, model_sc integer :: init_yr, init_mo, init_dy, init_hr, init_mn, init_sc integer :: kk ! time loop integer for cover type data module_is_initialized = .TRUE. ! get the size of our domain call mpp_get_compute_domain ( domain, is,ie,js,je ) call mpp_get_global_domain ( domain, xsize=gnlon, ysize=gnlat) !---- read and write namelist and version ----------------- if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr = 1; do while (ierr /= 0) read (unit, nml=land_properties_nml, iostat=io, end=10) ierr = check_nml_error (io, 'land_properties_nml') enddo 10 continue call close_file (unit) endif call write_version_number( version, tagname) if (mpp_pe() == mpp_root_pe()) then unit = stdlog() write (unit, nml=land_properties_nml) call close_file (unit) endif if (do_all_mcm) then veg_to_use = 'cons_ssl' soil_to_use = 'cons_ssl' do_mcm_masking = .true. use_single_geo = .true. geo_res_time = res_time_ssl factor_root = 1. factor_rough = 1. z_root_min = 0. max_snow = max_snow_ssl use_climap_mcm = .true. endif if (veg_to_use .eq. 'cons_ssl') then use_climap = .true. endif ! allocate private arrays nlon = size(lonb,1) - 1 nlat = size(latb,2) - 1 npart = size(land,3) allocate (ground_type (nlon, nlat, npart)) allocate (albedo_min_no_snow(nlon, nlat, npart)) allocate (albedo_max_no_snow(nlon, nlat, npart)) allocate (topo_stdev(nlon,nlat)) ! Read and re-grid each land input data field sb = -pi/2. ! !---- Cover (vegetation) type. There are 12 possible cases: veg_to_use can ! take 4 values, and glaciers can be (1) absent, (2) based on input cover ! field, or (3) based on climap albedo field. (use of climap albedo forces ! use of climap to locate glaciers, if glaciers are used.) the code ! ensures consistency between cover_type field and glacier mask, with ! the latter taking precedence: where glacier based on input cover field ! must be removed in deference to climap or if no glaciers are used, ! such points are assigned tundra (or global constant) type. ! !--- first get input cover field and implied glacier field, if either of ! these will be needed. then re-set ice cover type to tundra; it may be ! changed back to ice according to actual glacier specification later. if ( veg_to_use.eq.'variable' & .or. (use_glaciers.and..not.use_climap) ) then ! check that dynamic_cover_type and read_old_ascii_cover both are not true if (dynamic_cover_type .and. read_old_ascii_cover) & call error_mesg('land_properties_init', & 'The namelist variables dynamic_cover_type and read_old_ascii_cover cannot both have .true. values',FATAL) ! If dynamic_cover_type_netcdf is chosen, then the cover_type will be obtained ! from a netCDF file of annual data. The initial year of the data is ! specified in the namelist option cover_dataset_init_year. if (dynamic_cover_type) then if(.not.file_exist('INPUT/cover_type_field.nc')) & call error_mesg('land_properties_init','Cannot find file INPUT/cover_type_field.nc',FATAL) ! ! The cover type field file cannot be found. Provide this file or set up ! namelist parameters so it is not necessary. To do the latter, set ! veg_to_use to something other than 'variable' in the namelist ! land_properties_nml. ! call error_mesg('land_properties_init', 'Using dynamic cover type forcing data',NOTE) ! obtain the size on lon in the input data call field_size('INPUT/cover_type_field.nc', 'lon', siz) allocate (nlon_input(siz(1))) num_lon_input = size(nlon_input) ! obtain the size of lat in the input data call field_size('INPUT/cover_type_field.nc', 'lat', siz) allocate (nlat_input(siz(1))) num_lat_input = size(nlat_input) ! allocate arrays, read data as real and convert back to integer allocate (tmp(num_lon_input, num_lat_input, npart)) ! real temp array allocate (cover_type_input (num_lon_input, num_lat_input, npart)) ! input data ! check to make sure cover_dataset_entry has been entered if (cover_dataset_entry(1) == 1 .and. & cover_dataset_entry(2) == 1 .and. & cover_dataset_entry(3) == 1 .and. & cover_dataset_entry(4) == 0 .and. & cover_dataset_entry(5) == 0 .and. & cover_dataset_entry(6) == 0 ) & call error_mesg ('land_properties_mod', & 'must set cover_dataset_entry when using dynamic cover type inputs', FATAL) ! ! When specifying dynamic cover type forcing data inputs, the ! cover_dataset_entry time must be set. ! ! get the model initial time, current model time and cover type dataset time model_init_time = get_base_time() call get_date(model_init_time, init_yr, init_mo, init_dy, init_hr, init_mn, init_sc) call get_date(time, model_yr, model_mo, model_dy, model_hr, model_mn, model_sc) cover_entry = set_date (cover_dataset_entry(1), & cover_dataset_entry(2), & cover_dataset_entry(3), & cover_dataset_entry(4), & cover_dataset_entry(5), & cover_dataset_entry(6)) ! calculate the initial timelevel to be read from the netcdf file kk = ( model_yr - init_yr ) + & ( cover_dataset_entry(1) - cover_dataset_init_year ) if (mpp_pe() == mpp_root_pe()) then write (stdout(),'(a45,i4)') & 'The cover type dataset begins with year:', cover_dataset_init_year write (stdout(),'(a40,i4)') & 'The cover type data year being read is:', cover_dataset_entry(1) write (stdout(),'(a60,i4)') & 'This cover type data is mapped to current model year:', model_yr write (stdout(),'(a42,i4)') & 'The model init year:', init_yr endif ! read the initial cover type data from the netcdf file call read_data('INPUT/cover_type_field.nc','cover',tmp(:,:,1), & timelevel=kk+1, no_domain=.true.) cover_type_input(:,:,1) = int(tmp(:,:,1)) ! obtain the western boundary in the input data call field_size('INPUT/cover_type_field.nc', 'blon', siz) allocate (nblon_input(siz(1))) call read_data('INPUT/cover_type_field.nc','blon',nblon_input(:),timelevel=1,no_domain=.true.) ! regrid the data to the land model grid wb_degrees = (nblon_input(1)) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call regrid_discrete_field(cover_type_input(:,:,1), wb, sb, dx, dy, lonb,& latb, n_map_cover_types, land(:,:,1), cover_type(:,:,1) ) deallocate (cover_type_input) ! if using the read_old_ascii_cover option, read the ASCII data file for cover ! type elseif (read_old_ascii_cover) then call error_mesg('land_properties_init', 'Using static ASCII cover type dataset',NOTE) ! ! Using the read_old_ascii_cover option. Reading the ASCII cover type dataset. ! call mpp_open (unit, 'INPUT/cover_type_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( integer_input_1 (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (integer_input_1(ii,jj),ii=1,num_lon_input) end do close(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call regrid_discrete_field(integer_input_1, wb, sb, dx, dy, lonb, latb, & n_map_cover_types, land(:,:,1), cover_type(:,:,1) ) deallocate (integer_input_1) call close_file(unit) ! if not using dynamic_cover_type nor the read_old_ascii_cover option, use the ! static netcdf cover_type dataset if it exists elseif(file_exist('INPUT/cover_type_field.nc')) then call error_mesg('land_properties_init', 'Using the static netcdf cover type dataset',NOTE) ! ! Using the static NetCDF cover type dataset. ! ! obtain the size on lon in the input data call field_size('INPUT/cover_type_field.nc', 'lon', siz) allocate (nlon_input(siz(1))) num_lon_input = size(nlon_input) ! obtain the size of lat in the input data call field_size('INPUT/cover_type_field.nc', 'lat', siz) allocate (nlat_input(siz(1))) num_lat_input = size(nlat_input) ! allocate arrays, read data as real and convert back to integer allocate (tmp(num_lon_input, num_lat_input, npart)) ! real temp array allocate (cover_type_input (num_lon_input, num_lat_input, npart)) ! input data if (cover_dataset_entry(1) == 1 .and. & cover_dataset_entry(2) == 1 .and. & cover_dataset_entry(3) == 1 .and. & cover_dataset_entry(4) == 0 .and. & cover_dataset_entry(5) == 0 .and. & cover_dataset_entry(6) == 0 ) & call error_mesg ('land_properties_mod', & 'must set cover_dataset_entry when using static_cover_type_netcdf', FATAL) ! ! When specifying static_cover_type_netcdf forcing data inputs, the ! cover_dataset_entry time must be set. ! ! get the model initial time and cover type dataset time model_init_time = get_base_time() cover_entry = set_date (cover_dataset_entry(1), & cover_dataset_entry(2), & cover_dataset_entry(3), & cover_dataset_entry(4), & cover_dataset_entry(5), & cover_dataset_entry(6)) if (mpp_pe() == mpp_root_pe()) then write (stdout(),'(a45,i4)') & 'The cover type dataset begins with year:', cover_dataset_init_year write (stdout(),'(a40,i4)') & 'The cover type data year being read is:', cover_dataset_entry(1) endif ! calculate the static timelevel to read from the NetCDF file kk = cover_dataset_entry(1) - cover_dataset_init_year call read_data('INPUT/cover_type_field.nc','cover',tmp(:,:,1),timelevel=kk+1, no_domain=.true.) cover_type_input(:,:,1) = int(tmp(:,:,1)) ! obtain the western boundary in the input data call field_size('INPUT/cover_type_field.nc', 'blon', siz) allocate (nblon_input(siz(1))) call read_data('INPUT/cover_type_field.nc','blon',nblon_input(:),timelevel=1, no_domain=.true.) ! regrid the data onto the land grid wb_degrees = (nblon_input(1)) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call regrid_discrete_field(cover_type_input(:,:,1), wb, sb, dx, dy, lonb,& latb, n_map_cover_types, land(:,:,1), cover_type(:,:,1) ) deallocate (cover_type_input) ! if the static cover type forcing data option is chosen and the NetCDF file ! does not exist, read the GSWP2 cover type field in ASCII format else call error_mesg('land_properties_init','Cannot find the netcdf file INPUT/cover_type_field. Using the ASCII file.',NOTE) ! ! The netCDF cover type field file cannot be found. Using the ASCII file. ! call mpp_open (unit, 'INPUT/cover_type_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( integer_input_1 (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (integer_input_1(ii,jj),ii=1,num_lon_input) end do close(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call regrid_discrete_field(integer_input_1, wb, sb, dx, dy, lonb, latb, & n_map_cover_types, land(:,:,1), cover_type(:,:,1) ) deallocate (integer_input_1) call close_file(unit) endif endif !--- update cover type where (.not.land(:,:,1)) cover_type(:,:,1) = 0 glacier(:,:,1) = .false. where (cover_type(:,:,1).eq.veg_index_ice) glacier(:,:,1) = .true. cover_type(:,:,1) = veg_index_tundra endwhere !--- now re-define glacier mask, if necessary if (.not.use_glaciers) then glacier(:,:,1) = .false. else if (use_climap) then if(use_climap_mcm) then call get_climap_glacier_mcm(lonb, latb, gnlon, gnlat, is, is+size(land,1)-1, & js, js+size(land,2)-1, critical_glacier_albedo, glacier(:,:,1)) else call get_climap_glacier(lonb, latb, critical_glacier_albedo, glacier(:,:,1)) endif where (.not.land(:,:,1)) glacier(:,:,1) = .false. endif !--- next re-define cover type, if necessary, not yet worrying !--- about glacier type if (veg_to_use.eq.'constant') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_constant else if (veg_to_use.eq.'cons_tun') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_tuned else if (veg_to_use.eq.'cons_ssl') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_ssl_veg endif !--- finally, assign appropriate cover type to glacier cells (if any) if (veg_to_use.ne.'cons_ssl') then where (glacier(:,:,1)) cover_type(:,:,1) = veg_index_ice else where (glacier(:,:,1)) cover_type(:,:,1) = veg_index_ssl_ice endif ! call error_mesg ('land_properties_init', 'cover type defined', NOTE) ! !---- Ground (soil) type. To force consistency with glacier, already ! defined, any ice soil types are first re-set to coarse soil. ! ice types are then located on the basis of glacier. note that a distinct ! ice type is used only if soil_to_use='variable'; this is in contrast ! to the analogous treatment of vegetation types. ! if (soil_to_use .eq. 'variable') then if (.not.file_exist('INPUT/ground_type_field')) & call error_mesg('land_properties_init','Cannot find file INPUT/ground_type_field',FATAL) ! ! The ground (soil) type field file cannot be found. Provide this file or ! set up namelist parameters so it is not necessary. To do the latter, set ! soil_to_use to something other than 'variable' in the namelist ! land_properties_nml. ! call mpp_open (unit, 'INPUT/ground_type_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( integer_input_1 (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (integer_input_1(ii,jj),ii=1,num_lon_input) end do close(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) if (soil_index_ice_substitute.eq.0) then where (integer_input_1.eq.soil_index_ice) integer_input_1 = 0 endif call regrid_discrete_field(integer_input_1, wb, sb, dx, dy, lonb, latb, & n_map_ground_types, land(:,:,1), ground_type(:,:,1) ) deallocate (integer_input_1) where (.not.land(:,:,1)) ground_type(:,:,1) = 0 if (soil_index_ice_substitute.ne.0) then where (ground_type(:,:,1).eq.soil_index_ice) & ground_type(:,:,1) = soil_index_ice_substitute endif where (glacier(:,:,1)) ground_type(:,:,1) = soil_index_ice ! reconcile lake distribution in cover type and ground type fields if (reconcile_lakes) then where (ground_type(:,:,1).eq.soil_index_lake) & ground_type(:,:,1) = soil_index_lake_substitute where (cover_type(:,:,1).eq.veg_index_lake) & ground_type(:,:,1) = soil_index_lake endif call close_file(unit) else if (soil_to_use .eq. 'constant') then ground_type(:,:,1) = 0 where (land(:,:,1)) ground_type(:,:,1) = soil_index_constant else if (soil_to_use .eq. 'cons_ssl') then ground_type(:,:,1) = 0 where (land(:,:,1)) ground_type(:,:,1) = soil_index_ssl endif ! call error_mesg ('land_properties_init', 'ground type defined', NOTE) !---- groundwater residence time if (use_single_geo) then where (land(:,:,1)) tau_groundwater(:,:,1) = geo_res_time else if(.not.file_exist('INPUT/groundwater_residence_time_field')) & call error_mesg('land_properties_init','Cannot find file INPUT/groundwater_residence_time_field',FATAL) ! ! The groundwater residence time field file cannot be found. Provide this file ! or set up namelist parameters so it is not necessary. To do the latter, set ! use_single_geo to .true. in the namelist land_properties_nml. ! call mpp_open (unit,'INPUT/groundwater_residence_time_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( real_input (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (real_input(ii,jj),ii=1,num_lon_input) end do call close_file(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call create_horiz_interp_new (Interp_gw, num_lon_input, num_lat_input, & wb, sb, dx, dy, lonb, latb) call horiz_interp ( Interp_gw, real_input, tau_groundwater(:,:,1) ) deallocate (real_input) endif ! call error_mesg ('land_properties_init', 'groundwater residence time defined', NOTE) !---- river destination pointers lake(:,:,1) = .false. where (ground_type(:,:,1).eq.soil_index_lake) lake(:,:,1) = .true. do ii = 2, npart glacier (:,:,ii) = glacier (:,:,1) lake (:,:,ii) = lake (:,:,1) cover_type (:,:,ii) = cover_type (:,:,1) ground_type (:,:,ii) = ground_type (:,:,1) tau_groundwater(:,:,ii) = tau_groundwater(:,:,1) end do !---- assign properties that depend only on cover type ----------- veg_rs_min = 0 rough_momentum = 0 do i_cover_type = 1, n_dim_cover_types where (cover_type.eq.i_cover_type) albedo_min_no_snow = min_nosnow_alb_vec(i_cover_type) albedo_max_no_snow = max_nosnow_alb_vec(i_cover_type) veg_rs_min = factor_stomata * veg_rs_min_vec (i_cover_type) rough_momentum = factor_rough * rough_momentum_vec(i_cover_type) rough_heat = factor_rough * rough_momentum_vec(i_cover_type) & * exp(-k_over_B_vec(i_cover_type)) endwhere enddo !---- compute or read the topographic roughness ---- if (use_topo_rough) then call set_topo_rough (domain, lonb, latb, topo_stdev) stdev = min(topo_stdev*topo_rough_factor,max_topo_rough) do k = 1,size(rough_momentum,3) rough_scale(:,:,k) = max(stdev,rough_momentum(:,:,k)) end do else rough_scale = rough_momentum endif !---- assign properties that depend only on ground type ---------- do i_ground_type = 1, n_dim_ground_types where (ground_type.eq.i_ground_type) soil_therm_con(:,:,:,1) = soil_therm_dif_vec(i_ground_type) & *soil_therm_cap_vec(i_ground_type) soil_therm_cap(:,:,:,1) = soil_therm_cap_vec(i_ground_type) endwhere enddo do ii=2,size(soil_therm_con,4) where (land(:,:,:)) soil_therm_con(:,:,:,ii) = soil_therm_con(:,:,:,1) soil_therm_cap(:,:,:,ii) = soil_therm_cap(:,:,:,1) endwhere enddo do ii=1,num_sfc_layers where (land(:,:,:)) soil_therm_con(:,:,:,ii) = sfc_heat_factor * soil_therm_con(:,:,:,ii) soil_therm_cap(:,:,:,ii) = sfc_heat_factor * soil_therm_cap(:,:,:,ii) endwhere enddo !---- assign properties that depend on cover and ground types ---- max_water = 0 do i_cover_type = 1, n_dim_cover_types if ((i_cover_type .ne. veg_index_ice) .and. & (i_cover_type .ne. veg_index_ssl_ice)) then z_root = veg_zeta_vec(i_cover_type) & *log( veg_root_mass_vec(i_cover_type)/ & (veg_zeta_vec(i_cover_type)*factor_root*critical_root_density) ) z_root = max(z_root_min, z_root) do i_ground_type = 1, n_dim_ground_types where (cover_type.eq.i_cover_type.and.ground_type.eq.i_ground_type) max_water = z_root * soil_awc_vec(i_ground_type) endwhere enddo endif enddo !---- change snow-free albedo to CLIMAP array if requested ------- ! if (use_climap) then allocate (real_input(nlon, nlat)) if(use_climap_mcm) then call get_climap_albedo_mcm(lonb, latb, gnlon, gnlat, is, is+size(land,1)-1, & js, js+size(land,2)-1, 1, real_input) else call get_climap_albedo(lonb, latb, 1, real_input) endif do ii = 1, npart albedo_min_no_snow(:,:,ii) = real_input albedo_max_no_snow(:,:,ii) = real_input if(use_climap_mcm) then where(glacier(:,:,ii)) albedo_min_no_snow(:,:,ii) = min_nosnow_alb_vec(veg_index_ssl_ice) albedo_max_no_snow(:,:,ii) = max_nosnow_alb_vec(veg_index_ssl_ice) endwhere endif enddo deallocate (real_input) elseif (use_desert_albedo_map) then allocate (real_input(nlon, nlat)) call get_climap_albedo(lonb, latb, 5, real_input) do ii = 1, npart where (cover_type(:,:,ii).eq.veg_index_desert) albedo_min_no_snow(:,:,ii) = real_input albedo_max_no_snow(:,:,ii) = real_input endwhere enddo deallocate (real_input) endif max_snow_out = max_snow ! initialize land properties diagnostics call init_land_properties_diag(id_lon, id_lat, time) ! send static fields for diagnostic manager for output call diag_static(time, land) end subroutine land_properties_init ! !####################################################################### ! ! ! Updates slowly changing parameters, such as cover type. ! ! ! Updates slowly changing parameters, such as cover type. ! ! subroutine update_land_properties_slow ( & lonb, latb, land, time, domain, & glacier, & lake, & rough_momentum, & rough_heat, & rough_scale, & ! topographic roughness for drag scaling soil_therm_con, & soil_therm_cap, & veg_rs_min, & max_water, & tau_groundwater, & max_snow_out, & cover_type ) real, intent(in) :: lonb(:,:), latb(:,:) ! corners of the cells, ! radian logical, intent(in) :: land(:,:,:) ! land mask type(time_type), intent(in) :: time ! current time type(domain2d), intent(in) :: domain ! our domain logical, intent(out), dimension(:,:,:) :: glacier ! glacier mask logical, intent(out), dimension(:,:,:) :: lake ! lake mask real , intent(out), dimension(:,:,:) :: & rough_momentum, & ! roughness length for momentum rough_heat, & rough_scale, & veg_rs_min, & max_water, & tau_groundwater ! groundwater residence time real , intent(out), dimension(:,:,:,:) :: & soil_therm_con, & soil_therm_cap real , intent(out) :: max_snow_out integer, intent(inout), dimension(:,:,:) :: cover_type ! ! ---- local vars ---------------------------------------------------------- integer, allocatable :: integer_input_1(:,:) real, allocatable :: real_input (:,:) real, allocatable :: tmp(:,:,:) ! temp array for read data for cover_type integer, allocatable :: cover_type_input (:,:,:) ! cover_type input data integer, allocatable :: nlon_input (:) ! input data lon integer, allocatable :: nlat_input (:) ! input data lat real, allocatable :: nblon_input (:) ! input data western boundary real, dimension(size(land,1),size(land,2)):: stdev ! standard deviation of ! orography logical got_stdev character*80 input_format integer :: unit, ierr, i_cover_type, i_ground_type integer :: ii, jj, num_lon_input, num_lat_input, k real :: wb_degrees ! western boundary of the input grid cells real :: sb, wb, dx, dy, z_root integer :: is, js ! boundaries of our domain integer :: gnlon, gnlat integer :: siz(4) ! size of field in field_size call integer :: model_yr, model_mo, model_dy, model_hr, model_mn, model_sc integer :: init_yr, init_mo, init_dy, init_hr, init_mn, init_sc integer :: kk ! time loop integer for cover type data ! call only if using time-varying cover type if (dynamic_cover_type) then ! Read and re-grid each land input data field sb = -pi/2. ! !---- Cover (vegetation) type. There are 12 possible cases: veg_to_use can ! take 4 values, and glaciers can be (1) absent, (2) based on input cover ! field, or (3) based on climap albedo field. (use of climap albedo forces ! use of climap to locate glaciers, if glaciers are used.) the code ! ensures consistency between cover_type field and glacier mask, with ! the latter taking precedence: where glacier based on input cover field ! must be removed in deference to climap or if no glaciers are used, ! such points are assigned tundra (or global constant) type. ! !--- first get input cover field and implied glacier field, if either of ! these will be needed. then re-set ice cover type to tundra; it may be ! changed back to ice according to actual glacier specification later. ! If dynamic_cover_type_netcdf is chosen, then the cover_type will be obtained ! from a netCDF file of annual data. The initial year of the dataset is ! specified by the namelist option cover_dataset_init_year. if ( veg_to_use.eq.'variable' & .or. (use_glaciers.and..not.use_climap) ) then if(.not.file_exist('INPUT/cover_type_field.nc')) & call error_mesg('land_properties_slow','Cannot find file INPUT/cover_type_field.nc',FATAL) ! ! The cover type field file cannot be found. Provide this file or set up ! namelist parameters so it is not necessary. To do the latter, set ! veg_to_use to something other than 'variable' in the namelist ! land_properties_nml. ! ! obtain the size on lon in the input data call field_size('INPUT/cover_type_field.nc', 'lon', siz) allocate (nlon_input(siz(1))) num_lon_input = size(nlon_input) ! obtain the size of lat in the input data call field_size('INPUT/cover_type_field.nc', 'lat', siz) allocate (nlat_input(siz(1))) num_lat_input = size(nlat_input) ! allocate arrays, read data as real and convert back to integer allocate (tmp(num_lon_input, num_lat_input, npart)) ! real temp array allocate (cover_type_input (num_lon_input, num_lat_input, npart)) ! input data ! get the current year of the model call get_date(model_init_time, init_yr, init_mo, init_dy, init_hr, init_mn, init_sc) call get_date(time, model_yr, model_mo, model_dy, model_hr, model_mn, model_sc) ! calculate the timelevel read from the netcdf file kk = ( model_yr - init_yr ) + & ( cover_dataset_entry(1) - cover_dataset_init_year ) if (mpp_pe() == mpp_root_pe()) then write (stdout(),'(a45,i4)') & 'The cover type dataset begins with year:', cover_dataset_init_year write (stdout(),'(a40,i4)') & 'The cover type data year being read is:', cover_dataset_entry(1) write (stdout(),'(a60,i4)') & 'This cover type data is mapped to current model year:', model_yr write (stdout(),'(a42,i4)') & 'The model init year:', init_yr endif ! read the cover_type data from the netcdf file and covert it back to integer call read_data('INPUT/cover_type_field.nc','cover',tmp(:,:,1),timelevel=kk+1,no_domain=.true.) cover_type_input(:,:,1) = int(tmp(:,:,1)) ! obtain the western boundary in the input data call field_size('INPUT/cover_type_field.nc', 'blon', siz) allocate (nblon_input(siz(1))) call read_data('INPUT/cover_type_field.nc','blon',nblon_input(:),timelevel=1,no_domain=.true.) ! regrid the data to the land grid and update the cover type wb_degrees = (nblon_input(1)) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) call regrid_discrete_field(cover_type_input(:,:,1), wb, sb, dx, dy, lonb,& latb, n_map_cover_types, land(:,:,1), cover_type(:,:,1) ) deallocate (cover_type_input) where (.not.land(:,:,1)) cover_type(:,:,1) = 0 glacier(:,:,1) = .false. where (cover_type(:,:,1).eq.veg_index_ice) glacier(:,:,1) = .true. cover_type(:,:,1) = veg_index_tundra endwhere endif !--- next re-define cover type, if necessary, not yet worrying !--- about glacier type if (veg_to_use.eq.'constant') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_constant else if (veg_to_use.eq.'cons_tun') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_tuned else if (veg_to_use.eq.'cons_ssl') then cover_type(:,:,1) = 0 where (land(:,:,1)) cover_type(:,:,1) = veg_index_ssl_veg endif !--- finally, assign appropriate cover type to glacier cells (if any) if (veg_to_use.ne.'cons_ssl') then where (glacier(:,:,1)) cover_type(:,:,1) = veg_index_ice else where (glacier(:,:,1)) cover_type(:,:,1) = veg_index_ssl_ice endif ! call error_mesg ('land_properties_slow', 'cover type defined', NOTE) ! !---- Ground (soil) type. Any ice soil types are first re-set to coarse soil. ! ice types are then located on the basis of glacier. note that a distinct ! ice type is used only if soil_to_use='variable'; this is in contrast ! to the analogous treatment of vegetation types. ! if (soil_to_use .eq. 'variable') then if (.not.file_exist('INPUT/ground_type_field')) & call error_mesg('land_properties_slow','Cannot find file INPUT/ground_type_field',FATAL) ! ! The ground (soil) type field file cannot be found. Provide this file or ! set up namelist parameters so it is not necessary. To do the latter, set ! soil_to_use to something other than 'variable' in the namelist ! land_properties_nml. ! call mpp_open (unit, 'INPUT/ground_type_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( integer_input_1 (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (integer_input_1(ii,jj),ii=1,num_lon_input) end do close(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) if (soil_index_ice_substitute.eq.0) then where (integer_input_1.eq.soil_index_ice) integer_input_1 = 0 endif call regrid_discrete_field(integer_input_1, wb, sb, dx, dy, lonb, latb,& n_map_ground_types, land(:,:,1), ground_type(:,:,1) ) deallocate (integer_input_1) where (.not.land(:,:,1)) ground_type(:,:,1) = 0 if (soil_index_ice_substitute.ne.0) then where (ground_type(:,:,1).eq.soil_index_ice) & ground_type(:,:,1) = soil_index_ice_substitute endif where (glacier(:,:,1)) ground_type(:,:,1) = soil_index_ice call close_file(unit) else if (soil_to_use .eq. 'constant') then ground_type(:,:,1) = 0 where (land(:,:,1)) ground_type(:,:,1) = soil_index_constant else if (soil_to_use .eq. 'cons_ssl') then ground_type(:,:,1) = 0 where (land(:,:,1)) ground_type(:,:,1) = soil_index_ssl endif ! reconcile lake distribution in cover type and ground type fields if (reconcile_lakes) then where (ground_type(:,:,1).eq.soil_index_lake) & ground_type(:,:,1) = soil_index_lake_substitute where (cover_type(:,:,1).eq.veg_index_lake) & ground_type(:,:,1) = soil_index_lake endif lake(:,:,1) = .false. where (ground_type(:,:,1).eq.soil_index_lake) lake(:,:,1) = .true. ! call error_mesg ('land_properties_slow', 'ground type defined', NOTE) !---- groundwater residence time if (use_single_geo) then where (land(:,:,1)) tau_groundwater(:,:,1) = geo_res_time else if(.not.file_exist('INPUT/groundwater_residence_time_field')) & call error_mesg('land_properties_slow','Cannot find file INPUT/groundwater_residence_time_field',FATAL) ! ! The groundwater residence time field file cannot be found. Provide this ! file or set up namelist parameters so it is not necessary. To do the ! latter, set use_single_geo to .true. in the namelist land_properties_nml. ! call mpp_open (unit,'INPUT/groundwater_residence_time_field', action = MPP_RDONLY) read (unit,*) num_lon_input, num_lat_input, wb_degrees read (unit,*) input_format allocate ( real_input (num_lon_input, num_lat_input) ) do jj = 1, num_lat_input read(unit,input_format) (real_input(ii,jj),ii=1,num_lon_input) end do call close_file(unit) wb = pi*wb_degrees/180. dx = 2.*pi/float(num_lon_input) dy = pi/float(num_lat_input) ! initialized already !call create_horiz_interp_new (Interp_gw, num_lon_input, num_lat_input, & ! wb, sb, dx, dy, lonb, latb ) call horiz_interp ( Interp_gw, real_input, tau_groundwater(:,:,1) ) deallocate (real_input) endif ! call error_mesg ('land_properties_slow', 'groundwater residence time defined', NOTE) !---- river destination pointers do ii = 2, npart lake (:,:,ii) = lake (:,:,1) cover_type (:,:,ii) = cover_type (:,:,1) ground_type (:,:,ii) = ground_type (:,:,1) tau_groundwater(:,:,ii) = tau_groundwater(:,:,1) end do !---- assign properties that depend only on cover type ----------- veg_rs_min = 0 rough_momentum = 0 do i_cover_type = 1, n_dim_cover_types where (cover_type.eq.i_cover_type) albedo_min_no_snow = min_nosnow_alb_vec(i_cover_type) albedo_max_no_snow = max_nosnow_alb_vec(i_cover_type) veg_rs_min = factor_stomata * veg_rs_min_vec (i_cover_type) rough_momentum = factor_rough * rough_momentum_vec(i_cover_type) rough_heat = factor_rough * rough_momentum_vec(i_cover_type) & * exp(-k_over_B_vec(i_cover_type)) endwhere enddo !---- compute or read the topographic roughness ---- if (use_topo_rough) then call set_topo_rough (domain, lonb, latb, topo_stdev) stdev = min(topo_stdev*topo_rough_factor,max_topo_rough) do k = 1,size(rough_momentum,3) rough_scale(:,:,k) = max(stdev,rough_momentum(:,:,k)) end do else rough_scale = rough_momentum endif !---- assign properties that depend only on ground type ---------- do i_ground_type = 1, n_dim_ground_types where (ground_type.eq.i_ground_type) soil_therm_con(:,:,:,1) = soil_therm_dif_vec(i_ground_type) & *soil_therm_cap_vec(i_ground_type) soil_therm_cap(:,:,:,1) = soil_therm_cap_vec(i_ground_type) endwhere enddo do ii=2,size(soil_therm_con,4) where (land(:,:,:)) soil_therm_con(:,:,:,ii) = soil_therm_con(:,:,:,1) soil_therm_cap(:,:,:,ii) = soil_therm_cap(:,:,:,1) endwhere enddo do ii=1,num_sfc_layers where (land(:,:,:)) soil_therm_con(:,:,:,ii) = sfc_heat_factor * soil_therm_con(:,:,:,ii) soil_therm_cap(:,:,:,ii) = sfc_heat_factor * soil_therm_cap(:,:,:,ii) endwhere enddo !---- assign properties that depend on cover and ground types ---- max_water = 0 do i_cover_type = 1, n_dim_cover_types if ((i_cover_type .ne. veg_index_ice) .and. & (i_cover_type .ne. veg_index_ssl_ice)) then z_root = veg_zeta_vec(i_cover_type) & *log( veg_root_mass_vec(i_cover_type)/ & (veg_zeta_vec(i_cover_type)*factor_root*critical_root_density) ) z_root = max(z_root_min, z_root) do i_ground_type = 1, n_dim_ground_types where (cover_type.eq.i_cover_type.and.ground_type.eq.i_ground_type) max_water = z_root * soil_awc_vec(i_ground_type) endwhere enddo endif enddo !---- change snow-free albedo to CLIMAP array if requested ------- ! if (use_climap) then allocate (real_input(nlon, nlat)) if(use_climap_mcm) then call get_climap_albedo_mcm(lonb, latb, gnlon, gnlat, is, is+size(land,1)-1, & js, js+size(land,2)-1, 1, real_input) else call get_climap_albedo(lonb, latb, 1, real_input) endif do ii = 1, npart albedo_min_no_snow(:,:,ii) = real_input albedo_max_no_snow(:,:,ii) = real_input if(use_climap_mcm) then where(glacier(:,:,ii)) albedo_min_no_snow(:,:,ii) = min_nosnow_alb_vec(veg_index_ssl_ice) albedo_max_no_snow(:,:,ii) = max_nosnow_alb_vec(veg_index_ssl_ice) endwhere endif enddo deallocate (real_input) elseif (use_desert_albedo_map) then allocate (real_input(nlon, nlat)) call get_climap_albedo(lonb, latb, 5, real_input) do ii = 1, npart where (cover_type(:,:,ii).eq.veg_index_desert) albedo_min_no_snow(:,:,ii) = real_input albedo_max_no_snow(:,:,ii) = real_input endwhere enddo deallocate (real_input) endif max_snow_out = max_snow endif end subroutine update_land_properties_slow ! !####################################################################### ! ! ! Updates the rapidly changing parameters, such as albedo. ! ! ! Updates the rapidly changing parameters. Computes the albedo of ! hypothetical no-snow and deep-snow surfaces and uses snow mass to blend ! snow-free and deep-snow albedo values. ! ! ! subroutine update_land_properties_fast (snowmass, t_sfc, land, albedo, cover_type) !----------------------------------------------------------------------- ! ! INPUT ! snowmass = mass of snow on the ground (in kg/(m**2)) ! t_sfc = surface temperature (in degrees kelvin) ! land = logical land mask ! ! OUTPUT ! albedo = surface albedo real, intent(in), dimension(:,:,:) :: snowmass, t_sfc logical, intent(in) , dimension(:,:,:) :: land real, intent(out), dimension(:,:,:) :: albedo integer, intent(in), dimension(:,:,:) :: cover_type !----------------------------------------------------------------------- ! integer :: i_cover_type real :: t_crit real , dimension(size(t_sfc,1),size(t_sfc,2),size(t_sfc,3)) :: & albedo_no_snow, albedo_hi_snow, blend !---- albedo update ---------------------------------------------------- !--- compute albedo of hypothetical no-snow and deep-snow surfaces t_crit = tfreeze - t_range blend = (t_sfc - t_crit) / t_range where (blend .lt. 0.) blend = 0. where (blend .gt. 1.) blend = 1. where (land(:,:,:)) albedo_no_snow = albedo_max_no_snow & + blend * (albedo_min_no_snow - albedo_max_no_snow) endwhere do i_cover_type = 1, n_dim_cover_types where (cover_type.eq.i_cover_type) albedo_hi_snow = max_snow_alb_vec(i_cover_type) + blend & * (min_snow_alb_vec(i_cover_type) & - max_snow_alb_vec(i_cover_type) ) endwhere enddo !--- use snow mass to blend snow-free and deep-snow albedo values if (do_mcm_masking) then blend = 0. do i_cover_type = 1, n_dim_cover_types where (cover_type.eq.i_cover_type) & blend = 0.5 * sqrt(snowmass/crit_snowmass_vec(i_cover_type)) enddo where (blend .gt. 1.) blend = 1. else blend = 0. do i_cover_type = 1, n_dim_cover_types where (cover_type.eq.i_cover_type) & blend = snowmass/(snowmass+crit_snowmass_vec(i_cover_type)) enddo endif where (cover_type.gt. 0) albedo = albedo_no_snow & + (albedo_hi_snow - albedo_no_snow) * blend end subroutine update_land_properties_fast ! ! ! ! Regrids integer index data to any output grid from a uniformly spaced ! grid using a maximum-area voting scheme. ! ! ! Regrids integer index data to any output grid from a uniformly spaced ! grid using a maximum-area voting scheme. ! ! ! subroutine regrid_discrete_field (data_in, wb_in, sb_in, dlon_in, dlat_in, & lon_out, lat_out, ntype, & mask_out, data_out, mask_in) !----------------------------------------------------------------------- ! input: ! ----- ! data_in input data; dimensioned by mdim x ndim ! stored from south to north ! wb_in longitude corresponding to western boundary of box i=1 ! sb_in latitude corresponding to southern boundary of box j=1 ! dlon_in x axis grid spacing in degrees of longitude ! dlat_in y axis grid spacing in degrees of latitude ! ! lon_out longitudes of output data at grid box corners ! dimensioned by size(data_out,1)+1 by size(data_out,2)+1 ! lat_out latitudes of output data at grid box corners ! dimensioned by size(data_out,1)+1 by size(data_out,2)+1 ! ntype number of land cover types specified ! mask_out output mask that specifies where the data are defined ! ! output: ! ------ ! data_out output number of land cover type ! ! optional ! -------- ! mask_in input mask; =0.0 for data points that should not ! be used or have no data; has the same size as data_in ! !----------------------------------------------------------------------- integer, intent(in) :: ntype integer, intent(in) :: data_in(:,:) real, intent(in) :: sb_in,wb_in, dlat_in, dlon_in real, intent(in) :: lon_out(:,:), lat_out(:,:) logical, intent(in) :: mask_out(:,:) integer, intent(out) :: data_out(:,:) real, intent(in), optional :: mask_in(:,:) ! if (is_latlon(lon_out,lat_out)) then call regrid_discrete_field_latlon (data_in, wb_in, sb_in, dlon_in, dlat_in, & lon_out(:,1), lat_out(1,:), ntype, mask_out, data_out, mask_in) else call regrid_discrete_field_cube (data_in, wb_in, sb_in, dlon_in, dlat_in, & lon_out(:,:), lat_out(:,:), ntype, mask_out, data_out, mask_in) endif end subroutine regrid_discrete_field ! ! subroutine regrid_discrete_field_latlon (data_in, wb_in, sb_in, dlon_in, dlat_in, & lon_out, lat_out, ntype, & mask_out, data_out, mask_in) integer, intent(in) :: ntype integer, intent(in) :: data_in(:,:) real, intent(in) :: sb_in,wb_in, dlat_in, dlon_in real, intent(in) :: lon_out(:), lat_out(:) logical, intent(in) :: mask_out(:,:) integer, intent(out) :: data_out(:,:) real, intent(in), optional :: mask_in(:,:) ! --- local vars ----------------------------------------------------- real, dimension(size(data_in,2)) :: area_in real, dimension(size(lat_out(:))) :: ph ! real, dimension(ntype) :: sumtype ! real, dimension(size(data_out,1),size(data_out,2)) :: totarea integer, dimension(size(data_out,1)) :: is,ie integer, dimension(size(data_out,2)) :: js,je,jsave real, dimension(size(data_out,1)) :: fis,fie,facis,facie real, dimension(size(data_out,2)) :: fjs,fje,facjs,facje,fsave integer i,j,mdim,ndim,m360,nlon,nlat integer ismin,iemax, iis,iie,jjs,jje, kk real ratlat,ratlon,hpie,phs,phn, dsph real ffacis,ffacie,ffacjs,ffacje logical flip, enlarge integer, allocatable :: data(:,:) real, allocatable :: area(:,:) hpie= pi/2.0 mdim=size(data_in,1); ndim=size(data_in,2) m360=int(4.*hpie/dlon_in + 0.001) if (mdim < m360) call error_mesg ('regrid_discrete_field', & 'inner dimension for input data is too small.', FATAL) ! ! ! --- input (hires) grid resolution --- ratlat=1.0/dlat_in ratlon=1.0/dlon_in ! --- area of input (hires) grid boxes --- do j=1,ndim phs=sb_in + float(j-1)*dlat_in phn=sb_in + float(j) *dlat_in phs=min(phs, hpie); phs=max(phs,-hpie) phn=min(phn, hpie); phn=max(phn,-hpie) dsph=sin(phn)-sin(phs) area_in(j) = dsph * dlon_in enddo !----------------------------------------------------------------------- nlon=size(data_out,1); nlat=size(data_out,2) !*********************************************************************** !------ set up latitudinal indexing ------ !------ make sure output grid goes south to north ------ if (lat_out(1) < lat_out(nlat+1)) then ph(1:nlat+1) = lat_out(1:nlat+1) flip = .false. else ph(1:nlat+1) = lat_out(nlat+1:1:-1) flip = .true. endif fjs(1:nlat)=(ph(1:nlat )-sb_in)*ratlat+1.0 fje(1:nlat)=(ph(2:nlat+1)-sb_in)*ratlat+1.0 js=fjs where (fjs < 0.) js=js-1 where (js < 1 ) js=1 je=fje where (fje < 0. ) je=je-1 where (je > ndim) je=ndim facjs=float(js+1)-fjs facje=fje-float(je) if (flip) then jsave=js; js(1:nlat)=jsave(nlat:1:-1) jsave=je; je(1:nlat)=jsave(nlat:1:-1) fsave=facjs; facjs(1:nlat)=fsave(nlat:1:-1) fsave=facje; facje(1:nlat)=fsave(nlat:1:-1) endif !------ set up longitudinal indexing ------ fis(1:nlon)=(lon_out(1:nlon )-wb_in)*ratlon+1.0 fie(1:nlon)=(lon_out(2:nlon+1)-wb_in)*ratlon+1.0 is=fis where (fis < 0.) is=is-1 ie=fie where (fie < 0.) ie=ie-1 facis=float(is+1)-fis facie=fie-float(ie) !----- allocate expanded data arrays ------ ismin=min(minval(is),1) iemax=max(maxval(ie),mdim) allocate (data(ismin:iemax,1:ndim), area(ismin:iemax,1:ndim)) data(1:mdim,1:ndim) = data_in(1:mdim,1:ndim) if (ismin < 1) then data(ismin:0, 1:ndim) = data_in(ismin+mdim:mdim, 1:ndim) endif if (iemax > mdim) then data(mdim+1:iemax, 1:ndim) = data_in(1:iemax-mdim, 1:ndim) endif do j = 1,ndim area(:,j) = area_in(j) enddo if (present(mask_in)) then area(1:mdim,1:ndim) = area(1:mdim,1:ndim)*mask_in(1:mdim,1:ndim) if (ismin < 1) & area(ismin:0, 1:ndim) = area(ismin+mdim:mdim, 1:ndim) if (iemax > mdim) & area(mdim+1:iemax, 1:ndim) = area(1:iemax-mdim, 1:ndim) endif !----------------------------------------------------------------------- ! totarea=0. do j=1,nlat do i=1,nlon iis = is(i) iie = ie(i) jjs = js(j) jje = je(j) ffacis = facis(i) ffacie = facie(i) ffacjs = facjs(j) ffacje = facje(j) kk = 1 522 enlarge = .false. if (mask_out(i,j)) call typemax (data(iis:iie,jjs:jje), & area(iis:iie,jjs:jje), & ffacis,ffacie,ffacjs,ffacje,& ntype,data_out(i,j),enlarge) ! totarea(i,j)=sum(area(iis:iie,jjs:jje)) if(is(i)-kk.le.1.and.ie(i)+kk.ge.mdim.and. & js(j)-kk.le.1.and.je(j)+kk.ge.ndim) goto 523 if(enlarge) then if(is(i)-kk.ge.1) iis = is(i)-kk if(ie(i)+kk.le.mdim) iie = ie(i)+kk if(js(j)-kk.ge.1) jjs = js(j)-kk if(je(j)+kk.le.ndim) jje = je(j)+kk ffacis = 1. ffacie = 1. ffacjs = 1. ffacje = 1. kk=kk+1 go to 522 end if 523 continue enddo enddo ! sumtype = 0. ! do j=1,nlat ! do i=1,nlon ! do kk = 1, ntype ! if(data_out(i,j).eq.kk) sumtype(kk)=sumtype(kk)+totarea(i,j) ! end do ! enddo ! enddo !write(*,*)'ntype=',ntype !do i = 1, ntype ! write(*,*)'type,areatype=',i,sumtype(i) !end do deallocate(data, area) end subroutine regrid_discrete_field_latlon ! ! subroutine create_horiz_interp_new (Interp, nx, ny, wb, sb, dx, dy, lonb_out, latb_out, is_latlon_out) type(horiz_interp_type), intent(inout) :: Interp integer, intent(in) :: nx, ny real, intent(in) :: wb, sb, dx, dy, lonb_out(:,:), latb_out(:,:) logical, intent(in), optional :: is_latlon_out real :: lonb_in(nx+1), latb_in(ny+1), tpi integer :: i, j tpi = 2.*PI ! longitude boundaries do i = 1, nx+1 lonb_in(i) = wb + float(i-1)*dx enddo if (abs(lonb_in(nx+1)-lonb_in(1)-tpi) < epsilon(lonb_in)) & lonb_in(nx+1)=lonb_in(1)+tpi ! latitude boundaries do j = 2, ny latb_in(j) = sb + float(j-1)*dy enddo latb_in(1) = -0.5*PI latb_in(ny+1) = 0.5*PI call horiz_interp_new (Interp, lonb_in, latb_in, & lonb_out, latb_out, is_latlon_out=is_latlon_out) end subroutine create_horiz_interp_new ! ! ! ! Calculation of total area of each land cover type within the territory ! and the determination of the type occupying the biggest area within the ! territory. ! ! ! Calculation of total area of each land cover type within the territory ! and the determination of the type occupying the biggest area within the ! territory. ! ! ! subroutine typemax(data,area,facis,facie,facjs,facje,ntype,ntypemax,enlarge) integer, intent(in) :: data(:,:) real, intent(in) :: facis,facie,facjs,facje real, intent(in) :: area(:,:) integer, intent(in) :: ntype integer, intent(out) :: ntypemax ! number of land cover type to be assigned ! to the entire grid cell; logical, intent(inout) :: enlarge ! says if it's necessary to repeat the ! process for a larger area ! ! --- local vars ----------------------------------------------------- integer itype,id,jd,i,j real, dimension(size(area,1),size(area,2)) :: wt real, dimension(ntype) :: sumarea id=size(area,1); jd=size(area,2) wt = area wt( 1,:)=wt( 1,:)*facis wt(id,:)=wt(id,:)*facie wt(:, 1)=wt(:, 1)*facjs wt(:,jd)=wt(:,jd)*facje ! Calculation of total area of each land cover type within the territory: sumarea = 0. do i = 1, id do j = 1, jd do itype = 1, ntype if(int(data(i,j)).eq.itype) sumarea(itype)=sumarea(itype)+wt(i,j) end do end do end do ! Determination of the type occupying the biggest area within the territory: ntypemax = 0 do itype = 1, ntype if(maxval(sumarea).eq.sumarea(itype).and.sumarea(itype).gt.0.) & ntypemax = itype end do ! Check if the land cover type is defined or not (ntypemax = 0 means ! ocean or undefined type) if(ntypemax.eq.0) enlarge = .true. end subroutine typemax ! ! ! ! Regrids integer index data to any output grid from a uniformly spaced ! grid using a maximum-area voting scheme. ! ! ! Regrids integer index data to any output grid from a uniformly spaced ! grid using a maximum-area voting scheme. ! ! ! subroutine regrid_discrete_field_cube (data_in, wb_in, sb_in, dlon_in, dlat_in, & lon_out, lat_out, ntype, mask_out, data_out, mask_in) !----------------------------------------------------------------------- ! input: ! ----- ! data_in input data; dimensioned by mdim x ndim ! stored from south to north ! wb_in longitude corresponding to western boundary of box i=1 ! sb_in latitude corresponding to southern boundary of box j=1 ! dlon_in x axis grid spacing in degrees of longitude ! dlat_in y axis grid spacing in degrees of latitude ! ! lon_out longitudes of output data at grid box boundaries ! dimensioned by size(data_out,1)+1 ! lat_out latitudes of output data at grid box boundaries ! dimensioned by size(data_out,2)+1 ! ntype number of land cover types specified ! mask_out output mask that specifies where the data are defined ! ! output: ! ------ ! data_out output number of land cover type ! ! optional ! -------- ! mask_in input mask; =0.0 for data points that should not ! be used or have no data; has the same size as data_in ! !----------------------------------------------------------------------- integer, intent(in) :: ntype integer, intent(in) :: data_in(:,:) real, intent(in) :: sb_in, wb_in, dlat_in, dlon_in real, intent(in) :: lon_out(:,:), lat_out(:,:) logical, intent(in) :: mask_out(:,:) integer, intent(out) :: data_out(:,:) real, intent(in), optional :: mask_in(:,:) ! ! --- local vars ----------------------------------------------------- real :: lonb(2,2), latb(2,2) integer :: data1(1,1) logical :: mask1(1,1) integer :: i, j, np integer :: n, mdim, ndim, m360 real :: hpie type(horiz_interp_type) :: Interp ! setup interpolation mdim = size(data_in,1) ndim = size(data_in,2) hpie = pi/2.0 m360 = int(4.*hpie/dlon_in + 0.001) if (mdim < m360) call error_mesg ('regrid_discrete_field_cube', & 'inner dimension for input data is too small.', FATAL) ! ! call create_horiz_interp_new (Interp, mdim, ndim, wb_in, sb_in, dlon_in, dlat_in, lon_out, lat_out) !----------------------------------------------------------------------- ! loop thru number of types ! at each grid box find the type with the largest weight (i.e., area) call regrid_discrete_field_base (Interp, data_in, ntype, mask_out, data_out, mask_in) call horiz_interp_del (Interp) !----------------------------------------------------------------------- !if (count(mask_out.and.data_out==0) > 0) then ! print *, 'Bad discrete points, pe, nbad = ',mpp_pe(),count(mask_out.and.data_out==0) !endif ! search for missing/bad points do j = 1, size(data_out,2) do i = 1, size(data_out,1) if (mask_out(i,j) .and. data_out(i,j) == 0) then ! Initial grid box lonb = lon_out(i:i+1,j:j+1) latb = lat_out(i:i+1,j:j+1) mask1 = mask_out(i,j) np = 0 do ! grow the grid box call expand_cell (lonb,latb,1.5) ! create a new interp type call create_horiz_interp_new (Interp, mdim, ndim, wb_in, sb_in, dlon_in, dlat_in, lonb, latb, is_latlon_out=.false.) call regrid_discrete_field_base (Interp, data_in, ntype, mask1, data1, mask_in) call horiz_interp_del (Interp) np = np+1 if (data1(1,1) /= 0) then data_out(i,j)=data1(1,1) !print *, 'Fixed discrete point: pe, np = ',mpp_pe(),np exit endif if (np > 10) then print *, 'orig lon = ', lon_out(i:i+1,j:j+1) print *, 'orig lat = ', lat_out(i:i+1,j:j+1) print *, 'curr lon = ', lonb print *, 'curr lat = ', latb call error_mesg ('regrid_discrete_field_cube', 'failed to fix discrete point', FATAL) endif enddo endif enddo enddo !----------------------------------------------------------------------- end subroutine regrid_discrete_field_cube ! ! subroutine regrid_discrete_field_base (Intrp, data_in, ntype, mask_out, data_out, mask_in) type(horiz_interp_type), intent(inout) :: Intrp integer, intent(in) :: ntype integer, intent(in) :: data_in(:,:) logical, intent(in) :: mask_out(:,:) integer, intent(out) :: data_out(:,:) real, intent(in), optional :: mask_in(:,:) ! --- local vars --- real, dimension(size(data_in,1), size(data_in,2)) :: wt_in real, dimension(size(data_out,1),size(data_out,2)) :: wt_max, wt_out integer :: n ! loop thru number of types ! at each grid box find the type with the largest weight (i.e., area) data_out = 0 wt_max = 0. do n = 1, ntype where (data_in == n) wt_in = 1. elsewhere wt_in = 0. endwhere if (present(mask_in)) then wt_in = wt_in * mask_in endif call horiz_interp (Intrp, wt_in, wt_out) where (mask_out .and. wt_out > wt_max) wt_max = wt_out data_out = n endwhere enddo end subroutine regrid_discrete_field_base ! ! subroutine set_topo_rough (domain, lonb, latb, topo_stdev) type(domain2d), intent(in) :: domain real, intent(in), dimension(:,:) :: lonb, latb real, intent(out), dimension(:,:) :: topo_stdev logical :: got_stdev integer :: unit, nc character(len=128) :: topo_file_base ! ! get_topog_stdev failed to provide topography variance data. ! ! ! topo_rough_source is set to 'input', but input file name either ! not specified or specified incorrectly, so the program cannot ! find it. ! ! ! specified value of namelist parameter topo_rough_source is invalid; ! valid values are 'computed', 'input' or 'input/computed'. ! if (trim(topo_rough_source) == 'computed') then got_stdev = get_topog_stdev(lonb,latb,topo_stdev) if (.not.got_stdev) & call error_mesg ('land_properties_mod', & 'could not read topography data', FATAL) else if (trim(topo_rough_source)=='input' .or. trim(topo_rough_source)=='input/computed') then nc = len_trim(topo_rough_file) topo_file_base = topo_rough_file(1:nc-3) if (file_exist(topo_rough_file) .or. file_exist(trim(topo_file_base)//'.tile1.nc')) then call set_domain(domain) call read_data (topo_rough_file,topo_rough_var,topo_stdev) else if (file_exist(topo_file_base)) then unit = open_restart_file(topo_file_base,'read') call read_data(unit,topo_stdev) call close_file(unit) else if (trim(topo_rough_source)=='input') then call error_mesg('land_properties_mod', & 'input file for topography standard deviation "'// & trim(topo_rough_file)//'" does not exist', FATAL) else got_stdev = get_topog_stdev(lonb,latb,topo_stdev) if (.not.got_stdev) & call error_mesg ('land_properties_mod', & 'could not read topography data', FATAL) endif endif else call error_mesg('land_properties_mod','"'//trim(topo_rough_source)//& '" is not a valid value for topo_rough_source', FATAL) endif end subroutine set_topo_rough ! ! ! ! Initializes land properties diagnostics. ! ! ! Initializes land properties diagnostics. ! ! ! subroutine init_land_properties_diag (id_lon, id_lat, Time) integer, intent(in) :: id_lon ! ID of land longitude (X) axis integer, intent(in) :: id_lat ! ID of land longitude (X) axis type(time_type), intent(in) :: Time ! current time ! ! ---- local vars ---------------------------------------------------------- integer :: axes(2) axes = (/id_lon, id_lat/) id_alb_max = register_static_field(diag_mod_name, "albedo_max_no_snow", axes, & "max. snow-free land albedo","dimensionless", missing_value=-1.0) id_alb_min = register_static_field(diag_mod_name, "albedo_min_no_snow", axes, & "min. snow-free land albedo","dimensionless", missing_value=-1.0) id_ground_type = register_static_field(diag_mod_name,"ground_type", axes, & "land surface ground type", "dimensionless", missing_value = -1.0) end subroutine init_land_properties_diag ! ! ! ! Maximum snow-free land albedo ! ! ! Minimum snow-free land albedo ! ! ! Land surface cover type ! ! ! Land surface ground type ! ! ! ! ! Sends static fields to diagnostic output. ! ! ! Sends static fields to diagnostic output. ! ! ! subroutine diag_static ( Time,mask ) type(time_type), intent(in) :: Time logical, intent(in) :: mask(:,:,:) ! ! ---- local vars ---------------------------------------------------------- logical :: used real :: dummy (size(mask,1), size(mask,2)) if (id_alb_max > 0) & used = send_data ( id_alb_max, albedo_max_no_snow(:,:,1), Time, mask=mask(:,:,1) ) if (id_alb_min > 0) & used = send_data ( id_alb_min, albedo_min_no_snow(:,:,1), Time, mask=mask(:,:,1) ) if (id_ground_type > 0) then dummy = ground_type(:,:,1) used = send_data (id_ground_type, dummy, Time, mask=mask(:,:,1)) endif end subroutine diag_static ! !============================================================================= ! ! ! Deallocates the land property data. ! ! ! Deallocates the land property data. ! ! ! subroutine land_properties_end() ! module_is_initialized =.FALSE. deallocate(ground_type, albedo_min_no_snow, albedo_max_no_snow) end subroutine land_properties_end ! end module land_properties_mod