! ============================================================================ ! top-level core of the Land Dynamics (LaD) model code ! ============================================================================ module land_model_mod #include "shared/debug.inc" #include "shared/concat.inc" use time_manager_mod, only : time_type, get_time, increment_time, time_type_to_real, & operator(+) use mpp_domains_mod, only : domain2d, mpp_get_ntile_count #ifdef INTERNAL_FILE_NML use mpp_mod, only: input_nml_file #else use fms_mod, only: open_namelist_file #endif use mpp_mod, only : mpp_max, mpp_sum use fms_mod, only : write_version_number, error_mesg, FATAL, WARNING, NOTE, mpp_pe, & mpp_root_pe, file_exist, check_nml_error, close_file, & stdlog, stderr, mpp_clock_id, mpp_clock_begin, mpp_clock_end, string, & stdout, CLOCK_FLAG_DEFAULT, CLOCK_COMPONENT, CLOCK_ROUTINE use field_manager_mod, only : MODEL_LAND use data_override_mod, only : data_override use diag_manager_mod, only : diag_axis_init, register_static_field, & register_diag_field, send_data use constants_mod, only : radius, hlf, hlv, hls, tfreeze, pi, rdgas, rvgas, cp_air, & stefan use astronomy_mod, only : diurnal_solar use sphum_mod, only : qscomp use tracer_manager_mod, only : NO_TRACER use land_constants_mod, only : NBANDS, BAND_VIS, BAND_NIR, mol_air, mol_C, mol_co2 use glacier_mod, only : read_glac_namelist, glac_init, glac_end, glac_get_sfc_temp, & glac_radiation, glac_diffusion, glac_step_1, glac_step_2, save_glac_restart use lake_mod, only : read_lake_namelist, lake_init, lake_end, lake_get_sfc_temp, & lake_radiation, lake_diffusion, lake_step_1, lake_step_2, save_lake_restart use soil_mod, only : read_soil_namelist, soil_init, soil_end, soil_get_sfc_temp, & soil_radiation, soil_diffusion, soil_step_1, soil_step_2, soil_step_3, & save_soil_restart use snow_mod, only : read_snow_namelist, snow_init, snow_end, snow_get_sfc_temp, & snow_radiation, snow_diffusion, snow_get_depth_area, snow_step_1, snow_step_2, & save_snow_restart use vegetation_mod, only : read_vegn_namelist, vegn_init, vegn_end, vegn_get_cover, & vegn_radiation, vegn_diffusion, vegn_step_1, vegn_step_2, vegn_step_3, & update_vegn_slow, save_vegn_restart use cana_tile_mod, only : canopy_air_mass, canopy_air_mass_for_tracers, cana_tile_heat use canopy_air_mod, only : read_cana_namelist, cana_init, cana_end, cana_state,& cana_step_1, cana_step_2, cana_radiation, cana_roughness, & save_cana_restart use river_mod, only : river_init, river_end, update_river, river_stock_pe, & save_river_restart use topo_rough_mod, only : topo_rough_init, topo_rough_end, update_topo_rough use soil_tile_mod, only : soil_cover_cold_start, soil_tile_stock_pe, & soil_tile_heat use vegn_tile_mod, only : vegn_cover_cold_start, vegn_data_rs_min, & update_derived_vegn_data, vegn_tile_stock_pe, & vegn_tile_heat use lake_tile_mod, only : lake_cover_cold_start, lake_tile_stock_pe, & lake_tile_heat use glac_tile_mod, only : glac_pars_type, glac_cover_cold_start, & glac_tile_stock_pe, glac_tile_heat use snow_tile_mod, only : snow_tile_stock_pe, snow_tile_heat use land_numerics_mod, only : ludcmp, lubksb, nearest, & horiz_remap_type, horiz_remap_new, horiz_remap, horiz_remap_del, & horiz_remap_print use land_io_mod, only : read_land_io_namelist, input_buf_size use land_tile_mod, only : land_tile_type, land_tile_list_type, & land_tile_enum_type, new_land_tile, insert, nitems, & first_elmt, tail_elmt, next_elmt, current_tile, operator(/=), & get_elmt_indices, get_tile_tags, land_tile_carbon use land_data_mod, only : land_data_type, atmos_land_boundary_type, & land_state_type, land_data_init, land_data_end, lnd, & dealloc_land2cplr, realloc_land2cplr, & dealloc_cplr2land, realloc_cplr2land, & land_data_type_chksum, atm_lnd_bnd_type_chksum use nf_utils_mod, only : nfu_inq_var, nfu_inq_dim, nfu_get_var use land_utils_mod, only : put_to_tiles_r0d_fptr use land_tile_io_mod, only : print_netcdf_error, create_tile_out_file, & read_tile_data_r0d_fptr, write_tile_data_r0d_fptr, & write_tile_data_i0d_fptr, get_input_restart_name use land_tile_diag_mod, only : tile_diag_init, tile_diag_end, & register_tiled_diag_field, send_tile_data, dump_tile_diag_fields, & add_tiled_diag_field_alias, & OP_AVERAGE, OP_SUM use land_debug_mod, only : land_debug_init, land_debug_end, set_current_point, & is_watch_point, get_watch_point, check_temp_range, current_face, & get_current_point use static_vegn_mod, only : write_static_vegn use land_transitions_mod, only : & land_transitions_init, land_transitions_end, land_transitions, & save_land_transitions_restart use stock_constants_mod, only: ISTOCK_WATER, ISTOCK_HEAT, ISTOCK_SALT implicit none private ! ==== public interfaces ===================================================== public land_model_init ! initialize the land model public land_model_end ! finish land model calculations public land_model_restart ! saves the land model restart(s) public update_land_model_fast ! time-step integration public update_land_model_slow ! time-step integration public atmos_land_boundary_type ! data from coupler to land public land_data_type ! data from land to coupler public land_data_type_chksum ! routine to print checksums for land_data_type public atm_lnd_bnd_type_chksum ! routine to print checksums for atmos_land_boundary_type public :: Lnd_stock_pe ! return stocks of conservative quantities ! ==== end of public interfaces ============================================== ! ==== module constants ====================================================== character(len=*), parameter :: & module_name = 'land', & version = '$Id: land_model.F90,v 20.0 2013/12/13 23:29:26 fms Exp $', & tagname = '$Name: tikal $' ! ==== module variables ====================================================== ! ---- namelist -------------------------------------------------------------- logical :: use_old_conservation_equations = .false. logical :: lm2 = .false. logical :: do_age = .false. logical :: give_stock_details = .false. logical :: use_tfreeze_in_grnd_latent = .false. logical :: use_atmos_T_for_precip_T = .false. logical :: use_atmos_T_for_evap_T = .false. real :: cpw = 1952. ! specific heat of water vapor at constant pressure real :: clw = 4218. ! specific heat of water (liquid) real :: csw = 2106. ! specific heat of water (ice) real :: min_sum_lake_frac = 1.e-8 real :: min_frac = 0.0 ! minimum fraction of soil, lake, and glacier that is not discarded on cold start real :: gfrac_tol = 1.e-6 real :: discharge_tol = -1.e20 real :: con_fac_large = 1.e6 real :: con_fac_small = 1.e-6 integer :: num_c = 0 real :: tau_snow_T_adj = -1.0 ! time scale of snow temperature adjustment ! for the snow-free surface (s); negative means no adjustment logical :: prohibit_negative_canopy_water = .FALSE. ! if true, then in case of negative canopy ! water the evaporation is fixed and the equations are re-solved. ! Default retrievs old behavior. character(16) :: nearest_point_search = 'global' ! specifies where to look for ! nearest points for missing data, "global" or "face" logical :: print_remapping = .FALSE. ! if true, full land cover remapping ! information is printed on the cold start integer :: layout(2) = (/0,0/) integer :: io_layout(2) = (/0,0/) namelist /land_model_nml/ use_old_conservation_equations, & lm2, do_age, give_stock_details, & use_tfreeze_in_grnd_latent, & use_atmos_T_for_precip_T, & use_atmos_T_for_evap_T, & cpw, clw, csw, min_sum_lake_frac, min_frac, & gfrac_tol, discharge_tol, & con_fac_large, con_fac_small, num_c, & tau_snow_T_adj, prohibit_negative_canopy_water, & nearest_point_search, print_remapping, & layout, io_layout ! ---- end of namelist ------------------------------------------------------- logical :: module_is_initialized = .FALSE. logical :: stock_warning_issued = .FALSE. logical :: update_cana_co2 ! if false, cana_co2 is not updated during the model run. character(len=256) :: grid_spec_file="INPUT/grid_spec.nc" real :: delta_time ! duration of main land time step (s) integer :: num_species integer :: num_phys = 2 real, allocatable :: frac (:,:) ! fraction of land in cells logical, allocatable :: river_land_mask(:,:), missing_rivers(:,:) real, allocatable :: no_riv(:,:) ! ---- diag field IDs -------------------------------------------------------- integer :: & ! COLUMN VEGN SNOW GLAC/LAKE/SOIL CANOPY-AIR RIVER id_VWS, id_VWSc, & id_LWS, id_LWSv, id_LWSs, id_LWSg, & id_FWS, id_FWSv, id_FWSs, id_FWSg, & id_HS, id_HSv, id_HSs, id_HSg, id_HSc, & id_precip, & id_hprec, & id_lprec, id_lprecv, id_lprecs, id_lprecg, & id_hlprec, id_hlprecv, id_hlprecs, id_hlprecg, & id_fprec, id_fprecv, id_fprecs, & id_hfprec, id_hfprecv, id_hfprecs, & id_evap, & id_hevap, & id_levap, id_levapv, id_levaps, id_levapg, & id_hlevap, id_hlevapv, id_hlevaps, id_hlevapg, & id_fevap, id_fevapv, id_fevaps, id_fevapg, & id_hfevap, id_hfevapv, id_hfevaps, id_hfevapg, & id_runf, & id_hrunf, & id_lrunf, id_lrunfs, id_lrunfg, & id_hlrunf, id_hlrunfs, id_hlrunfg, & id_frunf, id_frunfs, & id_hfrunf, id_hfrunfs, & id_melt, id_meltv, id_melts, id_meltg, & id_fsw, id_fswv, id_fsws, id_fswg, & id_flw, id_flwv, id_flws, id_flwg, & id_sens, id_sensv, id_senss, id_sensg, & ! id_e_res_1, id_e_res_2, id_cd_m, id_cd_t, & id_cellarea, id_landarea, id_landfrac, id_no_riv, & id_geolon_t, id_geolat_t, & id_frac, id_area, id_ntiles, & id_dis_liq, id_dis_ice, id_dis_heat, id_dis_sink, & id_z0m, id_z0s, id_con_g_h, & id_transp, id_wroff, id_sroff, & id_htransp, id_huptake, id_hroff, id_gsnow, id_gequil, & id_grnd_flux, & id_soil_water_supply, id_levapg_max, & id_water, id_snow, & id_Trad, id_Tca, id_qca, id_qco2, id_qco2_dvmr, & id_swdn_dir, id_swdn_dif, id_swup_dir, id_swup_dif, id_lwdn, & id_fco2, & id_vegn_cover, id_cosz, & id_albedo_dir, id_albedo_dif, & id_vegn_refl_dir, id_vegn_refl_dif, id_vegn_refl_lw, & id_vegn_tran_dir, id_vegn_tran_dif, id_vegn_tran_lw, & id_vegn_sctr_dir, & id_subs_refl_dir, id_subs_refl_dif, id_subs_emis, id_grnd_T, id_total_C ! ---- global clock IDs integer :: landClock, landFastClock, landSlowClock ! ==== NetCDF declarations =================================================== include 'netcdf.inc' #define __NF_ASRT__(x) call print_netcdf_error((x),__FILE__,__LINE__) contains ! ============================================================================ subroutine land_model_init & (cplr2land, land2cplr, time_init, time, dt_fast, dt_slow) ! initialize land model using grid description file as an input. This routine ! reads land grid boundaries and area of land from a grid description file ! NOTES: theoretically, the grid description file can specify any regular ! rectangular grid for land, not just lon/lat grid. Therefore the variables ! "xbl" and "ybl" in NetCDF grid spec file are not necessarily lon and lat ! boundaries of the grid. ! However, at this time the module land_properties assumes that grid _is_ ! lon/lat and therefore the entire module also have to assume that the land ! grid is lon/lat. ! lon/lat grid is also assumed for the diagnostics, but this is probably not ! so critical. type(atmos_land_boundary_type), intent(inout) :: cplr2land ! boundary data type(land_data_type) , intent(inout) :: land2cplr ! boundary data type(time_type), intent(in) :: time_init ! initial time of simulation (?) type(time_type), intent(in) :: time ! current time type(time_type), intent(in) :: dt_fast ! fast time step type(time_type), intent(in) :: dt_slow ! slow time step ! ---- local vars ---------------------------------------------------------- integer :: ncid, varid integer :: unit, ierr, io integer :: id_lon, id_lat, id_band ! IDs of land diagnostic axes logical :: used ! return value of send_data diagnostics routine integer :: i,j,k type(land_tile_type), pointer :: tile type(land_tile_enum_type) :: ce, te character(len=256) :: restart_file_name logical :: restart_exists ! IDs of local clocks integer :: landInitClock module_is_initialized = .TRUE. ! [1] print out version number call write_version_number(version, tagname) ! initialize land model clocks landClock = mpp_clock_id('Land' ,CLOCK_FLAG_DEFAULT,CLOCK_COMPONENT) landFastClock = mpp_clock_id('Update-Land-Fast' ,CLOCK_FLAG_DEFAULT,CLOCK_ROUTINE) landSlowClock = mpp_clock_id('Update-Land-Slow' ,CLOCK_FLAG_DEFAULT,CLOCK_ROUTINE) landInitClock = mpp_clock_id('Land init' ,CLOCK_FLAG_DEFAULT,CLOCK_ROUTINE) call mpp_clock_begin(landInitClock) ! [ ] initialize land debug output call land_debug_init() ! [ ] initialize tile-specific diagnostics internals call tile_diag_init() ! [2] read namelists ! [2.1] read land model namelist #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=land_model_nml, iostat=io) ierr = check_nml_error(io, 'land_model_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr = 1; do while (ierr /= 0) read (unit, nml=land_model_nml, iostat=io, end=10) ierr = check_nml_error (io, 'land_model_nml') enddo 10 continue call close_file (unit) endif #endif if (mpp_pe() == mpp_root_pe()) then unit = stdlog() write (unit, nml=land_model_nml) call close_file (unit) endif ! [2.2] read sub-model namelists: then need to be read before initialization ! because they can affect the way cover and tiling is initialized on cold start. call read_land_io_namelist() call read_soil_namelist() call read_vegn_namelist() call read_lake_namelist() call read_glac_namelist() call read_snow_namelist() call read_cana_namelist() ! [ ] initialize land state data, including grid geometry and processor decomposition call land_data_init(layout, io_layout, time, dt_fast, dt_slow) delta_time = time_type_to_real(lnd%dt_fast) ! store in a module variable for convenience ! calculate land fraction allocate(frac(lnd%is:lnd%ie,lnd%js:lnd%je)) frac = lnd%area/lnd%cellarea ! [5] initialize tiling call get_input_restart_name('INPUT/land.res.nc',restart_exists,restart_file_name) if(restart_exists) then call error_mesg('land_model_init',& 'reading NetCDF restart "'//trim(restart_file_name)//'"',& NOTE) ! read map of tiles -- retrieve information from call land_cover_warm_start(restart_file_name,lnd) ! initialize land model data __NF_ASRT__(nf_open(restart_file_name,NF_NOWRITE,ncid)) if (nf_inq_varid(ncid,'lwup',varid)==NF_NOERR) & call read_tile_data_r0d_fptr(ncid,'lwup',land_lwup_ptr) if (nf_inq_varid(ncid,'e_res_1',varid)==NF_NOERR) & call read_tile_data_r0d_fptr(ncid,'e_res_1',land_e_res_1_ptr) if (nf_inq_varid(ncid,'e_res_2',varid)==NF_NOERR) & call read_tile_data_r0d_fptr(ncid,'e_res_2',land_e_res_2_ptr) __NF_ASRT__(nf_close(ncid)) else ! initialize map of tiles -- construct it by combining tiles ! from component models call error_mesg('land_model_init',& 'cold-starting land cover map',& NOTE) call land_cover_cold_start(lnd) endif ! [6] initialize land model diagnostics -- must be before *_data_init so that ! *_data_init can write static fields if necessary call land_diag_init( lnd%coord_glonb, lnd%coord_glatb, lnd%coord_glon, lnd%coord_glat, time, lnd%domain, & id_lon, id_lat, id_band ) ! set the land diagnostic axes ids for the flux exchange land2cplr%axes = (/id_lon,id_lat/) ! send some static diagnostic fields to output if ( id_cellarea > 0 ) used = send_data ( id_cellarea, lnd%cellarea, lnd%time ) if ( id_landarea > 0 ) used = send_data ( id_landarea, lnd%area, lnd%time ) if ( id_landfrac > 0 ) used = send_data ( id_landfrac, frac, lnd%time ) if ( id_geolon_t > 0 ) used = send_data ( id_geolon_t, lnd%lon*180.0/PI, lnd%time ) if ( id_geolat_t > 0 ) used = send_data ( id_geolat_t, lnd%lat*180.0/PI, lnd%time ) ! [7] initialize individual sub-models num_species = num_phys + num_c if (do_age) num_species = num_species + 1 call soil_init ( id_lon, id_lat, id_band ) call vegn_init ( id_lon, id_lat, id_band ) call lake_init ( id_lon, id_lat ) call glac_init ( id_lon, id_lat ) call snow_init ( id_lon, id_lat ) call cana_init ( id_lon, id_lat ) call topo_rough_init( lnd%time, lnd%lonb, lnd%latb, & lnd%domain, id_lon, id_lat) allocate (river_land_mask(lnd%is:lnd%ie,lnd%js:lnd%je)) allocate ( missing_rivers(lnd%is:lnd%ie,lnd%js:lnd%je)) allocate ( no_riv (lnd%is:lnd%ie,lnd%js:lnd%je)) call river_init( lnd%lon, lnd%lat, & lnd%time, lnd%dt_fast, lnd%domain, & frac, & id_lon, id_lat, & river_land_mask ) missing_rivers = frac.gt.0. .and. .not.river_land_mask no_riv = 0. where (missing_rivers) no_riv = 1. if ( id_no_riv > 0 ) used = send_data( id_no_riv, no_riv, lnd%time ) call land_transitions_init (id_lon, id_lat) ! [8] initialize boundary data ! [8.1] allocate storage for the boundary data call realloc_land2cplr ( land2cplr ) call realloc_cplr2land ( cplr2land ) ! [8.2] set the land mask to FALSE everywhere -- update_land_bc_fast ! will set it to true where necessary land2cplr%mask = .FALSE. land2cplr%tile_size = 0.0 ! [8.3] get the current state of the land boundary for the coupler ce = first_elmt(lnd%tile_map, & is=lbound(cplr2land%t_flux,1), & js=lbound(cplr2land%t_flux,2) ) te = tail_elmt(lnd%tile_map) do while(ce /= te) ! calculate indices of the current tile in the input arrays; ! assume all the cplr2land components have the same lbounds call get_elmt_indices(ce,i,j,k) ! set this point coordinates as current for debug output call set_current_point(i,j,k) ! get pointer to current tile tile => current_tile(ce) ! advance enumerator to the next tile ce=next_elmt(ce) call update_land_bc_fast (tile, i,j,k, land2cplr, is_init=.true.) enddo ! [8.4] update topographic roughness scaling call update_land_bc_slow( land2cplr ) ! mask error checking do j=lnd%js,lnd%je do i=lnd%is,lnd%ie if(frac(i,j)>0.neqv.ANY(land2cplr%mask(i,j,:))) then call error_mesg('land_model_init','masks are not equal',FATAL) endif enddo enddo ! [9] check the properties of co2 exchange with the atmosphere and set appropriate ! flags if (canopy_air_mass_for_tracers==0.and.lnd%ico2==NO_TRACER) then call error_mesg('land_model_init', & 'canopy_air_mass_for_tracers is set to zero, and CO2 exchange with the atmosphere is not set up: '// & 'canopy air CO2 concentration will not be updated',NOTE) update_cana_co2 = .FALSE. else update_cana_co2 = .TRUE. end if call mpp_clock_end(landInitClock) end subroutine land_model_init ! ============================================================================ subroutine land_model_end (cplr2land, land2cplr) type(atmos_land_boundary_type), intent(inout) :: cplr2land type(land_data_type) , intent(inout) :: land2cplr ! ---- local vars module_is_initialized = .FALSE. call error_mesg('land_model_end','writing NetCDF restart',NOTE) call land_model_restart() ! we still want to call the *_end procedures for component models, even ! if the number of tiles in this domain is zero, in case they are doing ! something else besides saving the restart, of if they want to save ! restart anyway call land_transitions_end() call glac_end () call lake_end () call soil_end () call snow_end () call vegn_end () call cana_end () call topo_rough_end() call river_end() deallocate(frac) deallocate(river_land_mask, missing_rivers, no_riv) call dealloc_land2cplr(land2cplr, dealloc_discharges=.TRUE.) call dealloc_cplr2land(cplr2land) call tile_diag_end() ! deallocate tiles call land_data_end() ! finish up the land debugging diagnostics call land_debug_end end subroutine land_model_end ! ============================================================================ ! write land model restarts subroutine land_model_restart(timestamp) character(*), intent(in), optional :: timestamp ! timestamp to add to the file name ! ---- local vars integer :: tile_dim_length ! length of tile dimension in output files ! global max of number of tiles per gridcell integer :: i,j,k integer :: unit ! netcdf id of the restart file character(256) :: timestamp_ ! [1] count all land tiles and determine the length of tile dimension ! sufficient for the current domain tile_dim_length = 0 do j = lnd%js, lnd%je do i = lnd%is, lnd%ie k = nitems(lnd%tile_map(i,j)) tile_dim_length = max(tile_dim_length,k) enddo enddo ! [2] calculate the tile dimension length by taking the max across all domains call mpp_max(tile_dim_length) if (tile_dim_length==0) then call error_mesg('land_model_restart',& 'No land points exist (tile_dim_length=0), therefore no land restarts will be saved',& WARNING) return endif ! [3] create tile output file timestamp_='' if (present(timestamp)) then if(trim(timestamp)/='') timestamp_=trim(timestamp)//'.' endif call create_tile_out_file(unit,'RESTART/'//trim(timestamp_)//'land.res.nc', & lnd%coord_glon, lnd%coord_glat, land_tile_exists, tile_dim_length) ! [4] write data fields ! write fractions and tile tags call write_tile_data_r0d_fptr(unit,'frac',land_frac_ptr,'fractional area of tile') call write_tile_data_i0d_fptr(unit,'glac',glac_tag_ptr,'tag of glacier tiles') call write_tile_data_i0d_fptr(unit,'lake',lake_tag_ptr,'tag of lake tiles') call write_tile_data_i0d_fptr(unit,'soil',soil_tag_ptr,'tag of soil tiles') call write_tile_data_i0d_fptr(unit,'vegn',vegn_tag_ptr,'tag of vegetation tiles') ! write the upward long-wave flux call write_tile_data_r0d_fptr(unit,'lwup',land_lwup_ptr,'upward long-wave flux') ! write energy residuals call write_tile_data_r0d_fptr(unit,'e_res_1',land_e_res_1_ptr,& 'energy residual in canopy air energy balance equation', 'W/m2') call write_tile_data_r0d_fptr(unit,'e_res_2',land_e_res_2_ptr,& 'energy residual in canopy energy balance equation', 'W/m2') ! [5] close file __NF_ASRT__(nf_close(unit)) ! [6] save component models' restarts call save_land_transitions_restart(timestamp_) call save_glac_restart(tile_dim_length,timestamp_) call save_lake_restart(tile_dim_length,timestamp_) call save_soil_restart(tile_dim_length,timestamp_) call save_snow_restart(tile_dim_length,timestamp_) call save_vegn_restart(tile_dim_length,timestamp_) call save_cana_restart(tile_dim_length,timestamp_) call save_river_restart(timestamp_) end subroutine land_model_restart ! ============================================================================ subroutine land_cover_cold_start(lnd) type(land_state_type), intent(inout) :: lnd ! ---- local vars real, dimension(:,:,:), pointer :: & glac, soil, lake, vegn ! arrays of fractions for respective sub-models logical, dimension(lnd%ie-lnd%is+1,lnd%je-lnd%js+1) :: & land_mask, valid_data, invalid_data integer :: iwatch,jwatch,kwatch,face integer :: i,j integer :: ps,pe ! boundaries of PE list for remapping type(horiz_remap_type) :: map ! calculate the global land mask land_mask = lnd%area > 0 ! get the global maps of fractional covers for each of the sub-models glac=>glac_cover_cold_start(land_mask,lnd%lonb,lnd%latb) lake=>lake_cover_cold_start(land_mask,lnd%lonb,lnd%latb,lnd%domain) soil=>soil_cover_cold_start(land_mask,lnd%lonb,lnd%latb) vegn=>vegn_cover_cold_start(land_mask,lnd%lonb,lnd%latb) ! remove any input lake fraction in coastal cells where (frac.lt. 1.-gfrac_tol) lake(:,:,1) = 0. ! NOTE that the lake area in the coastal cells can be set to non-zero ! again by the "ground fraction reconciliation code" below. Strictly ! speaking the above line of code should be replaced with the section ! commented out with "!-zero" below, but we preserve the old way to avoid ! backward incompatibility with older runs. This needs updating in the ! future when the decision about what to do with lakes in coastal cells is ! made. ! reconcile ground fractions with the land mask within compute domain valid_data = land_mask.and.(sum(glac,3)+sum(lake,3)+sum(soil,3)>0) invalid_data = land_mask.and..not.valid_data call get_watch_point(iwatch,jwatch,kwatch,face) if (face==lnd%face.and.(lnd%is<=iwatch.and.iwatch<=lnd%ie).and.(lnd%js<=jwatch.and.jwatch<=lnd%je)) then write(*,*)'###### land_cover_cold_start: input data #####' write(*,'(99(a,i4.2,x))')'iwatch=',iwatch,'jwatch=',jwatch,'face=',lnd%face write(*,'(99(a,g23.16,x))')'lon=',lnd%lon(iwatch,jwatch)*180/PI,'lat=',lnd%lat(iwatch,jwatch)*180/PI ! calculate local compute domain indices; we assume glac,lake,soil,vegn all ! have the same lbounds i = iwatch-lnd%is+lbound(glac,1); j = jwatch-lnd%js+lbound(glac,2) __DEBUG2__(lnd%is,lnd%js) write(*,'(a,99(a,i4.2,x))')'local indices:','i=',i,'j=',j __DEBUG3__(frac(iwatch,jwatch),land_mask(i,j),valid_data(i,j)) __DEBUG1__(glac(i,j,:)) __DEBUG1__(lake(i,j,:)) __DEBUG1__(soil(i,j,:)) __DEBUG1__(vegn(i,j,:)) endif if (trim(nearest_point_search)=='global') then ps=0 ; pe=size(lnd%pelist)-1 else if (trim(nearest_point_search)=='face') then ! this assumes that the number of PEs is divisible by the number of ! mosaic faces. lnd%pelist starts with 0 ps = size(lnd%pelist)/lnd%nfaces*(lnd%face-1) pe = size(lnd%pelist)/lnd%nfaces*lnd%face - 1 else call error_mesg('land_cover_cold_start',& 'option nearest_point_search="'//trim(nearest_point_search)//& '" is illegal, use "global" or "face"',& FATAL) endif call horiz_remap_new(invalid_data,valid_data,lnd%lon,lnd%lat,lnd%domain,& lnd%pelist(ps:pe),map) if (print_remapping) call horiz_remap_print(map,'land cover remap:') call horiz_remap(map,lnd%domain,glac) call horiz_remap(map,lnd%domain,lake) call horiz_remap(map,lnd%domain,soil) call horiz_remap_del(map) !-zero ! remove any input lake fraction in coastal cells !-zero do j = lnd%js,lnd%je !-zero do i = lnd%is,lnd%ie !-zero call set_current_point(i,j,1) !-zero if (frac(i,j) < 1-gfrac_tol) then !-zero lake(i,j,:) = 0.0 !-zero if(is_watch_point())then !-zero write(*,*)'###### land_cover_cold_start: lake fraction is set to zero #####' !-zero endif !-zero endif !-zero enddo !-zero enddo ! reconcile vegetation fractions with the land mask within compute domain valid_data = sum(vegn,3) > 0 invalid_data = .FALSE. do j = 1,size(land_mask,2) do i = 1,size(land_mask,1) if(.not.land_mask(i,j)) cycle ! skip ocean points if(valid_data(i,j)) cycle ! don't need to do anything with valid points if(sum(glac(i,j,:))+sum(lake(i,j,:))>=1) & cycle ! skip points fully covered by glaciers or lakes invalid_data(i,j)=.TRUE. enddo enddo call horiz_remap_new(invalid_data,valid_data,lnd%lon,lnd%lat,lnd%domain,& lnd%pelist(ps:pe),map) if (print_remapping) call horiz_remap_print(map,'vegetation cover remap:') call horiz_remap(map,lnd%domain,vegn) call horiz_remap_del(map) ! create tiles do j = 1,size(land_mask,2) do i = 1,size(land_mask,1) if(.not.land_mask(i,j)) cycle ! skip ocean points call set_current_point(i+lnd%is-1,j+lnd%js-1,1) call land_cover_cold_start_0d & (lnd%tile_map(i+lnd%is-1,j+lnd%js-1),glac(i,j,:),lake(i,j,:),soil(i,j,:),vegn(i,j,:)) if(nitems(lnd%tile_map(i+lnd%is-1,j+lnd%js-1))==0) then call error_mesg('land_cover_cold_start',& 'No tiles were created for a valid land point at i='& //trim(string(lnd%is+i-1))//' j='//trim(string(lnd%js+j-1))//' face='//trim(string(lnd%face)), FATAL) endif enddo enddo deallocate(glac,lake,soil,vegn) end subroutine land_cover_cold_start ! ============================================================================ subroutine land_cover_cold_start_0d (set,glac0,lake0,soil0,vegn0) type(land_tile_list_type), intent(inout) :: set real, dimension(:) , intent(in) :: & glac0,lake0,soil0,vegn0 ! fractions of area ! ---- local vars real :: glac(size(glac0(:))), lake(size(lake0(:))), & soil(size(soil0(:))), vegn(size(vegn0(:))) type(land_tile_type), pointer :: tile integer :: i,j,k real :: factor ! normalizing factor for the tile areas real :: frac type(land_tile_enum_type) :: first_non_vegn ! position of first non-vegetated tile in the list glac = glac0; lake = lake0; soil = soil0; vegn = vegn0 if (sum(glac)>1) & glac=glac/sum(glac) if (sum(lake)+sum(glac)>1)& lake = lake*(1-sum(glac))/sum(lake) if (sum(lake)1)& soil = soil*(1-sum(lake)-sum(glac))/sum(soil) ! make sure that the sum of the fractions of the soil, lake, and glaciers are ! either one or zero factor = sum(soil)+sum(glac)+sum(lake) if(factor>0)then glac = glac/factor lake = lake/factor soil = soil/factor endif ! remove soil/glac/lake fractions that are too small if (min_frac>0) then where (glac0)then glac = glac/factor lake = lake/factor soil = soil/factor endif endif if(is_watch_point()) then write(*,*)'#### land_cover_cold_start_0d input data ####' __DEBUG1__(glac0) __DEBUG1__(lake0) __DEBUG1__(soil0) __DEBUG1__(vegn0) __DEBUG1__(factor) write(*,*)'#### land_cover_cold_start_0d renormlaized fractions ####' __DEBUG1__(glac) __DEBUG1__(lake) __DEBUG1__(soil) __DEBUG1__(vegn) endif do i = 1,size(glac) if (glac(i)>0) then tile => new_land_tile(frac=glac(i),glac=i) call insert(tile,set) if(is_watch_point()) then write(*,*)'created glac tile: frac=',glac(i),' tag=',i endif endif enddo do i = 1,size(lake) if (lake(i)>0) then tile => new_land_tile(frac=lake(i),lake=i) call insert(tile,set) if(is_watch_point()) then write(*,*)'created lake tile: frac=',lake(i),' tag=',i endif endif enddo factor = sum(soil)*sum(vegn) if (factor/=0) factor = 1/factor factor = factor*(1-sum(glac)-sum(lake)) ! vegetation tiles, if any, are inserted in front of non-vegetated tiles; ! this really doesn't matter except for the static vegetation override ! case with the data saved by lm3v -- there the vegetation tiles are ! in front, so it works more consistently where lad2 has more than ! one tile (e.g. glac/soil or lake/soil), if lad2 vegetation tiles are ! also in front of the list. first_non_vegn=first_elmt(set) do i = 1,size(soil) do j = 1,size(vegn) frac = soil(i)*vegn(j)*factor if(frac>0) then tile => new_land_tile(frac=frac,soil=i,vegn=j) call insert(tile,first_non_vegn) if(is_watch_point()) then write(*,*)'created soil tile: frac=', frac, ' soil tag=',i, ' veg tag=',j endif endif enddo enddo end subroutine land_cover_cold_start_0d ! ============================================================================ ! reads the land restart file and restores the tiling structure from this file subroutine land_cover_warm_start ( restart_file_name, lnd ) character(len=*), intent(in) :: restart_file_name type(land_state_type), intent(inout) :: lnd ! ---- local vars integer, allocatable :: idx(:) ! compressed tile index integer, allocatable :: glac(:), lake(:), soil(:), snow(:), cana(:), vegn(:) ! tile tags real, allocatable :: frac(:) ! fraction of land covered by tile integer :: ncid ! unit number of the input file integer :: ntiles ! total number of land tiles in the input file integer :: bufsize ! size of the input buffer integer :: dimids(1) ! id of tile dimension character(NF_MAX_NAME) :: tile_dim_name ! name of the tile dimension and respective variable integer :: i,j,k,it type(land_tile_type), pointer :: tile; integer :: start, count ! slab for reading ! netcdf variable IDs integer :: id_idx, id_frac, id_glac, id_lake, id_soil, id_vegn __NF_ASRT__(nf_open(restart_file_name,NF_NOWRITE,ncid)) ! allocate the input data __NF_ASRT__(nfu_inq_var(ncid,'frac',id=id_frac,varsize=ntiles,dimids=dimids)) ! allocate input buffers for compression index and the variable bufsize=min(input_buf_size,ntiles) allocate(idx (bufsize), glac(bufsize), lake(bufsize), soil(bufsize), & snow(bufsize), cana(bufsize), vegn(bufsize), frac(bufsize) ) ! get the name of the fist (and only) dimension of the variable 'frac' -- this ! is supposed to be the compressed dimension, and associated variable will ! hold the compressed indices __NF_ASRT__(nfu_inq_dim(ncid,dimids(1),name=tile_dim_name)) __NF_ASRT__(nfu_inq_var(ncid,tile_dim_name,id=id_idx)) ! get the IDs of the varables to read __NF_ASRT__(nfu_inq_var(ncid,'glac',id=id_glac)) __NF_ASRT__(nfu_inq_var(ncid,'lake',id=id_lake)) __NF_ASRT__(nfu_inq_var(ncid,'soil',id=id_soil)) __NF_ASRT__(nfu_inq_var(ncid,'vegn',id=id_vegn)) do start = 1,ntiles,bufsize count = min(bufsize,ntiles-start+1) ! read the compressed tile indices __NF_ASRT__(nf_get_vara_int(ncid,id_idx,(/start/),(/count/),idx)) ! read input data -- fractions and tags __NF_ASRT__(nf_get_vara_double(ncid,id_frac,(/start/),(/count/),frac)) __NF_ASRT__(nf_get_vara_int(ncid,id_glac,(/start/),(/count/),glac)) __NF_ASRT__(nf_get_vara_int(ncid,id_lake,(/start/),(/count/),lake)) __NF_ASRT__(nf_get_vara_int(ncid,id_soil,(/start/),(/count/),soil)) __NF_ASRT__(nf_get_vara_int(ncid,id_vegn,(/start/),(/count/),vegn)) ! create tiles do it = 1,count k = idx(it) if (k<0) cycle ! skip negative indices i = modulo(k,lnd%nlon)+1; k = k/lnd%nlon j = modulo(k,lnd%nlat)+1; k = k/lnd%nlat k = k + 1 if (ilnd%ie) cycle if (jlnd%je) cycle ! the size of the tile set at the point (i,j) must be equal to k tile=>new_land_tile(frac=frac(it),& glac=glac(it),lake=lake(it),soil=soil(it),vegn=vegn(it)) call insert(tile,lnd%tile_map(i,j)) enddo enddo __NF_ASRT__(nf_close(ncid)) deallocate(idx, glac, lake, soil, snow, cana, vegn, frac) end subroutine ! ============================================================================ subroutine update_land_model_fast ( cplr2land, land2cplr ) type(atmos_land_boundary_type), intent(in) :: cplr2land type(land_data_type) , intent(inout) :: land2cplr ! ---- local vars real :: fco2_0,Dfco2Dq , & ! co2 flux from canopy air to the atmosphere ISa_dn_dir(NBANDS), & ! downward direct sw radiation at the top of the canopy ISa_dn_dif(NBANDS), & ! downward diffuse sw radiation at the top of the canopy cana_q ! variables for stock calculations real :: & cana_VMASS, cana_HEAT, & vegn_LMASS, vegn_FMASS, vegn_HEAT, & snow_LMASS, snow_FMASS, snow_HEAT, & subs_LMASS, subs_FMASS, subs_HEAT, & glac_LMASS, glac_FMASS, glac_HEAT, & lake_LMASS, lake_FMASS, lake_HEAT, & soil_LMASS, soil_FMASS, soil_HEAT real, dimension(lnd%is:lnd%ie,lnd%js:lnd%je) :: & runoff, & ! total (liquid+snow) runoff accumulated over tiles in cell runoff_snow, & ! runoff snow accumulated over tiles in cell runoff_heat, & ! runoff heat accumulated over tiles in cell heat_frac_liq, & ! fraction of runoff heat in liquid discharge_l, & ! discharge of liquid water to ocean discharge_sink ! container to collect small/negative values for later accounting real, dimension(lnd%is:lnd%ie,lnd%js:lnd%je,num_species) :: & runoff_c, & ! runoff of tracers accumulated over tiles in cell discharge_c ! discharge of tracers to ocean logical :: used ! return value of send_data diagnostics routine real, allocatable :: runoff_1d(:),runoff_snow_1d(:),runoff_heat_1d(:) integer :: i,j,k ! lon, lat, and tile indices integer :: i_species ! river tracer iterator integer :: i1 ! index used to iterate over grid cells efficiently integer :: is,ie,js,je ! horizontal bounds of the override buffer type(land_tile_enum_type) :: ce, te ! tile enumarator type(land_tile_type), pointer :: tile ! pointer to current tile ! variables for data override real, allocatable :: phot_co2_data(:,:) ! buffer for data logical :: phot_co2_overridden ! flag indicating successful override ! start clocks call mpp_clock_begin(landClock) call mpp_clock_begin(landFastClock) ! to avoid output of static vegetation after the transitions worked and ! changed the tiling structure, static vegetation output is done here. call write_static_vegn() ! override data at the beginning of the time step is=lbound(cplr2land%t_flux,1) ; ie = is+size(cplr2land%t_flux,1)-1 js=lbound(cplr2land%t_flux,2) ; je = js+size(cplr2land%t_flux,2)-1 allocate(phot_co2_data(is:ie,js:je)) call data_override('LND','phot_co2',phot_co2_data,lnd%time, & override=phot_co2_overridden) ! clear the runoff values, for accumulation over the tiles runoff = 0 ; runoff_snow = 0 ; runoff_heat = 0 ; runoff_c = 0 ! main tile loop !$OMP parallel do schedule(dynamic) default(shared) private(i1,i,j,k,ce,te,tile,fco2_0,Dfco2Dq,ISa_dn_dir,ISa_dn_dif) do i1 = 0,(ie-is+1)*(je-js+1)-1 i = mod(i1,ie-is+1)+is j = i1/(ie-is+1)+js ! __DEBUG4__(is,js,i-is+lnd%is,j-js+lnd%js) ce = first_elmt(lnd%tile_map(i-is+lnd%is,j-js+lnd%js)) te = tail_elmt (lnd%tile_map(i-is+lnd%is,j-js+lnd%js)) k = 0 do while (ce/=te) k = k+1 ; tile=>current_tile(ce) ; ce = next_elmt(ce) ! set this point coordinates as current for debug output call set_current_point(i-is+lnd%is,j-js+lnd%js,k) if (lnd%ico2/=NO_TRACER) then fco2_0 = cplr2land%tr_flux(i,j,k, lnd%ico2) Dfco2Dq = cplr2land%dfdtr (i,j,k, lnd%ico2) else fco2_0 = 0 Dfco2Dq = 0 endif ISa_dn_dir(BAND_VIS) = cplr2land%sw_flux_down_vis_dir(i,j,k) ISa_dn_dir(BAND_NIR) = cplr2land%sw_flux_down_total_dir(i,j,k)& -cplr2land%sw_flux_down_vis_dir(i,j,k) ISa_dn_dif(BAND_VIS) = cplr2land%sw_flux_down_vis_dif(i,j,k) ISa_dn_dif(BAND_NIR) = cplr2land%sw_flux_down_total_dif(i,j,k)& -cplr2land%sw_flux_down_vis_dif(i,j,k) call update_land_model_fast_0d(tile, i,j,k, land2cplr, & cplr2land%lprec(i,j,k), cplr2land%fprec(i,j,k), cplr2land%tprec(i,j,k), & cplr2land%t_flux(i,j,k), cplr2land%dhdt(i,j,k), & cplr2land%tr_flux(i,j,k, lnd%isphum), cplr2land%dfdtr(i,j,k, lnd%isphum), & fco2_0, Dfco2Dq, & ISa_dn_dir, ISa_dn_dif, cplr2land%lwdn_flux(i,j,k), & cplr2land%ustar(i,j,k), cplr2land%p_surf(i,j,k), cplr2land%drag_q(i,j,k), & phot_co2_overridden, phot_co2_data(i,j),& runoff(i,j), runoff_heat(i,j), runoff_snow(i,j) & ) ! some of the diagnostic variables are sent from here, purely for coding ! convenince: the compute domain-level 2d and 3d vars are generaly not ! available inside update_land_model_fast_0d, so the diagnostics for those ! was left here. call send_tile_data(id_area, tile%frac*lnd%area(i,j), tile%diag) call send_tile_data(id_z0m, land2cplr%rough_mom(i,j,k), tile%diag) call send_tile_data(id_z0s, land2cplr%rough_heat(i,j,k), tile%diag) call send_tile_data(id_Trad, land2cplr%t_surf(i,j,k), tile%diag) call send_tile_data(id_Tca, land2cplr%t_ca(i,j,k), tile%diag) call send_tile_data(id_qca, land2cplr%tr(i,j,k,lnd%isphum), tile%diag) call send_tile_data(id_cd_m, cplr2land%cd_m(i,j,k), tile%diag) call send_tile_data(id_cd_t, cplr2land%cd_t(i,j,k), tile%diag) enddo enddo ! set values of tracer fluxes runoff_c(:,:,1) = runoff_snow runoff_c(:,:,2) = runoff_heat do i_species = num_phys+1, num_species runoff_c(:,:,i_species) = 0 ! age, species enddo !================================================================================= ! update river state call update_river(runoff, runoff_c, discharge_l, discharge_c) !================================================================================= discharge_l = discharge_l/lnd%cellarea do i_species = 1, num_species discharge_c(:,:,i_species) = discharge_c(:,:,i_species)/lnd%cellarea enddo ! pass through to ocean the runoff that was not seen by river module because of land_frac diffs. ! need to multiply by gfrac to spread over whole cell where (missing_rivers) discharge_l = (runoff-runoff_c(:,:,1))*frac do i_species = 1, num_species where (missing_rivers) & discharge_c(:,:,i_species) = runoff_c(:,:,i_species)*frac enddo ! don't send negatives or insignificant values to ocean. put them in the sink instead. ! this code does not seem necessary, and default discharge_tol value should be used. discharge_sink = 0. where (discharge_l.le.discharge_tol) discharge_sink = discharge_sink + discharge_l discharge_l = 0. endwhere where (discharge_c(:,:,1).le.discharge_tol) discharge_sink = discharge_sink + discharge_c(:,:,1) discharge_c(:,:,1) = 0. endwhere ! find phase partitioning ratio for discharge sensible heat flux where (discharge_l.gt.0. .or. discharge_c(:,:,1).gt.0.) heat_frac_liq = clw*discharge_l / (clw*discharge_l+csw*discharge_c(:,:,1)) elsewhere heat_frac_liq = 1. endwhere ! scale up fluxes sent to ocean to compensate for non-ocean fraction of discharge cell. ! split heat into liquid and solid streams where (frac.lt.1.) land2cplr%discharge = discharge_l / (1-frac) land2cplr%discharge_snow = discharge_c(:,:,1) / (1-frac) land2cplr%discharge_heat = heat_frac_liq*discharge_c(:,:,2) / (1-frac) land2cplr%discharge_snow_heat = discharge_c(:,:,2) / (1-frac) & - land2cplr%discharge_heat endwhere ce = first_elmt(lnd%tile_map, & is=lbound(cplr2land%t_flux,1), & js=lbound(cplr2land%t_flux,2) ) te = tail_elmt(lnd%tile_map) do while(ce /= te) call get_elmt_indices(ce,i,j,k) tile => current_tile(ce) ce=next_elmt(ce) cana_VMASS = 0. ; cana_HEAT = 0. vegn_LMASS = 0. ; vegn_FMASS = 0. ; vegn_HEAT = 0. snow_LMASS = 0. ; snow_FMASS = 0. ; snow_HEAT = 0. subs_LMASS = 0. ; subs_FMASS = 0. ; subs_HEAT = 0. glac_LMASS = 0. ; glac_FMASS = 0. ; glac_HEAT = 0. lake_LMASS = 0. ; lake_FMASS = 0. ; lake_HEAT = 0. soil_LMASS = 0. ; soil_FMASS = 0. ; soil_HEAT = 0. if (associated(tile%cana)) then call cana_state( tile%cana, cana_q=cana_q ) cana_VMASS = canopy_air_mass*cana_q cana_HEAT = cana_tile_heat(tile%cana) endif if (associated(tile%vegn)) then call vegn_tile_stock_pe(tile%vegn, vegn_LMASS, vegn_FMASS) vegn_HEAT = vegn_tile_heat(tile%vegn) endif if(associated(tile%snow)) then call snow_tile_stock_pe(tile%snow, snow_LMASS, snow_FMASS) snow_HEAT = snow_tile_heat(tile%snow) endif if (associated(tile%glac)) then call glac_tile_stock_pe(tile%glac, subs_LMASS, subs_FMASS) subs_HEAT = glac_tile_heat(tile%glac) glac_LMASS = subs_LMASS glac_FMASS = subs_FMASS glac_HEAT = subs_HEAT else if (associated(tile%lake)) then call lake_tile_stock_pe(tile%lake, subs_LMASS, subs_FMASS) subs_HEAT = lake_tile_heat(tile%lake) lake_LMASS = subs_LMASS lake_FMASS = subs_FMASS lake_HEAT = subs_HEAT else if (associated(tile%soil)) then call soil_tile_stock_pe(tile%soil, subs_LMASS, subs_FMASS) subs_HEAT = soil_tile_heat(tile%soil) soil_LMASS = subs_LMASS soil_FMASS = subs_FMASS soil_HEAT = subs_HEAT endif call send_tile_data(id_VWS, cana_VMASS, tile%diag) call send_tile_data(id_VWSc, cana_VMASS, tile%diag) call send_tile_data(id_LWS, vegn_LMASS+snow_LMASS+subs_LMASS, tile%diag) call send_tile_data(id_LWSv, vegn_LMASS, tile%diag) call send_tile_data(id_LWSs, snow_LMASS, tile%diag) call send_tile_data(id_LWSg, subs_LMASS, tile%diag) call send_tile_data(id_FWS, vegn_FMASS+snow_FMASS+subs_FMASS, tile%diag) call send_tile_data(id_FWSv, vegn_FMASS, tile%diag) call send_tile_data(id_FWSs, snow_FMASS, tile%diag) call send_tile_data(id_FWSg, subs_FMASS, tile%diag) call send_tile_data(id_HS, vegn_HEAT+snow_HEAT+subs_HEAT+cana_HEAT, tile%diag) call send_tile_data(id_HSv, vegn_HEAT, tile%diag) call send_tile_data(id_HSs, snow_HEAT, tile%diag) call send_tile_data(id_HSg, subs_HEAT, tile%diag) call send_tile_data(id_HSc, cana_HEAT, tile%diag) call send_tile_data(id_water, subs_LMASS+subs_FMASS, tile%diag) call send_tile_data(id_snow, snow_LMASS+snow_FMASS, tile%diag) enddo ! advance land model time lnd%time = lnd%time + lnd%dt_fast ! send the accumulated diagnostics to the output call dump_tile_diag_fields(lnd%tile_map, lnd%time) if (id_dis_liq > 0) used = send_data (id_dis_liq, discharge_l, lnd%time) if (id_dis_ice > 0) used = send_data (id_dis_ice, discharge_c(:,:,1), lnd%time) if (id_dis_heat > 0) used = send_data (id_dis_heat, discharge_c(:,:,2), lnd%time) if (id_dis_sink > 0) used = send_data (id_dis_sink, discharge_sink, lnd%time) ! deallocate override buffer deallocate(phot_co2_data) call mpp_clock_end(landFastClock) call mpp_clock_end(landClock) end subroutine update_land_model_fast ! ============================================================================ subroutine update_land_model_fast_0d(tile, i,j,k, land2cplr, & precip_l, precip_s, atmos_T, & Ha0, DHaDTc, Ea0, DEaDqc, fco2_0, Dfco2Dq,& ISa_dn_dir, ISa_dn_dif, ILa_dn, & ustar, p_surf, drag_q, & phot_co2_overridden, phot_co2_data, & runoff, runoff_heat, runoff_snow & ) type (land_tile_type), pointer :: tile type(land_data_type), intent(inout) :: land2cplr integer, intent(in) :: i,j,k ! coordinates real, intent(in) :: & precip_l, precip_s, & ! liquid and solid precipitation, kg/(m2 s) atmos_T, & ! incoming precipitation temperature (despite its name), deg K Ha0, DHaDTc, & ! sensible heat flux from the canopy air to the atmosphere Ea0, DEaDqc, & ! water vapor flux from canopy air to the atmosphere fco2_0,Dfco2Dq,& ! co2 flux from canopy air to the atmosphere ISa_dn_dir(NBANDS), & ! downward direct sw radiation at the top of the canopy ISa_dn_dif(NBANDS), & ! downward diffuse sw radiation at the top of the canopy ILa_dn, & ! downward lw radiation at the top of the canopy ustar, & ! friction velocity, m/s p_surf, & ! surface pressure, Pa drag_q, & ! phot_co2_data ! data input for the CO2 for photosynthesis logical, intent(in):: phot_co2_overridden real, intent(inout) :: & runoff, runoff_heat, runoff_snow ! ---- local constants ! indices of variables and equations for implicit time stepping solution : integer, parameter :: iqc=1, iTc=2, iTv=3, iwl=4, iwf=5 ! ---- local vars real :: A(5,5),B0(5),B1(5),B2(5) ! implicit equation matrix and right-hand side vectors real :: A00(5,5),B10(5),B00(5) ! copy of the above, only for debugging integer :: indx(5) ! permutation vector ! linearization coefficients of various fluxes between components of land ! surface scheme real :: & G0, DGDTg, & ! ground heat flux Hv0, DHvDTv, DHvDTc, & ! sens heat flux from vegetation Et0, DEtDTv, DEtDqc, DEtDwl, DEtDwf, & ! transpiration Eli0, DEliDTv, DEliDqc, DEliDwl, DEliDwf, & ! evaporation of intercepted water Esi0, DEsiDTv, DEsiDqc, DEsiDwl, DEsiDwf, & ! sublimation of intercepted snow Hg0, DHgDTg, DHgDTc, & ! linearization of the sensible heat flux from ground Eg0, DEgDTg, DEgDqc, DEgDpsig, & ! linearization of evaporation from ground flwv0, DflwvDTg, DflwvDTv,& ! linearization of net LW radiation to the canopy flwg0, DflwgDTg, DflwgDTv,& ! linearization of net LW radiation to the canopy vegn_drip_l, vegn_drip_s, & ! drip rate of water and snow, respectively, kg/(m2 s) vegn_lai ! increments of respective variables over time step, results of the implicit ! time step: real :: delta_qc, delta_Tc, delta_Tv, delta_wl, delta_ws, delta_Tg, delta_psig real :: flwg ! updated value of long-wave ground energy balance real :: vegn_emis_lw, surf_emis_lw ! emissivities of ground and surface real :: vegn_emsn, surf_emsn ! emission by vegetation and surface, respectively real :: denom ! denominator in the LW radiative balance calculations real :: sum0, sum1 real :: & grnd_T, gT, & ! ground temperature and its value used for sensible heat advection vegn_T, vT, & ! vegetation (canopy) temperature cana_T, cT, & ! canopy air temperature evap_T, eT, & ! temperature assigned to vapor going between land and atmosphere soil_uptake_T, & ! average temperature of water taken up by the vegetation vegn_Wl, vegn_Ws, & ! water and snow mass of the canopy vegn_ifrac, & ! intercepted fraction of liquid or frozen precipitation vegn_hcap, & ! vegetation heat capacity, including intercepted water and snow vegn_fco2, & ! co2 flux from the vegetation, kg CO2/(m2 s) hlv_Tv, hlv_Tu, & ! latent heat of vaporization at vegn and uptake temperatures, respectively hls_Tv, & ! latent heat of sublimation at vegn temperature grnd_rh, & ! explicit relative humidity at ground surface grnd_rh_psi, & ! psi derivative of relative humidity at ground surface grnd_liq, grnd_ice, grnd_subl, & grnd_tf, & ! temperature of freezing on the ground grnd_latent, & grnd_flux, & grnd_E_min, & grnd_E_max, & soil_E_min, & soil_E_max, & soil_beta, & RSv(NBANDS), & ! net short-wave radiation balance of the canopy, W/m2 con_g_h, con_g_v, & ! turbulent cond. between ground and canopy air, for heat and vapor respectively snow_area, & cana_q, & ! specific humidity of canopy air cana_co2, & ! co2 moist mixing ratio in canopy air, kg CO2/kg wet air cana_co2_mol, & ! co2 dry mixing ratio in canopy air, mol CO2/mol dry air fswg, evapg, sensg, & subs_G, subs_G2, Mg_imp, snow_G_Z, snow_G_TZ, & snow_avrg_T, delta_T_snow, & ! vertically-average snow temperature and it's change due to s vegn_ovfl_l, vegn_ovfl_s, & ! overflow of liquid and solid water from the canopy vegn_ovfl_Hl, vegn_ovfl_Hs, & ! heat flux from canopy due to overflow delta_fprec, & ! correction of below-canopy solid precip in case it's average T > tfreeze hprec, & ! sensible heat flux carried by precipitation hevap, & ! sensible heat flux carried by total evapotranspiration land_evap, & ! total vapor flux from land to atmosphere land_sens, & ! turbulent sensible heat flux from land to atmosphere vegn_flw,vegn_sens,snow_sens,snow_levap,snow_fevap,snow_melt,& snow_lprec, snow_hlprec,snow_lrunf,vegn_levap,vegn_fevap,vegn_uptk,& vegn_fsw, vegn_melt,vegn_lprec,vegn_fprec,vegn_hlprec,vegn_hfprec,& precip_T,pT,snow_fsw,snow_flw,snow_frunf,snow_hlrunf,& snow_hfrunf,subs_fsw,subs_flw,subs_sens,& subs_DT, subs_M_imp, subs_evap, snow_Tbot, snow_Cbot, snow_C, subs_levap,& subs_fevap,subs_melt,subs_lrunf,subs_hlrunf,& subs_Ttop,subs_Ctop, subs_subl, new_T real :: soil_water_supply ! supply of water to roots, per unit active root biomass, kg/m2 real :: snow_T, snow_rh, snow_liq, snow_ice, snow_subl integer :: ii, jj ! indices for debug output integer :: ierr logical :: conserve_glacier_mass, snow_active, redo_leaf_water integer :: canopy_water_step real :: subs_z0m, subs_z0s, snow_z0m, snow_z0s, grnd_z0s soil_uptake_T = tfreeze ! just to avoid using un-initialized values soil_water_supply = 0.0 if (associated(tile%glac)) then call glac_step_1 ( tile%glac, & grnd_T, grnd_rh, grnd_liq, grnd_ice, grnd_subl, grnd_tf, & snow_G_Z, snow_G_TZ, conserve_glacier_mass ) grnd_E_min = -HUGE(grnd_E_min) grnd_E_max = HUGE(grnd_E_max) grnd_rh_psi = 0 else if (associated(tile%lake)) then call lake_step_1 ( ustar, p_surf, & lnd%lat(i,j), tile%lake, & grnd_T, grnd_rh, grnd_liq, grnd_ice, grnd_subl, grnd_tf, & snow_G_Z, snow_G_TZ) grnd_E_min = -HUGE(grnd_E_min) grnd_E_max = HUGE(grnd_E_max) grnd_rh_psi = 0 else if (associated(tile%soil)) then call soil_step_1 ( tile%soil, tile%vegn, tile%diag, & grnd_T, soil_uptake_T, soil_beta, soil_water_supply, soil_E_min, soil_E_max, & grnd_rh, grnd_rh_psi, grnd_liq, grnd_ice, grnd_subl, grnd_tf, & snow_G_Z, snow_G_TZ) grnd_E_min = soil_E_min grnd_E_max = soil_E_max grnd_liq = 0 ! sorry, but solver cannot handle implicit melt anymore grnd_ice = 0 ! sorry, but solver cannot handle implicit melt anymore ! no big loss, it's just the surface layer anyway else call get_current_point(face=ii) call error_mesg('update_land_model_fast','none of the surface tiles exist at ('//& trim(string(i))//','//trim(string(j))//','//trim(string(k))//& ', face='//trim(string(ii))//')',FATAL) endif subs_subl = grnd_subl call snow_step_1 ( tile%snow, snow_G_Z, snow_G_TZ, & snow_active, snow_T, snow_rh, snow_liq, snow_ice, & snow_subl, snow_area, G0, DGDTg ) if (snow_active) then grnd_T = snow_T; grnd_rh = snow_rh; grnd_liq = snow_liq grnd_rh_psi = 0 grnd_ice = snow_ice; grnd_subl = snow_subl; grnd_tf = tfreeze grnd_E_min = -HUGE(grnd_E_min) grnd_E_max = HUGE(grnd_E_max) endif call cana_state(tile%cana, cana_T, cana_q, cana_co2) if (associated(tile%vegn)) then ! Calculate net short-wave radiation input to the vegetation RSv = tile%Sv_dir*ISa_dn_dir + tile%Sv_dif*ISa_dn_dif call soil_diffusion(tile%soil, subs_z0s, subs_z0m) call snow_diffusion(tile%snow, snow_z0s, snow_z0m) grnd_z0s = exp( (1-snow_area)*log(subs_z0s) + snow_area*log(snow_z0s)) ! cana_co2 is moist mass mixing ratio [kg CO2/kg wet air], convert it to dry ! volumetric mixing ratio [mol CO2/mol dry air] cana_co2_mol = cana_co2*mol_air/mol_CO2/(1-cana_q) if (phot_co2_overridden) cana_co2_mol = phot_co2_data call vegn_step_1 ( tile%vegn, tile%soil, tile%diag, & p_surf, & ustar, & drag_q, & ISa_dn_dir+ISa_dn_dif, RSv, precip_l, precip_s, & tile%land_d, tile%land_z0s, tile%land_z0m, grnd_z0s, & soil_beta, soil_water_supply,& cana_T, cana_q, cana_co2_mol, & ! output con_g_h, con_g_v, & vegn_T, vegn_Wl, vegn_Ws, & ! temperature, water and snow mass on the canopy vegn_ifrac, vegn_lai, & vegn_drip_l, vegn_drip_s,& vegn_hcap, & ! total vegetation heat capacity (including intercepted water/snow) Hv0, DHvDTv, DHvDTc, & Et0, DEtDTv, DEtDqc, DEtDwl, DEtDwf, & Eli0, DEliDTv, DEliDqc, DEliDwl, DEliDwf, & Esi0, DEsiDTv, DEsiDqc, DEsiDwl, DEsiDwf ) if (LM2) then con_g_h = con_g_h * con_fac_large if (snow_active) then con_g_v = con_g_v * con_fac_large else con_g_v = con_g_v * con_fac_small endif endif else RSv = 0 con_g_h = con_fac_large ; con_g_v = con_fac_large if(associated(tile%glac).and.conserve_glacier_mass.and..not.snow_active) & con_g_v = con_fac_small vegn_T = cana_T ; vegn_Wl = 0 ; vegn_Ws = 0 vegn_ifrac = 0 ; vegn_lai = 0 vegn_drip_l = 0 ; vegn_drip_s = 0 vegn_hcap = 1.0 Hv0 =0; DHvDTv =0; DHvDTc=0; Et0 =0; DEtDTv =0; DEtDqc=0; DEtDwl=0; DEtDwf=0 Eli0=0; DEliDTv=0; DEliDqc=0; DEliDwl=0; DEliDwf=0 Esi0=0; DEsiDTv=0; DEsiDqc=0; DEsiDwl=0; DEsiDwf=0 endif ! calculate net shortwave for ground and canopy fswg = SUM(tile%Sg_dir*ISa_dn_dir + tile%Sg_dif*ISa_dn_dif) vegn_fsw = SUM(RSv) call cana_step_1 (tile%cana, p_surf, con_g_h, con_g_v, & grnd_t, grnd_rh, grnd_rh_psi, & Hg0, DHgDTg, DHgDTc, Eg0, DEgDTg, DEgDqc, DEgDpsig) ! [X.X] using long-wave optical properties, calculate the explicit long-wave ! radiative balances and their derivatives w.r.t. temperatures vegn_emis_lw = 1 - tile%vegn_refl_lw - tile%vegn_tran_lw surf_emis_lw = 1 - tile%surf_refl_lw denom = 1-tile%vegn_refl_lw*tile%surf_refl_lw vegn_emsn = vegn_emis_lw * stefan * vegn_T**4 surf_emsn = surf_emis_lw * stefan * grnd_T**4 flwv0 = ILa_dn * vegn_emis_lw*(1+tile%vegn_refl_lw*tile%surf_refl_lw/denom) & + vegn_emsn * (tile%surf_refl_lw*vegn_emis_lw/denom-2) & + surf_emsn * vegn_emis_lw/denom DflwvDTg = vegn_emis_lw/denom * surf_emis_lw * stefan * 4 * grnd_T**3 DflwvDTv = (tile%surf_refl_lw*vegn_emis_lw/denom-2)* vegn_emis_lw * stefan * 4 * vegn_T**3 flwg0 = (ILa_dn*tile%vegn_tran_lw + vegn_emsn)*(1-tile%surf_refl_lw)/denom & - surf_emsn*(1-tile%vegn_refl_lw)/denom DflwgDTg = -(1-tile%vegn_refl_lw)/denom * surf_emis_lw * stefan * 4 * grnd_T**3 DflwgDTv = (1-tile%surf_refl_lw)/denom * vegn_emis_lw * stefan * 4 * vegn_T**3 ! [X.0] calculate the latent heats of vaporization at appropriate temperatures if (use_tfreeze_in_grnd_latent) then grnd_latent = hlv + hlf*grnd_subl else grnd_latent = hlv + (cpw-clw)*(grnd_T-tfreeze) & + (hlf + (clw-csw)*(grnd_T-tfreeze)) * grnd_subl endif if (use_atmos_T_for_precip_T) then precip_T = atmos_T else precip_T = cana_T endif if (use_atmos_T_for_evap_T) then evap_T = atmos_T else evap_T = cana_T endif if (use_old_conservation_equations) then hlv_Tv = hlv - (cpw-clw)*tfreeze + cpw*vegn_T hls_Tv = hlv + hlf - (cpw-csw)*tfreeze + cpw*vegn_T hlv_Tu = hlv - (cpw-clw)*tfreeze + cpw*vegn_T - clw*soil_uptake_T pT = precip_T cT = cana_T eT = evap_T gT = grnd_T vT = vegn_T else hlv_Tv = hlv + cpw*(vegn_T-tfreeze) hls_Tv = hlf + hlv_Tv hlv_Tu = hlv_Tv - clw*(soil_uptake_T-tfreeze) pT = precip_T-tfreeze cT = cana_T-tfreeze eT = evap_T-tfreeze gT = grnd_T-tfreeze vT = vegn_T-tfreeze endif do canopy_water_step = 1,2 if(is_watch_point()) then write(*,*)'#### input data for the matrix ####' __DEBUG1__(delta_time) __DEBUG4__(vegn_T,vT,vegn_Wl,vegn_Ws) __DEBUG3__(grnd_T,gT,grnd_rh) __DEBUG3__(cana_T,cT,cana_q) __DEBUG2__(evap_T,eT) __DEBUG2__(vegn_emis_lw,surf_emis_lw) __DEBUG2__(vegn_emsn,surf_emsn) __DEBUG4__(precip_l, vegn_drip_l, pT, precip_T) __DEBUG2__(precip_s, vegn_drip_s) __DEBUG2__(vegn_ifrac, vegn_lai) __DEBUG1__(ILa_dn) __DEBUG2__(ISa_dn_dir(1),ISa_dn_dir(2)) __DEBUG2__(ISa_dn_dif(1),ISa_dn_dif(2)) __DEBUG2__(fswg, vegn_fsw) __DEBUG1__(vegn_hcap) __DEBUG3__(hlv_Tv, hlv_Tu, hls_Tv) __DEBUG2__(G0, DGDTg) __DEBUG2__(Ha0, DHaDTc) __DEBUG2__(Ea0, DEaDqc) __DEBUG3__(Hv0, DHvDTv, DHvDTc) __DEBUG5__(Et0, DEtDTv, DEtDqc, DEtDwl, DEtDwf) __DEBUG5__(Eli0, DEliDTv, DEliDqc, DEliDwl, DEliDwf) __DEBUG5__(Esi0, DEsiDTv, DEsiDqc, DEsiDwl, DEsiDwf) __DEBUG3__(Hg0, DHgDTg, DHgDTc) __DEBUG3__(Eg0, DEgDTg, DEgDqc) __DEBUG3__(flwv0, DflwvDTg, DflwvDTv) __DEBUG3__(flwg0, DflwgDTg, DflwgDTv) __DEBUG2__(tile%e_res_1,tile%e_res_2) endif ! [X.1] form the system of equations for implicit scheme, such that A*X = B1*delta_Tg+B2*delta_psig+B0 ! [X.1.1] equation of canopy air mass balance A(iqc,iqc) = canopy_air_mass/delta_time-DEtDqc-DEliDqc-DEsiDqc-DEgDqc+DEaDqc A(iqc,iTc) = 0 A(iqc,iTv) = -DEtDTv-DEliDTv-DEsiDTv A(iqc,iwl) = -DEtDwl-DEliDwl-DEsiDwl A(iqc,iwf) = -DEtDwf-DEliDwf-DEsiDwf B0(iqc) = Esi0+Eli0+Et0+Eg0-Ea0 B1(iqc) = DEgDTg B2(iqc) = DEgDpsig ! [X.1.2] equation of canopy air energy balance #ifdef USE_DRY_CANA_MASS A(iTc,iqc) = canopy_air_mass*cpw*cT/delta_time & #else A(iTc,iqc) = canopy_air_mass*(cpw-cp_air)*cT/delta_time & #endif - cpw*vT*(DEtDqc+DEliDqc+DEsiDqc) - cpw*gT*DEgDqc + cpw*eT*DEaDqc #ifdef USE_DRY_CANA_MASS A(iTc,iTc) = canopy_air_mass*cp_air/delta_time-DHvDTc-DHgDTc+DHaDTc #else A(iTc,iTc) = canopy_air_mass*(cp_air+cana_q*(cpw-cp_air))/delta_time-DHvDTc-DHgDTc+DHaDTc #endif A(iTc,iTv) = -DHvDTv-cpw*vT*(DEtDTv+DEliDTv+DEsiDTv) A(iTc,iwl) = -cpw*vT*(DEtDwl+DEliDwl+DEsiDwl) A(iTc,iwf) = -cpw*vT*(DEtDwf+DEliDwf+DEsiDwf) B0(iTc) = Hv0 + Hg0 - Ha0 + cpw*(vT*(Et0+Eli0+Esi0)+gT*Eg0-eT*Ea0) - tile%e_res_1 B1(iTc) = DHgDTg + cpw*gT*DEgDTg B2(iTc) = cpw*gT*DEgDpsig ! [X.1.3] equation of canopy energy balance A(iTv,iqc) = hlv_Tu*DEtDqc + hlv_Tv*DEliDqc + hls_Tv*DEsiDqc A(iTv,iTc) = DHvDTc A(iTv,iTv) = vegn_hcap/delta_time-DflwvDTv + DHvDTv + & hlv_Tu*DEtDTv + hlv_Tv*DEliDTv + hls_Tv*DEsiDTv + clw*vegn_drip_l + csw*vegn_drip_s A(iTv,iwl) = clw*vT/delta_time + hlv_Tu*DEtDwl + hlv_Tv*DEliDwl + hls_Tv*DEsiDwl A(iTv,iwf) = csw*vT/delta_time + hlv_Tu*DEtDwf + hlv_Tv*DEliDwf + hls_Tv*DEsiDwf B0(iTv) = vegn_fsw + flwv0 - Hv0 - hlv_Tu*Et0 - Hlv_Tv*Eli0 - hls_Tv*Esi0 & + clw*precip_l*vegn_ifrac*pT + csw*precip_s*vegn_ifrac*pT & - clw*vegn_drip_l*vT - csw*vegn_drip_s*vT - tile%e_res_2 B1(iTv) = DflwvDTg B2(iTv) = 0 ! [X.1.4] equation of intercepted liquid water mass balance A(iwl,iqc) = DEliDqc A(iwl,iTc) = 0 A(iwl,iTv) = DEliDTv A(iwl,iwl) = 1.0/delta_time + DEliDwl A(iwl,iwf) = DEliDwf B0(iwl) = -Eli0 + precip_l*vegn_ifrac - vegn_drip_l B1(iwl) = 0 B2(iwl) = 0 ! [X.1.5] equation of intercepted frozen water mass balance A(iwf,iqc) = DEsiDqc A(iwf,iTc) = 0 A(iwf,iTv) = DEsiDTv A(iwf,iwl) = DEsiDwl A(iwf,iwf) = 1.0/delta_time + DEsiDwf B0(iwf) = -Esi0 + precip_s*vegn_ifrac - vegn_drip_s B1(iwf) = 0 B2(iwf) = 0 ! [X.1.6] if LAI becomes zero (and, therefore, all fluxes from vegetation and their ! derivatives must be zero too) we get a degenerate case. Still, the drip may be non-zero ! because some water may remain from before leaf drop, and non-zero energy residual can be ! carried over from the previous time step. ! To prevent temperature from going haywire in those cases, we simply replace the equations ! of canopy energy and mass balance with the following: ! vegn_T + delta_Tv = cana_T + delta_Tc ! delta_Wl = -vegn_drip_l*delta_time ! delta_Ws = -vegn_drip_s*delta_time ! the residual vegn_Wl and vegn_Ws, if any, are taken care of by the overflow calculations if(vegn_hcap==0) then ! vegn_T + delta_Tv = cana_T + delta_Tc A(iTv,:) = 0 A(iTv,iTc) = -1 A(iTv,iTv) = +1 B0(iTv) = cana_T - vegn_T B1(iTv) = 0 ! delta_Wl = -vegn_drip_l*delta_time A(iwl,:) = 0 A(iwl,iwl) = 1 B0(iwl) = -vegn_drip_l*delta_time B1(iwl) = 0 ! delta_Ws = -vegn_drip_s*delta_time A(iwf,:) = 0 A(iwf,iwf) = 1 B0(iwf) = -vegn_drip_s*delta_time B1(iwf) = 0 endif if(is_watch_point()) then write(*,*)'#### A, B0, B1, B2 ####' do ii = 1, size(A,1) write(*,'(99g23.16)')(A(ii,jj),jj=1,size(A,2)),B0(ii),B1(ii),B2(ii) enddo endif A00 = A B00 = B0 B10 = B1 ! [X.2] solve the system for free terms and delta_Tg and delta_psig terms, getting ! linear equation for delta_Tg and delta_psig call ludcmp(A,indx, ierr) if (ierr/=0)& write(*,*) 'Matrix is singular',i,j,k call lubksb(A,indx,B0) call lubksb(A,indx,B1) call lubksb(A,indx,B2) if(is_watch_point()) then write(*,*)'#### solution: B0, B1, B2 ####' do ii = 1, size(A,1) __DEBUG3__(B0(ii),B1(ii),B2(ii)) enddo !!$ write(*,*)'#### solution check ####' !!$ do ii = 1, size(A,1) !!$ sum0 = 0; sum1 = 0; !!$ do jj = 1, size(A,2) !!$ sum0 = sum0 + A00(ii,jj)*B0(jj) !!$ sum1 = sum1 + A00(ii,jj)*B1(jj) !!$ enddo !!$ write(*,'(99g)')sum0-B00(ii),sum1-B10(ii) !!$ enddo endif ! the result of this solution is a set of expressions for delta_xx in terms ! of delta_Tg and delta_psig: ! delta_xx(i) = B0(i) + B1(i)*delta_Tg + B2(i)*delta_psig. Note that A, B0, B1 and B2 ! are destroyed in the process: A is replaced with LU-decomposition, and ! B0, B1, B2 are replaced with solutions ! solve the non-linear equation for energy balance at the surface. call land_surface_energy_balance( & grnd_T, grnd_liq, grnd_ice, grnd_latent, grnd_Tf, grnd_E_min, & grnd_E_max, fswg, & flwg0 + b0(iTv)*DflwgDTv, DflwgDTg + b1(iTv)*DflwgDTv, b2(iTv)*DflwgDTv, & Hg0 + b0(iTc)*DHgDTc, DHgDTg + b1(iTc)*DHgDTc, b2(iTc)*DHgDTc, & Eg0 + b0(iqc)*DEgDqc, DEgDTg + b1(iqc)*DEgDqc, DEgDpsig + b2(iqc)*DEgDqc, & G0, DGDTg, & ! output delta_Tg, delta_psig, Mg_imp ) ! [X.5] calculate final value of other tendencies delta_qc = B0(iqc) + B1(iqc)*delta_Tg + B2(iqc)*delta_psig delta_Tc = B0(iTc) + B1(iTc)*delta_Tg + B2(iTc)*delta_psig delta_Tv = B0(iTv) + B1(iTv)*delta_Tg + B2(iTv)*delta_psig delta_wl = B0(iwl) + B1(iwl)*delta_Tg + B2(iwl)*delta_psig delta_ws = B0(iwf) + B1(iwf)*delta_Tg + B2(iwf)*delta_psig ! [X.6] calculate updated values of energy balance components used in further ! calculations flwg = flwg0 + DflwgDTg*delta_Tg + DflwgDTv*delta_Tv evapg = Eg0 + DEgDTg*delta_Tg + DEgDpsig*delta_psig + DEgDqc*delta_qc sensg = Hg0 + DHgDTg*delta_Tg + DHgDTc*delta_Tc grnd_flux = G0 + DGDTg*delta_Tg vegn_sens = Hv0 + DHvDTv*delta_Tv + DHvDTc*delta_Tc vegn_levap = Eli0 + DEliDTv*delta_Tv + DEliDqc*delta_qc + DEliDwl*delta_wl + DEliDwf*delta_ws vegn_fevap = Esi0 + DEsiDTv*delta_Tv + DEsiDqc*delta_qc + DEsiDwl*delta_wl + DEsiDwf*delta_ws vegn_uptk = Et0 + DEtDTv*delta_Tv + DEtDqc*delta_qc + DEtDwl*delta_wl + DEtDwf*delta_ws vegn_flw = flwv0 + DflwvDTv*delta_Tv + DflwvDTg*delta_Tg land_evap = Ea0 + DEaDqc*delta_qc land_sens = Ha0 + DHaDTc*delta_Tc ! [X.7] calculate energy residuals due to cross-product of time tendencies #ifdef USE_DRY_CANA_MASS tile%e_res_1 = canopy_air_mass*cpw*delta_qc*delta_Tc/delta_time #else tile%e_res_1 = canopy_air_mass*(cpw-cp_air)*delta_qc*delta_Tc/delta_time #endif tile%e_res_2 = delta_Tv*(clw*delta_Wl+csw*delta_Ws)/delta_time ! calculate the final value upward long-wave radiation flux from the land, to be ! returned to the flux exchange. tile%lwup = ILa_dn - vegn_flw - flwg if(is_watch_point())then write(*,*)'#### ground balance' __DEBUG2__(fswg,flwg) __DEBUG2__(sensg,evapg*grnd_latent) __DEBUG1__(grnd_flux) __DEBUG1__(Mg_imp) write(*,*)'#### implicit time steps' __DEBUG3__(delta_Tg, grnd_T, grnd_T+delta_Tg ) __DEBUG1__(delta_psig ) __DEBUG3__(delta_qc, cana_q, cana_q+delta_qc ) __DEBUG3__(delta_Tc, cana_T, cana_T+delta_Tc ) __DEBUG3__(delta_Tv, vegn_T, vegn_T+delta_Tv ) __DEBUG3__(delta_wl, vegn_Wl, vegn_Wl+delta_wl) __DEBUG3__(delta_ws, vegn_Ws, vegn_Ws+delta_ws) __DEBUG2__(tile%e_res_1, tile%e_res_2) write(*,*)'#### resulting fluxes' __DEBUG4__(flwg, evapg, sensg, grnd_flux) __DEBUG3__(vegn_levap,vegn_fevap,vegn_uptk) __DEBUG2__(vegn_sens,vegn_flw) __DEBUG1__(Ea0+DEaDqc*delta_qc) __DEBUG2__(tile%cana%prog%q,cana_q) endif if (.not.prohibit_negative_canopy_water) exit ! do no corrections redo_leaf_water = .FALSE. if (vegn_Wl+delta_wl<0) then redo_leaf_water = .TRUE. Eli0 = vegn_Wl/delta_time + precip_l*vegn_ifrac - vegn_drip_l DEliDTv = 0.0; DEliDqc = 0.0 DEliDwl = 0.0; DEliDwf = 0.0 endif if (vegn_Ws+delta_ws<0) then redo_leaf_water = .TRUE. Esi0 = vegn_Ws/delta_time + precip_s*vegn_ifrac - vegn_drip_s DEsiDTv = 0.0; DEsiDqc = 0.0 DEsiDwl = 0.0; DEsiDwf = 0.0 endif if (.not.redo_leaf_water) exit ! from loop enddo ! canopy_water_step call cana_step_2 ( tile%cana, delta_Tc, delta_qc ) if(associated(tile%vegn)) then call vegn_step_2 ( tile%vegn, tile%diag, & delta_Tv, delta_wl, delta_ws, & vegn_melt, & vegn_ovfl_l, vegn_ovfl_s, & vegn_ovfl_Hl, vegn_ovfl_Hs ) ! calculate total amount of liquid and solid precipitation below the canopy vegn_lprec = (1-vegn_ifrac)*precip_l + vegn_drip_l + vegn_ovfl_l vegn_fprec = (1-vegn_ifrac)*precip_s + vegn_drip_s + vegn_ovfl_s ! calculate heat carried by liquid and solid precipitation below the canopy vegn_hlprec = clw*((1-vegn_ifrac)*precip_l*(precip_T-tfreeze) & + vegn_drip_l*(vegn_T+delta_Tv-tfreeze)) & + vegn_ovfl_Hl vegn_hfprec = csw*((1-vegn_ifrac)*precip_s*(precip_T-tfreeze) & + vegn_drip_s*(vegn_T+delta_Tv-tfreeze)) & + vegn_ovfl_Hs ! make sure the temperature of the snow falling below canopy is below freezing ! this correction was introduced in an attempt to fix the problem with fictitious ! heat accumulating in near-zero-mass snow; however it does not seem to make a ! difference. if(vegn_hfprec>0)then ! solid precipitation from vegetation carries positive energy -- we can't have ! that, because that would bring snow T above tfreeze, so convert excess to ! liquid delta_fprec = min(vegn_fprec,vegn_hfprec/hlf) vegn_fprec = vegn_fprec - delta_fprec vegn_lprec = vegn_lprec + delta_fprec vegn_hfprec = vegn_hfprec - hlf*delta_fprec ! we don't need to correct the vegn_hlprec since the temperature of additional ! liquid precip is tfreeze, and therefore its contribution to vegn_hlprec is ! exactly zero endif ! possibly we need to correct for the opposite situation: negative energy carried ! by liquid precipitation. else vegn_lprec = precip_l vegn_fprec = precip_s vegn_hlprec = precip_l*clw*(precip_T-tfreeze) vegn_hfprec = precip_s*csw*(precip_T-tfreeze) ! the fields below are only used in diagnostics vegn_melt = 0 vegn_fsw = 0 endif call snow_step_2 ( tile%snow, & snow_subl, vegn_lprec, vegn_fprec, vegn_hlprec, vegn_hfprec, & delta_Tg, Mg_imp, evapg, fswg, flwg, sensg, & use_tfreeze_in_grnd_latent, & ! output: subs_DT, subs_M_imp, subs_evap, subs_fsw, subs_flw, subs_sens, & snow_fsw, snow_flw, snow_sens, & snow_levap, snow_fevap, snow_melt, & snow_lprec, snow_hlprec, snow_lrunf, snow_frunf, & snow_hlrunf, snow_hfrunf, snow_Tbot, snow_Cbot, snow_C, snow_avrg_T ) if(is_watch_point()) then write(*,*) 'subs_M_imp', subs_M_imp endif if (snow_active) then subs_G = snow_G_Z+snow_G_TZ*subs_DT else subs_G = 0 endif if (associated(tile%glac)) then call glac_step_2 & ( tile%glac, tile%diag, subs_subl, snow_lprec, snow_hlprec, & subs_DT, subs_M_imp, subs_evap, & subs_levap, subs_fevap, & subs_melt, subs_lrunf, subs_hlrunf, subs_Ttop, subs_Ctop ) else if (associated(tile%lake)) then call lake_step_2 & ( tile%lake, tile%diag, subs_subl, snow_lprec, snow_hlprec, & subs_DT, subs_M_imp, subs_evap, & use_tfreeze_in_grnd_latent, subs_levap, subs_fevap, & subs_melt, subs_Ttop, subs_Ctop ) subs_lrunf = 0. subs_hlrunf = 0. else if (associated(tile%soil)) then call soil_step_2 & ( tile%soil, tile%vegn, tile%diag, subs_subl, snow_lprec, snow_hlprec, & vegn_uptk, subs_DT, subs_M_imp, subs_evap, & use_tfreeze_in_grnd_latent, & ! output: subs_levap, subs_fevap, & subs_melt, subs_lrunf, subs_hlrunf, subs_Ttop, subs_Ctop ) endif ! TEMP FIX: MAIN PROG SHOULD NOT TOUCH CONTENTS OF PROG VARS. ****** ! ALSO, DIAGNOSTICS IN COMPONENT MODULES SHOULD _FOLLOW_ THIS ADJUSTMENT****** if (LM2) then tile%snow%T = subs_Ttop subs_G2 = 0. else if (sum(tile%snow%ws(:))>0)then new_T = (subs_Ctop*subs_Ttop +snow_Cbot*snow_Tbot) & / (subs_Ctop+snow_Cbot) tile%snow%T(size(tile%snow%T)) = new_T if(associated(tile%glac)) tile%glac%prog(1)%T = new_T if(associated(tile%lake)) tile%lake%prog(1)%T = new_T if(associated(tile%soil)) tile%soil%prog(1)%T = new_T subs_G2 = subs_Ctop*(new_T-subs_Ttop)/delta_time else if(tau_snow_T_adj>=0) then delta_T_snow = subs_Ctop*(subs_Ttop-snow_avrg_T)/& (subs_Ctop*tau_snow_T_adj/delta_time+subs_Ctop+snow_C) tile%snow%T(:) = snow_avrg_T + delta_T_snow new_T = subs_Ttop-snow_C/subs_Ctop*delta_T_snow if(associated(tile%glac)) tile%glac%prog(1)%T = new_T if(associated(tile%lake)) tile%lake%prog(1)%T = new_T if(associated(tile%soil)) tile%soil%prog(1)%T = new_T subs_G2 = subs_Ctop*(new_T-subs_Ttop)/delta_time else subs_G2 = 0. endif endif endif vegn_fco2 = 0 if (associated(tile%vegn)) then ! do the calculations that require updated land surface prognostic variables call vegn_step_3 (tile%vegn, tile%soil, tile%cana%prog%T, precip_l+precip_s, & vegn_fco2, tile%diag) ! if vegn is present, then soil must be too call soil_step_3(tile%soil, tile%diag) endif ! update co2 concentration in the canopy air. It would be more consistent to do that ! in the same place and fashion as the rest of prognostic variables: that is, have the ! vegn_step_1 (and perhaps other *_step_1 procedures) calculate fluxes and their ! derivatives, then solve the linear equation(s), and finally have cana_step_2 update ! the concentration. if(update_cana_co2) then tile%cana%prog%co2 = tile%cana%prog%co2 + & (vegn_fco2 - fco2_0)/(canopy_air_mass_for_tracers/delta_time+Dfco2Dq) endif if(is_watch_point())then __DEBUG1__(tile%cana%prog%co2) __DEBUG3__(fco2_0,Dfco2Dq,vegn_fco2) endif call update_land_bc_fast (tile, i,j,k, land2cplr) runoff = runoff + (snow_frunf + subs_lrunf + snow_lrunf )*tile%frac runoff_heat = runoff_heat + (snow_hfrunf + subs_hlrunf + snow_hlrunf)*tile%frac runoff_snow = runoff_snow + snow_frunf*tile%frac hprec = (clw*precip_l+csw*precip_s)*(precip_T-tfreeze) hevap = cpw*land_evap*(evap_T-tfreeze) ! ---- diagnostic section ---------------------------------------------- call send_tile_data(id_frac, tile%frac, tile%diag) call send_tile_data(id_ntiles, 1.0, tile%diag) call send_tile_data(id_precip, precip_l+precip_s, tile%diag) call send_tile_data(id_hprec, hprec, tile%diag) call send_tile_data(id_lprec, precip_l, tile%diag) call send_tile_data(id_lprecv, precip_l-vegn_lprec, tile%diag) call send_tile_data(id_lprecs, vegn_lprec-snow_lprec, tile%diag) call send_tile_data(id_lprecg, snow_lprec, tile%diag) call send_tile_data(id_hlprec, clw*precip_l*(precip_T-tfreeze), tile%diag) call send_tile_data(id_hlprecv, clw*precip_l*(precip_T-tfreeze)-vegn_hlprec, & tile%diag) call send_tile_data(id_hlprecs, vegn_hlprec-snow_hlprec, tile%diag) call send_tile_data(id_hlprecg, snow_hlprec, tile%diag) call send_tile_data(id_fprec, precip_s, tile%diag) call send_tile_data(id_fprecv, precip_s-vegn_fprec, tile%diag) call send_tile_data(id_fprecs, vegn_fprec, tile%diag) call send_tile_data(id_hfprec, csw*precip_s*(precip_T-tfreeze), tile%diag) call send_tile_data(id_hfprecv, csw*precip_s*(precip_T-tfreeze)-vegn_hfprec, & tile%diag) call send_tile_data(id_hfprecs, vegn_hfprec, tile%diag) call send_tile_data(id_evap, land_evap, tile%diag) call send_tile_data(id_hevap, hevap, tile%diag) call send_tile_data(id_levap, vegn_levap+snow_levap+subs_levap+vegn_uptk, & tile%diag) call send_tile_data(id_levapv, vegn_levap, tile%diag) call send_tile_data(id_levaps, snow_levap, tile%diag) call send_tile_data(id_levapg, subs_levap, tile%diag) call send_tile_data(id_hlevap, cpw*vegn_levap*(vegn_T-tfreeze) & +cpw*snow_levap*(snow_T-tfreeze) & +cpw*subs_levap*(grnd_T-tfreeze), tile%diag) call send_tile_data(id_hlevapv, cpw*vegn_levap*(vegn_T-tfreeze), tile%diag) call send_tile_data(id_hlevaps, cpw*snow_levap*(snow_T-tfreeze), tile%diag) call send_tile_data(id_hlevapg, cpw*subs_levap*(grnd_T-tfreeze), tile%diag) call send_tile_data(id_fevap, vegn_fevap+snow_fevap+subs_fevap, tile%diag) call send_tile_data(id_fevapv, vegn_fevap, tile%diag) call send_tile_data(id_fevaps, snow_fevap, tile%diag) call send_tile_data(id_fevapg, subs_fevap, tile%diag) call send_tile_data(id_hfevap, cpw*vegn_fevap*(vegn_T-tfreeze) & +cpw*snow_fevap*(snow_T-tfreeze) & +cpw*subs_fevap*(grnd_T-tfreeze), tile%diag) call send_tile_data(id_hfevapv, cpw*vegn_fevap*(vegn_T-tfreeze), tile%diag) call send_tile_data(id_hfevaps, cpw*snow_fevap*(snow_T-tfreeze), tile%diag) call send_tile_data(id_hfevapg, cpw*subs_fevap*(grnd_T-tfreeze), tile%diag) call send_tile_data(id_runf, snow_lrunf+snow_frunf+subs_lrunf, tile%diag) call send_tile_data(id_hrunf, snow_hlrunf+snow_hfrunf+subs_hlrunf,tile%diag) call send_tile_data(id_lrunf, snow_lrunf+subs_lrunf, tile%diag) call send_tile_data(id_lrunfs, snow_lrunf, tile%diag) call send_tile_data(id_lrunfg, subs_lrunf, tile%diag) call send_tile_data(id_hlrunf, snow_hlrunf+subs_hlrunf, tile%diag) call send_tile_data(id_hlrunfs, snow_hlrunf, tile%diag) call send_tile_data(id_hlrunfg, subs_hlrunf, tile%diag) call send_tile_data(id_frunf, snow_frunf, tile%diag) call send_tile_data(id_frunfs, snow_frunf, tile%diag) call send_tile_data(id_hfrunf, snow_hfrunf, tile%diag) call send_tile_data(id_hfrunfs, snow_hfrunf, tile%diag) call send_tile_data(id_melt, vegn_melt+snow_melt+subs_melt, tile%diag) call send_tile_data(id_meltv, vegn_melt, tile%diag) call send_tile_data(id_melts, snow_melt, tile%diag) call send_tile_data(id_meltg, subs_melt, tile%diag) call send_tile_data(id_fsw, vegn_fsw+snow_fsw+subs_fsw, tile%diag) call send_tile_data(id_fswv, vegn_fsw, tile%diag) call send_tile_data(id_fsws, snow_fsw, tile%diag) call send_tile_data(id_fswg, subs_fsw, tile%diag) call send_tile_data(id_flw, vegn_flw+snow_flw+subs_flw, tile%diag) call send_tile_data(id_flwv, vegn_flw, tile%diag) call send_tile_data(id_flws, snow_flw, tile%diag) call send_tile_data(id_flwg, subs_flw, tile%diag) call send_tile_data(id_sens, land_sens, tile%diag) call send_tile_data(id_sensv, vegn_sens, tile%diag) call send_tile_data(id_senss, snow_sens, tile%diag) call send_tile_data(id_sensg, subs_sens, tile%diag) call send_tile_data(id_e_res_1, tile%e_res_1, tile%diag) call send_tile_data(id_e_res_2, tile%e_res_2, tile%diag) call send_tile_data(id_con_g_h, con_g_h, tile%diag) call send_tile_data(id_transp, vegn_uptk, tile%diag) call send_tile_data(id_wroff, snow_lrunf+subs_lrunf, tile%diag) call send_tile_data(id_sroff, snow_frunf, tile%diag) call send_tile_data(id_htransp, cpw*vegn_uptk*(vegn_T-tfreeze), tile%diag) call send_tile_data(id_huptake, clw*vegn_uptk*(soil_uptake_T-tfreeze), & tile%diag) call send_tile_data(id_hroff, snow_hlrunf+subs_hlrunf+snow_hfrunf, & tile%diag) call send_tile_data(id_gsnow, subs_G, tile%diag) call send_tile_data(id_gequil, subs_G2, tile%diag) call send_tile_data(id_grnd_flux, grnd_flux, tile%diag) call send_tile_data(id_soil_water_supply, soil_water_supply, tile%diag) if(grnd_E_max.lt.0.5*HUGE(grnd_E_Max)) & call send_tile_data(id_levapg_max, grnd_E_max, tile%diag) call send_tile_data(id_qco2, tile%cana%prog%co2, tile%diag) call send_tile_data(id_qco2_dvmr,& tile%cana%prog%co2*mol_air/mol_co2/(1-tile%cana%prog%q), tile%diag) call send_tile_data(id_fco2, vegn_fco2*mol_C/mol_co2, tile%diag) call send_tile_data(id_swdn_dir, ISa_dn_dir, tile%diag) call send_tile_data(id_swdn_dif, ISa_dn_dif, tile%diag) call send_tile_data(id_swup_dir, ISa_dn_dir*tile%land_refl_dir, tile%diag) call send_tile_data(id_swup_dif, ISa_dn_dif*tile%land_refl_dif, tile%diag) call send_tile_data(id_lwdn, ILa_dn, tile%diag) call send_tile_data(id_subs_emis,surf_emis_lw, tile%diag) if (id_total_C > 0) & call send_tile_data(id_total_C, land_tile_carbon(tile), tile%diag) end subroutine update_land_model_fast_0d ! ============================================================================ subroutine update_land_model_slow ( cplr2land, land2cplr ) type(atmos_land_boundary_type), intent(inout) :: cplr2land type(land_data_type) , intent(inout) :: land2cplr ! ---- local vars integer :: i,j,k type(land_tile_type), pointer :: tile type(land_tile_enum_type) :: ce, te call mpp_clock_begin(landClock) call mpp_clock_begin(landSlowClock) call land_transitions( lnd%time ) call update_vegn_slow( ) ! send the accumulated diagnostics to the output call dump_tile_diag_fields(lnd%tile_map, lnd%time) ! land_transitions may have changed the number of tiles per grid cell: reallocate ! boundary conditions, if necessary call realloc_cplr2land( cplr2land ) call realloc_land2cplr( land2cplr ) ! set the land mask to FALSE everywhere -- update_land_bc_fast will set it to ! true where necessary land2cplr%mask = .FALSE. land2cplr%tile_size = 0.0 ! get the current state of the land boundary for the coupler ce = first_elmt(lnd%tile_map, & is=lbound(cplr2land%t_flux,1), & js=lbound(cplr2land%t_flux,2) ) te = tail_elmt(lnd%tile_map) do while(ce /= te) ! calculate indices of the current tile in the input arrays; ! assume all the cplr2land components have the same lbounds call get_elmt_indices(ce,i,j,k) ! set this point coordinates as current for debug output call set_current_point(i,j,k) ! get pointer to current tile tile => current_tile(ce) ! advance enumerator to the next tile ce=next_elmt(ce) call update_land_bc_fast (tile, i,j,k, land2cplr, is_init=.true.) enddo call update_land_bc_slow( land2cplr ) call mpp_clock_end(landClock) call mpp_clock_end(landSlowClock) end subroutine update_land_model_slow ! ============================================================================ ! solve for surface temperature. ensure that melt does not exceed available ! snow or soil ice (which would create energy-balance problems). also ensure ! that melt and temperature are consistent and that evaporation from liquid ! soil water does not exceed exfiltration rate limit, if one is supplied. ! because the possible combinations of active constraints has multiplied ! greatly, we do not allow phase change for any surface (i.e., soil) at which ! we might apply a constraint on Eg. subroutine land_surface_energy_balance ( & ! surface parameters grnd_T, & ! ground temperature grnd_liq, grnd_ice, & ! amount of water available for freeze or melt on the surface, kg/m2 grnd_latent, & ! specific heat of vaporization for the ground grnd_Tf, & ! ground freezing temperature grnd_E_min, & ! Eg floor of 0 if condensation is prohibited grnd_E_max, & ! exfiltration rate limit, kg/(m2 s) ! components of the ground energy balance linearization. Note that those ! are full derivatives, which include the response of the other surface ! scheme parameters to the change in ground temperature. fswg, & ! net short-wave flwg0, DflwgDTg, DflwgDpsig, & ! net long-wave Hg0, DHgDTg, DHgDpsig, & ! sensible heat Eg0, DEgDTg, DEgDpsig, & ! latent heat G0, DGDTg, & ! sub-surface heat delta_Tg, & ! surface temperature change for the time step delta_psig, & Mg_imp ) ! implicit melt, kg/m2 real, intent(in) :: & grnd_T, & ! ground temperature grnd_liq, grnd_ice, & ! amount of water available for freeze or melt on the surface grnd_latent, & ! specific heat of vaporization for the ground grnd_Tf, & ! ground freezing temperature grnd_E_min, & ! Eg floor of 0 if condensation is prohibited grnd_E_max, & ! exfiltration rate limit, kg/(m2 s) fswg, & ! net short-wave flwg0, DflwgDTg, DflwgDpsig, & ! net long-wave Hg0, DHgDTg, DHgDpsig, & ! sensible heat Eg0, DEgDTg, DEgDpsig, & ! latent heat G0, DGDTg ! sub-surface heat real, intent(out) :: & delta_Tg, & ! change in surface temperature delta_psig, & ! change in surface soil-water matric head Mg_imp ! mass of surface ice melted (or water frozen) during the ! time step, kg/m2 real :: grnd_B ! surface energy balance real :: grnd_DBDTg ! full derivative of grnd_B w.r.t. surface temperature real :: grnd_DBDpsig ! full derivative of grnd_B w.r.t. surface soil-water matric head real :: grnd_E_force, Eg_trial, Eg_check, determinant grnd_B = fswg + flwg0 - Hg0 - grnd_latent*Eg0 - G0 grnd_DBDTg = DflwgDTg - DHgDTg - grnd_latent*DEgDTg - DGDTg grnd_DBDpsig = DflwgDpsig - DHgDpsig - grnd_latent*DEgDpsig if(is_watch_point())then write(*,*)'#### ground balance input' __DEBUG1__(grnd_T) __DEBUG2__(grnd_liq, grnd_ice) __DEBUG1__(grnd_latent) __DEBUG1__(grnd_Tf) __DEBUG1__(grnd_E_min) __DEBUG1__(grnd_E_max) __DEBUG1__(fswg) __DEBUG3__(flwg0, DflwgDTg,DflwgDpsig) __DEBUG3__(Hg0, DHgDTg, DHgDpsig ) __DEBUG3__(Eg0, DEgDTg, DEgDpsig ) __DEBUG2__(G0, DGDTg) write(*,*)'#### end of ground balance input' __DEBUG3__(grnd_B, grnd_DBDTg, grnd_DBDpsig) endif ! determine the ground temperature change under the assumptions that ! (1) no phase change occurs at surface (always true for soil now), and ! (2) delta_psig is zero (always true for snow, lake, glacier now) delta_Tg = - grnd_B/grnd_DBDTg delta_psig = 0 ! calculate phase change on the ground, if necessary if (grnd_ice>0.and.grnd_T+delta_Tg>grnd_Tf) then ! melt > 0 Mg_imp = min(grnd_ice, grnd_DBDTg*(grnd_Tf-grnd_T-delta_Tg)*delta_time/hlf) elseif (grnd_liq>0.and.grnd_T+delta_Tg0.01) snow_area_rad = 1 endif call cana_radiation( lm2, & subs_refl_dir, subs_refl_dif, subs_refl_lw, & snow_refl_dir, snow_refl_dif, snow_refl_lw, & snow_area_rad, & vegn_refl_dir, vegn_refl_dif, vegn_tran_dir, vegn_tran_dif, & vegn_tran_dir_dir, vegn_refl_lw, vegn_tran_lw, & vegn_cover, & tile%Sg_dir, tile%Sg_dif, tile%Sv_dir, tile%Sv_dif, & tile%land_refl_dir, tile%land_refl_dif ) call cana_roughness( lm2, & subs_z0m, subs_z0s, & snow_z0m, snow_z0s, snow_area, & vegn_cover, vegn_height, vegn_lai, vegn_sai, & tile%land_d, tile%land_z0m, tile%land_z0s, & associated(tile%lake).or.associated(tile%glac)) if(is_watch_point()) then write(*,*) '#### update_land_bc_fast ### checkpoint 2 ####' write(*,*) 'Sg_dir', tile%Sg_dir write(*,*) 'Sg_dif', tile%Sg_dif write(*,*) 'Sv_dir', tile%Sv_dir write(*,*) 'Sv_dif', tile%Sv_dif write(*,*) 'land_albedo_dir', tile%land_refl_dir write(*,*) 'land_albedo_dif', tile%land_refl_dif write(*,*) 'land_z0m', tile%land_z0m write(*,*) '#### update_land_bc_fast ### end of checkpoint 2 ####' endif land2cplr%t_surf (i,j,k) = tfreeze land2cplr%t_ca (i,j,k) = tfreeze land2cplr%tr (i,j,k, lnd%isphum) = 0.0 land2cplr%albedo (i,j,k) = 0.0 land2cplr%albedo_vis_dir (i,j,k) = 0.0 land2cplr%albedo_nir_dir (i,j,k) = 0.0 land2cplr%albedo_vis_dif (i,j,k) = 0.0 land2cplr%albedo_nir_dif (i,j,k) = 0.0 land2cplr%rough_mom (i,j,k) = 0.1 land2cplr%rough_heat (i,j,k) = 0.1 ! Calculate radiative surface temperature. lwup can't be calculated here ! based on the available temperatures because it's a result of the implicit ! time step: lwup = lwup0 + DlwupDTg*delta_Tg + ..., so we have to carry it ! from the update_land_fast ! Consequence: since update_landbc_fast is called once at the very beginning of ! every run (before update_land_fast is called) lwup from the previous step ! must be stored in the in the restart for reproducibility land2cplr%t_surf(i,j,k) = ( tile%lwup/stefan ) ** 0.25 if (associated(tile%glac)) call glac_get_sfc_temp(tile%glac, grnd_T) if (associated(tile%lake)) call lake_get_sfc_temp(tile%lake, grnd_T) if (associated(tile%soil)) call soil_get_sfc_temp(tile%soil, grnd_T) if (snow_area > 0) call snow_get_sfc_temp(tile%snow, grnd_T) ! set the boundary conditions for the flux exchange land2cplr%mask (i,j,k) = .TRUE. land2cplr%tile_size (i,j,k) = tile%frac call cana_state ( tile%cana, land2cplr%t_ca(i,j,k), & land2cplr%tr(i,j,k,lnd%isphum), cana_co2) ! land2cplr%t_ca (i,j,k) = tile%cana%prog%T ! land2cplr%tr (i,j,k,lnd%isphum) = tile%cana%prog%q if(lnd%ico2/=NO_TRACER) then land2cplr%tr(i,j,k,lnd%ico2) = cana_co2 endif land2cplr%albedo_vis_dir (i,j,k) = tile%land_refl_dir(BAND_VIS) land2cplr%albedo_nir_dir (i,j,k) = tile%land_refl_dir(BAND_NIR) land2cplr%albedo_vis_dif (i,j,k) = tile%land_refl_dif(BAND_VIS) land2cplr%albedo_nir_dif (i,j,k) = tile%land_refl_dif(BAND_NIR) land2cplr%albedo (i,j,k) = SUM(tile%land_refl_dir + tile%land_refl_dif)/4 ! incorrect, replace with proper weighting later land2cplr%rough_mom (i,j,k) = tile%land_z0m land2cplr%rough_heat (i,j,k) = tile%land_z0s if(is_watch_point()) then write(*,*)'#### update_land_bc_fast ### output ####' write(*,*)'land2cplr%mask',land2cplr%mask(i,j,k) write(*,*)'land2cplr%tile_size',land2cplr%tile_size(i,j,k) write(*,*)'land2cplr%t_surf',land2cplr%t_surf(i,j,k) write(*,*)'land2cplr%t_ca',land2cplr%t_ca(i,j,k) write(*,*)'land2cplr%albedo',land2cplr%albedo(i,j,k) write(*,*)'land2cplr%rough_mom',land2cplr%rough_mom(i,j,k) write(*,*)'land2cplr%rough_heat',land2cplr%rough_heat(i,j,k) write(*,*)'land2cplr%tr',land2cplr%tr(i,j,k,:) write(*,*)'#### update_land_bc_fast ### end of output ####' endif ! ---- diagnostic section call send_tile_data(id_vegn_cover, vegn_cover, tile%diag) call send_tile_data(id_cosz, cosz, tile%diag) call send_tile_data(id_albedo_dir, tile%land_refl_dir, tile%diag) call send_tile_data(id_albedo_dif, tile%land_refl_dif, tile%diag) call send_tile_data(id_vegn_refl_dir, vegn_refl_dir, tile%diag) call send_tile_data(id_vegn_refl_dif, vegn_refl_dif, tile%diag) call send_tile_data(id_vegn_refl_lw, vegn_refl_lw, tile%diag) call send_tile_data(id_vegn_tran_dir, vegn_tran_dir_dir, tile%diag) call send_tile_data(id_vegn_tran_dif, vegn_tran_dif, tile%diag) call send_tile_data(id_vegn_tran_lw, vegn_tran_lw, tile%diag) call send_tile_data(id_vegn_sctr_dir, vegn_tran_dir, tile%diag) call send_tile_data(id_subs_refl_dir, subs_refl_dir, tile%diag) call send_tile_data(id_subs_refl_dif, subs_refl_dif, tile%diag) call send_tile_data(id_grnd_T, grnd_T, tile%diag) ! --- debug section call check_temp_range(land2cplr%t_ca(i,j,k),'update_land_bc_fast','T_ca',lnd%time) end subroutine update_land_bc_fast ! ============================================================================ subroutine update_land_bc_slow (land2cplr) type(land_data_type), intent(inout) :: land2cplr ! ---- local vars integer :: i,j,k,face ! coordinates of the watch point, for debug printout call update_topo_rough(land2cplr%rough_scale) where (land2cplr%mask) & land2cplr%rough_scale = max(land2cplr%rough_mom,land2cplr%rough_scale) call get_watch_point(i,j,k,face) if ( lnd%face==face.and. & lnd%is<=i.and.i<=lnd%ie.and. & lnd%js<=j.and.j<=lnd%je.and. & k<=size(land2cplr%rough_scale,3)) then write(*,*)'#### update_land_bc_slow ### output ####' write(*,*)'land2cplr%rough_scale',land2cplr%rough_scale(i,j,k) write(*,*)'#### update_land_bc_slow ### end of output ####' endif end subroutine update_land_bc_slow ! ============================================================================ subroutine Lnd_stock_pe(bnd,index,value) type(land_data_type), intent(in) :: bnd integer , intent(in) :: index real , intent(out) :: value ! Domain water (Kg) or heat (Joules) integer :: i,j,n type(land_tile_enum_type) :: ce,te type(land_tile_type), pointer :: tile character(len=128) :: message integer :: is,ie,js,je real :: area_factor, river_value ! *_twd are tile water densities (kg water per m2 of tile) real twd_gas_cana, twd_liq_glac, twd_sol_glac, & twd_liq_lake, twd_sol_lake, twd_liq_soil, twd_sol_soil, & twd_liq_snow, twd_sol_snow, twd_liq_vegn, twd_sol_vegn ! *_gcwd are grid-cell water densities (kg water per m2 of land in grid cell) real gcwd_cana, gcwd_glac, gcwd_lake, gcwd_soil, gcwd_snow, gcwd_vegn ! v_* are global masses of water real v_cana, v_glac, v_lake, v_soil, v_snow, v_vegn real cana_q, a_globe value = 0.0 v_cana = 0. v_glac = 0. v_lake = 0. v_soil = 0. v_snow = 0. v_vegn = 0. if(.not.bnd%pe) return is = lnd%is ie = lnd%ie js = lnd%js je = lnd%je ! The following is a dirty getaround if(lnd%cellarea(is,js) < 1.0) then area_factor = 4*pi*radius**2 ! lnd%area is fraction of globe else area_factor = 1.0 ! lnd%area is actual area (m**2) endif select case(index) case(ISTOCK_WATER) do j = js, je do i = is, ie ce = first_elmt(lnd%tile_map(i,j)) te = tail_elmt (lnd%tile_map(i,j)) gcwd_cana = 0.0; gcwd_glac = 0.0; gcwd_lake = 0.0 gcwd_soil = 0.0; gcwd_snow = 0.0; gcwd_vegn = 0.0 do while(ce /= te) tile => current_tile(ce) twd_gas_cana = 0.0 twd_liq_glac = 0.0 ; twd_sol_glac = 0.0 twd_liq_lake = 0.0 ; twd_sol_lake = 0.0 twd_liq_soil = 0.0 ; twd_sol_soil = 0.0 twd_liq_snow = 0.0 ; twd_sol_snow = 0.0 twd_liq_vegn = 0.0 ; twd_sol_vegn = 0.0 if(associated(tile%cana)) then call cana_state ( tile%cana, cana_q=cana_q ) twd_gas_cana = canopy_air_mass*cana_q endif if(associated(tile%glac)) & call glac_tile_stock_pe(tile%glac, twd_liq_glac, twd_sol_glac) if(associated(tile%lake)) & call lake_tile_stock_pe(tile%lake, twd_liq_lake, twd_sol_lake) if(associated(tile%soil)) & call soil_tile_stock_pe(tile%soil, twd_liq_soil, twd_sol_soil) if(associated(tile%snow)) & call snow_tile_stock_pe(tile%snow, twd_liq_snow, twd_sol_snow) if(associated(tile%vegn)) & call vegn_tile_stock_pe(tile%vegn, twd_liq_vegn, twd_sol_vegn) gcwd_cana = gcwd_cana + twd_gas_cana * tile%frac gcwd_glac = gcwd_glac + (twd_liq_glac + twd_sol_glac) * tile%frac gcwd_lake = gcwd_lake + (twd_liq_lake + twd_sol_lake) * tile%frac gcwd_soil = gcwd_soil + (twd_liq_soil + twd_sol_soil) * tile%frac gcwd_snow = gcwd_snow + (twd_liq_snow + twd_sol_snow) * tile%frac gcwd_vegn = gcwd_vegn + (twd_liq_vegn + twd_sol_vegn) * tile%frac ce=next_elmt(ce) enddo v_cana = v_cana + gcwd_cana * lnd%area(i,j)*area_factor v_glac = v_glac + gcwd_glac * lnd%area(i,j)*area_factor v_lake = v_lake + gcwd_lake * lnd%area(i,j)*area_factor v_soil = v_soil + gcwd_soil * lnd%area(i,j)*area_factor v_snow = v_snow + gcwd_snow * lnd%area(i,j)*area_factor v_vegn = v_vegn + gcwd_vegn * lnd%area(i,j)*area_factor enddo enddo value = v_cana + v_glac + v_lake + v_soil + v_snow + v_vegn a_globe = 4. * pi * radius**2 case(ISTOCK_HEAT) if(.not.stock_warning_issued) then call error_mesg('Lnd_stock_pe','Heat stock not yet implemented',NOTE) stock_warning_issued = .true. endif ! do j = js, je ! do i = is, ie ! ce = first_elmt(lnd%tile_map(i,j)) ! te = tail_elmt (lnd%tile_map(i,j)) ! grid_cell_heat_density = 0.0 ! do while(ce /= te) ! tile => current_tile(ce) ! tile_heat_density = 0.0 ! if(associated(tile%soil)) then ! do n=1, size(tile%soil%prog) ! tile_heat_density = tile_heat_density + (tile%soil%prog(n)%T-tfreeze)* & ! (tile%soil%heat_capacity_dry(n)*dz(n) + & ! clw*tile%soil%prog(n)%wl + & ! csw*tile%soil%prog(n)%ws) ! enddo ! tile_heat_density = tile_heat_density + clw*soil%prog(1)%groundwater*(soil%prog(1)%groundwater_T-tfreeze) ! Why is this outside n loop? ! endif ! grid_cell_heat_density = grid_cell_heat_density + tile_heat_density * tile%frac ! ce=next_elmt(ce) ! enddo ! grid_cell_heat = grid_cell_heat_density * lnd%area(i,j)*area_factor ! value = value + grid_cell_heat ! enddo ! enddo case(ISTOCK_SALT) if(.not.stock_warning_issued) then call error_mesg('Lnd_stock_pe','Salt stock not yet implemented',NOTE) stock_warning_issued = .true. endif case default write(message,'(i2,a,i2,a,i2,a,i2,a)') & index,' is an invalid stock index. Must be ISTOCK_WATER or ISTOCK_HEAT or ISTOCK_SALT (',ISTOCK_WATER,' or ',ISTOCK_HEAT,' or ', ISTOCK_SALT,')' call error_mesg('Lnd_stock_pe',message,FATAL) end select call river_stock_pe(index, river_value) value = value + river_value if (index.eq.ISTOCK_WATER.and.give_stock_details) then call mpp_sum(river_value, pelist=lnd%pelist) call mpp_sum(v_cana, pelist=lnd%pelist) call mpp_sum(v_glac, pelist=lnd%pelist) call mpp_sum(v_lake, pelist=lnd%pelist) call mpp_sum(v_soil, pelist=lnd%pelist) call mpp_sum(v_snow, pelist=lnd%pelist) call mpp_sum(v_vegn, pelist=lnd%pelist) write (message,'(a,f10.5)') 'total land storage:',v_cana/a_globe+v_glac/a_globe+ & v_lake/a_globe+v_soil/a_globe+v_snow/a_globe+v_vegn/a_globe+river_value/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '...canopy air:',v_cana/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '......glacier:',v_glac/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '.........lake:',v_lake/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '.........soil:',v_soil/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '.........snow:',v_snow/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '.........vegn:',v_vegn/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) write (message,'(a,7f10.5)') '.......rivers:',river_value/a_globe call error_mesg('Lnd_stock_pe',message,NOTE) endif end subroutine Lnd_stock_pe ! ============================================================================ ! initialize horizontal axes for land grid so that all sub-modules can use them, ! instead of creating their own subroutine land_diag_init(clonb, clatb, clon, clat, time, domain, & id_lon, id_lat, id_band) real, intent(in) :: & clonb(:), clatb(:), & ! longitudes and latitudes of grid cells vertices, ! specified for the global grid clon(:), clat(:) ! lon and lat of respective grid cell centers type(time_type), intent(in) :: time ! initial time for diagnostic fields type(domain2d), intent(in) :: domain integer, intent(out) :: & id_lon, id_lat, id_band ! IDs of respective diag. manager axes ! ---- local vars ---------------------------------------------------------- integer :: id_lonb, id_latb ! IDs for cell boundaries integer :: nlon, nlat ! sizes of respective axes integer :: axes(2) ! array of axes for 2-D fields integer :: i nlon = size(clon) nlat = size(clat) if(mpp_get_ntile_count(lnd%domain)==1) then ! grid has just one tile, so we assume that the grid is regular lat-lon ! define longitude axes and its edges id_lonb = diag_axis_init ( & 'lonb', clonb, 'degrees_E', 'X', 'longitude edges', & set_name='land', domain2=domain ) id_lon = diag_axis_init ( & 'lon', clon, 'degrees_E', 'X', & 'longitude', set_name='land', edges=id_lonb, domain2=domain ) ! define latitude axes and its edges id_latb = diag_axis_init ( & 'latb', clatb, 'degrees_N', 'Y', 'latitude edges', & set_name='land', domain2=domain ) id_lat = diag_axis_init ( & 'lat', clat, 'degrees_N', 'Y', & 'latitude', set_name='land', edges=id_latb, domain2=domain ) else id_lon = diag_axis_init ( 'grid_xt', (/(real(i),i=1,nlon)/), 'degrees_E', 'X', & 'T-cell longitude', set_name='land', domain2=domain, aux='geolon_t' ) id_lat = diag_axis_init ( 'grid_yt', (/(real(i),i=1,nlat)/), 'degrees_N', 'Y', & 'T-cell latitude', set_name='land', domain2=domain, aux='geolat_t' ) endif id_band = diag_axis_init ( & 'band', (/1.0,2.0/), 'unitless', 'Z', & 'spectral band', set_name='land' ) ! set up an array of axes, for convenience axes = (/id_lon, id_lat/) ! register auxilary coordinate variables id_geolon_t = register_static_field ( module_name, 'geolon_t', axes, & 'longitude of grid cell centers', 'degrees_E', missing_value = -1.0e+20 ) id_geolat_t = register_static_field ( module_name, 'geolat_t', axes, & 'latitude of grid cell centers', 'degrees_N', missing_value = -1.0e+20 ) ! register static diagnostic fields id_cellarea = register_static_field ( module_name, 'cell_area', axes, & 'total area in grid cell', 'm2', missing_value=-1.0 ) id_landarea = register_static_field ( module_name, 'land_area', axes, & 'land area in grid cell', 'm2', missing_value=-1.0 ) id_landfrac = register_static_field ( module_name, 'land_frac', axes, & 'fraction of land in grid cell','unitless', missing_value=-1.0 ) id_no_riv = register_static_field ( module_name, 'no_riv', axes, & 'indicator of land without rivers','unitless', missing_value=-1.0 ) ! register regular (dynamic) diagnostic fields id_VWS = register_tiled_diag_field ( module_name, 'VWS', axes, time, & 'vapor storage on land', 'kg/m2', missing_value=-1.0e+20 ) id_VWSc = register_tiled_diag_field ( module_name, 'VWSc', axes, time, & 'vapor mass in canopy air', 'kg/m2', missing_value=-1.0e+20 ) id_LWS = register_tiled_diag_field ( module_name, 'LWS', axes, time, & 'liquid storage on land', 'kg/m2', missing_value=-1.0e+20 ) id_LWSv = register_tiled_diag_field ( module_name, 'LWSv', axes, time, & 'liquid interception storage', 'kg/m2', missing_value=-1.0e+20 ) id_LWSs = register_tiled_diag_field ( module_name, 'LWSs', axes, time, & 'liquid storage in snowpack', 'kg/m2', missing_value=-1.0e+20 ) id_LWSg = register_tiled_diag_field ( module_name, 'LWSg', axes, time, & 'liquid ground storage', 'kg/m2', missing_value=-1.0e+20 ) id_FWS = register_tiled_diag_field ( module_name, 'FWS', axes, time, & 'frozen storage on land', 'kg/m2', missing_value=-1.0e+20 ) id_FWSv = register_tiled_diag_field ( module_name, 'FWSv', axes, time, & 'frozen interception storage', 'kg/m2', missing_value=-1.0e+20 ) id_FWSs = register_tiled_diag_field ( module_name, 'FWSs', axes, time, & 'frozen storage in snowpack', 'kg/m2', missing_value=-1.0e+20 ) id_FWSg = register_tiled_diag_field ( module_name, 'FWSg', axes, time, & 'frozen ground storage', 'kg/m2', missing_value=-1.0e+20 ) id_HS = register_tiled_diag_field ( module_name, 'HS', axes, time, & 'land heat storage', 'J/m2', missing_value=-1.0e+20 ) id_HSv = register_tiled_diag_field ( module_name, 'HSv', axes, time, & 'interception heat storage', 'J/m2', missing_value=-1.0e+20 ) id_HSs = register_tiled_diag_field ( module_name, 'HSs', axes, time, & 'heat storage in snowpack', 'J/m2', missing_value=-1.0e+20 ) id_HSg = register_tiled_diag_field ( module_name, 'HSg', axes, time, & 'ground heat storage', 'J/m2', missing_value=-1.0e+20 ) id_HSc = register_tiled_diag_field ( module_name, 'HSc', axes, time, & 'canopy-air heat storage', 'J/m2', missing_value=-1.0e+20 ) id_dis_liq = register_diag_field ( module_name, 'dis_liq', axes, & time, 'liquid discharge to ocean', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_dis_ice = register_diag_field ( module_name, 'dis_ice', axes, & time, 'ice discharge to ocean', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_dis_heat = register_diag_field ( module_name, 'dis_heat', axes, & time, 'heat of mass discharge to ocean', 'W/m2', missing_value=-1.0e+20 ) id_dis_sink = register_diag_field ( module_name, 'dis_sink', axes, & time, 'burial rate of small/negative discharge', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_precip = register_tiled_diag_field ( module_name, 'precip', axes, time, & 'precipitation rate', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_hprec = register_tiled_diag_field ( module_name, 'hprec', axes, time, & 'sensible heat of precipitation', 'W/m2', missing_value=-1.0e+20 ) id_lprec = register_tiled_diag_field ( module_name, 'lprec_l', axes, time, & 'rainfall to land', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_lprecv = register_tiled_diag_field ( module_name, 'lprecv', axes, time, & 'net rainfall to vegetation', 'kg/(m2 s)', missing_value=-1.0e+20) id_lprecs = register_tiled_diag_field ( module_name, 'lprecs', axes, time, & 'rainfall to snow, minus drainage', 'kg/(m2 s)', missing_value=-1.0e+20) id_lprecg = register_tiled_diag_field ( module_name, 'lprecg', axes, time, & 'effective rainfall to ground sfc', 'kg/(m2 s)', missing_value=-1.0e+20) id_hlprec = register_tiled_diag_field ( module_name, 'hlprec', axes, time, & 'total liq precipitation heat', 'W/m2', missing_value=-1.0e+20) id_hlprecv = register_tiled_diag_field ( module_name, 'hlprecv', axes, time, & 'net liq heat to vegetation', 'W/m2', missing_value=-1.0e+20) id_hlprecs = register_tiled_diag_field ( module_name, 'hlprecs', axes, time, & 'net liq heat to snow', 'W/m2', missing_value=-1.0e+20) id_hlprecg = register_tiled_diag_field ( module_name, 'hlprecg', axes, time, & 'net liq heat to ground sfc', 'W/m2', missing_value=-1.0e+20) id_fprec = register_tiled_diag_field ( module_name, 'fprec_l', axes, time, & 'snowfall to land', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_fprecv = register_tiled_diag_field ( module_name, 'fprecv', axes, time, & 'net snowfall to vegetation', 'kg/(m2 s)', missing_value=-1.0e+20) id_fprecs = register_tiled_diag_field ( module_name, 'fprecs', axes, time, & 'effective snowfall to snowpack', 'kg/(m2 s)', missing_value=-1.0e+20) id_hfprec = register_tiled_diag_field ( module_name, 'hfprec', axes, time, & 'sens heat of snowfall', 'W/m2', missing_value=-1.0e+20) id_hfprecv = register_tiled_diag_field ( module_name, 'hfprecv', axes, time, & 'net sol heat to vegetation', 'W/m2', missing_value=-1.0e+20) id_hfprecs = register_tiled_diag_field ( module_name, 'hfprecs', axes, time, & 'net sol heat to snow', 'W/m2', missing_value=-1.0e+20) id_evap = register_tiled_diag_field ( module_name, 'evap', axes, time, & 'vapor flux up from land', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_hevap = register_tiled_diag_field ( module_name, 'hevap', axes, time, & 'sensible heat of evap', 'W/m2', missing_value=-1.0e+20 ) id_levap = register_tiled_diag_field ( module_name, 'levap', axes, time, & 'vapor flux from all liq (inc Tr)', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_levapv = register_tiled_diag_field ( module_name, 'levapv', axes, time, & 'vapor flux leaving intercepted liquid', 'kg/(m2 s)', missing_value=-1.0e+20) id_levaps = register_tiled_diag_field ( module_name, 'levaps', axes, time, & 'vapor flux leaving snow liquid', 'kg/(m2 s)', missing_value=-1.0e+20) id_levapg = register_tiled_diag_field ( module_name, 'levapg', axes, time, & 'vapor flux from ground liquid', 'kg/(m2 s)', missing_value=-1.0e+20) id_hlevap = register_tiled_diag_field ( module_name, 'hlevap', axes, time, & 'vapor flux heat from liq source', 'W/m2', missing_value=-1.0e+20) id_hlevapv = register_tiled_diag_field ( module_name, 'hlevapv', axes, time, & 'vapor heat from liq interc', 'W/m2', missing_value=-1.0e+20) id_hlevaps = register_tiled_diag_field ( module_name, 'hlevaps', axes, time, & 'vapor heat from snow liq', 'W/m2', missing_value=-1.0e+20) id_hlevapg = register_tiled_diag_field ( module_name, 'hlevapg', axes, time, & 'vapor heat from ground liq', 'W/m2', missing_value=-1.0e+20) id_fevap = register_tiled_diag_field ( module_name, 'fevap', axes, time, & 'vapor flux from all ice', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_fevapv = register_tiled_diag_field ( module_name, 'fevapv', axes, time, & 'vapor flux leaving vegn ice', 'kg/(m2 s)', missing_value=-1.0e+20) id_fevaps = register_tiled_diag_field ( module_name, 'fevaps', axes, time, & 'vapor flux leaving snow ice', 'kg/(m2 s)', missing_value=-1.0e+20) id_fevapg = register_tiled_diag_field ( module_name, 'fevapg', axes, time, & 'vapor flux leaving ground ice', 'kg/(m2 s)', missing_value=-1.0e+20) id_hfevap = register_tiled_diag_field ( module_name, 'hfevap', axes, time, & 'vapor flux heat from solid source', 'W/m2', missing_value=-1.0e+20) id_hfevapv = register_tiled_diag_field ( module_name, 'hfevapv', axes, time, & 'vapor heat from sol interc', 'W/m2', missing_value=-1.0e+20) id_hfevaps = register_tiled_diag_field ( module_name, 'hfevaps', axes, time, & 'vapor heat from snow sol', 'W/m2', missing_value=-1.0e+20) id_hfevapg = register_tiled_diag_field ( module_name, 'hfevapg', axes, time, & 'vapor heat from ground sol', 'W/m2', missing_value=-1.0e+20) id_runf = register_tiled_diag_field ( module_name, 'runf', axes, time, & 'total runoff', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_hrunf = register_tiled_diag_field ( module_name, 'hrunf', axes, time, & 'sensible heat of total runoff', 'W/m2', missing_value=-1.0e+20 ) id_lrunf = register_tiled_diag_field ( module_name, 'lrunf', axes, time, & 'total rate of liq runoff', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_lrunfs = register_tiled_diag_field ( module_name, 'lrunfs', axes, time, & 'rate of liq runoff via calving', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_lrunfg = register_tiled_diag_field ( module_name, 'lrunfg', axes, time, & 'rate of liq runoff, ground', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_hlrunf = register_tiled_diag_field ( module_name, 'hlrunf', axes, time, & 'heat of liq runoff', 'W/m2', missing_value=-1.0e+20 ) id_hlrunfs = register_tiled_diag_field ( module_name, 'hlrunfs', axes, time, & 'heat of liq runoff from snow pack', 'W/m2', missing_value=-1.0e+20 ) id_hlrunfg = register_tiled_diag_field ( module_name, 'hlrunfg', axes, time, & 'heat of liq surface runoff', 'W/m2', missing_value=-1.0e+20 ) id_frunf = register_tiled_diag_field ( module_name, 'frunf', axes, time, & 'total rate of solid runoff', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_frunfs = register_tiled_diag_field ( module_name, 'frunfs', axes, time, & 'rate of solid calving', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_hfrunf = register_tiled_diag_field ( module_name, 'hfrunf', axes, time, & 'heat of total ice runoff', 'W/m2', missing_value=-1.0e+20 ) id_hfrunfs = register_tiled_diag_field ( module_name, 'hfrunfs', axes, time, & 'heat of sol snow runoff', 'W/m2', missing_value=-1.0e+20 ) id_melt = register_tiled_diag_field ( module_name, 'melt', axes, time, & 'total rate of melt', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_meltv = register_tiled_diag_field ( module_name, 'meltv', axes, time, & 'rate of melt, interception', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_melts = register_tiled_diag_field ( module_name, 'melts', axes, time, & 'rate of snow melt', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_meltg = register_tiled_diag_field ( module_name, 'meltg', axes, time, & 'rate of substrate thaw', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_fsw = register_tiled_diag_field ( module_name, 'fsw', axes, time, & 'net sw rad to land', 'W/m2', missing_value=-1.0e+20 ) id_fswv = register_tiled_diag_field ( module_name, 'fswv', axes, time, & 'net sw rad to vegetation', 'W/m2', missing_value=-1.0e+20 ) id_fsws = register_tiled_diag_field ( module_name, 'fsws', axes, time, & 'net sw rad to snow', 'W/m2', missing_value=-1.0e+20 ) id_fswg = register_tiled_diag_field ( module_name, 'fswg', axes, time, & 'net sw rad to ground', 'W/m2', missing_value=-1.0e+20 ) id_flw = register_tiled_diag_field ( module_name, 'flw', axes, time, & 'net lw rad to land', 'W/m2', missing_value=-1.0e+20 ) id_flwv = register_tiled_diag_field ( module_name, 'flwv', axes, time, & 'net lw rad to vegetation', 'W/m2', missing_value=-1.0e+20 ) id_flws = register_tiled_diag_field ( module_name, 'flws', axes, time, & 'net lw rad to snow', 'W/m2', missing_value=-1.0e+20 ) id_flwg = register_tiled_diag_field ( module_name, 'flwg', axes, time, & 'net lw rad to ground', 'W/m2', missing_value=-1.0e+20 ) id_cd_m = register_tiled_diag_field ( module_name, 'cd_m', axes, time, & 'drag coefficient for momentum', missing_value=-1e20) id_cd_t = register_tiled_diag_field ( module_name, 'cd_t', axes, time, & 'drag coefficient for heat and tracers', missing_value=-1e20) id_sens = register_tiled_diag_field ( module_name, 'sens', axes, time, & 'sens heat flux from land', 'W/m2', missing_value=-1.0e+20 ) id_sensv = register_tiled_diag_field ( module_name, 'sensv', axes, time, & 'sens heat flux from vegn', 'W/m2', missing_value=-1.0e+20 ) id_senss = register_tiled_diag_field ( module_name, 'senss', axes, time, & 'sens heat flux from snow', 'W/m2', missing_value=-1.0e+20 ) id_sensg = register_tiled_diag_field ( module_name, 'sensg', axes, time, & 'sens heat flux from ground', 'W/m2', missing_value=-1.0e+20 ) id_e_res_1 = register_tiled_diag_field ( module_name, 'e_res_1', axes, time, & 'canopy air energy residual due to nonlinearities', 'W/m2', missing_value=-1e20) id_e_res_2 = register_tiled_diag_field ( module_name, 'e_res_2', axes, time, & 'canopy energy residual due to nonlinearities', 'W/m2', missing_value=-1e20) id_frac = register_tiled_diag_field(module_name,'frac', axes,& time, 'fraction of land area', 'unitless', missing_value=-1.0, op=OP_SUM ) id_area = register_tiled_diag_field(module_name,'area', axes,& time, 'area in the grid cell', 'm2', missing_value=-1.0, op=OP_SUM ) id_ntiles = register_tiled_diag_field(module_name,'ntiles',axes, & time, 'number of tiles', 'unitless', missing_value=-1.0, op=OP_SUM) id_z0m = register_tiled_diag_field ( module_name, 'z0m', axes, time, & 'momentum roughness of land', 'm', missing_value=-1.0e+20 ) id_z0s = register_tiled_diag_field ( module_name, 'z0s', axes, time, & 'scalar roughness of land', 'm', missing_value=-1.0e+20 ) id_con_g_h = register_tiled_diag_field ( module_name, 'con_g_h', axes, time, & 'conductance for sensible heat between ground surface and canopy air', & 'm/s', missing_value=-1.0 ) id_transp = register_tiled_diag_field ( module_name, 'transp', axes, time, & 'transpiration; = uptake by roots', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_wroff = register_tiled_diag_field ( module_name, 'wroff', axes, time, & 'rate of liquid runoff to rivers', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_sroff = register_tiled_diag_field ( module_name, 'sroff', axes, time, & 'rate of solid runoff to rivers', 'kg/(m2 s)', missing_value=-1.0e+20 ) id_htransp = register_tiled_diag_field ( module_name, 'htransp', axes, time, & 'heat of transpired vapior', 'W/m2', missing_value=-1.0e+20 ) id_huptake = register_tiled_diag_field ( module_name, 'huptk', axes, time, & 'heat of soil water uptake', 'W/m2', missing_value=-1.0e+20 ) id_hroff = register_tiled_diag_field ( module_name, 'hroff', axes, time, & 'sensible heat of runoff', 'W/m2', missing_value=-1.0e+20 ) id_gsnow = register_tiled_diag_field ( module_name, 'gsnow', axes, time, & 'sens heat into ground from snow', 'W/m2', missing_value=-1.0e+20 ) call add_tiled_diag_field_alias ( id_gsnow, module_name, 'gflux', axes, time, & 'obsolete, please use "gsnow" instead', 'W/m2', missing_value=-1.0e+20 ) id_gequil = register_tiled_diag_field ( module_name, 'gequil', axes, time, & 'snow-subs equilibration flux', 'W/m2', missing_value=-1.0e+20 ) id_grnd_flux = register_tiled_diag_field ( module_name, 'grnd_flux', axes, time, & 'sensible heat into ground from surface', 'W/m2', missing_value=-1.0e+20 ) id_soil_water_supply = register_tiled_diag_field ( module_name, 'soil_water_supply', axes, time, & 'maximum rate of soil water supply to vegetation', 'kg/(m2 s)', missing_value=-1e20) id_levapg_max = register_tiled_diag_field ( module_name, 'Eg_max', axes, time, & 'soil_water limit on vapor flux from ground liquid', 'kg/(m2 s)', missing_value=-1.0e+20) id_water = register_tiled_diag_field ( module_name, 'water', axes, time, & 'column-integrated soil water', 'kg/m2', missing_value=-1.0e+20 ) id_snow = register_tiled_diag_field ( module_name, 'snow', axes, time, & 'column-integrated snow water', 'kg/m2', missing_value=-1.0e+20 ) id_Trad = register_tiled_diag_field ( module_name, 'Trad', axes, time, & 'radiative sfc temperature', 'degK', missing_value=-1.0e+20 ) id_Tca = register_tiled_diag_field ( module_name, 'Tca', axes, time, & 'canopy-air temperature', 'degK', missing_value=-1.0e+20 ) id_qca = register_tiled_diag_field ( module_name, 'qca', axes, time, & 'canopy-air specific humidity', 'kg/kg', missing_value=-1.0 ) id_qco2 = register_tiled_diag_field ( module_name, 'qco2', axes, time, & 'canopy-air CO2 moist mass mixing ratio', 'kg/kg', missing_value=-1.0 ) id_qco2_dvmr = register_tiled_diag_field ( module_name, 'qco2_dvmr', axes, time, & 'canopy-air CO2 dry volumetric mixing ratio', 'mol CO2/mol air', missing_value=-1.0 ) id_fco2 = register_tiled_diag_field ( module_name, 'fco2', axes, time, & 'flux of CO2 to canopy air', 'kg C/(m2 s)', missing_value=-1.0 ) id_swdn_dir = register_tiled_diag_field ( module_name, 'swdn_dir', (/id_lon,id_lat,id_band/), time, & 'downward direct short-wave radiation flux to the land surface', 'W/m2', missing_value=-999.0) id_swdn_dif = register_tiled_diag_field ( module_name, 'swdn_dif', (/id_lon,id_lat,id_band/), time, & 'downward diffuse short-wave radiation flux to the land surface', 'W/m2', missing_value=-999.0) id_swup_dir = register_tiled_diag_field ( module_name, 'swup_dir', (/id_lon,id_lat,id_band/), time, & 'direct short-wave radiation flux reflected by the land surface', 'W/m2', missing_value=-999.0) id_swup_dif = register_tiled_diag_field ( module_name, 'swup_dif', (/id_lon,id_lat,id_band/), time, & 'diffuse short-wave radiation flux reflected by the land surface', 'W/m2', missing_value=-999.0) id_lwdn = register_tiled_diag_field ( module_name, 'lwdn', axes, time, & 'downward long-wave radiation flux to the land surface', 'W/m2', missing_value=-999.0) id_vegn_cover = register_tiled_diag_field ( module_name, 'vegn_cover', axes, time, & 'fraction covered by vegetation', missing_value=-1.0 ) id_cosz = register_tiled_diag_field ( module_name, 'coszen', axes, time, & 'cosine of zenith angle', missing_value=-2.0 ) id_albedo_dir = register_tiled_diag_field ( module_name, 'albedo_dir', & (/id_lon,id_lat,id_band/), time, & 'land surface albedo for direct light', missing_value=-1.0 ) id_albedo_dif = register_tiled_diag_field ( module_name, 'albedo_dif', & (/id_lon,id_lat,id_band/), time, & 'land surface albedo for diffuse light', missing_value=-1.0 ) id_vegn_refl_dir = register_tiled_diag_field(module_name, 'vegn_refl_dir', & (/id_lon, id_lat, id_band/), time, & 'black-background canopy reflectivity for direct light',missing_value=-1.0) id_vegn_refl_dif = register_tiled_diag_field(module_name, 'vegn_refl_dif', & (/id_lon, id_lat, id_band/), time, & 'black-background canopy reflectivity for diffuse light',missing_value=-1.0) id_vegn_refl_lw = register_tiled_diag_field ( module_name, 'vegn_refl_lw', axes, time, & 'canopy reflectivity for thermal radiation', missing_value=-1.0) id_vegn_tran_dir = register_tiled_diag_field(module_name, 'vegn_tran_dir', & (/id_lon, id_lat, id_band/), time, & 'part of direct light that passes through canopy unscattered',missing_value=-1.0) id_vegn_tran_dif = register_tiled_diag_field(module_name, 'vegn_tran_dif', & (/id_lon, id_lat, id_band/), time, & 'black-background canopy transmittance for diffuse light',missing_value=-1.0) id_vegn_tran_lw = register_tiled_diag_field ( module_name, 'vegn_tran_lw', axes, time, & 'canopy transmittance for thermal radiation', missing_value=-1.0) id_vegn_sctr_dir = register_tiled_diag_field(module_name, 'vegn_sctr_dir', & (/id_lon, id_lat, id_band/), time, & 'part of direct light scattered downward by canopy',missing_value=-1.0) id_subs_refl_dir = register_tiled_diag_field(module_name, 'subs_refl_dir', & (/id_lon, id_lat, id_band/), time, & 'substrate reflectivity for direct light',missing_value=-1.0) id_subs_refl_dif = register_tiled_diag_field(module_name, 'subs_refl_dif', & (/id_lon, id_lat, id_band/), time, & 'substrate reflectivity for diffuse light',missing_value=-1.0) id_subs_emis = register_tiled_diag_field(module_name, 'subs_emis', & (/id_lon, id_lat/), time, & 'substrate emissivity for long-wave radiation',missing_value=-1.0) id_grnd_T = register_tiled_diag_field ( module_name, 'Tgrnd', axes, time, & 'ground surface temperature', 'degK', missing_value=-1.0 ) id_total_C = register_tiled_diag_field ( module_name, 'Ctot', axes, time, & 'total land carbon', 'kg C/m2', missing_value=-1.0 ) end subroutine land_diag_init ! the code below defines the accessor routines that are used to access fields of the ! tile data structure in collective operations, like restart i/o. Fore example, a statement ! DEFINE_LAND_ACCESSOR_0D(real,lwup) ! defines a function "land_lwup_ptr" that, given a tile, returns a pointer to the field ! called "lwup" in this tile. The procedure implementing a collective operation would ! enumerate all the tiles within the domain, call accessor routine for each of them, and ! get or set the value pointed to by the accessor routine. #define DEFINE_LAND_ACCESSOR_0D(xtype,x) subroutine CONCAT3(land_,x,_ptr(t,p));\ type(land_tile_type),pointer::t;xtype,pointer::p;p=>NULL();if(associated(t))p=>t%x;end subroutine DEFINE_LAND_ACCESSOR_0D(real,frac) DEFINE_LAND_ACCESSOR_0D(real,lwup) DEFINE_LAND_ACCESSOR_0D(real,e_res_1) DEFINE_LAND_ACCESSOR_0D(real,e_res_2) ! ============================================================================ ! tile existence detector: returns a logical value indicating wether component ! model tile exists or not logical function land_tile_exists(tile) type(land_tile_type), pointer :: tile land_tile_exists = associated(tile) end function land_tile_exists #define DEFINE_TAG_ACCESSOR(x) subroutine CONCAT2(x,_tag_ptr(t,p));\ type(land_tile_type),pointer::t;integer,pointer::p;p=>NULL();if(associated(t))\ then;if (associated(t%x)) p=>t%x%tag;endif;end subroutine DEFINE_TAG_ACCESSOR(glac) DEFINE_TAG_ACCESSOR(lake) DEFINE_TAG_ACCESSOR(soil) DEFINE_TAG_ACCESSOR(vegn) end module land_model_mod