#include module land_transitions_mod #include "../shared/debug.inc" use constants_mod, only : PI #ifdef INTERNAL_FILE_NML use mpp_mod, only: input_nml_file #else use fms_mod, only: open_namelist_file #endif use fms_mod, only : write_version_number, string, error_mesg, FATAL, WARNING, NOTE, & mpp_pe, write_version_number, file_exist, close_file, & check_nml_error, stdlog, mpp_root_pe, fms_error_handler use mpp_io_mod, only : mpp_open, mpp_close, MPP_RDONLY, MPP_ASCII use time_manager_mod, only : time_type, set_date, get_date, set_time, & operator(+), operator(-), operator(>), operator(<), operator(<=), operator(/), & operator(//), days_in_year, print_date, increment_date, get_time, & valid_calendar_types, get_calendar_type use get_cal_time_mod, only : get_cal_time use horiz_interp_mod, only : horiz_interp_type, horiz_interp_init, & horiz_interp_new, horiz_interp_del, horiz_interp use time_interp_mod, only : time_interp use diag_manager_mod, only : register_diag_field, send_data use nfu_mod, only : nfu_validtype, nfu_inq_var, nfu_get_dim_bounds, nfu_get_rec, & nfu_get_dim, nfu_get_valid_range, nfu_is_valid use vegn_data_mod, only : & N_LU_TYPES, LU_NTRL, LU_SCND, landuse_name, landuse_longname use cana_tile_mod, only : cana_tile_heat use snow_tile_mod, only : snow_tile_heat use vegn_tile_mod, only : vegn_tile_heat use soil_tile_mod, only : soil_tile_heat use land_tile_mod, only : & land_tile_type, land_tile_list_type, land_tile_enum_type, new_land_tile, delete_land_tile, & first_elmt, tail_elmt, next_elmt, operator(/=), operator(==), current_tile, & land_tile_list_init, land_tile_list_end, & empty, erase, remove, insert, land_tiles_can_be_merged, merge_land_tiles, & get_tile_water, land_tile_carbon, land_tile_heat use land_tile_io_mod, only : print_netcdf_error use land_data_mod, only : & land_data_type, lnd use vegn_tile_mod, only : & vegn_tile_type, vegn_tran_priority use vegn_harvesting_mod, only : & vegn_cut_forest use land_debug_mod, only : set_current_point, is_watch_point, get_current_point implicit none private ! ==== public interface ===================================================== public :: land_transitions_init public :: land_transitions_end public :: save_land_transitions_restart public :: land_transitions ! ==== end of public interface ============================================== ! ==== module constants ===================================================== character(len=*), parameter :: & version = '$Id: transitions.F90,v 20.0 2013/12/13 23:30:55 fms Exp $', & tagname = '$Name: tikal $', & module_name = 'land_transitions_mod', & diag_mod_name = 'landuse' ! selectors for overshoot handling options, for efficiency integer, parameter :: & OPT_IGNORE = 0, & OPT_STOP = 1, & OPT_REPORT = 2 ! ==== NetCDF declarations =================================================== include 'netcdf.inc' #define __NF_ASRT__(x) call print_netcdf_error((x),module_name,__LINE__) ! ==== data types =========================================================== type :: tran_type integer :: donor = 0 ! kind of donor tile integer :: acceptor = 0 ! kind of acceptor tile real :: frac = 0 ! area of transition end type tran_type ! ==== module data ========================================================== logical :: module_is_initialized = .FALSE. integer :: ncid ! netcd id of the input file integer :: input_ids (N_LU_TYPES,N_LU_TYPES) ! id's of input transition rate fields integer :: diag_ids (N_LU_TYPES,N_LU_TYPES) real, allocatable :: buffer_in(:,:) ! buffer for input data type(horiz_interp_type), save :: interp type(time_type), allocatable :: time_in(:) ! time axis in input data type(time_type) :: time0 ! time of previous transition calculations integer :: overshoot_opt ! selector for overshoot handling options, for efficiency integer :: conservation_opt ! selector for non-conservation handling options, for efficiency ! ---- namelist variables --------------------------------------------------- logical :: do_landuse_change = .FALSE. ! if true, then the landuse changes with time character(len=512) :: input_file = '' ! sets how to handle transition overshoot: that is, the situation when transition ! is larger than available area of the given land use type. character(len=16) :: overshoot_handling = 'report' ! or 'stop', or 'ignore' real :: overshoot_tolerance = 1e-4 ! tolerance interval for overshoots ! specifies how to handle non-conservation character(len=16) :: conservation_handling = 'stop' ! or 'report', or 'ignore' namelist/landuse_nml/input_file, do_landuse_change, & overshoot_handling, overshoot_tolerance, & conservation_handling contains ! ################################################################### ! ============================================================================ subroutine land_transitions_init(id_lon, id_lat) integer, intent(in) :: id_lon, id_lat ! the IDs of land diagnostic axes ! ---- local vars logical :: grid_initialized = .false. integer :: len, unit, ierr, io integer :: year,month,day,hour,min,sec integer :: k1,k2,i real,allocatable :: lon_in(:,:),lat_in(:,:) character(len=12) :: fieldname integer :: dimids(NF_MAX_VAR_DIMS), dimlens(NF_MAX_VAR_DIMS) integer :: nrec ! number of records in the input file real, allocatable :: time(:) ! real values of time coordinate real, allocatable :: mask_in(:,:) ! valid data mask on the input data grid type(nfu_validtype) :: v ! valid values range integer :: timedim ! id of the record (time) dimension integer :: timevar ! id of the time variable character(len=NF_MAX_NAME) :: timename ! name of the time variable character(len=256) :: timeunits ! units ot time in the file character(len=24) :: calendar ! model calendar if(module_is_initialized) return call horiz_interp_init call write_version_number(version, tagname) #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=landuse_nml, iostat=io) ierr = check_nml_error(io, 'landuse_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr = 1; do while (ierr /= 0) read (unit, nml=landuse_nml, iostat=io, end=10) ierr = check_nml_error (io, 'landuse_nml') enddo 10 continue call close_file (unit) endif #endif if (mpp_pe() == mpp_root_pe()) then unit=stdlog() write(unit, nml=landuse_nml) endif ! read restart file, if any if (file_exist('INPUT/landuse.res')) then call error_mesg('land_transitions_init',& 'reading restart "INPUT/landuse.res"',& NOTE) call mpp_open(unit,'INPUT/landuse.res', action=MPP_RDONLY, form=MPP_ASCII) read(unit,*) year,month,day,hour,min,sec time0 = set_date(year,month,day,hour,min,sec) call mpp_close(unit) else call error_mesg('land_transitions_init',& 'cold-starting land transitions',& NOTE) time0 = set_date(0001,01,01); endif ! parse the overshoot handling option if (trim(overshoot_handling)=='stop') then overshoot_opt = OPT_STOP else if (trim(overshoot_handling)=='ignore') then overshoot_opt = OPT_IGNORE else if (trim(overshoot_handling)=='report') then overshoot_opt = OPT_REPORT else call error_mesg('land_transitions_init','overshoot_handling value "'//& trim(overshoot_handling)//'" is illegal, use "stop", "report", or "ignore"',& FATAL) endif ! parse the non-conservation handling option if (trim(conservation_handling)=='stop') then conservation_opt = OPT_STOP else if (trim(conservation_handling)=='ignore') then conservation_opt = OPT_IGNORE else if (trim(conservation_handling)=='report') then conservation_opt = OPT_REPORT else call error_mesg('land_transitions_init','conservation_handling value "'//& trim(conservation_handling)//'" is illegal, use "stop", "report", or "ignore"',& FATAL) endif if (do_landuse_change) then if (trim(input_file)=='') call error_mesg('landuse_init', & 'do_landuse_change is requested, but landuse transition file is not specified', & FATAL) ierr=nf_open(input_file,NF_NOWRITE,ncid) if(ierr/=NF_NOERR) call error_mesg('landuse_init', & 'do_landuse_change is requested, but landuse transition file "'// & trim(input_file)//'" could not be opened because '//nf_strerror(ierr), FATAL) ! initialize array of input field ids input_ids(:,:) = 0 do k1 = 1,size(input_ids,1) do k2 = 1,size(input_ids,2) ! construct a name of input field and register the field fieldname = trim(landuse_name(k1))//'2'//trim(landuse_name(k2)) if(trim(fieldname)=='2') cycle ! skip unspecified tiles ierr = nfu_inq_var(ncid,fieldname, id=input_ids(k1,k2)) if (ierr/=NF_NOERR) then if (ierr==NF_ENOTVAR) then input_ids(k1,k2)=-1 ! to indicate that input field is not present in the input data else call error_mesg('landuse_init',& 'error initializing field "'//trim(fieldname)//& '" from file "'//trim(input_file)//'" : '//& nf_strerror(ierr), & FATAL) endif endif ! initialize the input data grid and horizontal interpolator if ((.not.grid_initialized).and.(input_ids(k1,k2)>0)) then ! we assume that all transition rate fields are specified on the same grid, ! in both horizontal and time "directions". Therefore there is a single grid ! for all fields, initialized only once. __NF_ASRT__(nfu_inq_var(ncid,fieldname,dimids=dimids,dimlens=dimlens,nrec=nrec)) ! allocate temporary variables allocate(lon_in(dimlens(1)+1,1), & lat_in(1,dimlens(2)+1), & time(nrec), mask_in(dimlens(1),dimlens(2)) ) ! allocate module data allocate(buffer_in(dimlens(1),dimlens(2)),time_in(nrec)) ! get the boundaries of the horizontal axes and initialize horizontal ! interpolator __NF_ASRT__(nfu_get_dim_bounds(ncid, dimids(1), lon_in(:,1))) __NF_ASRT__(nfu_get_dim_bounds(ncid, dimids(2), lat_in(1,:))) ! get the first record from variable and obtain the mask of valid data ! assume that valid mask does not change with time __NF_ASRT__(nfu_get_rec(ncid,fieldname,1,buffer_in)) ! get the valid range for the variable __NF_ASRT__(nfu_get_valid_range(ncid,fieldname,v)) ! get the mask where (nfu_is_valid(buffer_in,v)) mask_in = 1 elsewhere mask_in = 0 end where ! add mask_in and mask_out to this call call horiz_interp_new(interp, lon_in*PI/180,lat_in*PI/180, & lnd%lonb, lnd%latb, & interp_method='conservative',& mask_in=mask_in, is_latlon_in=.TRUE. ) ! get the time axis __NF_ASRT__(nf_inq_unlimdim(ncid, timedim)) __NF_ASRT__(nfu_get_dim(ncid, timedim, time)) ! get units of time __NF_ASRT__(nf_inq_dimname(ncid, timedim, timename)) __NF_ASRT__(nf_inq_varid(ncid,timename, timevar)) timeunits = ' ' __NF_ASRT__(nf_get_att_text(ncid,timevar,'units',timeunits)) ! get model calendar calendar=valid_calendar_types(get_calendar_type()) ! loop through the time axis and get time_type values in time_in do i = 1,size(time) time_in(i) = get_cal_time(time(i),timeunits,calendar) end do grid_initialized = .true. ! get rid of allocated data deallocate(lon_in,lat_in,time,mask_in) endif enddo enddo endif ! do_landuse_changes ! initialize diagnostics diag_ids(:,:) = 0 do k1 = 1,size(diag_ids,1) do k2 = 1,size(diag_ids,2) ! skip unnamed tiles if(landuse_name(k1)=='')cycle if(landuse_name(k2)=='')cycle ! construct a name of input field and register the field fieldname = trim(landuse_name(k1))//'2'//trim(landuse_name(k2)) diag_ids(k1,k2) = register_diag_field(diag_mod_name,fieldname,(/id_lon,id_lat/), lnd%time, & 'rate of transition from '//trim(landuse_longname(k1))//' to '//trim(landuse_longname(k2)),& units='1/year', missing_value=-1.0) enddo enddo module_is_initialized=.TRUE. end subroutine ! ============================================================================ subroutine land_transitions_end() if (do_landuse_change) & call horiz_interp_del(interp) if(allocated(time_in)) & deallocate(time_in,buffer_in) module_is_initialized=.FALSE. end subroutine ! ============================================================================ subroutine save_land_transitions_restart(timestamp) character(*), intent(in) :: timestamp ! timestamp to add to the file name integer :: unit,year,month,day,hour,min,sec call mpp_open( unit, 'RESTART/'//trim(timestamp)//'landuse.res', nohdrs=.TRUE. ) if (mpp_pe() == mpp_root_pe()) then call get_date(time0, year,month,day,hour,min,sec) write(unit,'(6i6,8x,a)') year,month,day,hour,min,sec, & 'Time of previous landuse transition calculation' endif call mpp_close(unit) end subroutine save_land_transitions_restart ! ============================================================================ ! performs transitions between tiles, e.g. conversion of forests to crop, etc. subroutine land_transitions (time) type(time_type), intent(in) :: time ! ---- local vars. integer :: i,j,k1,k2 type(tran_type), pointer :: transitions(:,:,:) => NULL() integer :: second, minute, hour, day0, day1, month0, month1, year0, year1 if (.not.do_landuse_change) & return ! do nothing if landuse change not requested ! NB: in this case file/interp/data are not initialized, so it is ! not even possible to use the code below call get_date(time, year0,month0,day0,hour,minute,second) call get_date(time-lnd%dt_slow, year1,month1,day1,hour,minute,second) if(year0 == year1) & !!$ if(day0 == day1) & return ! do nothing during a year ! get transition rates for current time: read map of transitions, and accumulate ! as many layers in array of transitions as necessary. Note that "transitions" ! array gets reallocated inside get_transitions as necessary, it has only as many ! layers as the max number of transitions occuring at a point at the time. do k1 = 1,N_LU_TYPES do k2 = 1,N_LU_TYPES call get_transitions(time0,time,k1,k2,transitions) enddo enddo ! perform the transitions do j = lnd%js,lnd%je do i = lnd%is,lnd%ie if(empty(lnd%tile_map(i,j))) cycle ! skip cells where there is no land ! set current point for debugging call set_current_point(i,j,1) ! transiton land area between different tile types call land_transitions_0d(lnd%tile_map(i,j), & transitions(i,j,:)%donor, & transitions(i,j,:)%acceptor,& transitions(i,j,:)%frac ) enddo enddo ! deallocate array of transitions if (associated(transitions)) deallocate(transitions) ! store current time for future reference time0=time end subroutine land_transitions ! ============================================================================= ! performs tile transitions in a given grid cell subroutine land_transitions_0d(d_list,d_kinds,a_kinds,area) type(land_tile_list_type), intent(inout) :: d_list ! list of tiles integer, intent(in) :: d_kinds(:) ! array of donor tile kinds integer, intent(in) :: a_kinds(:) ! array of acceptor tile kinds real , intent(in) :: area(:) ! array of areas changing from donor tiles to acceptor tiles ! ---- local vars integer :: i type(land_tile_type), pointer :: ptr type(land_tile_list_type) :: a_list type(land_tile_enum_type) :: ts, te real :: atot ! total fraction of tiles that can be involved in transitions ! variable used for conservation check: real :: lmass0, fmass0, cmass0, heat0, & soil_heat0, vegn_heat0, cana_heat0, snow_heat0 ! pre-transition values real :: lmass1, fmass1, cmass1, heat1, & soil_heat1, vegn_heat1, cana_heat1, snow_heat1 ! post-transition values real :: lm, fm ! buffers for transition calulations ! conservation check code, part 1: calculate the pre-transition grid ! cell totals lmass0 = 0 ; fmass0 = 0 ; cmass0 = 0 ; heat0 = 0 soil_heat0 = 0 ; vegn_heat0 = 0 ; cana_heat0 = 0 ; snow_heat0 = 0 ts = first_elmt(d_list) ; te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) call get_tile_water(ptr,lm,fm) lmass0 = lmass0 + lm*ptr%frac ; fmass0 = fmass0 + fm*ptr%frac heat0 = heat0 + land_tile_heat (ptr)*ptr%frac cmass0 = cmass0 + land_tile_carbon(ptr)*ptr%frac if(associated(ptr%soil)) soil_heat0 = soil_heat0 + soil_tile_heat(ptr%soil)*ptr%frac if(associated(ptr%vegn)) vegn_heat0 = vegn_heat0 + vegn_tile_heat(ptr%vegn)*ptr%frac if(associated(ptr%cana)) cana_heat0 = cana_heat0 + cana_tile_heat(ptr%cana)*ptr%frac if(associated(ptr%snow)) snow_heat0 = snow_heat0 + snow_tile_heat(ptr%snow)*ptr%frac enddo ! calculate the area that can participate in land transitions atot = 0 ; ts = first_elmt(d_list) ; te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) if (.not.associated(ptr%vegn)) cycle atot = atot + ptr%frac enddo if (is_watch_point()) then write(*,*)'### land_transitions_0d: input parameters ###' do i = 1, size(d_kinds) __DEBUG4__(i,d_kinds(i),a_kinds(i),area(i)) ! write(*,'(a,i2.2,100(2x,a,g))')'i='i,d_kinds(i)d_kinds(i),a_kinds(i),) enddo write(*,*)'### land_transitions_0d: land fractions before transitions (initial state) ###' ts = first_elmt(d_list) ; te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) if (.not.associated(ptr%vegn)) cycle write(*,*)'landuse=',ptr%vegn%landuse,' area=',ptr%frac enddo write(*,'(a,g23.16)')'total area=',atot endif ! split each donor tile and gather the parts that undergo a ! transition into a separate list. Note that the kind of the landuse is ! changed during this transition, including forest harvesting if necessary. ! This has to occur at some time before the tiles are merged, and it seems ! to be the most convenient place as both original and final landuse kind ! is known for each part. call land_tile_list_init(a_list) do i = 1,size(d_kinds) call split_changing_tile_parts(d_list,d_kinds(i),a_kinds(i),area(i)*atot,a_list) ! the factor atot normalizes the transitions to the total area in the grid cell ! available for the land use, that is, the area of land excluding lakes and glaciers enddo if (is_watch_point()) then write(*,*)'### land_transitions_0d: land fractions after splitting changing parts ###' atot = 0 ; ts = first_elmt(d_list) ; te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) if (.not.associated(ptr%vegn)) cycle write(*,'(2(a,g23.16,2x))')' donor: landuse=',ptr%vegn%landuse,' area=',ptr%frac atot = atot + ptr%frac enddo ts = first_elmt(a_list); te=tail_elmt(a_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) if (.not.associated(ptr%vegn)) cycle write(*,'(2(a,g23.16,2x))')'acceptor: landuse=',ptr%vegn%landuse,' area=',ptr%frac atot = atot + ptr%frac enddo write(*,'(a,g23.16)')'total area=',atot endif ! move all tiles from the donor list to the acceptor list -- this will ensure ! that all the tiles that can be merged at this time will be te = tail_elmt(d_list) do ts=first_elmt(d_list) if(ts==te) exit ! reached the end of the list ptr=>current_tile(ts) if(ptr%frac <= 0.0) then call erase(ts) ! if area of the tile is zero, free it else ! othervise, move it to a_list call remove(ts) call insert(ptr,a_list) endif enddo ! d_list is empty at this point ! merge all generated tiles into the source (donor) list te = tail_elmt(a_list) do ts=first_elmt(a_list) if(ts==te) exit ! break out of loop ptr=>current_tile(ts) call remove(ts) call land_tile_merge(ptr,d_list) enddo ! a_list is empty at this point call land_tile_list_end(a_list) if (is_watch_point()) then write(*,*)'### land_transitions_0d: land fractions final state ###' atot = 0 ts = first_elmt(d_list); te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) if (.not.associated(ptr%vegn)) cycle write(*,'(2(a,g23.16,2x))')'landuse=',ptr%vegn%landuse,' area=',ptr%frac atot = atot + ptr%frac enddo write(*,'(a,g23.16)')'total area=',atot endif ! conservation check part 2: calculate grid cell totals in final state, and ! compare them with pre-transition totals lmass1 = 0 ; fmass1 = 0 ; cmass1 = 0 ; heat1 = 0 soil_heat1 = 0 ; vegn_heat1 = 0 ; cana_heat1 = 0 ; snow_heat1 = 0 ts = first_elmt(d_list) ; te=tail_elmt(d_list) do while (ts /= te) ptr=>current_tile(ts); ts=next_elmt(ts) call get_tile_water(ptr,lm,fm) lmass1 = lmass1 + lm*ptr%frac ; fmass1 = fmass1 + fm*ptr%frac heat1 = heat1 + land_tile_heat (ptr)*ptr%frac cmass1 = cmass1 + land_tile_carbon(ptr)*ptr%frac if(associated(ptr%soil)) soil_heat1 = soil_heat1 + soil_tile_heat(ptr%soil)*ptr%frac if(associated(ptr%vegn)) vegn_heat1 = vegn_heat1 + vegn_tile_heat(ptr%vegn)*ptr%frac if(associated(ptr%cana)) cana_heat1 = cana_heat1 + cana_tile_heat(ptr%cana)*ptr%frac if(associated(ptr%snow)) snow_heat1 = snow_heat1 + snow_tile_heat(ptr%snow)*ptr%frac enddo call check_conservation ('liquid water', lmass0, lmass1, 1e-6) call check_conservation ('frozen water', fmass0, fmass1, 1e-6) call check_conservation ('carbon' , cmass0, cmass1, 1e-6) call check_conservation ('canopy air heat content', cana_heat0 , cana_heat1 , 1e-6) call check_conservation ('vegetation heat content', vegn_heat0 , vegn_heat1 , 1e-6) call check_conservation ('snow heat content', snow_heat0 , snow_heat1 , 1e-6) call check_conservation ('soil heat content', soil_heat0 , soil_heat1 , 1e-6) call check_conservation ('heat content', heat0 , heat1 , 1e-6) end subroutine ! ============================================================================== ! given a pointer to a tile and a tile list, insert the tile into the list so that ! if tile can be merged with any one already present, it is merged; otherwise ! the tile is inserted into the list subroutine land_tile_merge(tile, list) type(land_tile_type), pointer :: tile type(land_tile_list_type), intent(inout) :: list ! ---- local vars type(land_tile_type), pointer :: ptr type(land_tile_enum_type) :: ct,et ! try to find a tile that we can merge to ct = first_elmt(list) ; et = tail_elmt(list) do while(ct/=et) ptr=>current_tile(ct) ; ct = next_elmt(ct) if (land_tiles_can_be_merged(tile,ptr)) then call merge_land_tiles(tile,ptr) call delete_land_tile(tile) return ! break out of the subroutine endif enddo ! we reach here only if no suitable files was found in the list ! if no suitable tile was found, just insert given tile into the list call insert(tile,list) end subroutine land_tile_merge ! ============================================================================= ! splits changing parts of donor tiles into a separate tile list, performing ! land use changes in the process subroutine split_changing_tile_parts(d_list,d_kind,a_kind,dfrac,a_list) type(land_tile_list_type), intent(in) :: d_list ! list of donor tiles integer, intent(in) :: d_kind ! kind of donor tiles integer, intent(in) :: a_kind ! kind of acceptor tiles real, intent(in) :: dfrac ! fraction of land area that changes kind type(land_tile_list_type), intent(inout) :: a_list ! list of acceptors ! ---- local vars type(land_tile_enum_type) :: ct, et type(land_tile_type), pointer :: temp type(land_tile_type), pointer :: tile real :: area, darea, area0, area1 real :: x0,x1,x2 ! values of transition intensity real, parameter :: eps = 1e-6 ! area calculation precision integer :: iter real :: factor = 1.6 integer :: severity ! severity of overshoot errors integer :: i,j,k,face ! coordinates of current point, for overshoot diagnostics ! calculate total area of the tiles that should be transitioned to another kind area = 0 ct = first_elmt(d_list); et=tail_elmt(d_list) do while (ct /= et) tile=>current_tile(ct); ct=next_elmt(ct) if (.not.associated(tile%vegn)) cycle if (tile%vegn%landuse == d_kind) & area = area + tile%frac enddo ! check for overshoot situtaion: that is, a case where the transition area is ! larger than the available area if(overshoot_opt /= OPT_IGNORE.and.dfrac>area+overshoot_tolerance) then severity = WARNING if (overshoot_opt==OPT_STOP) severity = FATAL call get_current_point(i,j,k,face) call error_mesg('she_landuse',& 'transition at ('//trim(string(i))//','//trim(string(j))//& ',face='//trim(string(face))//& ') from "'//trim(landuse_name(d_kind))// & '" to "' //trim(landuse_name(a_kind))//& '" ('//trim(string(dfrac))//') is larger than area of "'& //trim(landuse_name(d_kind))//'" ('//trim(string(area))//')', & severity) endif ! if area of the tiles of requested kind is zero we cannot transition ! anything, so just return if (area==0) return ! transition cannot be more than current total area of specified kind darea = min(dfrac, area) ! solve equation to get transition intensity ! (1) bracket transition intensity interval so that requested area is within it x0=0.0; area0 = total_transition_area(d_list, d_kind, a_kind, x0) x1=1.0; area1 = total_transition_area(d_list, d_kind, a_kind, x1) iter = 0 do if ((area0<=darea).and.(area1>=darea)) exit if (area0>darea) then x0 = x0-(x1-x0)*factor area0 = total_transition_area(d_list, d_kind, a_kind, x0) else x1 = x1+(x1-x0)*factor area1 = total_transition_area(d_list, d_kind, a_kind, x1) endif iter = iter+1 if (iter>50) then call error_mesg('veg_tile_transitions',& 'cannot braket transition intensity interval after 50 iterations',& FATAL) endif enddo ! find solution for transition intensity by binary search do iter = 1,50 x2 = (x0+x1)/2 area = total_transition_area(d_list, d_kind, a_kind, x2) if (abs(x1-x2)darea) then x1=x2 else x0=x2 endif enddo ! do tile transitions to destination list ct = first_elmt(d_list); et=tail_elmt(d_list) do while (ct /= et) tile=>current_tile(ct) ; ct=next_elmt(ct) if(.not.associated(tile%vegn)) cycle if(tile%vegn%landuse == d_kind) then darea = vegn_tran_priority(tile%vegn, a_kind, x2) if(darea > 0) then ! make a copy temp => new_land_tile(tile) temp%frac = tile%frac*darea tile%frac = tile%frac*(1.0-darea) ! convert land use type of the tile: ! cut the forest, if necessary if(temp%vegn%landuse==LU_NTRL.or.temp%vegn%landuse==LU_SCND) & call vegn_cut_forest(temp%vegn, a_kind) ! change landuse type of the tile temp%vegn%landuse = a_kind ! add the new tile to the resulting list call insert(temp, a_list) ! insert tile into output list endif endif enddo end subroutine split_changing_tile_parts ! ============================================================================ ! calculates total area (fraction of grid cell area) participating in ! vegetation transition from src_kind to dst_kind for given transition ! intensity tau function total_transition_area(list,src_kind,dst_kind,tau) result (total_area) real :: total_area type(land_tile_list_type), intent(in) :: list ! list of tiles integer , intent(in) :: src_kind, dst_kind ! source and destination kinds real , intent(in) :: tau ! transition intensity ! ---- local vars type(land_tile_enum_type) :: ct, et type(land_tile_type), pointer :: tile total_area = 0 ct = first_elmt(list) ; et = tail_elmt(list) do while (ct/=et) tile=>current_tile(ct); ct = next_elmt(ct) if (.not.associated(tile%vegn)) cycle ! skip non-vegetated tiles if(tile%vegn%landuse == src_kind) & total_area = total_area + tile%frac*vegn_tran_priority(tile%vegn,dst_kind,tau) enddo end function ! ============================================================================ subroutine get_transitions(time0,time1,k1,k2,tran) type(time_type), intent(in) :: time0 ! time of previous calculation of ! transitions (the integral transitinos will be calculated between time0 ! and time) type(time_type), intent(in) :: time1 ! current time integer, intent(in) :: k1,k2 ! kinds of tiles type(tran_type), pointer :: tran(:,:,:) ! transition info ! ---- local vars integer :: i,j,k,sec,days type(tran_type), pointer :: ptr(:,:,:) => NULL() real :: frac(lnd%is:lnd%ie,lnd%js:lnd%je) real :: part_of_year logical :: used ! allocate array of transitions, if necessary if (.not.associated(tran)) then allocate(tran(lnd%is:lnd%ie,lnd%js:lnd%je,1) ) end if ! get transition rate for this specific transition frac(:,:) = 0.0 if(input_ids(k1,k2)>0) then call integral_transition(time0,time1,input_ids(k1,k2),frac) do j = lnd%js,lnd%je do i = lnd%is,lnd%ie if(frac(i,j) == 0) cycle ! skip points where transition rate is zero ! find the first empty transition element for the current indices k = 1 do while ( k <= size(tran,3) ) if(tran(i,j,k)%donor == 0) exit k = k+1 enddo if (k>size(tran,3)) then ! if there is no room, make the array of transitions larger allocate(ptr(lnd%is:lnd%ie,lnd%js:lnd%je,size(tran,3)*2)) ptr(:,:,1:size(tran,3)) = tran deallocate(tran) tran => ptr nullify(ptr) end if ! store the transition element tran(i,j,k) = tran_type(k1,k2,frac(i,j)) enddo enddo endif ! send transition data to diagnostics if(diag_ids(k1,k2)>0) then call get_time(time1-time0, sec,days) part_of_year = (days+sec/86400.0)/days_in_year(time0) used = send_data(diag_ids(k1,k2),frac/part_of_year,time1) endif end subroutine get_transitions ! =========================================================================== ! given beginnig and end of the time period, and the id of the external ! tile transition field, calculates total transition during the specified period. ! The transition rate data are assumed to be in fraction of land area per year, ! timestamped at the beginning of the year subroutine integral_transition(t1, t2, id, frac, err_msg) type(time_type), intent(in) :: t1,t2 ! time boundaries integer , intent(in) :: id ! id of the field real , intent(out) :: frac(:,:) character(len=*),intent(out), optional :: err_msg ! ---- local vars integer :: n ! size of time axis type(time_type) :: ts,te integer :: i1,i2 real :: w ! time interpolation weight real :: dt ! current time interval, in years real :: sum(size(frac,1),size(frac,2)) integer :: i,j character(len=256) :: msg msg = '' ! adjust the integration limits, in case they are out of range n = size(time_in) ts = t1; if (tstime_in(n)) ts = time_in(n) te = t2 if (tetime_in(n)) te = time_in(n) call time_interp(ts, time_in, w, i1,i2, err_msg=msg) if(msg /= '') then if(fms_error_handler('integral_transition','Message from time_interp: '//trim(msg),err_msg)) return endif __NF_ASRT__(nfu_get_rec(ncid,id,i1,buffer_in)) frac = 0; call horiz_interp(interp,buffer_in,frac) dt = (time_in(i2)-time_in(i1))//set_time(0,days_in_year((time_in(i2)+time_in(i1))/2)) sum = -frac*w*dt do while(time_in(i2)<=te) __NF_ASRT__(nfu_get_rec(ncid,id,i1,buffer_in)) call horiz_interp(interp,buffer_in,frac) dt = (time_in(i2)-time_in(i1))//set_time(0,days_in_year((time_in(i2)+time_in(i1))/2)) sum = sum+frac*dt i2 = i2+1 i1 = i2-1 if(i2>size(time_in)) exit ! from loop enddo call time_interp(te,time_in,w,i1,i2, err_msg=msg) if(msg /= '') then if(fms_error_handler('integral_transition','Message from time_interp: '//trim(msg),err_msg)) return endif __NF_ASRT__(nfu_get_rec(ncid,id,i1,buffer_in)) call horiz_interp(interp,buffer_in,frac) dt = (time_in(i2)-time_in(i1))//set_time(0,days_in_year((time_in(i2)+time_in(i1))/2)) frac = sum+frac*w*dt do i = 1,size(frac,1) do j = 1,size(frac,2) if(frac(i,j)<0) then call error_mesg('get','transition rate is below zero',FATAL) endif enddo enddo end subroutine ! ============================================================================== ! checks conservation and aborts with fatal error if tolerance is exceeded subroutine check_conservation(name, d1, d2, tolerance) character(*), intent(in) :: name ! name of the component real, intent(in) :: d1,d2 ! values to check real, intent(in) :: tolerance ! tolerance of the test integer :: curr_i, curr_j, face integer :: severity ! severity of the generated message character(256) :: message if (conservation_opt == OPT_IGNORE) return ! do nothing severity = WARNING if (overshoot_opt==OPT_STOP) severity = FATAL if (abs(d1-d2)>tolerance) then call get_current_point(i=curr_i,j=curr_j,face=face) write(message,'(a,3(x,a,i4), 2(x,a,g23.16))')& 'conservation of '//trim(name)//' is violated', & 'at i=',curr_i,'j=',curr_j,'face=',face, & 'value before=', d1, 'after=', d2 call error_mesg('land_transitions',message,severity) endif end subroutine end module