diff --git a/documentation/docs/user_guide/inputs/cable_nml.md b/documentation/docs/user_guide/inputs/cable_nml.md index 5a8616781..6e51c7757 100644 --- a/documentation/docs/user_guide/inputs/cable_nml.md +++ b/documentation/docs/user_guide/inputs/cable_nml.md @@ -25,7 +25,8 @@ The cable.nml file includes some settings that are common across all CABLE appli | filename%lai | character(len=500) | any string of max. 500 characters | uninitialised | Default leaf area index file name. | | filename%soilcolor | character(len=500) | any string of max. 500 characters | uninitialised | Soil color file name. E.g. "CABLE-AUX/offline/soilcolor_global_1x1.nc" | | filename%gw_elev | character(len=500) | any string of max. 500 characters | uninitialised | Elevation data filename. This data is used by the groundwater option only. | -| filename%fxpft | character(len=500) | any string of max. 500 characters | uninitialised | Plant functional type fraction, wood harvest and secondary harvest file. | +| filename%fxpft1 | character(len=500) | any string of max. 500 characters | uninitialised | Plant functional type fraction, wood harvest and secondary harvest file at current year. | +| filename%fxpft0 | character(len=500) | any string of max. 500 characters | uninitialised | Plant functional type fraction, wood harvest and secondary harvest file at next year. | | filename%fxluh2cable | character(len=500) | any string of max. 500 characters | uninitialised | 12 land-use states into 17 CABLE plant functional types mapping file name. | | filename%gridnew | character(len=500) | any string of max. 500 characters | uninitialised | Updated gridinfo file name. | | filename%trunk_sumbal | character(len=500) | any string of max. 500 characters | '.trunk_sumbal' | Input filename to read combined energy and water balance at each timestep (control run). Used when `consistency_check` is TRUE | diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 6ae6ffb2c..23ab2c3e3 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -326,6 +326,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) real(r_2), dimension(:,:,:), allocatable, save :: luc_atransit real(r_2), dimension(:,:), allocatable, save :: luc_fharvw real(r_2), dimension(:,:,:), allocatable, save :: luc_xluh2cable + real(r_2), dimension(:,:), allocatable, save :: luc_delarea real(r_2), dimension(:), allocatable, save :: arealand integer, dimension(:,:), allocatable, save :: landmask integer, dimension(:), allocatable, save :: cstart,cend,nap @@ -1294,6 +1295,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) allocate(luc_atransit(mland,mvmax,mvmax)) allocate(luc_fharvw(mland,mharvw)) allocate(luc_xluh2cable(mland,mvmax,mstate)) + allocate(luc_delarea(mland,mvmax)) allocate(landmask(mlon,mlat)) allocate(arealand(mland)) allocate(patchfrac_new(mlon,mlat,mvmax)) @@ -1305,14 +1307,10 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) nap(m) = landpt(m)%nap enddo - call landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable) - - call landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, & - phen,casapool,casabal,casamet,casabiome,casaflux,bgc,rad, & - cstart,cend,nap,lucmp) - - print *, 'writing new gridinfo: landuse' - print *, 'new patch information. mland= ',mland + call landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) + call landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, & + phen,casapool,casabal,casamet,casabiome,casaflux,bgc,rad, & + cstart,cend,nap,lucmp,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) do m=1,mland do np=cstart(m),cend(m) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 2faec5942..85039e49f 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -263,6 +263,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) real(r_2), dimension(:,:,:), allocatable, save :: luc_atransit real(r_2), dimension(:,:), allocatable, save :: luc_fharvw real(r_2), dimension(:,:,:), allocatable, save :: luc_xluh2cable + real(r_2), dimension(:,:), allocatable, save :: luc_delarea real(r_2), dimension(:), allocatable, save :: arealand integer, dimension(:,:), allocatable, save :: landmask integer, dimension(:), allocatable, save :: cstart,cend,nap @@ -953,6 +954,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) allocate(luc_atransit(mland,mvmax,mvmax)) allocate(luc_fharvw(mland,mharvw)) allocate(luc_xluh2cable(mland,mvmax,mstate)) + allocate(luc_delarea(mland,mvmax)) allocate(landmask(mlon,mlat)) allocate(arealand(mland)) allocate(patchfrac_new(mlon,mlat,mvmax)) @@ -964,10 +966,10 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) nap(m) = landpt(m)%nap enddo - call landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable) - call landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, & - phen,casapool,casabal,casamet,casabiome,casaflux,bgc,rad, & - cstart,cend,nap,lucmp) + call landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) + call landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, & + phen,casapool,casabal,casamet,casabiome,casaflux,bgc,rad, & + cstart,cend,nap,lucmp,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) do m=1,mland do np=cstart(m),cend(m) diff --git a/src/offline/landuse_inout.F90 b/src/offline/landuse_inout.F90 index 57adf02d3..9e088ec5a 100644 --- a/src/offline/landuse_inout.F90 +++ b/src/offline/landuse_inout.F90 @@ -1,4 +1,4 @@ - subroutine landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable) +subroutine landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) use netcdf use cable_abort_module, ONLY: nc_abort use cable_common_module, ONLY: filename @@ -10,6 +10,7 @@ subroutine landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_ real(r_2), dimension(mland,mvmax,mvmax) :: luc_atransit real(r_2), dimension(mland,mharvw) :: luc_fharvw real(r_2), dimension(mland,mvmax,mstate) :: luc_xluh2cable + real(r_2), dimension(mland,mvmax) :: luc_delarea integer, dimension(mlon,mlat) :: landmask real(r_2), dimension(mland) :: arealand ! "mland" variables @@ -64,7 +65,7 @@ subroutine landuse_data(mlon,mlat,landmask,arealand,luc_atransit,luc_fharvw,luc_ ! get the mapping matrix (landuse type to PFT) call landuse_getxluh2(mlat,mlon,landmask,filename%fxluh2cable,luc_xluh2cable) !"xluh2cable" - call landuse_getdata(mlat,mlon,landmask,filename%fxpft,luc_atransit,luc_fharvw) + call landuse_getdata(mlat,mlon,landmask,filename%fxpft0,filename%fxpft1,luc_atransit,luc_fharvw,luc_delarea) end subroutine landuse_data @@ -85,7 +86,7 @@ SUBROUTINE landuse_getxluh2(mlat,mlon,landmask,fxluh2cable,luc_xluh2cable) allocate(xluh2cable(mlon,mlat,21,mstate)) ok = nf90_open(fxluh2cable,nf90_nowrite,ncid2) - ok = nf90_inq_varid(ncid2,"xluh2cable",varxid) + ok = nf90_inq_varid(ncid2,"XLUH2CABLE_ACCESS",varxid) ok = nf90_get_var(ncid2,varxid,xluh2cable) ok = nf90_close(ncid2) ! assig the values of luc%variables @@ -120,48 +121,101 @@ SUBROUTINE landuse_getxluh2(mlat,mlon,landmask,fxluh2cable,luc_xluh2cable) END SUBROUTINE landuse_getxluh2 -SUBROUTINE landuse_getdata(mlat,mlon,landmask,fxpft,luc_atransit,luc_fharvw) +SUBROUTINE landuse_getdata(mlat,mlon,landmask,fxpft0,fxpft1,luc_atransit,luc_fharvw,luc_delarea) ! get LUC data USE netcdf USE cable_def_types_mod, ONLY: mland,r_2 + USE casa_ncdf_module, ONLY: handle_err use landuse_constant, ONLY: mstate,mvmax,mharvw IMPLICIT NONE - character*500 fxpft + character*500 fxpft0,fxpft1 + character*200 msg integer mlat,mlon integer, dimension(mlon,mlat) :: landmask real(r_2), dimension(mland,mvmax,mvmax) :: luc_atransit real(r_2), dimension(mland,mharvw) :: luc_fharvw + real(r_2), dimension(mland,mvmax) :: luc_delarea ! local variables - real(r_2), dimension(:,:,:), allocatable :: fracharvw - real(r_2), dimension(:,:,:,:), allocatable :: transitx - integer ok,ncid1,varxid - integer i,j,m,k,ivt + real, dimension(:,:,:), allocatable :: fracharvw + real, dimension(:,:,:,:), allocatable :: transitx + real, dimension(:,:,:), allocatable :: fracpatch0,fracpatch1 + real, dimension(:,:), allocatable :: fracp0,fracp1 + integer,dimension(:,:,:), allocatable :: cablepft0,cablepft1 + integer ok,ncid0,ncid1,varxid + integer i,j,m,k,ivt,pft0,pft1 allocate(fracharvw(mlon,mlat,mharvw)) allocate(transitx(mlon,mlat,mvmax,mvmax)) + allocate(fracpatch0(mlon,mlat,mvmax),fracpatch1(mlon,mlat,mvmax)) + allocate(cablepft0(mlon,mlat,mvmax),cablepft1(mlon,mlat,mvmax)) + allocate(fracp0(mland,mvmax),fracp1(mland,mvmax)) + + ok = nf90_open(fxpft0,nf90_nowrite,ncid0) + ok = nf90_inq_varid(ncid0,"CABLEpft",varxid) + ok = nf90_get_var(ncid0,varxid,cablepft0) + msg = 'cablepft0 was not read correctly from file= '//trim(fxpft0) + call handle_err(ok,msg) + + ok = nf90_inq_varid(ncid0,"patchfrac",varxid) + ok = nf90_get_var(ncid0,varxid,fracpatch0) + msg = 'patchfrac0 was not read correctly from file= '//trim(fxpft0) + call handle_err(ok,msg) + + ok = nf90_inq_varid(ncid0,"harvest",varxid) + ok = nf90_get_var(ncid0,varxid,fracharvw) + msg = 'harvest was not read correctly from file= '//trim(fxpft0) + call handle_err(ok,msg) + + ok = nf90_inq_varid(ncid0,"transition",varxid) + ok = nf90_get_var(ncid0,varxid,transitx) + msg = 'transition matrix was not read correctly from file= '//trim(fxpft0) + call handle_err(ok,msg) + + ok = nf90_close(ncid0) + + ok = nf90_open(fxpft1,nf90_nowrite,ncid1) + ok = nf90_inq_varid(ncid1,"CABLEpft",varxid) + ok = nf90_get_var(ncid1,varxid,cablepft1) + msg = 'cablepft was not read correctly from file= '//trim(fxpft1) + call handle_err(ok,msg) + + ok = nf90_inq_varid(ncid1,"patchfrac",varxid) + ok = nf90_get_var(ncid1,varxid,fracpatch1) + msg = 'patchfrac was not read correctly from file= '//trim(fxpft1) + call handle_err(ok,msg) - ok = nf90_open(fxpft,nf90_nowrite,ncid1) - ok = nf90_inq_varid(ncid1,"harvest",varxid) - ok = nf90_get_var(ncid1,varxid,fracharvw) - ok = nf90_inq_varid(ncid1,"transition",varxid) - ok = nf90_get_var(ncid1,varxid,transitx) ok = nf90_close(ncid1) - ! assig the values of luc%variables - luc_fharvw(:,:) =0.0; luc_atransit(:,:,:)=0.0 + ! assig the values of luc_variables + luc_fharvw(:,:) =0.0; luc_atransit(:,:,:)=0.0;luc_delarea(:,:) = 0.0; fracp0(:,:)=0.0; fracp1(:,:)=0.0 + ! remove some missing values + transitx = min(1.0,max(0.0,transitx)) m = 0 - do i=1,mlon do j=1,mlat + do i=1,mlon if(landmask(i,j) ==1) then m= m +1 luc_atransit(m,:,:) = transitx(i,j,:,:) luc_fharvw(m,:) = fracharvw(i,j,:) + do k=1,mvmax + pft0=cablepft0(i,j,k) + pft1=cablepft1(i,j,k) + fracp0(m,pft0) = min(1.0,max(0.0,fracpatch0(i,j,k))) + fracp1(m,pft1) = min(1.0,max(0.0,fracpatch1(i,j,k))) + enddo + do k=1,mvmax + luc_delarea(m,k) = min(1.0,max(-1.0,fracp1(m,k)-fracp0(m,k))) + enddo endif enddo enddo deallocate(fracharvw) deallocate(transitx) + deallocate(fracpatch0,fracpatch1) + deallocate(fracp0,fracp1) + deallocate(cablepft0,cablepft1) + END SUBROUTINE landuse_getdata subroutine create_new_gridinfo(fgridold,fgridnew,mlon,mlat,landmask,patchfrac_new) @@ -913,8 +967,8 @@ SUBROUTINE WRITE_LANDUSE_CASA_RESTART_NC(mpx, lucmp, CASAONLY ) ! CHARACTER(len=20),DIMENSION(3), PARAMETER :: A4 = (/ 'csoil', 'nsoil', 'psoil' /) ! 1 dim arrays (npt ) - CHARACTER(len=20),DIMENSION(12) :: A1 - CHARACTER(len=20),DIMENSION(2) :: AI1 + CHARACTER(len=20),DIMENSION(13) :: A1 + CHARACTER(len=20),DIMENSION(3) :: AI1 ! 2 dim arrays (npt,mplant) CHARACTER(len=20),DIMENSION(3) :: A2 ! 2 dim arrays (npt,mlitter) @@ -942,9 +996,11 @@ SUBROUTINE WRITE_LANDUSE_CASA_RESTART_NC(mpx, lucmp, CASAONLY ) A1(10) = 'phen' A1(11) = 'aphen' A1(12) = 'nsoilmin' + A1(13) = 'patchfrac' AI1(1) = 'phase' AI1(2) = 'doyphase3' + AI1(3) = 'iveg' A2(1) = 'cplant' @@ -1074,12 +1130,18 @@ SUBROUTINE WRITE_LANDUSE_CASA_RESTART_NC(mpx, lucmp, CASAONLY ) STATUS = NF90_PUT_VAR(FILE_ID, VID1(12), lucmp%Nsoilmin ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) + STATUS = NF90_PUT_VAR(FILE_ID, VID1(13), lucmp%patchfrac ) + IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) + STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(1), lucmp%phase ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(2), lucmp%doyphase3 ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) + STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(3), lucmp%iveg ) + IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) + STATUS = NF90_PUT_VAR(FILE_ID, VID2(1), lucmp%cplant ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) diff --git a/src/science/landuse/landuse3.F90 b/src/science/landuse/landuse3.F90 index faccd6a58..b3bda5ae6 100644 --- a/src/science/landuse/landuse3.F90 +++ b/src/science/landuse/landuse3.F90 @@ -1,3 +1,16 @@ +! landuse change module developed by YP Wang, CSIRO Oceans and Atmosphere +! some notation conventions +! (1) arrange index sequence: +! for any two-dimensional variable: x(mlon,mlat) +! +! (2) translation from 2d field inot 1d land point +! nland=0 +! do j=1,mlat +! do i=1,mlon +! if(landmask(i,j)==1) nland=nland+1 +! enddo +! enddo +! MODULE landuse_variable use landuse_constant IMPLICIT NONE @@ -195,6 +208,7 @@ MODULE landuse_variable REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: fharvw REAL(r_2), DIMENSION(:,:,:), ALLOCATABLE :: xluh2cable REAL(r_2), DIMENSION(:,:,:), ALLOCATABLE :: atransit + REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: delarea END TYPE landuse_mland TYPE landuse_mp @@ -380,7 +394,8 @@ SUBROUTINE landuse_allocate_mland(imland,luc) ALLOCATE(luc%pftfrac(imland,mvtype), & luc%fharvw(imland,mharvw), & luc%xluh2cable(imland,mvmax,mstate), & - luc%atransit(imland,mvmax,mvmax)) + luc%atransit(mland,mvmax,mvmax), & + luc%delarea(mland,mvmax)) ! luc%phen_y(imland,mvmax), & ! luc%aphen_y(imland,mvmax), & @@ -483,7 +498,8 @@ SUBROUTINE landuse_allocate_mland(imland,luc) luc%psoilsorb_y = 0.0; luc%psoilocc_y = 0.0 luc%cwoodprod_y=0.0; luc%nwoodprod_y=0.0; luc%pwoodprod_y=0.0 - luc%pftfrac = 0.0; luc%fharvw=0.0; luc%xluh2cable=0.0; luc%atransit=0.0 + luc%pftfrac = 0.0; luc%fharvw=0.0; luc%xluh2cable=0.0 + luc%atransit=0.0; luc%delarea=0.0 END SUBROUTINE landuse_allocate_mland SUBROUTINE landuse_deallocate_mland(luc) @@ -613,7 +629,8 @@ SUBROUTINE landuse_deallocate_mland(luc) DEALLOCATE(luc%pftfrac, & luc%fharvw, & luc%xluh2cable, & - luc%atransit) + luc%atransit, & + luc%delarea) END SUBROUTINE landuse_deallocate_mland @@ -733,7 +750,7 @@ END MODULE landuse_variable subroutine landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, & phen,casapool,casabal,casamet,casabiome,casaflux,bgc,rad, & - cstart,cend,nap,lucmp) + cstart,cend,nap,lucmp,luc_atransit,luc_fharvw,luc_xluh2cable,luc_delarea) !! Main driver for the land-use change ! USE cable_IO_vars_module, ONLY: mask,patch,landpt, latitude, longitude @@ -768,9 +785,13 @@ subroutine landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, ! output ! "mland" variables integer, dimension(mland) :: cstart,cend,nap + real(r_2) luc_atransit(mland,mvmax,mvmax) + real(r_2) luc_fharvw(mland,mharvw) + real(r_2) luc_xluh2cable(mland,mvmax,mstate) + real(r_2) luc_delarea(mland,mvmax) - character*500 fxpft,fxluh2cable - integer ivt,ee,hh,np,p,q,np1 + character*500 fxpft1,fxluh2cable + integer ivt,ee,hh,np,p,q,np1,k1,k2 integer ncid,ok,xID,yID,varID,i,j,m,mpx print *, 'calling allocate mp: landuse' @@ -857,6 +878,23 @@ subroutine landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, luc%pftfrac(m,ivt) = patch(np)%frac endif enddo + do k1=1,mvmax + do k2=1,mvmax + luc%atransit(m,k1,k2) = luc_atransit(m,k1,k2) + enddo + enddo + + do k1=1,mharvw + luc%fharvw(m,k1) = luc_fharvw(m,k1) + enddo + + do k1=1,mvmax + luc%delarea(m,k1) = luc_delarea(m,k1) + do k2=1,mstate + luc%xluh2cable(m,k1,k2) = luc_xluh2cable(m,k1,k2) + enddo + enddo + enddo print *, 'point F: landuse' @@ -956,17 +994,11 @@ subroutine landuse_driver(mlon,mlat,landmask,arealand,ssnow,soil,veg,bal,canopy, enddo enddo -! print *, 'calling land2mpx: landuse' - call landuse_land2mpx(luc,lucmp,mpx) - ! call landuse_land2mpx(luc,lucmp,mpx,cstart,cend,nap) + ! sort all variables (mland,mvmax,:) by descending patchfrac orde + call landuse_land2mpx(luc,lucmp,mpx,cstart,cend,nap) -! print *, 'calling deallocate mland: landuse' call landuse_deallocate_mland(luc) -! print *, 'landuse: exit landuse_driver mpx', mpx - - close(21) -211 format(i4,a120) end subroutine landuse_driver SUBROUTINE landuse_mp2land(luc,lucmp,mp,cstart,cend) @@ -1252,7 +1284,7 @@ SUBROUTINE landuse_transitx(luc,casabiome) if(d<11.or.d>13) then delfwhpri(d) = luc%fharvw(p,1) * luc%xluh2cable(p,d,1) ! donor (positive) else - delfwhpri(d) = -luc%fharvw(p,1) * luc%xluh2cable(p,r,3) ! receiver (negative) + delfwhpri(d) = -luc%fharvw(p,1) * luc%xluh2cable(p,d,3) ! receiver (negative) endif enddo ! of "d" call landuse_redistribution(p,mvmax,delfwhpri,afwhpri) @@ -1795,7 +1827,10 @@ SUBROUTINE landuse_update_mland(luc) END SUBROUTINE landuse_update_mland - SUBROUTINE landuse_land2mpx(luc,lucmp,mpx) + SUBROUTINE landuse_land2mpx(luc,lucmp,mpx,cstart,cend,nap) + ! assign luc%varx(mland,mvmax,:) to lucmp%varx(mp,:) + ! within each land cell, order the variable by descending patchfracions + ! !! Maps `luc%var_y(mland,mvmax)` to `lucmp%var(mp)` USE landuse_constant, ONLY: mvmax USE landuse_variable @@ -1803,19 +1838,35 @@ SUBROUTINE landuse_land2mpx(luc,lucmp,mpx) IMPLICIT NONE TYPE(landuse_mland) :: luc TYPE(landuse_mp) :: lucmp -! integer, dimension(mland) :: cstart,cend,nap + integer, dimension(mland) :: cstart,cend,nap integer mpx - integer np,np1,p,q,n,npnew,npold + integer np,np1,p,q,n,npnew,napx + integer, dimension(mvmax) :: tmpint + real(r_2),dimension(mvmax) :: tmpx - npnew=0; npold=0 do p=1,mland do q=1,mvmax - if(luc%patchfrac_x(p,q)>thresh_frac) then - npold=npold +1 + tmpx(q) = luc%patchfrac_y(p,q) + tmpint(q)= q + enddo + ! look through "mland" and order the patches within the cell in descending "patchfrac" + call sortx(mvmax,tmpx,tmpint) + napx = 0 + do q=1,mvmax + if(tmpx(q)> thresh_frac) then + napx = napx +1 endif - if(luc%patchfrac_y(p,q)>thresh_frac) then - npnew = npnew +1 - lucmp%iveg(npnew) = q + enddo + + if(napx/=nap(p)) then + print *, 'number of patches not equal at land cell! Stop ', p,napx,nap(p) + stop + endif + + do npnew=cstart(p),cend(p) + lucmp%iveg(npnew) = tmpint(npnew-cstart(p)+1) + lucmp%patchfrac(npnew) = tmpx(npnew-cstart(p)+1) + q = lucmp%iveg(npnew) lucmp%isoil(npnew) = luc%isoil_y(p,q) lucmp%soilorder(npnew) = luc%soilorder_y(p,q) lucmp%phase(npnew) = luc%phase_y(p,q) @@ -1890,15 +1941,38 @@ SUBROUTINE landuse_land2mpx(luc,lucmp,mpx) lucmp%psoilocc(npnew) = luc%psoilocc_y(p,q) lucmp%pwoodprod(npnew,:) = luc%pwoodprod_y(p,q,:) endif - endif - !update patch_type - enddo ! end of "q" - enddo ! end of "p" - print *, 'npnew npold', npnew,npold + enddo ! end of "npnew" + enddo ! end of "p" END SUBROUTINE landuse_land2mpx + subroutine sortx(mvmax,tmpx,tmpint) + ! based on numerical recipes, straight insertion method p322 + USE cable_def_types_mod, ONLY: r_2 + implicit none + integer, intent(in) :: mvmax + integer, intent(inout) :: tmpint(mvmax) + real(r_2), intent(inout) :: tmpx(mvmax) + integer :: i, j, na + real(r_2) :: xa + + do j=2,mvmax + xa = tmpx(j) + na = tmpint(j) + do i=j-1,1,-1 + if(tmpx(i)>=xa) exit + tmpx(i+1) = tmpx(i) + tmpint(i+1) = tmpint(i) + enddo + if(tmpx(i)