diff --git a/.circleci/config.yml b/.circleci/config.yml index 11bd27865..f3eeea071 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,5 +1,9 @@ version: 2.1 +# Anchors to prevent forgetting to update a version +baselibs_version: &baselibs_version v7.5.0 +bcs_version: &bcs_version v10.22.3 + orbs: ci: geos-esm/circleci-tools@1 @@ -14,6 +18,7 @@ workflows: matrix: parameters: compiler: [gfortran, ifort] + baselibs_version: *baselibs_version repo: GEOSgcm checkout_fixture: true mepodevelop: true diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 2d4e599e5..7a4871936 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -17,7 +17,7 @@ jobs: if: "!contains(github.event.pull_request.labels.*.name, '0 diff trivial')" runs-on: ubuntu-latest container: - image: gmao/ubuntu20-geos-env-mkl:v6.2.13-openmpi_4.1.2-gcc_11.2.0 + image: gmao/ubuntu20-geos-env-mkl:v7.5.0-openmpi_4.1.2-gcc_11.2.0 credentials: username: ${{ secrets.DOCKERHUB_USERNAME }} password: ${{ secrets.DOCKERHUB_TOKEN }} @@ -35,10 +35,13 @@ jobs: with: access_token: ${{ github.token }} - name: Checkout GCM - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: repository: GEOS-ESM/GEOSgcm fetch-depth: 1 + - name: Set all directories as git safe + run: | + git config --global --add safe.directory '*' - name: Versions etc. run: | gfortran --version diff --git a/GEOSagcm_GridComp/GEOS_AgcmGridComp.F90 b/GEOSagcm_GridComp/GEOS_AgcmGridComp.F90 index fb4a56992..c9cd0aa9f 100644 --- a/GEOSagcm_GridComp/GEOS_AgcmGridComp.F90 +++ b/GEOSagcm_GridComp/GEOS_AgcmGridComp.F90 @@ -1026,6 +1026,8 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) VERIFY_(STATUS) + call MAPL_TimerAdd(GC, name="AGCM_BARRIER" ,RC=STATUS) + VERIFY_(STATUS) ! All done !--------- @@ -1361,6 +1363,7 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Local derived type aliases + type (ESMF_VM) :: VMG type (MAPL_MetaComp), pointer :: STATE type (ESMF_GridComp), pointer :: GCS(:) type (ESMF_State), pointer :: GIM(:) @@ -1554,7 +1557,7 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! ----------------------------------------------------------- Iam = "Run" - call ESMF_GridCompGet ( GC, name=COMP_NAME, RC=STATUS ) + call ESMF_GridCompGet ( GC, VM=VMG, name=COMP_NAME, RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // trim(Iam) @@ -2399,6 +2402,9 @@ subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) call DO_UPDATE_PHY ('DPEDT') call DO_UPDATE_PHY ('DTDT' ) + call MAPL_TimerOn (STATE,"AGCM_BARRIER" ) + call ESMF_VMBarrier(VMG, rc=status); VERIFY_(STATUS) + call MAPL_TimerOff(STATE,"AGCM_BARRIER" ) ! Run RUN2 of SuperDynamics (ADD_INCS) to add Physics Diabatic Tendencies !------------------------------------------------------------------------ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/CMakeLists.txt index 7abedf76c..24e12a4ee 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/CMakeLists.txt @@ -10,6 +10,7 @@ install( FILES ${resource_files} ) set (srcs + getids.F90 StieglitzSnow.F90 SurfParams.F90) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/getids.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/getids.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/getids.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/getids.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt index 21d1777ea..b4b947833 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt @@ -17,7 +17,7 @@ if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) list(APPEND srcs findloc.F90) endif () -esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_LandShared esmf NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran OpenMP::OpenMP_C) +esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared esmf NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran OpenMP::OpenMP_C) if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) target_compile_definitions(${this} PRIVATE USE_EXTERNAL_FINDLOC) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CMakeLists.txt index 78a37bd8d..320bc3430 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CMakeLists.txt @@ -1,7 +1,8 @@ esma_set_this () set(srcs - getids.F90 + CatchmentRst.F90 + CatchmentCNRst.F90 ) set (exe_srcs @@ -16,14 +17,15 @@ set (exe_srcs mk_LakeLandiceSaltRestarts.F90 mk_RouteRestarts.F90 mk_GEOSldasRestarts.F90 + mk_catchANDcnRestarts.F90 ) esma_add_library (${this} SRCS ${srcs} - DEPENDENCIES MAPL esmf NetCDF::NetCDF_Fortran) + DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared esmf NetCDF::NetCDF_Fortran) foreach (src ${exe_srcs}) - string (REGEX REPLACE ".F90" "" exe ${src}) + string (REGEX REPLACE ".F90" ".x" exe ${src}) ecbuild_add_executable ( TARGET ${exe} SOURCES ${src} @@ -31,3 +33,12 @@ foreach (src ${exe_srcs}) endforeach () install(PROGRAMS mk_Restarts DESTINATION bin) +foreach (src ${exe_srcs}) + string (REGEX REPLACE ".F90" ".x" exe ${src}) + string (REGEX REPLACE ".F90" "" lname ${src}) + install(CODE "execute_process( \ + COMMAND ${CMAKE_COMMAND} -E create_symlink \ + ${CMAKE_INSTALL_PREFIX}/bin/${exe} \ + ${CMAKE_INSTALL_PREFIX}/bin/${lname} \ + )") +endforeach () diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 new file mode 100644 index 000000000..0ec4d1692 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 @@ -0,0 +1,1363 @@ +#include "MAPL_Generic.h" + +module CatchmentCNRstMod + use mk_restarts_getidsMod + use mpi + use MAPL + use CatchmentRstMod, only : CatchmentRst + implicit none + + real, parameter :: ECCENTRICITY = 0.0167 + real, parameter :: PERIHELION = 102.0 + real, parameter :: OBLIQUITY = 23.45 + integer, parameter :: EQUINOX = 80 + + integer, parameter :: nveg = 4 + integer, parameter :: nzone = 3 + integer, parameter :: VAR_COL_CLM40 = 40 ! number of CN column restart variables + integer, parameter :: VAR_PFT_CLM40 = 74 ! number of CN PFT variables per column + integer, parameter :: npft = 19 + integer, parameter :: npft_clm45 = 19 + integer, parameter :: VAR_COL_CLM45 = 35 ! number of CN column restart variables + integer, parameter :: VAR_PFT_CLM45 = 75 ! number of CN PFT variables per column + real, parameter :: nan = O'17760000000' + real, parameter :: fmin= 1.e-4 ! ignore vegetation fractions at or below this value + integer :: iclass(npft) = (/1,1,2,3,3,4,5,5,6,7,8,9,10,11,12,11,12,11,12/) + + type, extends(CatchmentRst) :: CatchmentCNRst + logical :: isCLM45 + integer :: VAR_COL + integer :: VAR_PFT + real, allocatable :: cnity(:,:) + real, allocatable :: fvg(:,:) + real, allocatable :: tg(:,:) + real, allocatable :: tgwm(:,:) + real, allocatable :: rzmm(:,:) + real, allocatable :: sfmm(:,:) + real, allocatable :: TILE_ID(:) + real, allocatable :: ndep(:) + real, allocatable :: t2(:) + real, allocatable :: BGALBVR(:) + real, allocatable :: BGALBVF(:) + real, allocatable :: BGALBNR(:) + real, allocatable :: BGALBNF(:) + real, allocatable :: CNCOL(:,:) + real, allocatable :: CNPFT(:,:) + real, allocatable :: ABM (:) + real, allocatable :: FIELDCAP(:) + real, allocatable :: HDM (:) + real, allocatable :: GDP (:) + real, allocatable :: PEATF (:) + ! below is not necessary. It is not read. It is set to 0 during writing + !real, allocatable :: bflowm(:) + !real, allocatable :: totwatm(:) + !real, allocatable :: tairm(:) + !real, allocatable :: tpm(:) + !real, allocatable :: cnsum(:) + !real, allocatable :: sndzm(:) + !real, allocatable :: asnowm(:) + !real, allocatable :: ar1m(:) + !real, allocatable :: rainfm(:) + !real, allocatable :: rhm(:) + !real, allocatable :: runsrfm(:) + !real, allocatable :: snowfm(:) + !real, allocatable :: windm(:) + !real, allocatable :: tprec10d(:) + !real, allocatable :: tprec60d(:) + !real, allocatable :: t2m10d(:) + !real, allocatable :: sfmcm(:) + !real, allocatable :: psnsunm(:,:,:) + !real, allocatable :: psnsham(:,:,:) + + contains + procedure :: write_nc4 + procedure :: allocate_cn + procedure :: add_bcs_to_cnrst + procedure :: re_tile + endtype CatchmentCNRst + + interface CatchmentCNRst + module procedure CatchmentCNRst_Create + end interface + +contains + + function CatchmentCNRst_create(filename, cnclm, time, rc) result (catch) + type(CatchmentCNRst) :: catch + character(*), intent(in) :: filename + character(*), intent(in) :: cnclm + character(*), intent(in) :: time + integer, optional, intent(out) :: rc + integer :: status + type(Netcdf4_fileformatter) :: formatter + integer :: filetype, ntiles, unit + integer :: j, dim1,dim2, myid, mpierr + type(Variable), pointer :: myVariable + character(len=:), pointer :: dname + type(FileMetadata) :: meta + character(len=256) :: Iam = "CatchmentCNRst_create" + + call MAPL_NCIOGetFileType(filename, filetype, __RC__) + if (filetype /= 0) then + _ASSERT( .false., "CatchmentCN only support nc4 file restart") + endif + + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + + catch%isCLM45 = .false. + call formatter%open(filename, pFIO_READ, __RC__) + meta = formatter%read(__RC__) + ntiles = meta%get_dimension('tile', __RC__) + catch%ntiles = ntiles + catch%meta = meta + catch%time = time + if (index(cnclm, '40') /=0) then + catch%VAR_COL = VAR_COL_CLM40 + catch%VAR_PFT = VAR_PFT_CLM40 + endif + if (index(cnclm, '45') /=0) then + catch%VAR_COL = VAR_COL_CLM45 + catch%VAR_PFT = VAR_PFT_CLM45 + catch%isCLM45 = .true. + endif + + if (myid == 0) then + call catch%allocate_cn(__RC__) + call catch%read_shared_nc4(formatter, __RC__) + + myVariable => meta%get_variable("ITY") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + do j=1,dim1 + call MAPL_VarRead(formatter,"ITY",catch%cnity(:,j),offset1=j, __RC__) + call MAPL_VarRead(formatter,"FVG",catch%fvg(:,j),offset1=j, __RC__) + enddo + + call MAPL_VarRead(formatter,"TG",catch%tg, __RC__) + call MAPL_VarRead(formatter,"TILE_ID",catch%TILE_ID, __RC__) + call MAPL_VarRead(formatter,"NDEP",catch%ndep, __RC__) + call MAPL_VarRead(formatter,"CLI_T2M",catch%t2, __RC__) + call MAPL_VarRead(formatter,"BGALBVR",catch%BGALBVR, __RC__) + call MAPL_VarRead(formatter,"BGALBVF",catch%BGALBVF, __RC__) + call MAPL_VarRead(formatter,"BGALBNR",catch%BGALBNR, __RC__) + call MAPL_VarRead(formatter,"BGALBNF",catch%BGALBNF, __RC__) + + myVariable => meta%get_variable("CNCOL") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + if( catch%isCLM45) then + call MAPL_VarRead(formatter,"ABM", catch%ABM, __RC__) + call MAPL_VarRead(formatter,"FIELDCAP",catch%FIELDCAP, __RC__) + call MAPL_VarRead(formatter,"HDM", catch%HDM , __RC__) + call MAPL_VarRead(formatter,"GDP", catch%GDP , __RC__) + call MAPL_VarRead(formatter,"PEATF", catch%PEATF , __RC__) + endif + do j=1,dim1 + call MAPL_VarRead(formatter,"CNCOL",catch%CNCOL(:,j),offset1=j, __RC__) + enddo + ! The following three lines were added as a bug fix by smahanam on 5 Oct 2020 + ! (to be merged into the "develop" branch in late 2020): + ! The length of the 2nd dim of CNPFT differs from that of CNCOL. Prior to this fix, + ! CNPFT was not read in its entirety and some elements remained uninitialized (or zero), + ! resulting in bad values in the "regridded" (re-tiled) restart file. + ! This impacted re-tiled restarts for both CNCLM40 and CLCLM45. + ! - reichle, 23 Nov 2020 + myVariable => meta%get_variable("CNPFT") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + do j=1,dim1 + call MAPL_VarRead(formatter,"CNPFT",catch%CNPFT(:,j),offset1=j, __RC__) + enddo + endif + + call formatter%close() + + if (present(rc)) rc =0 + end function CatchmentCNRst_Create + + function CatchmentCNRst_empty(meta, cnclm, time, rc) result (catch) + type(CatchmentCNRst) :: catch + type(FileMetadata), intent(in) :: meta + character(*), intent(in) :: cnclm + character(*), intent(in) :: time + integer, optional, intent(out) :: rc + integer :: status, myid, mpierr + character(len=256) :: Iam = "CatchmentCNRst_empty" + + catch%isCLM45 = .false. + catch%ntiles = meta%get_dimension('tile', __RC__) + catch%time = time + catch%meta = meta + if (index(cnclm, '40') /=0) then + catch%VAR_COL = VAR_COL_CLM40 + catch%VAR_PFT = VAR_PFT_CLM40 + endif + if (index(cnclm, '45') /=0) then + catch%VAR_COL = VAR_COL_CLM45 + catch%VAR_PFT = VAR_PFT_CLM45 + catch%isCLM45 = .true. + endif + + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + if (myid ==0) call catch%allocate_cn(__RC__) + if(present(rc)) rc = 0 + end function CatchmentCNRst_empty + + subroutine write_nc4(this, filename, rc) + class(CatchmentCNRst), intent(inout):: this + character(*), intent(in) :: filename + integer, optional, intent(out):: rc + + type(Netcdf4_fileformatter) :: formatter + integer :: status + character(256) :: Iam = "write_nc4" + integer :: i,j, dim1,dim2 + real, dimension (:), allocatable :: var + type(Variable), pointer :: myVariable + character(len=:), pointer :: dname + type(FileMetadata) :: meta + + meta = this%meta + call formatter%create(filename, __RC__) + call formatter%write(meta, __RC__) + + call this%write_shared_nc4(formatter, __RC__) + + myVariable => meta%get_variable("ITY") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + do j=1,dim1 + call MAPL_VarWrite(formatter,"ITY",this%cnity(:,j),offset1=j) + call MAPL_VarWrite(formatter,"FVG",this%fvg(:,j),offset1=j) + call MAPL_VarWrite(formatter,"TG",this%tg(:,j),offset1=j) + enddo + + call MAPL_VarWrite(formatter,"TILE_ID",this%TILE_ID) + call MAPL_VarWrite(formatter,"NDEP",this%NDEP) + call MAPL_VarWrite(formatter,"CLI_T2M",this%t2) + call MAPL_VarWrite(formatter,"BGALBVR",this%BGALBVR) + call MAPL_VarWrite(formatter,"BGALBVF",this%BGALBVF) + call MAPL_VarWrite(formatter,"BGALBNR",this%BGALBNR) + call MAPL_VarWrite(formatter,"BGALBNF",this%BGALBNF) + myVariable => meta%get_variable("CNCOL") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + + do j=1,dim1 + call MAPL_VarWrite(formatter,"CNCOL",this%CNCOL(:,j),offset1=j) + enddo + myVariable => meta%get_variable("CNPFT") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + do j=1,dim1 + call MAPL_VarWrite(formatter,"CNPFT",this%CNPFT(:,j),offset1=j) + enddo + + dim1 = meta%get_dimension('tile') + + allocate (var(dim1),source = 0.) + + call MAPL_VarWrite(formatter,"BFLOWM", var) + call MAPL_VarWrite(formatter,"TOTWATM",var) + call MAPL_VarWrite(formatter,"TAIRM", var) + call MAPL_VarWrite(formatter,"TPM", var) + call MAPL_VarWrite(formatter,"CNSUM", var) + call MAPL_VarWrite(formatter,"SNDZM", var) + call MAPL_VarWrite(formatter,"ASNOWM", var) + + myVariable => meta%get_variable("TGWM") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + do j=1,dim1 + call MAPL_VarWrite(formatter,"TGWM",var,offset1=j) + call MAPL_VarWrite(formatter,"RZMM",var,offset1=j) + end do + + if (this%isCLM45) then + do j=1,dim1 + call MAPL_VarWrite(formatter,"SFMM", var,offset1=j) + enddo + + call MAPL_VarWrite(formatter,"ABM", this%ABM, rc =rc ) + call MAPL_VarWrite(formatter,"FIELDCAP",this%FIELDCAP) + call MAPL_VarWrite(formatter,"HDM", this%HDM ) + call MAPL_VarWrite(formatter,"GDP", this%GDP ) + call MAPL_VarWrite(formatter,"PEATF", this%PEATF ) + call MAPL_VarWrite(formatter,"RHM", var) + call MAPL_VarWrite(formatter,"WINDM", var) + call MAPL_VarWrite(formatter,"RAINFM", var) + call MAPL_VarWrite(formatter,"SNOWFM", var) + call MAPL_VarWrite(formatter,"RUNSRFM", var) + call MAPL_VarWrite(formatter,"AR1M", var) + call MAPL_VarWrite(formatter,"T2M10D", var) + call MAPL_VarWrite(formatter,"TPREC10D",var) + call MAPL_VarWrite(formatter,"TPREC60D",var) + else + call MAPL_VarWrite(formatter,"SFMCM", var) + endif + myVariable => meta%get_variable("PSNSUNM") + dname => myVariable%get_ith_dimension(2) + dim1 = meta%get_dimension(dname) + dname => myVariable%get_ith_dimension(3) + dim2 = meta%get_dimension(dname) + do i=1,dim2 + do j=1,dim1 + call MAPL_VarWrite(formatter,"PSNSUNM",var,offset1=j,offset2=i) + call MAPL_VarWrite(formatter,"PSNSHAM",var,offset1=j,offset2=i) + end do + end do + call formatter%close() + + _RETURN(_SUCCESS) + end subroutine write_nc4 + + subroutine allocate_cn(this,rc) + class(CatchmentCNRst), intent(inout) :: this + integer, optional, intent(out):: rc + integer :: status + integer :: ncol,npft, ntiles + + ntiles = this%ntiles + ncol = nzone* this%VAR_COL + npft = nzone*nveg*this%VAR_PFT + + call this%CatchmentRst%allocate_catch(__RC__) + + ! W.Jiang notes : some varaiables are not allocated because they are set to zero directly during write + allocate(this%cnity(ntiles,nveg)) + allocate(this%fvg(ntiles,nveg)) + allocate(this%tg(ntiles,nveg)) + allocate(this%tgwm(ntiles,nzone)) + allocate(this%rzmm(ntiles,nzone)) + allocate(this%TILE_ID(ntiles)) + allocate(this%ndep(ntiles)) + allocate(this%t2(ntiles)) + allocate(this%BGALBVR(ntiles)) + allocate(this%BGALBVF(ntiles)) + allocate(this%BGALBNR(ntiles)) + allocate(this%BGALBNF(ntiles)) + allocate(this%CNCOL(ntiles,ncol)) + allocate(this%CNPFT(ntiles,npft)) + allocate(this%ABM(ntiles)) + allocate(this%FIELDCAP(ntiles)) + allocate(this%HDM(ntiles)) + allocate(this%GDP(ntiles)) + allocate(this%PEATF(ntiles)) + _RETURN(_SUCCESS) + end subroutine allocate_cn + + SUBROUTINE add_bcs_to_cnrst (this, SURFLAY, OutBcsDir,rc) + class(CatchmentCNRst), intent(inout) :: this + real, intent (in) :: SURFLAY + character(*), intent (in) :: OutBcsDir + integer, optional, intent(out) :: rc + real, allocatable :: CLMC_pf1(:), CLMC_pf2(:), CLMC_sf1(:), CLMC_sf2(:) + real, allocatable :: CLMC_pt1(:), CLMC_pt2(:), CLMC_st1(:), CLMC_st2(:) + real, allocatable :: NDEP(:), BVISDR(:), BVISDF(:), BNIRDR(:), BNIRDF(:) + real, allocatable :: T2(:), var1(:), hdm(:), fc(:), gdp(:), peatf(:) + integer, allocatable :: ity(:), abm (:) + integer :: STATUS, ntiles, unit27, unit28, unit29, unit30 + integer :: idum, i,j,n, ib, nv + real :: rdum, zdep1, zdep2, zdep3, zmet, term1, term2, bare,fvg(4) + logical :: NEWLAND + logical :: file_exists + + type(NetCDF4_Fileformatter) :: CatchCNFmt + character*256 :: Iam = "add_bcs" + + open (10,file =trim(OutBcsDir)//"/clsm/catchment.def",status='old',form='formatted') + read (10,*) ntiles + close (10, status = 'keep') + + !ntiles = this%ntiles + !call this%CatchmentRst%add_bcs_to_rst(surflay, OutBcsDir, __RC__) + + allocate (BVISDR(ntiles), BVISDF(ntiles), BNIRDR(ntiles) ) + allocate (BNIRDF(ntiles), T2(ntiles), NDEP(ntiles) ) + allocate (CLMC_pf1(ntiles), CLMC_pf2(ntiles), CLMC_sf1(ntiles)) + allocate (CLMC_sf2(ntiles), CLMC_pt1(ntiles), CLMC_pt2(ntiles)) + allocate (CLMC_st1(ntiles), CLMC_st2(ntiles)) + allocate (hdm(ntiles), fc(ntiles), gdp(ntiles)) + allocate (peatf(ntiles), abm(ntiles), var1(ntiles)) + + inquire(file = trim(OutBcsDir)//'/clsm/catchcn_params.nc4', exist=file_exists) + inquire(file = trim(OutBcsDir)//'/clsm/CLM_veg_typs_fracs', exist=NewLand ) + _ASSERT(Newland, "catchcn should get bc from newland") + + if(file_exists) then + call CatchCNFmt%Open(trim(OutBcsDir)//'/clsm/catchcn_params.nc4', pFIO_READ, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'BGALBNF', BNIRDF, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'BGALBNR', BNIRDR, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'BGALBVF', BVISDF, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'BGALBVR', BVISDR, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'NDEP', NDEP, __RC__) + call MAPL_VarRead ( CatchCNFmt ,'T2_M', T2, __RC__) + call MAPL_VarRead(CatchCNFmt,'ITY',CLMC_pt1,offset1=1, __RC__) ! 30 + call MAPL_VarRead(CatchCNFmt,'ITY',CLMC_pt2,offset1=2, __RC__) ! 31 + call MAPL_VarRead(CatchCNFmt,'ITY',CLMC_st1,offset1=3, __RC__) ! 32 + call MAPL_VarRead(CatchCNFmt,'ITY',CLMC_st2,offset1=4, __RC__) ! 33 + call MAPL_VarRead(CatchCNFmt,'FVG',CLMC_pf1,offset1=1, __RC__) ! 34 + call MAPL_VarRead(CatchCNFmt,'FVG',CLMC_pf2,offset1=2, __RC__) ! 35 + call MAPL_VarRead(CatchCNFmt,'FVG',CLMC_sf1,offset1=3, __RC__) ! 36 + call MAPL_VarRead(CatchCNFmt,'FVG',CLMC_sf2,offset1=4, __RC__) ! 37 + call CatchCNFmt%close() + else + + open(newunit=unit27, file=trim(OutBcsDir)//'/clsm/CLM_veg_typs_fracs' ,form='formatted') + open(newunit=unit28, file=trim(OutBcsDir)//'/clsm/CLM_NDep_SoilAlb_T2m' ,form='formatted') + + do n=1,ntiles + read (unit27, *) i,j, CLMC_pt1(n), CLMC_pt2(n), CLMC_st1(n), CLMC_st2(n), & + CLMC_pf1(n), CLMC_pf2(n), CLMC_sf1(n), CLMC_sf2(n) + + read (unit28, *) NDEP(n), BVISDR(n), BVISDF(n), BNIRDR(n), BNIRDF(n), T2(n) ! MERRA-2 Annual Mean Temp is default. + end do + + CLOSE (unit27, STATUS = 'KEEP') + CLOSE (unit28, STATUS = 'KEEP') + + endif + + if (this%isCLM45 ) then + + open(newunit=unit30, file=trim(OutBcsDir)//'/clsm/CLM4.5_abm_peatf_gdp_hdm_fc' ,form='formatted') + do n=1,ntiles + read (unit30, *) i, j, abm(n), peatf(n), & + gdp(n), hdm(n), fc(n) + end do + CLOSE (unit30, STATUS = 'KEEP') + endif + + do n=1,ntiles + BVISDR(n) = amax1(1.e-6, BVISDR(n)) + BVISDF(n) = amax1(1.e-6, BVISDF(n)) + BNIRDR(n) = amax1(1.e-6, BNIRDR(n)) + BNIRDF(n) = amax1(1.e-6, BNIRDF(n)) + + ! convert % to fractions + + CLMC_pf1(n) = CLMC_pf1(n) / 100. + CLMC_pf2(n) = CLMC_pf2(n) / 100. + CLMC_sf1(n) = CLMC_sf1(n) / 100. + CLMC_sf2(n) = CLMC_sf2(n) / 100. + + fvg(1) = CLMC_pf1(n) + fvg(2) = CLMC_pf2(n) + fvg(3) = CLMC_sf1(n) + fvg(4) = CLMC_sf2(n) + + BARE = 1. + + DO NV = 1, NVEG + BARE = BARE - FVG(NV)! subtract vegetated fractions + END DO + + if (BARE /= 0.) THEN + IB = MAXLOC(FVG(:),1) + FVG (IB) = FVG(IB) + BARE ! This also corrects all cases sum ne 0. + ENDIF + + CLMC_pf1(n) = fvg(1) + CLMC_pf2(n) = fvg(2) + CLMC_sf1(n) = fvg(3) + CLMC_sf2(n) = fvg(4) + enddo + + NDEP = NDEP * 1.e-9 + + ! prevent trivial fractions + ! ------------------------- + do n = 1,ntiles + if(CLMC_pf1(n) <= 1.e-4) then + CLMC_pf2(n) = CLMC_pf2(n) + CLMC_pf1(n) + CLMC_pf1(n) = 0. + endif + + if(CLMC_pf2(n) <= 1.e-4) then + CLMC_pf1(n) = CLMC_pf1(n) + CLMC_pf2(n) + CLMC_pf2(n) = 0. + endif + + if(CLMC_sf1(n) <= 1.e-4) then + if(CLMC_sf2(n) > 1.e-4) then + CLMC_sf2(n) = CLMC_sf2(n) + CLMC_sf1(n) + else if(CLMC_pf2(n) > 1.e-4) then + CLMC_pf2(n) = CLMC_pf2(n) + CLMC_sf1(n) + else if(CLMC_pf1(n) > 1.e-4) then + CLMC_pf1(n) = CLMC_pf1(n) + CLMC_sf1(n) + else + stop 'fveg3' + endif + CLMC_sf1(n) = 0. + endif + + if(CLMC_sf2(n) <= 1.e-4) then + if(CLMC_sf1(n) > 1.e-4) then + CLMC_sf1(n) = CLMC_sf1(n) + CLMC_sf2(n) + else if(CLMC_pf2(n) > 1.e-4) then + CLMC_pf2(n) = CLMC_pf2(n) + CLMC_sf2(n) + else if(CLMC_pf1(n) > 1.e-4) then + CLMC_pf1(n) = CLMC_pf1(n) + CLMC_sf2(n) + else + stop 'fveg4' + endif + CLMC_sf2(n) = 0. + endif + enddo + + this%cnity = reshape([CLMC_pt1,CLMC_pt2,CLMC_st1,CLMC_st2],[ntiles,4]) + this%fvg = reshape([CLMC_pf1,CLMC_pf2,CLMC_sf1,CLMC_sf2],[ntiles,4]) + + this%ndep = ndep + this%t2 = t2 + this%BGALBVR = BVISDR + this%BGALBVF = BVISDF + this%BGALBNR = BNIRDR + this%BGALBNF = BNIRDF + + if (this%isCLM45) then + this%abm = real(abm) + this%fieldcap = fc + this%hdm = hdm + this%gdp = gdp + this%peatf = peatf + endif + + deallocate (BVISDR, BVISDF, BNIRDR ) + deallocate (BNIRDF, T2, NDEP ) + deallocate (CLMC_pf1, CLMC_pf2, CLMC_sf1) + deallocate (CLMC_sf2, CLMC_pt1, CLMC_pt2) + deallocate (CLMC_st1,CLMC_st2) + + _RETURN(_SUCCESS) + end subroutine add_bcs_to_cnrst + + subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) + class(CatchmentCNRst), intent(inout) :: this + character(*), intent(in) :: InTileFile + character(*), intent(in) :: OutBcsDir + character(*), intent(in) :: OutTileFile + real, intent(in) :: surflay + integer, optional, intent(out) :: rc + + real , allocatable, dimension (:) :: DAYX + integer, allocatable, dimension (:) :: low_ind, upp_ind, nt_local + integer, allocatable, dimension (:,:) :: Id_glb_cn, id_loc_cn + integer, allocatable, dimension (:) :: tid_offl, id_loc + real, allocatable, dimension (:) :: CLMC_pf1, CLMC_pf2, CLMC_sf1, CLMC_sf2, & + CLMC_pt1, CLMC_pt2,CLMC_st1,CLMC_st2, var_dum2, var_dum3 + integer :: AGCM_YY,AGCM_MM,AGCM_DD,AGCM_HR=0,AGCM_DATE + real, allocatable, dimension(:,:) :: fveg_offl, ityp_offl, tg_tmp + real, allocatable :: var_off_col (:,:,:), var_off_pft (:,:,:,:) + integer :: status, in_ntiles, out_ntiles, numprocs + logical :: root_proc + integer :: mpierr, n, i, k, tag, req, st, ed, myid, L, iv, nv,nz, var_col, var_pft + character(*), parameter :: Iam = "CatchmentCN::Re_tile" + + + in_ntiles = this%ntiles + var_pft = this%var_pft + var_col = this%var_col + call this%CatchmentRst%re_tile(InTileFile, OutBcsDir, OutTileFile, surflay, _RC) + + out_ntiles = this%ntiles + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, mpierr ) + + root_proc = .false. + if (myid == 0) root_proc = .true. + + allocate(low_ind (numprocs)) + allocate(upp_ind (numprocs)) + allocate(nt_local(numprocs)) + low_ind (:) = 1 + upp_ind (:) = out_ntiles + nt_local(:) = out_ntiles + + if (numprocs > 1) then + do i = 1, numprocs - 1 + upp_ind(i) = low_ind(i) + (out_ntiles/numprocs) - 1 + low_ind(i+1) = upp_ind(i) + 1 + nt_local(i) = upp_ind(i) - low_ind(i) + 1 + end do + nt_local(numprocs) = upp_ind(numprocs) - low_ind(numprocs) + 1 + endif + + allocate (CLMC_pf1(nt_local (myid + 1))) + allocate (CLMC_pf2(nt_local (myid + 1))) + allocate (CLMC_sf1(nt_local (myid + 1))) + allocate (CLMC_sf2(nt_local (myid + 1))) + allocate (CLMC_pt1(nt_local (myid + 1))) + allocate (CLMC_pt2(nt_local (myid + 1))) + allocate (CLMC_st1(nt_local (myid + 1))) + allocate (CLMC_st2(nt_local (myid + 1))) + allocate (ityp_offl (in_ntiles,nveg)) + allocate (fveg_offl (in_ntiles,nveg)) + allocate (tid_offl(in_ntiles)) + allocate (id_loc_cn (nt_local (myid + 1),nveg)) + + do n = 1, in_ntiles + tid_offl(n) = n + enddo + + if (root_proc) then + + allocate (DAYX (out_ntiles)) + + READ(this%time(1:8),'(I8)') AGCM_DATE + AGCM_YY = AGCM_DATE / 10000 + AGCM_MM = (AGCM_DATE - AGCM_YY*10000) / 100 + AGCM_DD = (AGCM_DATE - AGCM_YY*10000 - AGCM_MM*100) + + call compute_dayx ( & + out_NTILES, AGCM_YY, AGCM_MM, AGCM_DD, AGCM_HR, & + this%LATG, DAYX) + + ! save the old vaues dimension (in_ntiles, nv) + ityp_offl = this%cnity + fveg_offl = this%fvg + + do n = 1, in_ntiles + do nv = 1,nveg + if(ityp_offl(n,nv)<0 .or. ityp_offl(n,nv)>npft) stop 'ityp' + if(fveg_offl(n,nv)<0..or. fveg_offl(n,nv)>1.00001) stop 'fveg' + end do + + if (nint(this%tile_id(n)) /= n) stop ("cannot assign ity_offl to cnity and fvg_offl to fvg") + + if((ityp_offl(N,3) == 0).and.(ityp_offl(N,4) == 0)) then + if(ityp_offl(N,1) /= 0) then + ityp_offl(N,3) = ityp_offl(N,1) + else + ityp_offl(N,3) = ityp_offl(N,2) + endif + endif + + if((ityp_offl(N,1) == 0).and.(ityp_offl(N,2) /= 0)) ityp_offl(N,1) = ityp_offl(N,2) + if((ityp_offl(N,2) == 0).and.(ityp_offl(N,1) /= 0)) ityp_offl(N,2) = ityp_offl(N,1) + if((ityp_offl(N,3) == 0).and.(ityp_offl(N,4) /= 0)) ityp_offl(N,3) = ityp_offl(N,4) + if((ityp_offl(N,4) == 0).and.(ityp_offl(N,3) /= 0)) ityp_offl(N,4) = ityp_offl(N,3) + end do + endif + + call MPI_BCAST(ityp_offl,size(ityp_offl),MPI_REAL ,0,MPI_COMM_WORLD,mpierr) + call MPI_BCAST(fveg_offl,size(fveg_offl),MPI_REAL ,0,MPI_COMM_WORLD,mpierr) + + if (root_proc ) then + + ! after this call, the cnity and fvg is the dimension of (out_ntiles, nveg) + call this%add_bcs_to_cnrst(surflay, OutBcsDir, __RC__) + + do i = 1, numprocs -1 + st = low_ind(i+1) + l = nt_local(i+1) + tag = i*numprocs + call MPI_send(this%cnity(st,1),l, MPI_REAL, i, tag, MPI_COMM_WORLD, mpierr) + call MPI_send(this%cnity(st,2),l, MPI_REAL, i, tag+1, MPI_COMM_WORLD, mpierr) + call MPI_send(this%cnity(st,3),l, MPI_REAL, i, tag+2, MPI_COMM_WORLD, mpierr) + call MPI_send(this%cnity(st,4),l, MPI_REAL, i, tag+3, MPI_COMM_WORLD, mpierr) + call MPI_send(this%fvg(st,1),l, MPI_REAL, i, tag+4, MPI_COMM_WORLD, mpierr) + call MPI_send(this%fvg(st,2),l, MPI_REAL, i, tag+5, MPI_COMM_WORLD, mpierr) + call MPI_send(this%fvg(st,3),l, MPI_REAL, i, tag+6, MPI_COMM_WORLD, mpierr) + call MPI_send(this%fvg(st,4),l, MPI_REAL, i, tag+7, MPI_COMM_WORLD, mpierr) + enddo + st = low_ind(1) + l = nt_local(1) + ed = st + l -1 + CLMC_pt1 = this%cnity(st:ed,1) + CLMC_pt2 = this%cnity(st:ed,2) + CLMC_st1 = this%cnity(st:ed,3) + CLMC_st2 = this%cnity(st:ed,4) + CLMC_pf1 = this%fvg(st:ed,1) + CLMC_pf2 = this%fvg(st:ed,2) + CLMC_sf1 = this%fvg(st:ed,3) + CLMC_sf2 = this%fvg(st:ed,4) + else + tag = myid*numprocs + call MPI_RECV(CLMC_pt1,nt_local(myid+1) , MPI_REAL, 0, tag, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_pt2,nt_local(myid+1) , MPI_REAL, 0, tag+1, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_st1,nt_local(myid+1) , MPI_REAL, 0, tag+2, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_st2,nt_local(myid+1) , MPI_REAL, 0, tag+3, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_pf1,nt_local(myid+1) , MPI_REAL, 0, tag+4, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_pf2,nt_local(myid+1) , MPI_REAL, 0, tag+5, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_sf1,nt_local(myid+1) , MPI_REAL, 0, tag+6, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(CLMC_sf2,nt_local(myid+1) , MPI_REAL, 0, tag+7, MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + endif + + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + + if(root_proc) print*, "GetIDs...." + + call GetIds(this%lonc,this%latc,this%lonn,this%latt,id_loc_cn, tid_offl, & + CLMC_pf1, CLMC_pf2, CLMC_sf1, CLMC_sf2, CLMC_pt1, CLMC_pt2,CLMC_st1,CLMC_st2, & + fveg_offl, ityp_offl) + + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + + if(root_proc) allocate (id_glb_cn (out_ntiles,nveg)) + + allocate (id_loc (out_ntiles)) + deallocate (CLMC_pf1, CLMC_pf2, CLMC_sf1, CLMC_sf2) + deallocate (CLMC_pt1, CLMC_pt2, CLMC_st1, CLMC_st2) + + do nv = 1, nveg + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + do i = 1, numprocs + if((I == 1).and.(myid == 0)) then + id_loc(low_ind(i) : upp_ind(i)) = Id_loc_cn(:,nv) + else if (I > 1) then + if(I-1 == myid) then + ! send to root + call MPI_ISend(id_loc_cn(:,nv),nt_local(i),MPI_INTEGER,0,994,MPI_COMM_WORLD,req,mpierr) + call MPI_WAIT (req,MPI_STATUS_IGNORE,mpierr) + else if (myid == 0) then + ! root receives + call MPI_RECV(id_loc(low_ind(i) : upp_ind(i)),nt_local(i) , MPI_INTEGER, i-1,994,MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + endif + endif + end do + + if(root_proc) id_glb_cn (:,nv) = id_loc + + end do + + if(root_proc) then + + allocate (var_off_col (1: in_ntiles, 1 : nzone,1 : var_col)) + allocate (var_off_pft (1: in_ntiles, 1 : nzone,1 : nveg, 1 : var_pft)) + allocate (var_dum2 (1:in_ntiles)) + + this%tile_id = [(i*1.0, i=1, out_ntiles)] + + allocate (tg_tmp(out_ntiles, 4),source = 0.) + do i = 1, 3 + tg_tmp(:,i) = this%tg(this%id_glb(:),i) + enddo + this%tg = tg_tmp + + i = 1 + do nv = 1,VAR_COL + do nz = 1,nzone + var_off_col(:,nz,nv) = this%cncol(:,i) + i = i + 1 + end do + end do + + i = 1 + do iv = 1,VAR_PFT + do nv = 1,nveg + do nz = 1,nzone + var_off_pft(:, nz,nv,iv) = this%cnpft(:,i) + i = i + 1 + end do + end do + end do + + where(isnan(var_off_pft)) var_off_pft = 0. + where(var_off_pft /= var_off_pft) var_off_pft = 0. + + print *, 'calculating regridded carbn' + + call regrid_carbon (out_NTILES, in_ntiles,id_glb_cn, & + DAYX, var_off_col,var_off_pft, ityp_offl, fveg_offl) + deallocate (var_off_col,var_off_pft) + endif + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + + _RETURN(_SUCCESS) + + contains + SUBROUTINE regrid_carbon (NTILES, in_ntiles, id_glb, & + DAYX, var_off_col, var_off_pft, ityp_offl, fveg_offl) + + ! write out regridded carbon variables + implicit none + integer, intent (in) :: NTILES, in_ntiles,id_glb (ntiles,nveg) + real, intent (in) :: DAYX (NTILES), var_off_col(in_ntiles,NZONE,var_col), var_off_pft(in_ntiles,NZONE, NVEG, var_pft) + real, intent (in), dimension(in_ntiles,nveg) :: fveg_offl, ityp_offl + real, allocatable, dimension (:) :: CLMC_pf1, CLMC_pf2, CLMC_sf1, CLMC_sf2, & + CLMC_pt1, CLMC_pt2,CLMC_st1,CLMC_st2, var_dum + real, allocatable :: var_col_out (:,:,:), var_pft_out (:,:,:,:) + integer :: N, STATUS, nv, nx, offl_cell, ityp_new, i, j, nz, iv + real :: fveg_new + character(256) :: Iam = "write_regridded_carbon" + + + allocate (CLMC_pf1(NTILES)) + allocate (CLMC_pf2(NTILES)) + allocate (CLMC_sf1(NTILES)) + allocate (CLMC_sf2(NTILES)) + allocate (CLMC_pt1(NTILES)) + allocate (CLMC_pt2(NTILES)) + allocate (CLMC_st1(NTILES)) + allocate (CLMC_st2(NTILES)) + allocate (VAR_DUM (NTILES)) + + CLMC_pt1 = this%cnity(:,1) + CLMC_pt2 = this%cnity(:,2) + CLMC_st1 = this%cnity(:,3) + CLMC_st2 = this%cnity(:,4) + CLMC_pf1 = this%fvg(:,1) + CLMC_pf2 = this%fvg(:,2) + CLMC_sf1 = this%fvg(:,3) + CLMC_sf2 = this%fvg(:,4) + + allocate (var_col_out (1: NTILES, 1 : nzone,1 : var_col)) + allocate (var_pft_out (1: NTILES, 1 : nzone,1 : nveg, 1 : var_pft)) + + var_col_out = 0. + var_pft_out = NaN + + OUT_TILE : DO N = 1, NTILES + + ! if(mod (n,1000) == 0) print *, myid +1, n, Id_glb(n,:) + + NVLOOP2 : do nv = 1, nveg + + if(nv <= 2) then ! index for secondary PFT index if primary or primary if secondary + nx = nv + 2 + else + nx = nv - 2 + endif + + if (nv == 1) ityp_new = CLMC_pt1(n) + if (nv == 1) fveg_new = CLMC_pf1(n) + if (nv == 2) ityp_new = CLMC_pt2(n) + if (nv == 2) fveg_new = CLMC_pf2(n) + if (nv == 3) ityp_new = CLMC_st1(n) + if (nv == 3) fveg_new = CLMC_sf1(n) + if (nv == 4) ityp_new = CLMC_st2(n) + if (nv == 4) fveg_new = CLMC_sf2(n) + if (fveg_new > fmin) then + + offl_cell = Id_glb(n,nv) + + if(ityp_new == ityp_offl (offl_cell,nv) .and. fveg_offl (offl_cell,nv)> fmin) then + iv = nv ! same type fraction (primary of secondary) + else if(ityp_new == ityp_offl (offl_cell,nx) .and. fveg_offl (offl_cell,nx)> fmin) then + iv = nx ! not same fraction + else if(iclass(ityp_new)==iclass(ityp_offl(offl_cell,nv)) .and. fveg_offl (offl_cell,nv)> fmin) then + iv = nv ! primary, other type (same class) + else if(fveg_offl (offl_cell,nx)> fmin) then + iv = nx ! secondary, other type (same class) + endif + + ! Get col and pft variables for the Id_glb(nv) grid cell from offline catchcn_internal_rst + ! ---------------------------------------------------------------------------------------- + + ! call NCDF_reshape_getOput (NCFID,Id_glb(n,nv),var_off_col,var_off_pft,.true.) + + var_pft_out (n,:,nv,:) = var_off_pft(Id_glb(n,nv), :,iv,:) + var_col_out (n,:,:) = var_col_out(n,:,:) + fveg_new * var_off_col(Id_glb(n,nv), :,:) ! gkw: column state simple weighted mean; ! could use "woody" fraction? + + ! Check whether var_pft_out is realistic + do nz = 1, nzone + do j = 1, VAR_PFT + if (isnan(var_pft_out (n, nz,nv,j))) print *,j,nv,nz,n,var_pft_out (n, nz,nv,j),fveg_new + !if(isnan(var_pft_out (n, nz,nv,69))) var_pft_out (n, nz,nv,69) = 1.e-6 + !if(isnan(var_pft_out (n, nz,nv,70))) var_pft_out (n, nz,nv,70) = 1.e-6 + !if(isnan(var_pft_out (n, nz,nv,73))) var_pft_out (n, nz,nv,73) = 1.e-6 + !if(isnan(var_pft_out (n, nz,nv,74))) var_pft_out (n, nz,nv,74) = 1.e-6 + end do + end do + endif + + end do NVLOOP2 + + ! reset carbon if negative < 10g + ! ------------------------ + + NZLOOP : do nz = 1, nzone + + if(var_col_out (n, nz,14) < 10.) then + + var_col_out(n, nz, 1) = max(var_col_out(n, nz, 1), 0.) + var_col_out(n, nz, 2) = max(var_col_out(n, nz, 2), 0.) + var_col_out(n, nz, 3) = max(var_col_out(n, nz, 3), 0.) + var_col_out(n, nz, 4) = max(var_col_out(n, nz, 4), 0.) + var_col_out(n, nz, 5) = max(var_col_out(n, nz, 5), 0.) + var_col_out(n, nz,10) = max(var_col_out(n, nz,10), 0.) + var_col_out(n, nz,11) = max(var_col_out(n, nz,11), 0.) + var_col_out(n, nz,12) = max(var_col_out(n, nz,12), 0.) + var_col_out(n, nz,13) = max(var_col_out(n, nz,13),10.) ! soil4c + var_col_out(n, nz,14) = max(var_col_out(n, nz,14), 0.) + var_col_out(n, nz,15) = max(var_col_out(n, nz,15), 0.) + var_col_out(n, nz,16) = max(var_col_out(n, nz,16), 0.) + var_col_out(n, nz,17) = max(var_col_out(n, nz,17), 0.) + var_col_out(n, nz,18) = max(var_col_out(n, nz,18), 0.) + var_col_out(n, nz,19) = max(var_col_out(n, nz,19), 0.) + var_col_out(n, nz,20) = max(var_col_out(n, nz,20), 0.) + var_col_out(n, nz,24) = max(var_col_out(n, nz,24), 0.) + var_col_out(n, nz,25) = max(var_col_out(n, nz,25), 0.) + var_col_out(n, nz,26) = max(var_col_out(n, nz,26), 0.) + var_col_out(n, nz,27) = max(var_col_out(n, nz,27), 0.) + var_col_out(n, nz,28) = max(var_col_out(n, nz,28), 1.) + var_col_out(n, nz,29) = max(var_col_out(n, nz,29), 0.) + + NVLOOP3 : do nv = 1,nveg + + if (nv == 1) ityp_new = CLMC_pt1(n) + if (nv == 1) fveg_new = CLMC_pf1(n) + if (nv == 2) ityp_new = CLMC_pt2(n) + if (nv == 2) fveg_new = CLMC_pf2(n) + if (nv == 3) ityp_new = CLMC_st1(n) + if (nv == 3) fveg_new = CLMC_sf1(n) + if (nv == 4) ityp_new = CLMC_st2(n) + if (nv == 4) fveg_new = CLMC_sf2(n) + + if(fveg_new > fmin) then + var_pft_out(n, nz,nv, 1) = max(var_pft_out(n, nz,nv, 1),0.) + var_pft_out(n, nz,nv, 2) = max(var_pft_out(n, nz,nv, 2),0.) + var_pft_out(n, nz,nv, 3) = max(var_pft_out(n, nz,nv, 3),0.) + var_pft_out(n, nz,nv, 4) = max(var_pft_out(n, nz,nv, 4),0.) + + if(ityp_new <= 12) then ! tree or shrub deadstemc + var_pft_out(n, nz,nv, 5) = max(var_pft_out(n, nz,nv, 5),0.1) + else + var_pft_out(n, nz,nv, 5) = max(var_pft_out(n, nz,nv, 5),0.0) + endif + + var_pft_out(n, nz,nv, 6) = max(var_pft_out(n, nz,nv, 6),0.) + var_pft_out(n, nz,nv, 7) = max(var_pft_out(n, nz,nv, 7),0.) + var_pft_out(n, nz,nv, 8) = max(var_pft_out(n, nz,nv, 8),0.) + var_pft_out(n, nz,nv, 9) = max(var_pft_out(n, nz,nv, 9),0.) + var_pft_out(n, nz,nv,10) = max(var_pft_out(n, nz,nv,10),0.) + var_pft_out(n, nz,nv,11) = max(var_pft_out(n, nz,nv,11),0.) + var_pft_out(n, nz,nv,12) = max(var_pft_out(n, nz,nv,12),0.) + + if(ityp_new <=2 .or. ityp_new ==4 .or. ityp_new ==5 .or. ityp_new == 9) then + var_pft_out(n, nz,nv,13) = max(var_pft_out(n, nz,nv,13),1.) ! leaf carbon display for evergreen + var_pft_out(n, nz,nv,14) = max(var_pft_out(n, nz,nv,14),0.) + else + var_pft_out(n, nz,nv,13) = max(var_pft_out(n, nz,nv,13),0.) + var_pft_out(n, nz,nv,14) = max(var_pft_out(n, nz,nv,14),1.) ! leaf carbon storage for deciduous + endif + + var_pft_out(n, nz,nv,15) = max(var_pft_out(n, nz,nv,15),0.) + var_pft_out(n, nz,nv,16) = max(var_pft_out(n, nz,nv,16),0.) + var_pft_out(n, nz,nv,17) = max(var_pft_out(n, nz,nv,17),0.) + var_pft_out(n, nz,nv,18) = max(var_pft_out(n, nz,nv,18),0.) + var_pft_out(n, nz,nv,19) = max(var_pft_out(n, nz,nv,19),0.) + var_pft_out(n, nz,nv,20) = max(var_pft_out(n, nz,nv,20),0.) + var_pft_out(n, nz,nv,21) = max(var_pft_out(n, nz,nv,21),0.) + var_pft_out(n, nz,nv,22) = max(var_pft_out(n, nz,nv,22),0.) + var_pft_out(n, nz,nv,23) = max(var_pft_out(n, nz,nv,23),0.) + var_pft_out(n, nz,nv,25) = max(var_pft_out(n, nz,nv,25),0.) + var_pft_out(n, nz,nv,26) = max(var_pft_out(n, nz,nv,26),0.) + var_pft_out(n, nz,nv,27) = max(var_pft_out(n, nz,nv,27),0.) + var_pft_out(n, nz,nv,41) = max(var_pft_out(n, nz,nv,41),0.) + var_pft_out(n, nz,nv,42) = max(var_pft_out(n, nz,nv,42),0.) + var_pft_out(n, nz,nv,44) = max(var_pft_out(n, nz,nv,44),0.) + var_pft_out(n, nz,nv,45) = max(var_pft_out(n, nz,nv,45),0.) + var_pft_out(n, nz,nv,46) = max(var_pft_out(n, nz,nv,46),0.) + var_pft_out(n, nz,nv,47) = max(var_pft_out(n, nz,nv,47),0.) + var_pft_out(n, nz,nv,48) = max(var_pft_out(n, nz,nv,48),0.) + var_pft_out(n, nz,nv,49) = max(var_pft_out(n, nz,nv,49),0.) + var_pft_out(n, nz,nv,50) = max(var_pft_out(n, nz,nv,50),0.) + var_pft_out(n, nz,nv,51) = max(var_pft_out(n, nz,nv, 5)/500.,0.) + var_pft_out(n, nz,nv,52) = max(var_pft_out(n, nz,nv,52),0.) + var_pft_out(n, nz,nv,53) = max(var_pft_out(n, nz,nv,53),0.) + var_pft_out(n, nz,nv,54) = max(var_pft_out(n, nz,nv,54),0.) + var_pft_out(n, nz,nv,55) = max(var_pft_out(n, nz,nv,55),0.) + var_pft_out(n, nz,nv,56) = max(var_pft_out(n, nz,nv,56),0.) + var_pft_out(n, nz,nv,57) = max(var_pft_out(n, nz,nv,13)/25.,0.) + var_pft_out(n, nz,nv,58) = max(var_pft_out(n, nz,nv,14)/25.,0.) + var_pft_out(n, nz,nv,59) = max(var_pft_out(n, nz,nv,59),0.) + var_pft_out(n, nz,nv,60) = max(var_pft_out(n, nz,nv,60),0.) + var_pft_out(n, nz,nv,61) = max(var_pft_out(n, nz,nv,61),0.) + var_pft_out(n, nz,nv,62) = max(var_pft_out(n, nz,nv,62),0.) + var_pft_out(n, nz,nv,63) = max(var_pft_out(n, nz,nv,63),0.) + var_pft_out(n, nz,nv,64) = max(var_pft_out(n, nz,nv,64),0.) + var_pft_out(n, nz,nv,65) = max(var_pft_out(n, nz,nv,65),0.) + var_pft_out(n, nz,nv,66) = max(var_pft_out(n, nz,nv,66),0.) + var_pft_out(n, nz,nv,67) = max(var_pft_out(n, nz,nv,67),0.) + var_pft_out(n, nz,nv,68) = max(var_pft_out(n, nz,nv,68),0.) + var_pft_out(n, nz,nv,69) = max(var_pft_out(n, nz,nv,69),0.) + var_pft_out(n, nz,nv,70) = max(var_pft_out(n, nz,nv,70),0.) + var_pft_out(n, nz,nv,73) = max(var_pft_out(n, nz,nv,73),0.) + var_pft_out(n, nz,nv,74) = max(var_pft_out(n, nz,nv,74),0.) + if(this%isCLM45) var_pft_out(n, nz,nv,75) = max(var_pft_out(n, nz,nv,75),0.) + endif + end do NVLOOP3 ! end veg loop + endif ! end carbon check + end do NZLOOP ! end zone loop + + ! Update dayx variable var_pft_out (:,:,28) + + do j = 28, 28 ! 1,VAR_PFT var_pft_out (:,:,:,28) + do nv = 1,nveg + do nz = 1,nzone + var_pft_out (n, nz,nv,j) = dayx(n) + end do + end do + end do + + ! call NCDF_reshape_getOput (OutID,N,var_col_out,var_pft_out,.false.) + + ! column vars clm40 clm45 + ! ----------------- --------------------- + ! 1 clm3%g%l%c%ccs%col_ctrunc ! 1 ccs%col_ctrunc_vr (:,1) + ! 2 clm3%g%l%c%ccs%cwdc ! 2 ccs%decomp_cpools_vr(:,1,4) ! cwdc + ! 3 clm3%g%l%c%ccs%litr1c ! 3 ccs%decomp_cpools_vr(:,1,1) ! litr1c + ! 4 clm3%g%l%c%ccs%litr2c ! 4 ccs%decomp_cpools_vr(:,1,2) ! litr2c + ! 5 clm3%g%l%c%ccs%litr3c ! 5 ccs%decomp_cpools_vr(:,1,3) ! litr3c + ! 6 clm3%g%l%c%ccs%pcs_a%totvegc ! 6 ccs%totvegc_col + ! 7 clm3%g%l%c%ccs%prod100c ! 7 ccs%prod100c + ! 8 clm3%g%l%c%ccs%prod10c ! 8 ccs%prod10c + ! 9 clm3%g%l%c%ccs%seedc ! 9 ccs%seedc + ! 10 clm3%g%l%c%ccs%soil1c ! 10 ccs%decomp_cpools_vr(:,1,5) ! soil1c + ! 11 clm3%g%l%c%ccs%soil2c ! 11 ccs%decomp_cpools_vr(:,1,6) ! soil2c + ! 12 clm3%g%l%c%ccs%soil3c ! 12 ccs%decomp_cpools_vr(:,1,7) ! soil3c + ! 13 clm3%g%l%c%ccs%soil4c ! 13 ccs%decomp_cpools_vr(:,1,8) ! soil4c + ! 14 clm3%g%l%c%ccs%totcolc ! 14 ccs%totcolc + ! 15 clm3%g%l%c%ccs%totlitc ! 15 ccs%totlitc + ! 16 clm3%g%l%c%cns%col_ntrunc ! 16 cns%col_ntrunc_vr (:,1) + ! 17 clm3%g%l%c%cns%cwdn ! 17 cns%decomp_npools_vr(:,1,4) ! cwdn + ! 18 clm3%g%l%c%cns%litr1n ! 18 cns%decomp_npools_vr(:,1,1) ! litr1n + ! 19 clm3%g%l%c%cns%litr2n ! 19 cns%decomp_npools_vr(:,1,2) ! litr2n + ! 20 clm3%g%l%c%cns%litr3n ! 20 cns%decomp_npools_vr(:,1,3) ! litr3n + ! 21 clm3%g%l%c%cns%prod100n ! 21 cns%prod100n + ! 22 clm3%g%l%c%cns%prod10n ! 22 cns%prod10n + ! 23 clm3%g%l%c%cns%seedn ! 23 cns%seedn + ! 24 clm3%g%l%c%cns%sminn ! 24 cns%sminn_vr (:,1) + ! 25 clm3%g%l%c%cns%soil1n ! 25 cns%decomp_npools_vr(:,1,5) ! soil1n + ! 26 clm3%g%l%c%cns%soil2n ! 26 cns%decomp_npools_vr(:,1,6) ! soil2n + ! 27 clm3%g%l%c%cns%soil3n ! 27 cns%decomp_npools_vr(:,1,7) ! soil3n + ! 28 clm3%g%l%c%cns%soil4n ! 28 cns%decomp_npools_vr(:,1,8) ! soil4n + ! 29 clm3%g%l%c%cns%totcoln ! 29 cns%totcoln + ! 30 clm3%g%l%c%cps%ann_farea_burned ! 30 cps%fpg + ! 31 clm3%g%l%c%cps%annsum_counter ! 31 cps%annsum_counter + ! 32 clm3%g%l%c%cps%cannavg_t2m ! 32 cps%cannavg_t2m + ! 33 clm3%g%l%c%cps%cannsum_npp ! 33 cps%cannsum_npp + ! 34 clm3%g%l%c%cps%farea_burned ! 34 cps%farea_burned + ! 35 clm3%g%l%c%cps%fire_prob ! 35 cps%fpi_vr (:,1) + ! 36 clm3%g%l%c%cps%fireseasonl ! OLD ! 30 cps%altmax + ! 37 clm3%g%l%c%cps%fpg ! OLD ! 31 cps%annsum_counter + ! 38 clm3%g%l%c%cps%fpi ! OLD ! 32 cps%cannavg_t2m + ! 39 clm3%g%l%c%cps%me ! OLD ! 33 cps%cannsum_npp + ! 40 clm3%g%l%c%cps%mean_fire_prob ! OLD ! 34 cps%farea_burned + ! OLD ! 35 cps%altmax_lastyear + ! OLD ! 36 cps%altmax_indx + ! OLD ! 37 cps%fpg + ! OLD ! 38 cps%fpi_vr (:,1) + ! OLD ! 39 cps%altmax_lastyear_indx + + ! PFT vars CLM40 CLM45 + ! -------------- ----- + ! 1 clm3%g%l%c%p%pcs%cpool ! 1 pcs%cpool + ! 2 clm3%g%l%c%p%pcs%deadcrootc ! 2 pcs%deadcrootc + ! 3 clm3%g%l%c%p%pcs%deadcrootc_storage ! 3 pcs%deadcrootc_storage + ! 4 clm3%g%l%c%p%pcs%deadcrootc_xfer ! 4 pcs%deadcrootc_xfer + ! 5 clm3%g%l%c%p%pcs%deadstemc ! 5 pcs%deadstemc + ! 6 clm3%g%l%c%p%pcs%deadstemc_storage ! 6 pcs%deadstemc_storage + ! 7 clm3%g%l%c%p%pcs%deadstemc_xfer ! 7 pcs%deadstemc_xfer + ! 8 clm3%g%l%c%p%pcs%frootc ! 8 pcs%frootc + ! 9 clm3%g%l%c%p%pcs%frootc_storage ! 9 pcs%frootc_storage + ! 10 clm3%g%l%c%p%pcs%frootc_xfer ! 10 pcs%frootc_xfer + ! 11 clm3%g%l%c%p%pcs%gresp_storage ! 11 pcs%gresp_storage + ! 12 clm3%g%l%c%p%pcs%gresp_xfer ! 12 pcs%gresp_xfer + ! 13 clm3%g%l%c%p%pcs%leafc ! 13 pcs%leafc + ! 14 clm3%g%l%c%p%pcs%leafc_storage ! 14 pcs%leafc_storage + ! 15 clm3%g%l%c%p%pcs%leafc_xfer ! 15 pcs%leafc_xfer + ! 16 clm3%g%l%c%p%pcs%livecrootc ! 16 pcs%livecrootc + ! 17 clm3%g%l%c%p%pcs%livecrootc_storage ! 17 pcs%livecrootc_storage + ! 18 clm3%g%l%c%p%pcs%livecrootc_xfer ! 18 pcs%livecrootc_xfer + ! 19 clm3%g%l%c%p%pcs%livestemc ! 19 pcs%livestemc + ! 20 clm3%g%l%c%p%pcs%livestemc_storage ! 20 pcs%livestemc_storage + ! 21 clm3%g%l%c%p%pcs%livestemc_xfer ! 21 pcs%livestemc_xfer + ! 22 clm3%g%l%c%p%pcs%pft_ctrunc ! 22 pcs%pft_ctrunc + ! 23 clm3%g%l%c%p%pcs%xsmrpool ! 23 pcs%xsmrpool + ! 24 clm3%g%l%c%p%pepv%annavg_t2m ! 24 pepv%annavg_t2m + ! 25 clm3%g%l%c%p%pepv%annmax_retransn ! 25 pepv%annmax_retransn + ! 26 clm3%g%l%c%p%pepv%annsum_npp ! 26 pepv%annsum_npp + ! 27 clm3%g%l%c%p%pepv%annsum_potential_gpp ! 27 pepv%annsum_potential_gpp + ! 28 clm3%g%l%c%p%pepv%dayl ! 28 pepv%dayl + ! 29 clm3%g%l%c%p%pepv%days_active ! 29 pepv%days_active + ! 30 clm3%g%l%c%p%pepv%dormant_flag ! 30 pepv%dormant_flag + ! 31 clm3%g%l%c%p%pepv%offset_counter ! 31 pepv%offset_counter + ! 32 clm3%g%l%c%p%pepv%offset_fdd ! 32 pepv%offset_fdd + ! 33 clm3%g%l%c%p%pepv%offset_flag ! 33 pepv%offset_flag + ! 34 clm3%g%l%c%p%pepv%offset_swi ! 34 pepv%offset_swi + ! 35 clm3%g%l%c%p%pepv%onset_counter ! 35 pepv%onset_counter + ! 36 clm3%g%l%c%p%pepv%onset_fdd ! 36 pepv%onset_fdd + ! 37 clm3%g%l%c%p%pepv%onset_flag ! 37 pepv%onset_flag + ! 38 clm3%g%l%c%p%pepv%onset_gdd ! 38 pepv%onset_gdd + ! 39 clm3%g%l%c%p%pepv%onset_gddflag ! 39 pepv%onset_gddflag + ! 40 clm3%g%l%c%p%pepv%onset_swi ! 40 pepv%onset_swi + ! 41 clm3%g%l%c%p%pepv%prev_frootc_to_litter ! 41 pepv%prev_frootc_to_litter + ! 42 clm3%g%l%c%p%pepv%prev_leafc_to_litter ! 42 pepv%prev_leafc_to_litter + ! 43 clm3%g%l%c%p%pepv%tempavg_t2m ! 43 pepv%tempavg_t2m + ! 44 clm3%g%l%c%p%pepv%tempmax_retransn ! 44 pepv%tempmax_retransn + ! 45 clm3%g%l%c%p%pepv%tempsum_npp ! 45 pepv%tempsum_npp + ! 46 clm3%g%l%c%p%pepv%tempsum_potential_gpp ! 46 pepv%tempsum_potential_gpp + ! 47 clm3%g%l%c%p%pepv%xsmrpool_recover ! 47 pepv%xsmrpool_recover + ! 48 clm3%g%l%c%p%pns%deadcrootn ! 48 pns%deadcrootn + ! 49 clm3%g%l%c%p%pns%deadcrootn_storage ! 49 pns%deadcrootn_storage + ! 50 clm3%g%l%c%p%pns%deadcrootn_xfer ! 50 pns%deadcrootn_xfer + ! 51 clm3%g%l%c%p%pns%deadstemn ! 51 pns%deadstemn + ! 52 clm3%g%l%c%p%pns%deadstemn_storage ! 52 pns%deadstemn_storage + ! 53 clm3%g%l%c%p%pns%deadstemn_xfer ! 53 pns%deadstemn_xfer + ! 54 clm3%g%l%c%p%pns%frootn ! 54 pns%frootn + ! 55 clm3%g%l%c%p%pns%frootn_storage ! 55 pns%frootn_storage + ! 56 clm3%g%l%c%p%pns%frootn_xfer ! 56 pns%frootn_xfer + ! 57 clm3%g%l%c%p%pns%leafn ! 57 pns%leafn + ! 58 clm3%g%l%c%p%pns%leafn_storage ! 58 pns%leafn_storage + ! 59 clm3%g%l%c%p%pns%leafn_xfer ! 59 pns%leafn_xfer + ! 60 clm3%g%l%c%p%pns%livecrootn ! 60 pns%livecrootn + ! 61 clm3%g%l%c%p%pns%livecrootn_storage ! 61 pns%livecrootn_storage + ! 62 clm3%g%l%c%p%pns%livecrootn_xfer ! 62 pns%livecrootn_xfer + ! 63 clm3%g%l%c%p%pns%livestemn ! 63 pns%livestemn + ! 64 clm3%g%l%c%p%pns%livestemn_storage ! 64 pns%livestemn_storage + ! 65 clm3%g%l%c%p%pns%livestemn_xfer ! 65 pns%livestemn_xfer + ! 66 clm3%g%l%c%p%pns%npool ! 66 pns%npool + ! 67 clm3%g%l%c%p%pns%pft_ntrunc ! 67 pns%pft_ntrunc + ! 68 clm3%g%l%c%p%pns%retransn ! 68 pns%retransn + ! 69 clm3%g%l%c%p%pps%elai ! 69 pps%elai + ! 70 clm3%g%l%c%p%pps%esai ! 70 pps%esai + ! 71 clm3%g%l%c%p%pps%hbot ! 71 pps%hbot + ! 72 clm3%g%l%c%p%pps%htop ! 72 pps%htop + ! 73 clm3%g%l%c%p%pps%tlai ! 73 pps%tlai + ! 74 clm3%g%l%c%p%pps%tsai ! 74 pps%tsai + ! 75 pepv%plant_ndemand + ! OLD ! 75 pps%gddplant + ! OLD ! 76 pps%gddtsoi + ! OLD ! 77 pps%peaklai + ! OLD ! 78 pps%idop + ! OLD ! 79 pps%aleaf + ! OLD ! 80 pps%aleafi + ! OLD ! 81 pps%astem + ! OLD ! 82 pps%astemi + ! OLD ! 83 pps%htmx + ! OLD ! 84 pps%hdidx + ! OLD ! 85 pps%vf + ! OLD ! 86 pps%cumvd + ! OLD ! 87 pps%croplive + ! OLD ! 88 pps%cropplant + ! OLD ! 89 pps%harvdate + ! OLD ! 90 pps%gdd1020 + ! OLD ! 91 pps%gdd820 + ! OLD ! 92 pps%gdd020 + ! OLD ! 93 pps%gddmaturity + ! OLD ! 94 pps%huileaf + ! OLD ! 95 pps%huigrain + ! OLD ! 96 pcs%grainc + ! OLD ! 97 pcs%grainc_storage + ! OLD ! 98 pcs%grainc_xfer + ! OLD ! 99 pns%grainn + ! OLD !100 pns%grainn_storage + ! OLD !101 pns%grainn_xfer + ! OLD !102 pepv%fert_counter + ! OLD !103 pnf%fert + ! OLD !104 pepv%grain_flag + + end do OUT_TILE + + i = 1 + deallocate(this%cncol) + allocate(this%cncol(NTILES, nzone*VAR_COL)) + do nv = 1,VAR_COL + do nz = 1,nzone + this%cncol(:,i) = var_col_out(:, nz,nv) + !STATUS = NF_PUT_VARA_REAL(OutID,VarID(OutID,'CNCOL'), (/1,i/), (/NTILES,1 /),var_col_out(:, nz,nv)) ; VERIFY_(STATUS) + i = i + 1 + end do + end do + + i = 1 + deallocate(this%cnpft) + allocate(this%cnpft(NTILES,VAR_PFT*nveg*nzone)) + if(this%isclm45) then + do iv = 1,VAR_PFT + do nv = 1,nveg + do nz = 1,nzone + if(iv <= 74) then + this%cnpft(:,i) = var_pft_out(:, nz,nv,iv) + !STATUS = NF_PUT_VARA_REAL(OutID,VarID(OutID,'CNPFT'), (/1,i/), (/NTILES,1 /),var_pft_out(:, nz,nv,iv)) ; VERIFY_(STATUS) + else + if((iv == 78) .OR. (iv == 89)) then ! idop and harvdate + var_dum = 999 + this%cnpft(:,i) = var_dum + !STATUS = NF_PUT_VARA_REAL(OutID,VarID(OutID,'CNPFT'), (/1,i/), (/NTILES,1 /),var_dum) ; VERIFY_(STATUS) + else + var_dum = 0. + this%cnpft(:,i) = var_dum + !STATUS = NF_PUT_VARA_REAL(OutID,VarID(OutID,'CNPFT'), (/1,i/), (/NTILES,1 /),var_dum) ; VERIFY_(STATUS) + endif + endif + i = i + 1 + end do + end do + end do + else + do iv = 1,VAR_PFT + do nv = 1,nveg + do nz = 1,nzone + this%cnpft(:,i) = var_pft_out(:, nz,nv,iv) + !STATUS = NF_PUT_VARA_REAL(OutID,VarID(OutID,'CNPFT'), (/1,i/), (/NTILES,1 /),var_pft_out(:, nz,nv,iv)) ; VERIFY_(STATUS) + i = i + 1 + end do + end do + end do + endif + + deallocate (var_col_out,var_pft_out) + deallocate (CLMC_pf1, CLMC_pf2, CLMC_sf1, CLMC_sf2) + deallocate (CLMC_pt1, CLMC_pt2, CLMC_st1, CLMC_st2) + + end subroutine regrid_carbon + + subroutine compute_dayx ( & + NTILES, AGCM_YY, AGCM_MM, AGCM_DD, AGCM_HR, & + LATT, DAYX) + + implicit none + + integer, intent (in) :: NTILES,AGCM_YY,AGCM_MM,AGCM_DD,AGCM_HR + real, dimension (NTILES), intent (in) :: LATT + real, dimension (NTILES), intent (out) :: DAYX + integer, parameter :: DT = 900 + integer, parameter :: ncycle = 1461 ! number of days in a 4-year leap cycle (365*4 + 1) + real, dimension(ncycle) :: zc, zs + integer :: dofyr, sec,YEARS_PER_CYCLE, DAYS_PER_CYCLE, year, iday, idayp1, nn, n + real :: fac, YEARLEN, zsin, zcos, declin + + dofyr = AGCM_DD + if(AGCM_MM > 1) dofyr = dofyr + 31 + if(AGCM_MM > 2) then + dofyr = dofyr + 28 + if(mod(AGCM_YY,4) == 0) dofyr = dofyr + 1 + endif + if(AGCM_MM > 3) dofyr = dofyr + 31 + if(AGCM_MM > 4) dofyr = dofyr + 30 + if(AGCM_MM > 5) dofyr = dofyr + 31 + if(AGCM_MM > 6) dofyr = dofyr + 30 + if(AGCM_MM > 7) dofyr = dofyr + 31 + if(AGCM_MM > 8) dofyr = dofyr + 31 + if(AGCM_MM > 9) dofyr = dofyr + 30 + if(AGCM_MM > 10) dofyr = dofyr + 31 + if(AGCM_MM > 11) dofyr = dofyr + 30 + + sec = AGCM_HR * 3600 - DT ! subtract DT to get time of previous physics step + fac = real(sec) / 86400. + + + call orbit_create(zs,zc,ncycle) ! GEOS5 leap cycle routine + + YEARLEN = 365.25 + + ! Compute length of leap cycle + !------------------------------ + + if(YEARLEN-int(YEARLEN) > 0.) then + YEARS_PER_CYCLE = nint(1./(YEARLEN-int(YEARLEN))) + else + YEARS_PER_CYCLE = 1 + endif + + DAYS_PER_CYCLE=nint(YEARLEN*YEARS_PER_CYCLE) + + ! declination & daylength + ! ----------------------- + + YEAR = mod(AGCM_YY-1,YEARS_PER_CYCLE) + + IDAY = YEAR*int(YEARLEN)+dofyr + IDAYP1 = mod(IDAY,DAYS_PER_CYCLE) + 1 + + ZSin = ZS(IDAYP1)*FAC + ZS(IDAY)*(1.-FAC) ! sine of solar declination + ZCos = ZC(IDAYP1)*FAC + ZC(IDAY)*(1.-FAC) ! cosine of solar declination + + nn = 0 + do n = 1,days_per_cycle + nn = nn + 1 + if(nn > 365) nn = nn - 365 + ! print *, 'cycle:',n,nn,asin(ZS(n)) + end do + declin = asin(ZSin) + + ! compute daylength on input tile space (accounts for any change in physics time step) + ! do n = 1,ntiles_cn + ! fac = -(sin((latc(n)/zoom)*(MAPL_PI/180.))*zsin)/(cos((latc(n)/zoom)*(MAPL_PI/180.))*zcos) + ! fac = min(1.,max(-1.,fac)) + ! dayl(n) = (86400./MAPL_PI) * acos(fac) ! daylength (seconds) + ! end do + + ! compute daylength on output tile space (accounts for lat shift due to split & change in time step) + + do n = 1,ntiles + fac = -(sin(latt(n)*(MAPL_PI/180.))*zsin)/(cos(latt(n)*(MAPL_PI/180.))*zcos) + fac = min(1.,max(-1.,fac)) + dayx(n) = (86400./MAPL_PI) * acos(fac) ! daylength (seconds) + end do + + ! print *,'DAYX : ', minval(dayx),maxval(dayx), minval(latt), maxval(latt), zsin, zcos, dofyr, iday, idayp1, declin + + end subroutine compute_dayx + + ! ***************************************************************************** + + subroutine orbit_create(zs,zc,ncycle) + implicit none + + integer, intent(in) :: ncycle + real, intent(out), dimension(ncycle) :: zs, zc + + integer :: YEARS_PER_CYCLE, DAYS_PER_CYCLE + integer :: K, KP !, KM + real*8 :: T1, T2, T3, T4, FUN, Y, SOB, OMG, PRH, TT + real*8 :: YEARLEN + + ! STATEMENT FUNCTION + + FUN(Y) = OMG*(1.0-ECCENTRICITY*cos(Y-PRH))**2 + + YEARLEN = 365.25 + + ! Factors involving the orbital parameters + !------------------------------------------ + + OMG = (2.0*MAPL_PI/YEARLEN) / (sqrt(1.-ECCENTRICITY**2)**3) + PRH = PERIHELION*(MAPL_PI/180.) + SOB = sin(OBLIQUITY*(MAPL_PI/180.)) + + ! Compute length of leap cycle + !------------------------------ + + if(YEARLEN-int(YEARLEN) > 0.) then + YEARS_PER_CYCLE = nint(1./(YEARLEN-int(YEARLEN))) + else + YEARS_PER_CYCLE = 1 + endif + + + DAYS_PER_CYCLE=nint(YEARLEN*YEARS_PER_CYCLE) + + if(days_per_cycle /= ncycle) stop 'bad cycle' + + ! ZS: Sine of declination + ! ZC: Cosine of declination + + ! Begin integration at vernal equinox + + KP = EQUINOX + TT = 0.0 + ZS(KP) = sin(TT)*SOB + ZC(KP) = sqrt(1.0-ZS(KP)**2) + + ! Integrate orbit for entire leap cycle using Runge-Kutta + + do K=2,DAYS_PER_CYCLE + T1 = FUN(TT ) + T2 = FUN(TT+T1*0.5) + T3 = FUN(TT+T2*0.5) + T4 = FUN(TT+T3 ) + KP = mod(KP,DAYS_PER_CYCLE) + 1 + TT = TT + (T1 + 2.0*(T2 + T3) + T4) / 6.0 + ZS(KP) = sin(TT)*SOB + ZC(KP) = sqrt(1.0-ZS(KP)**2) + end do + end subroutine orbit_create + + end subroutine re_tile + +end module CatchmentCNRstMod diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 new file mode 100644 index 000000000..0ea043b54 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -0,0 +1,1283 @@ +#include "MAPL_Generic.h" + +module CatchmentRstMod + use mk_restarts_getidsMod + use MAPL + use mpi + use LSM_ROUTINES, ONLY: & + catch_calc_soil_moist, & + catch_calc_tp, & + catch_calc_ght + + USE CATCH_CONSTANTS, ONLY: & + N_GT => CATCH_N_GT, & + DZGT => CATCH_DZGT, & + PEATCLSM_POROS_THRESHOLD + + implicit none +#ifndef __GFORTRAN__ + integer :: ftell + external :: ftell +#endif + + type :: CatchmentRst + integer :: ntiles + type(FileMetadata) :: meta + character(len=:), allocatable :: time !yyyymmddhhmm + real, allocatable :: bf1(:) + real, allocatable :: bf2(:) + real, allocatable :: bf3(:) + real, allocatable :: vgwmax(:) + real, allocatable :: cdcr1(:) + real, allocatable :: cdcr2(:) + real, allocatable :: psis(:) + real, allocatable :: bee(:) + real, allocatable :: poros(:) + real, allocatable :: wpwet(:) + real, allocatable :: cond(:) + real, allocatable :: gnu(:) + real, allocatable :: ars1(:) + real, allocatable :: ars2(:) + real, allocatable :: ars3(:) + real, allocatable :: ara1(:) + real, allocatable :: ara2(:) + real, allocatable :: ara3(:) + real, allocatable :: ara4(:) + real, allocatable :: arw1(:) + real, allocatable :: arw2(:) + real, allocatable :: arw3(:) + real, allocatable :: arw4(:) + real, allocatable :: tsa1(:) + real, allocatable :: tsa2(:) + real, allocatable :: tsb1(:) + real, allocatable :: tsb2(:) + real, allocatable :: atau(:) + real, allocatable :: btau(:) + real, allocatable :: ity(:) + real, allocatable :: tc(:,:) + real, allocatable :: qc(:,:) + real, allocatable :: capac(:) + real, allocatable :: catdef(:) + real, allocatable :: rzexc(:) + real, allocatable :: srfexc(:) + real, allocatable :: ghtcnt1(:) + real, allocatable :: ghtcnt2(:) + real, allocatable :: ghtcnt3(:) + real, allocatable :: ghtcnt4(:) + real, allocatable :: ghtcnt5(:) + real, allocatable :: ghtcnt6(:) + real, allocatable :: tsurf(:) + real, allocatable :: wesnn1(:) + real, allocatable :: wesnn2(:) + real, allocatable :: wesnn3(:) + real, allocatable :: htsnnn1(:) + real, allocatable :: htsnnn2(:) + real, allocatable :: htsnnn3(:) + real, allocatable :: sndzn1(:) + real, allocatable :: sndzn2(:) + real, allocatable :: sndzn3(:) + real, allocatable :: ch(:,:) + real, allocatable :: cm(:,:) + real, allocatable :: cq(:,:) + real, allocatable :: fr(:,:) + real, allocatable :: ww(:,:) + ! save old values for scale + real, allocatable :: old_catdef(:) + real, allocatable :: old_cdcr1 (:) + real, allocatable :: old_cdcr2 (:) + real, allocatable :: old_rzexc (:) + real, allocatable :: old_vgwmax(:) + real, allocatable :: old_ghtcnt1(:) + real, allocatable :: old_ghtcnt2(:) + real, allocatable :: old_ghtcnt3(:) + real, allocatable :: old_ghtcnt4(:) + real, allocatable :: old_ghtcnt5(:) + real, allocatable :: old_ghtcnt6(:) + real, allocatable :: old_poros(:) + real, allocatable :: old_sndzn3(:) + ! intermediate + real, allocatable, dimension(:) :: lonc,latc,lonn,latt, latg + integer, allocatable, dimension(:) :: id_glb + contains + procedure :: read_GEOSldas_rst_bin + procedure :: write_nc4 + procedure :: read_shared_nc4 + procedure :: write_shared_nc4 + procedure :: add_bcs_to_rst + procedure :: allocate_catch + procedure :: re_tile + procedure :: re_scale + procedure :: set_scale_var + endtype CatchmentRst + + interface CatchmentRst + module procedure CatchmentRst_Create + module procedure CatchmentRst_empty + end interface + +contains + + function CatchmentRst_create(filename, time, rc) result (catch) + type(CatchmentRst) :: catch + character(*), intent(in) :: filename + character(*), intent(in) :: time ! yyyymmddhhmm format + integer, optional, intent(out) :: rc + integer :: status + character(len=256) :: Iam = "CatchmentRst_create" + type(Netcdf4_fileformatter) :: formatter + integer :: filetype, ntiles, unit + type(FileMetadata) :: meta + integer :: bpos, epos, n, myid, mpierr + + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + catch%time = time + call MAPL_NCIOGetFileType(filename, filetype, __RC__) + + if(filetype == 0) then + ! nc4 format + call formatter%open(filename, pFIO_READ, __RC__) + meta = formatter%read(__RC__) + ntiles = meta%get_dimension('tile', __RC__) + catch%meta = meta + catch%ntiles = ntiles + if (myid ==0) then + call catch%allocate_catch() + call catch%read_shared_nc4(formatter, __RC__) + call MAPL_VarRead(formatter,"OLD_ITY",catch%ity, __RC__) + endif + call formatter%close() + else + !if ( .not. present(time)) then + ! _ASSERT(.false., 'Please provide time for binary catch, format yyyymmddhhmm') + !endif + !GEOSldas binary + open(newunit=unit, file=filename, form='unformatted', action = 'read') + bpos=0 + read(unit) + epos = ftell(unit) ! ending position of file pointer + close(unit) + ntiles = (epos-bpos)/4-2 ! record size (in 4 byte words; + catch%ntiles = ntiles + catch%meta = create_meta(ntiles, time) + if (myid ==0) then + call catch%allocate_catch() + call catch%read_GEOSldas_rst_bin(filename, __RC__) + endif + endif + _RETURN(_SUCCESS) + end function CatchmentRst_Create + + function CatchmentRst_empty(meta, time, rc) result (catch) + type(CatchmentRst) :: catch + character(*), intent(in) :: time + type(FileMetadata), intent(in) :: meta + + integer, optional, intent(out) :: rc + integer :: status, myid, mpierr + character(len=256) :: Iam = "CatchmentRst_create" + type(Netcdf4_fileformatter) :: formatter + + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + ! nc4 format + catch%ntiles = meta%get_dimension('tile', __RC__) + catch%meta = meta + catch%time = time + if (myid ==0) then + call catch%allocate_catch() + endif + _RETURN(_SUCCESS) + end function CatchmentRst_empty + + subroutine read_GEOSldas_rst_bin(this, filename, rc) + class(CatchmentRst), intent(inout) :: this + character(*), intent(in) :: filename + integer, optional, intent(out) :: rc + integer :: unit + open(newunit=unit, file=filename, form='unformatted', action='read') + read(unit) this% bf1 + read(unit) this% bf2 + read(unit) this% bf3 + read(unit) this% vgwmax + read(unit) this% cdcr1 + read(unit) this% cdcr2 + read(unit) this% psis + read(unit) this% bee + read(unit) this% poros + read(unit) this% wpwet + read(unit) this% cond + read(unit) this% gnu + read(unit) this% ars1 + read(unit) this% ars2 + read(unit) this% ars3 + read(unit) this% ara1 + read(unit) this% ara2 + read(unit) this% ara3 + read(unit) this% ara4 + read(unit) this% arw1 + read(unit) this% arw2 + read(unit) this% arw3 + read(unit) this% arw4 + read(unit) this% tsa1 + read(unit) this% tsa2 + read(unit) this% tsb1 + read(unit) this% tsb2 + read(unit) this% atau + read(unit) this% btau + read(unit) this% ity + read(unit) this% tc + read(unit) this% qc + read(unit) this% capac + read(unit) this% catdef + read(unit) this% rzexc + read(unit) this% srfexc + read(unit) this% ghtcnt1 + read(unit) this% ghtcnt2 + read(unit) this% ghtcnt3 + read(unit) this% ghtcnt4 + read(unit) this% ghtcnt5 + read(unit) this% ghtcnt6 + read(unit) this% tsurf + read(unit) this% wesnn1 + read(unit) this% wesnn2 + read(unit) this% wesnn3 + read(unit) this% htsnnn1 + read(unit) this% htsnnn2 + read(unit) this% htsnnn3 + read(unit) this% sndzn1 + read(unit) this% sndzn2 + read(unit) this% sndzn3 + read(unit) this% ch + read(unit) this% cm + read(unit) this% cq + read(unit) this% fr + read(unit) this% ww + close(unit) + _RETURN(_SUCCESS) + end subroutine read_GEOSldas_rst_bin + + subroutine read_shared_nc4(this, formatter, rc) + class(CatchmentRst), intent(inout):: this + type(Netcdf4_fileformatter),intent(inout) :: formatter + integer, optional, intent(out):: rc + integer :: status + + call MAPL_VarRead(formatter,"BF1",this%bf1, __RC__) + call MAPL_VarRead(formatter,"BF2",this%bf2, __RC__) + call MAPL_VarRead(formatter,"BF3",this%bf3, __RC__) + call MAPL_VarRead(formatter,"VGWMAX",this%vgwmax, __RC__) + call MAPL_VarRead(formatter,"CDCR1",this%cdcr1, __RC__) + call MAPL_VarRead(formatter,"CDCR2",this%cdcr2, __RC__) + call MAPL_VarRead(formatter,"PSIS",this%psis, __RC__) + call MAPL_VarRead(formatter,"BEE",this%bee, __RC__) + call MAPL_VarRead(formatter,"POROS",this%poros, __RC__) + call MAPL_VarRead(formatter,"WPWET",this%wpwet, __RC__) + call MAPL_VarRead(formatter,"COND",this%cond, __RC__) + call MAPL_VarRead(formatter,"GNU",this%gnu, __RC__) + call MAPL_VarRead(formatter,"ARS1",this%ars1, __RC__) + call MAPL_VarRead(formatter,"ARS2",this%ars2, __RC__) + call MAPL_VarRead(formatter,"ARS3",this%ars3, __RC__) + call MAPL_VarRead(formatter,"ARA1",this%ara1, __RC__) + call MAPL_VarRead(formatter,"ARA2",this%ara2, __RC__) + call MAPL_VarRead(formatter,"ARA3",this%ara3, __RC__) + call MAPL_VarRead(formatter,"ARA4",this%ara4, __RC__) + call MAPL_VarRead(formatter,"ARW1",this%arw1, __RC__) + call MAPL_VarRead(formatter,"ARW2",this%arw2, __RC__) + call MAPL_VarRead(formatter,"ARW3",this%arw3, __RC__) + call MAPL_VarRead(formatter,"ARW4",this%arw4, __RC__) + call MAPL_VarRead(formatter,"TSA1",this%tsa1, __RC__) + call MAPL_VarRead(formatter,"TSA2",this%tsa2, __RC__) + call MAPL_VarRead(formatter,"TSB1",this%tsb1, __RC__) + call MAPL_VarRead(formatter,"TSB2",this%tsb2, __RC__) + call MAPL_VarRead(formatter,"ATAU",this%atau, __RC__) + call MAPL_VarRead(formatter,"BTAU",this%btau, __RC__) + call MAPL_VarRead(formatter,"TC",this%tc, __RC__) + call MAPL_VarRead(formatter,"QC",this%qc, __RC__) +! +! call MAPL_VarRead(formatter,"OLD_ITY",this%ity, __RC__) +! + call MAPL_VarRead(formatter,"CAPAC",this%capac, __RC__) + call MAPL_VarRead(formatter,"CATDEF",this%catdef, __RC__) + call MAPL_VarRead(formatter,"RZEXC",this%rzexc, __RC__) + call MAPL_VarRead(formatter,"SRFEXC",this%srfexc, __RC__) + call MAPL_VarRead(formatter,"GHTCNT1",this%ghtcnt1, __RC__) + call MAPL_VarRead(formatter,"GHTCNT2",this%ghtcnt2, __RC__) + call MAPL_VarRead(formatter,"GHTCNT3",this%ghtcnt3, __RC__) + call MAPL_VarRead(formatter,"GHTCNT4",this%ghtcnt4, __RC__) + call MAPL_VarRead(formatter,"GHTCNT5",this%ghtcnt5, __RC__) + call MAPL_VarRead(formatter,"GHTCNT6",this%ghtcnt6, __RC__) + if (this%meta%has_variable('TSURF')) then + call MAPL_VarRead(formatter,"TSURF",this%tsurf, __RC__) + endif + call MAPL_VarRead(formatter,"WESNN1",this%wesnn1, __RC__) + call MAPL_VarRead(formatter,"WESNN2",this%wesnn2, __RC__) + call MAPL_VarRead(formatter,"WESNN3",this%wesnn3, __RC__) + call MAPL_VarRead(formatter,"HTSNNN1",this%htsnnn1, __RC__) + call MAPL_VarRead(formatter,"HTSNNN2",this%htsnnn2, __RC__) + call MAPL_VarRead(formatter,"HTSNNN3",this%htsnnn3, __RC__) + call MAPL_VarRead(formatter,"SNDZN1",this%sndzn1, __RC__) + call MAPL_VarRead(formatter,"SNDZN2",this%sndzn2, __RC__) + call MAPL_VarRead(formatter,"SNDZN3",this%sndzn3, __RC__) + if (this%meta%has_variable('CH')) then + call MAPL_VarRead(formatter,"CH",this%ch, __RC__) + endif + if (this%meta%has_variable('CM')) then + call MAPL_VarRead(formatter,"CM",this%cm, __RC__) + endif + if (this%meta%has_variable('CQ')) then + call MAPL_VarRead(formatter,"CQ",this%cq, __RC__) + endif + if (this%meta%has_variable('FR')) then + call MAPL_VarRead(formatter,"FR",this%fr, __RC__) + endif + if (this%meta%has_variable('WW')) then + call MAPL_VarRead(formatter,"WW",this%ww, __RC__) + endif + _RETURN(_SUCCESS) + end subroutine + + subroutine write_nc4 (this, filename, rc) + class(CatchmentRst), intent(inout):: this + character(*), intent(in) :: filename + integer, optional, intent(out):: rc + + type(Netcdf4_fileformatter) :: formatter + integer :: status + character(256) :: Iam = "write_nc4" + + call formatter%create(filename, __RC__) + call formatter%write(this%meta, __RC__) + call this%write_shared_nc4(formatter, __RC__) + call MAPL_VarWrite(formatter,"OLD_ITY",this%ity) + call formatter%close() + _RETURN(_SUCCESS) + + end subroutine write_nc4 + + subroutine write_shared_nc4(this, formatter, rc) + class(CatchmentRst), intent(inout):: this + type(Netcdf4_fileformatter),intent(inout) :: formatter + integer, optional, intent(out):: rc + integer :: status + call MAPL_VarWrite(formatter,"BF1",this%bf1) + call MAPL_VarWrite(formatter,"BF2",this%bf2) + call MAPL_VarWrite(formatter,"BF3",this%bf3) + call MAPL_VarWrite(formatter,"VGWMAX",this%vgwmax) + call MAPL_VarWrite(formatter,"CDCR1",this%cdcr1) + call MAPL_VarWrite(formatter,"CDCR2",this%cdcr2) + call MAPL_VarWrite(formatter,"PSIS",this%psis) + call MAPL_VarWrite(formatter,"BEE",this%bee) + call MAPL_VarWrite(formatter,"POROS",this%poros) + call MAPL_VarWrite(formatter,"WPWET",this%wpwet) + call MAPL_VarWrite(formatter,"COND",this%cond) + call MAPL_VarWrite(formatter,"GNU",this%gnu) + call MAPL_VarWrite(formatter,"ARS1",this%ars1) + call MAPL_VarWrite(formatter,"ARS2",this%ars2) + call MAPL_VarWrite(formatter,"ARS3",this%ars3) + call MAPL_VarWrite(formatter,"ARA1",this%ara1) + call MAPL_VarWrite(formatter,"ARA2",this%ara2) + call MAPL_VarWrite(formatter,"ARA3",this%ara3) + call MAPL_VarWrite(formatter,"ARA4",this%ara4) + call MAPL_VarWrite(formatter,"ARW1",this%arw1) + call MAPL_VarWrite(formatter,"ARW2",this%arw2) + call MAPL_VarWrite(formatter,"ARW3",this%arw3) + call MAPL_VarWrite(formatter,"ARW4",this%arw4) + call MAPL_VarWrite(formatter,"TSA1",this%tsa1) + call MAPL_VarWrite(formatter,"TSA2",this%tsa2) + call MAPL_VarWrite(formatter,"TSB1",this%tsb1) + call MAPL_VarWrite(formatter,"TSB2",this%tsb2) + call MAPL_VarWrite(formatter,"ATAU",this%atau) + call MAPL_VarWrite(formatter,"BTAU",this%btau) + call MAPL_VarWrite(formatter,"TC",this%tc) + call MAPL_VarWrite(formatter,"QC",this%qc) + call MAPL_VarWrite(formatter,"CAPAC",this%capac) + call MAPL_VarWrite(formatter,"CATDEF",this%catdef) + call MAPL_VarWrite(formatter,"RZEXC",this%rzexc) + call MAPL_VarWrite(formatter,"SRFEXC",this%srfexc) + call MAPL_VarWrite(formatter,"GHTCNT1",this%ghtcnt1) + call MAPL_VarWrite(formatter,"GHTCNT2",this%ghtcnt2) + call MAPL_VarWrite(formatter,"GHTCNT3",this%ghtcnt3) + call MAPL_VarWrite(formatter,"GHTCNT4",this%ghtcnt4) + call MAPL_VarWrite(formatter,"GHTCNT5",this%ghtcnt5) + call MAPL_VarWrite(formatter,"GHTCNT6",this%ghtcnt6) + if (this%meta%has_variable('TSURF')) then + call MAPL_VarWrite(formatter,"TSURF",this%tsurf) + endif + call MAPL_VarWrite(formatter,"WESNN1",this%wesnn1) + call MAPL_VarWrite(formatter,"WESNN2",this%wesnn2) + call MAPL_VarWrite(formatter,"WESNN3",this%wesnn3) + call MAPL_VarWrite(formatter,"HTSNNN1",this%htsnnn1) + call MAPL_VarWrite(formatter,"HTSNNN2",this%htsnnn2) + call MAPL_VarWrite(formatter,"HTSNNN3",this%htsnnn3) + call MAPL_VarWrite(formatter,"SNDZN1",this%sndzn1) + call MAPL_VarWrite(formatter,"SNDZN2",this%sndzn2) + call MAPL_VarWrite(formatter,"SNDZN3",this%sndzn3) + + if (this%meta%has_variable('CH')) then + call MAPL_VarWrite(formatter,"CH",this%ch) + endif + if (this%meta%has_variable('CM')) then + call MAPL_VarWrite(formatter,"CM",this%cm) + endif + if (this%meta%has_variable('CQ')) then + call MAPL_VarWrite(formatter,"CQ",this%cq) + endif + if (this%meta%has_variable('FR')) then + call MAPL_VarWrite(formatter,"FR",this%fr) + endif + if (this%meta%has_variable('WW')) then + call MAPL_VarWrite(formatter,"WW",this%ww) + endif + + _RETURN(_SUCCESS) + + end subroutine write_shared_nc4 + + subroutine allocate_catch(this,rc) + class(CatchmentRst), intent(inout) :: this + integer, optional, intent(out):: rc + integer :: ntiles + ntiles = this%ntiles + allocate( this% bf1(ntiles) ) + allocate( this% bf2(ntiles) ) + allocate( this% bf3(ntiles) ) + allocate( this% vgwmax(ntiles) ) + allocate( this% cdcr1(ntiles) ) + allocate( this% cdcr2(ntiles) ) + allocate( this% psis(ntiles) ) + allocate( this% bee(ntiles) ) + allocate( this% poros(ntiles) ) + allocate( this% wpwet(ntiles) ) + allocate( this% cond(ntiles) ) + allocate( this% gnu(ntiles) ) + allocate( this% ars1(ntiles) ) + allocate( this% ars2(ntiles) ) + allocate( this% ars3(ntiles) ) + allocate( this% ara1(ntiles) ) + allocate( this% ara2(ntiles) ) + allocate( this% ara3(ntiles) ) + allocate( this% ara4(ntiles) ) + allocate( this% arw1(ntiles) ) + allocate( this% arw2(ntiles) ) + allocate( this% arw3(ntiles) ) + allocate( this% arw4(ntiles) ) + allocate( this% tsa1(ntiles) ) + allocate( this% tsa2(ntiles) ) + allocate( this% tsb1(ntiles) ) + allocate( this% tsb2(ntiles) ) + allocate( this% atau(ntiles) ) + allocate( this% btau(ntiles) ) + allocate( this% ity(ntiles) ) + allocate( this% tc(ntiles,4) ) + allocate( this% qc(ntiles,4) ) + allocate( this% capac(ntiles) ) + allocate( this% catdef(ntiles) ) + allocate( this% rzexc(ntiles) ) + allocate( this% srfexc(ntiles) ) + allocate( this% ghtcnt1(ntiles) ) + allocate( this% ghtcnt2(ntiles) ) + allocate( this% ghtcnt3(ntiles) ) + allocate( this% ghtcnt4(ntiles) ) + allocate( this% ghtcnt5(ntiles) ) + allocate( this% ghtcnt6(ntiles) ) + + if (this%meta%has_variable('TSURF')) then + allocate( this% tsurf(ntiles) ) + endif + + allocate( this% wesnn1(ntiles) ) + allocate( this% wesnn2(ntiles) ) + allocate( this% wesnn3(ntiles) ) + allocate( this% htsnnn1(ntiles) ) + allocate( this% htsnnn2(ntiles) ) + allocate( this% htsnnn3(ntiles) ) + allocate( this% sndzn1(ntiles) ) + allocate( this% sndzn2(ntiles) ) + allocate( this% sndzn3(ntiles) ) + + if (this%meta%has_variable('CH')) then + allocate( this% ch(ntiles,4) ) + endif + if (this%meta%has_variable('CM')) then + allocate( this% cm(ntiles,4) ) + endif + if (this%meta%has_variable('CQ')) then + allocate( this% cq(ntiles,4) ) + endif + if (this%meta%has_variable('FR')) then + allocate( this% fr(ntiles,4) ) + endif + if (this%meta%has_variable('WW')) then + allocate( this% ww(ntiles,4) ) + endif + _RETURN(_SUCCESS) + end subroutine allocate_catch + + ! This subroutine reads BCs from BCSDIR and hydrological varable + subroutine add_bcs_to_rst(this, surflay, DataDir, rc) + class(CatchmentRst), intent(inout) :: this + real, intent(in) :: surflay + character(*), intent(in) :: DataDir + integer, optional, intent(out) :: rc + + real, allocatable :: DP2BR(:) + real :: CanopH + integer, allocatable :: ity(:) + integer :: ntiles, STATUS + integer :: idum, i,j,n, ib, nv + real :: rdum, zdep1, zdep2, zdep3, zmet, term1, term2, bare,fvg(4) + logical :: NEWLAND + logical :: file_exists + + type(NetCDF4_Fileformatter) :: CatchFmt + + character*256 :: Iam = "add_bcs" + + ntiles = this%ntiles + + allocate (DP2BR(ntiles), ity(ntiles) ) + + inquire(file = trim(DataDir)//'/clsm/catch_params.nc4', exist=file_exists) + inquire(file = trim(DataDir)//"/clsm/CLM_veg_typs_fracs",exist=NewLand ) + + if (size(this%ara1) /= this%ntiles ) then + ! it is just re-allocate + this%ity = DP2BR + this%ARA1 = DP2BR + this%ARA2 = DP2BR + this%ARA3 = DP2BR + this%ARA4 = DP2BR + this%ARS1 = DP2BR + this%ARS2 = DP2BR + this%ARS3 = DP2BR + this%ARW1 = DP2BR + this%ARW2 = DP2BR + this%ARW3 = DP2BR + this%ARW4 = DP2BR + this%ATAU = DP2BR + this%BTAU = DP2BR + this%PSIS = DP2BR + this%BEE = DP2BR + this%BF1 = DP2BR + this%BF2 = DP2BR + this%BF3 = DP2BR + this%TSA1 = DP2BR + this%TSA2 = DP2BR + this%TSB1 = DP2BR + this%TSB2 = DP2BR + this%COND = DP2BR + this%WPWET = DP2BR + this%POROS = DP2BR + this%VGWMAX = DP2BR + this%cdcr1 = DP2BR + this%cdcr2 = DP2BR + endif + + if(file_exists) then + call CatchFmt%Open(trim(DataDir)//'/clsm/catch_params.nc4', pFIO_READ, __RC__) + call MAPL_VarRead ( CatchFmt ,'OLD_ITY', this%ITY, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARA1', this%ARA1, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARA2', this%ARA2, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARA3', this%ARA3, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARA4', this%ARA4, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARS1', this%ARS1, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARS2', this%ARS2, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARS3', this%ARS3, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARW1', this%ARW1, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARW2', this%ARW2, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARW3', this%ARW3, __RC__) + call MAPL_VarRead ( CatchFmt ,'ARW4', this%ARW4, __RC__) + + if( SURFLAY.eq.20.0 ) then + call MAPL_VarRead ( CatchFmt ,'ATAU2', this%ATAU, __RC__) + call MAPL_VarRead ( CatchFmt ,'BTAU2', this%BTAU, __RC__) + endif + + if( SURFLAY.eq.50.0 ) then + call MAPL_VarRead ( CatchFmt ,'ATAU5', this%ATAU, __RC__) + call MAPL_VarRead ( CatchFmt ,'BTAU5', this%BTAU, __RC__) + endif + + call MAPL_VarRead ( CatchFmt ,'PSIS', this%PSIS, __RC__) + call MAPL_VarRead ( CatchFmt ,'BEE', this%BEE, __RC__) + call MAPL_VarRead ( CatchFmt ,'BF1', this%BF1, __RC__) + call MAPL_VarRead ( CatchFmt ,'BF2', this%BF2, __RC__) + call MAPL_VarRead ( CatchFmt ,'BF3', this%BF3, __RC__) + call MAPL_VarRead ( CatchFmt ,'TSA1', this%TSA1, __RC__) + call MAPL_VarRead ( CatchFmt ,'TSA2', this%TSA2, __RC__) + call MAPL_VarRead ( CatchFmt ,'TSB1', this%TSB1, __RC__) + call MAPL_VarRead ( CatchFmt ,'TSB2', this%TSB2, __RC__) + call MAPL_VarRead ( CatchFmt ,'COND', this%COND, __RC__) + call MAPL_VarRead ( CatchFmt ,'GNU', this%GNU, __RC__) + call MAPL_VarRead ( CatchFmt ,'WPWET', this%WPWET, __RC__) + call MAPL_VarRead ( CatchFmt ,'DP2BR', DP2BR, __RC__) + call MAPL_VarRead ( CatchFmt ,'POROS', this%POROS, __RC__) + call CatchFmt%close() + else + open(unit=21, file=trim(DataDir)//'/clsm/mosaic_veg_typs_fracs',form='formatted') + open(unit=22, file=trim(DataDir)//'/clsm//bf.dat' ,form='formatted') + open(unit=23, file=trim(DataDir)//'/clsm/soil_param.dat' ,form='formatted') + open(unit=24, file=trim(DataDir)//'/clsm/ar.new' ,form='formatted') + open(unit=25, file=trim(DataDir)//'/clsm/ts.dat' ,form='formatted') + open(unit=26, file=trim(DataDir)//'/clsm/tau_param.dat' ,form='formatted') + + do n=1,ntiles + ! W.J notes: CanopH is not used. If CLM_veg_typs_fracs exists, the read some dummy ???? Ask Sarith + if (NewLand) then + read(21,*) I, j, ITY(N),idum, rdum, rdum, CanopH + else + read(21,*) I, j, ITY(N),idum, rdum, rdum + endif + + read (22, *) i,j, this%GNU(n), this%BF1(n), this%BF2(n), this%BF3(n) + + read (23, *) i,j, idum, idum, this%BEE(n), this%PSIS(n),& + this%POROS(n), this%COND(n), this%WPWET(n), DP2BR(n) + + read (24, *) i,j, rdum, this%ARS1(n), this%ARS2(n), this%ARS3(n), & + this%ARA1(n), this%ARA2(n), this%ARA3(n), this%ARA4(n), & + this%ARW1(n), this%ARW2(n), this%ARW3(n), this%ARW4(n) + + read (25, *) i,j, rdum, this%TSA1(n), this%TSA2(n), this%TSB1(n), this%TSB2(n) + + if( SURFLAY.eq.20.0 ) read (26, *) i,j, this%ATAU(n), this%BTAU(n), rdum, rdum ! for old soil params + if( SURFLAY.eq.50.0 ) read (26, *) i,j, rdum , rdum, this%ATAU(n), this%BTAU(n) ! for new soil params + + end do + this%ity = real(ity) + CLOSE (21, STATUS = 'KEEP') + CLOSE (22, STATUS = 'KEEP') + CLOSE (23, STATUS = 'KEEP') + CLOSE (24, STATUS = 'KEEP') + CLOSE (25, STATUS = 'KEEP') + CLOSE (26, STATUS = 'KEEP') + endif + + do n=1,ntiles + zdep2=1000. + zdep3=amax1(1000.,DP2BR(n)) + + if (zdep2 .gt.0.75*zdep3) then + zdep2 = 0.75*zdep3 + end if + + zdep1=20. + zmet=zdep3/1000. + + term1=-1.+((this%PSIS(n)-zmet)/this%PSIS(n))**((this%BEE(n)-1.)/this%BEE(n)) + term2=this%PSIS(n)*this%BEE(n)/(this%BEE(n)-1) + + this%VGWMAX(n) = this%POROS(n)*zdep2 + this%CDCR1(n) = 1000.*this%POROS(n)*(zmet-(-term2*term1)) + this%CDCR2(n) = (1.-this%WPWET(n))*this%POROS(n)*zdep3 + enddo + + _RETURN(_SUCCESS) + end subroutine add_bcs_to_rst + + type(FileMetadata) function create_meta(ntiles, t, rc) result(meta) + integer, intent(in) :: ntiles + character(*), intent(in) :: t ! yyyymmddmmhh format + integer, optional, intent(out) :: rc + character(64), dimension(58, 3) :: fields + integer :: n, status + character(:), allocatable :: s + type(Variable) :: var + + fields(1,:) = [character(len=64)::"ARA1" , "shape_param_1" , "m+2 kg-1"] + fields(2,:) = [character(len=64)::"ARA2" , "shape_param_2" , "1"] + fields(3,:) = [character(len=64)::"ARA3" , "shape_param_3" , "m+2 kg-1"] + fields(4,:) = [character(len=64)::"ARA4" , "shape_param_4" , "1"] + fields(5,:) = [character(len=64)::"ARS1" , "wetness_param_1" , "m+2 kg-1"] + fields(6,:) = [character(len=64)::"ARS2" , "wetness_param_2" , "m+2 kg-1"] + fields(7,:) = [character(len=64)::"ARS3" , "wetness_param_3" , "m+4 kg-2"] + fields(8,:) = [character(len=64)::"ARW1" , "min_theta_param_1" , "m+2 kg-1"] + fields(9,:) = [character(len=64)::"ARW2" , "min_theta_param_2" , "m+2 kg-1"] + fields(10,:) = [character(len=64)::"ARW3" , "min_theta_param_3" , "m+4 kg-2"] + fields(11,:) = [character(len=64)::"ARW4" , "min_theta_param_4" , "1"] + fields(12,:) = [character(len=64)::"ATAU" , "water_transfer_param_5" , "1"] + fields(13,:) = [character(len=64)::"BEE" , "clapp_hornberger_b" , "1"] + fields(14,:) = [character(len=64)::"BF1" , "topo_baseflow_param_1" , "kg m-4"] + fields(15,:) = [character(len=64)::"BF2" , "topo_baseflow_param_2" , "m"] + fields(16,:) = [character(len=64)::"BF3" , "topo_baseflow_param_3" , "log(m)"] + fields(17,:) = [character(len=64)::"BTAU" , "water_transfer_param_6" , "1"] + fields(18,:) = [character(len=64)::"CAPAC" , "interception_reservoir_capac" , "kg m-2"] + fields(19,:) = [character(len=64)::"CATDEF" , "catchment_deficit" , "kg m-2"] + fields(20,:) = [character(len=64)::"CDCR1" , "moisture_threshold" , "kg m-2"] + fields(21,:) = [character(len=64)::"CDCR2" , "max_water_content" , "kg m-2"] + fields(22,:) = [character(len=64)::"COND" , "sfc_sat_hydraulic_conduct" , "m s-1"] + fields(23,:) = [character(len=64)::"GHTCNT1" , "soil_heat_content_layer_1" , "J m-2"] + fields(24,:) = [character(len=64)::"GHTCNT2" , "soil_heat_content_layer_2" , "J_m-2"] + fields(25,:) = [character(len=64)::"GHTCNT3" , "soil_heat_content_layer_3" , "J m-2"] + fields(26,:) = [character(len=64)::"GHTCNT4" , "soil_heat_content_layer_4" , "J m-2"] + fields(27,:) = [character(len=64)::"GHTCNT5" , "soil_heat_content_layer_5" , "J m-2"] + fields(28,:) = [character(len=64)::"GHTCNT6" , "soil_heat_content_layer_6" , "J m-2"] + fields(29,:) = [character(len=64)::"GNU" , "vertical_transmissivity" , "m-1"] + fields(30,:) = [character(len=64)::"HTSNNN1" , "heat_content_snow_layer_1" , "J m-2"] + fields(31,:) = [character(len=64)::"HTSNNN2" , "heat_content_snow_layer_2" , "J m-2"] + fields(32,:) = [character(len=64)::"HTSNNN3" , "heat_content_snow_layer_3" , "J m-2"] + fields(33,:) = [character(len=64)::"OLD_ITY" , "Placeholder. Used to be vegetation_type." , "1"] + fields(34,:) = [character(len=64)::"POROS" , "soil_porosity" , "1"] + fields(35,:) = [character(len=64)::"PSIS" , "saturated_matric_potential" , "m"] + fields(36,:) = [character(len=64)::"RZEXC" , "root_zone_excess" , "kg m-2"] + fields(37,:) = [character(len=64)::"SNDZN1" , "snow_depth_layer_1" , "m"] + fields(38,:) = [character(len=64)::"SNDZN2" , "snow_depth_layer_2" , "m"] + fields(39,:) = [character(len=64)::"SNDZN3" , "snow_depth_layer_3" , "m"] + fields(40,:) = [character(len=64)::"SRFEXC" , "surface_excess" , "kg m-2"] + fields(41,:) = [character(len=64)::"TILE_ID" , "catchment_tile_id" , "1"] + fields(42,:) = [character(len=64)::"TSA1" , "water_transfer_param_1" , "1"] + fields(43,:) = [character(len=64)::"TSA2" , "water_transfer_param_2" , "1"] + fields(44,:) = [character(len=64)::"TSB1" , "water_transfer_param_3" , "1"] + fields(45,:) = [character(len=64)::"TSB2" , "water_transfer_param_4" , "1"] + fields(46,:) = [character(len=64)::"TSURF" , "mean_catchment_temp_incl_snw" , "K"] + fields(47,:) = [character(len=64)::"VGWMAX" , "max_rootzone_water_content" , "kg m-2"] + fields(48,:) = [character(len=64)::"WESNN1" , "snow_mass_layer_1" , "kg m-2"] + fields(49,:) = [character(len=64)::"WESNN2" , "snow_mass_layer_2" , "kg m-2"] + fields(50,:) = [character(len=64)::"WESNN3" , "snow_mass_layer_3" , "kg m-2"] + fields(51,:) = [character(len=64)::"WPWET" , "wetness_at_wilting_point" , "1"] + fields(52,:) = [character(len=64)::"CH" , "surface_heat_exchange_coefficient" , "kg m-2 s-1"] + fields(53,:) = [character(len=64)::"CM" , "surface_momentum_exchange_coefficient" , "kg m-2 s-1"] + fields(54,:) = [character(len=64)::"CQ" , "surface_moisture_exchange_coffiecient" , "kg m-2 s-1"] + fields(55,:) = [character(len=64)::"FR" , "subtile_fractions" , "1"] + fields(56,:) = [character(len=64)::"QC" , "canopy_specific_humidity" , "kg kg-1"] + fields(57,:) = [character(len=64)::"TC" , "canopy_temperature" , "K"] + fields(58,:) = [character(len=64)::"WW" , "vertical_velocity_scale_squared" , "m+2 s-2"] + + call meta%add_dimension('tile', ntiles) + call meta%add_dimension('subtile', 4) + call meta%add_dimension('time',1) + + do n = 1, 58 + if (n >=52) then + var = Variable(type=pFIO_REAL32, dimensions='tile,subtile') + else + var = Variable(type=pFIO_REAL32, dimensions='tile') + endif + call var%add_attribute('long_name', trim(fields(n,2))) + call var%add_attribute('units', trim(fields(n,3))) + call meta%add_variable(trim(fields(n,1)), var) + enddo + var = Variable(type=pFIO_REAL32, dimensions='time') + s = "minutes since "//t(1:4)//"-"//t(5:6)//"-"//t(7:8)//" "//t(9:10)//":00:00" + call var%add_attribute('units', s) + call meta%add_variable('time', var) + _RETURN(_SUCCESS) + end function create_meta + + subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) + class(CatchmentRst), intent(inout) :: this + character(*), intent(in) :: InTileFile + character(*), intent(in) :: OutBcsDir + character(*), intent(in) :: OutTileFile + real, intent(in) :: surflay + integer, optional, intent(out) :: rc + integer :: status, in_ntiles, out_ntiles, myid, numprocs + real, allocatable :: var_out(:), tmp2d(:,:) + real , allocatable , dimension (:) :: lonn,latt + real , pointer, dimension (:) :: long, latg, lonc, latc + integer, allocatable , dimension (:) :: low_ind, upp_ind, nt_local + integer, allocatable , dimension (:) :: Id_glb, id_loc, tid_offl + logical :: root_proc + integer :: mpierr, n, i, k, req + type(CatchmentRst) :: xgrid + type(FileMetadata) :: meta + character(*), parameter :: Iam = "Catchment::Re_tile" + + open (10,file =trim(OutBcsDir)//"/clsm/catchment.def",status='old',form='formatted') + read (10,*) out_ntiles + close (10, status = 'keep') + + in_ntiles = this%ntiles + + call this%meta%modify_dimension('tile', out_ntiles, __RC__) + + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, mpierr ) + root_proc = .false. + if (myid == 0) root_proc = .true. + + if(root_proc) then + print *,'ntiles in target BCs : ',out_ntiles + print *,'ntiles in restarts : ',in_ntiles + endif + + + ! Domain decomposition + ! -------------------- + + + allocate(low_ind ( numprocs)) + allocate(upp_ind ( numprocs)) + allocate(nt_local( numprocs)) + low_ind (:) = 1 + upp_ind (:) = out_ntiles + nt_local(:) = out_ntiles + + if (numprocs > 1) then + do i = 1, numprocs - 1 + upp_ind(i) = low_ind(i) + (out_ntiles/numprocs) - 1 + low_ind(i+1) = upp_ind(i) + 1 + nt_local(i) = upp_ind(i) - low_ind(i) + 1 + end do + nt_local(numprocs) = upp_ind(numprocs) - low_ind(numprocs) + 1 + endif + + allocate (id_loc (nt_local (myid + 1))) + allocate (lonn (nt_local (myid + 1))) + allocate (latt (nt_local (myid + 1))) + allocate (lonc (1:in_ntiles)) + allocate (latc (1:in_ntiles)) + allocate (tid_offl (in_ntiles)) + + if (root_proc) then + allocate (long (out_ntiles)) + allocate (latg (out_ntiles)) + call ReadTileFile_RealLatLon ( OutTileFile, n, long, latg) + _ASSERT( n == out_ntiles, "Out tile number should match") + this%latg = latg + call ReadTileFile_RealLatLon ( InTileFile, n, lonc, latc) + _ASSERT( n == in_ntiles, "In tile number should match") + endif + + do i = 1, in_ntiles + tid_offl(i) = i + end do + + ! create mapping, nearest + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + + do i = 1, numprocs + if((I == 1).and.(myid == 0)) then + lonn(:) = long(low_ind(i) : upp_ind(i)) + latt(:) = latg(low_ind(i) : upp_ind(i)) + else if (I > 1) then + if(I-1 == myid) then + ! receiving from root + call MPI_RECV(lonn,nt_local(i) , MPI_REAL, 0,995,MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + call MPI_RECV(latt,nt_local(i) , MPI_REAL, 0,994,MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + else if (myid == 0) then + ! root sends + call MPI_ISend(long(low_ind(i) : upp_ind(i)),nt_local(i),MPI_REAL,i-1,995,MPI_COMM_WORLD,req,mpierr) + call MPI_WAIT (req,MPI_STATUS_IGNORE,mpierr) + call MPI_ISend(latg(low_ind(i) : upp_ind(i)),nt_local(i),MPI_REAL,i-1,994,MPI_COMM_WORLD,req,mpierr) + call MPI_WAIT (req,MPI_STATUS_IGNORE,mpierr) + endif + endif + enddo + if(root_proc) deallocate (long) + + call MPI_BCAST(lonc,in_ntiles,MPI_REAL,0,MPI_COMM_WORLD,mpierr) + call MPI_BCAST(latc,in_ntiles,MPI_REAL,0,MPI_COMM_WORLD,mpierr) + + ! -------------------------------------------------------------------------------- + ! Here we create transfer index array to map offline restarts to output tile space + ! -------------------------------------------------------------------------------- + + ! id_glb for hydrologic variable + this%lonc = lonc + this%latc = latc + this%lonn = lonn + this%latt = latt + + call GetIds(lonc, latc, lonn, latt, id_loc, tid_offl) + if(root_proc) allocate (id_glb (out_ntiles)) + + call MPI_Barrier(MPI_COMM_WORLD, STATUS) + + do i = 1, numprocs + if((I == 1).and.(myid == 0)) then + id_glb(low_ind(i) : upp_ind(i)) = Id_loc(:) + else if (I > 1) then + if(I-1 == myid) then + ! send to root + call MPI_ISend(id_loc,nt_local(i),MPI_INTEGER,0,993,MPI_COMM_WORLD,req,mpierr) + call MPI_WAIT (req,MPI_STATUS_IGNORE,mpierr) + else if (myid == 0) then + ! root receives + call MPI_RECV(id_glb(low_ind(i) : upp_ind(i)),nt_local(i) , MPI_INTEGER, i-1,993,MPI_COMM_WORLD,MPI_STATUS_IGNORE,mpierr) + endif + endif + end do + + deallocate (id_loc) + + this%ntiles = out_ntiles + if (root_proc) then + ! regrid + this%id_glb = id_glb + var_out = this%poros(id_glb(:)) + this%poros = var_out + + var_out = this%cond(id_glb(:)) + this%cond = var_out + + var_out = this%psis(id_glb(:)) + this%psis = var_out + + var_out = this%bee(id_glb(:)) + this%bee = var_out + + var_out = this%wpwet(id_glb(:)) + this%wpwet = var_out + + var_out = this%gnu(id_glb(:)) + this%gnu = var_out + + var_out = this%vgwmax(id_glb(:)) + this%vgwmax = var_out + + var_out = this%bf1(id_glb(:)) + this%bf1 = var_out + + var_out = this%bf2(id_glb(:)) + this%bf2 = var_out + + var_out = this%bf3(id_glb(:)) + this%bf3 = var_out + + var_out = this%cdcr1(id_glb(:)) + this%cdcr1 = var_out + + var_out = this%cdcr2(id_glb(:)) + this%cdcr2 = var_out + + var_out = this%ars1(id_glb(:)) + this%ars1 = var_out + + var_out = this%ars2(id_glb(:)) + this%ars2 = var_out + + var_out = this%ars3(id_glb(:)) + this%ars3 = var_out + + var_out = this%ara1(id_glb(:)) + this%ara1 = var_out + + var_out = this%ara2(id_glb(:)) + this%ara2 = var_out + + var_out = this%ara3(id_glb(:)) + this%ara3 = var_out + + var_out = this%ara4(id_glb(:)) + this%ara4 = var_out + + var_out = this%arw1(id_glb(:)) + this%arw1 = var_out + + var_out = this%arw2(id_glb(:)) + this%arw2 = var_out + + var_out = this%arw3(id_glb(:)) + this%arw3 = var_out + + var_out = this%arw4(id_glb(:)) + this%arw4 = var_out + + var_out = this%tsa1(id_glb(:)) + this%tsa1 = var_out + + var_out = this%tsa2(id_glb(:)) + this%tsa2 = var_out + + var_out = this%tsb1(id_glb(:)) + this%tsb1 = var_out + + var_out = this%tsb2(id_glb(:)) + this%tsb2 = var_out + + var_out = this%atau(id_glb(:)) + this%atau = var_out + + var_out = this%btau(id_glb(:)) + this%btau = var_out + + this%ity = [(k*1., k=1, out_ntiles)] + + tmp2d = this%tc + deallocate(this%tc) + allocate(this%tc(out_ntiles, 4)) + do k = 1, 4 + this%tc(:,k) = tmp2d(id_glb(:),k) + enddo + + tmp2d = this%qc + deallocate(this%qc) + allocate(this%qc(out_ntiles, 4)) + do k = 1, 4 + this%qc(:,k) = tmp2d(id_glb(:),k) + enddo + + var_out = this%capac(id_glb(:)) + this%capac = var_out + + var_out = this%catdef(id_glb(:)) + this%catdef = var_out + + var_out = this%rzexc(id_glb(:)) + this%rzexc = var_out + + var_out = this%SRFEXC(id_glb(:)) + this%SRFEXC = var_out + + var_out = this%GHTCNT1(id_glb(:)) + this%GHTCNT1 = var_out + + var_out = this%GHTCNT2(id_glb(:)) + this%GHTCNT2 = var_out + + var_out = this%GHTCNT3(id_glb(:)) + this%GHTCNT3 = var_out + + var_out = this%GHTCNT4(id_glb(:)) + this%GHTCNT4 = var_out + + var_out = this%GHTCNT5(id_glb(:)) + this%GHTCNT5 = var_out + + var_out = this%GHTCNT6(id_glb(:)) + this%GHTCNT6 = var_out + + var_out = this%WESNN1(id_glb(:)) + this%WESNN1 = var_out + + var_out = this%WESNN2(id_glb(:)) + this%WESNN2 = var_out + + var_out = this%WESNN3(id_glb(:)) + this%WESNN3 = var_out + + var_out = this%HTSNNN1(id_glb(:)) + this%HTSNNN1 = var_out + + var_out = this%HTSNNN2(id_glb(:)) + this%HTSNNN2 = var_out + + var_out = this%HTSNNN3(id_glb(:)) + this%HTSNNN3 = var_out + + var_out = this%SNDZN1(id_glb(:)) + this%SNDZN1 = var_out + + var_out = this%SNDZN2(id_glb(:)) + this%SNDZN2 = var_out + + var_out = this%SNDZN3(id_glb(:)) + this%SNDZN3 = var_out + + !set tsurf to zero + if (this%meta%has_variable('TSURF')) then + var_out = this%tsurf(id_glb(:)) + this%tsurf = var_out + endif + + ! CH CM CQ FR WW + ! WW + if(allocated(tmp2d)) deallocate(tmp2d) + allocate(tmp2d(out_ntiles,4)) + !tmp2d = 0.001 + if (this%meta%has_variable('CH')) then + do k = 1,4 + tmp2d(:,k) = this%ch(id_glb(:),k) + enddo + this%ch = tmp2d + endif + if (this%meta%has_variable('CM')) then + do k = 1,4 + tmp2d(:,k) = this%cm(id_glb(:),k) + enddo + this%cm = tmp2d + endif + if (this%meta%has_variable('CQ')) then + do k = 1,4 + tmp2d(:,k) = this%cq(id_glb(:),k) + enddo + this%cq = tmp2d + endif + !tmp2d = 0.25 + if (this%meta%has_variable('FR')) then + do k = 1,4 + tmp2d(:,k) = this%fr(id_glb(:),k) + enddo + this%fr = tmp2d + endif + !tmp2d = 0.1 + if (this%meta%has_variable('WW')) then + do k = 1,4 + tmp2d(:,k) = this%ww(id_glb(:),k) + enddo + this%ww = tmp2d + endif + + call this%set_scale_var() + endif + + if(associated(long)) deallocate(long) + if(associated(latg)) deallocate(latg) + if(associated(lonc)) deallocate(lonc) + if(associated(latc)) deallocate(latc) + + _RETURN(_SUCCESS) + end subroutine re_tile + + subroutine re_scale(this, surflay, wemin_in, wemin_out, rc) + class(CatchmentRst), intent(inout) :: this + real, intent(in) :: surflay + real, intent(in) :: wemin_in + real, intent(in) :: wemin_out + integer, optional, intent(out) :: rc + integer :: ntiles + + real, allocatable, dimension(:) :: dzsf, ar1, ar2, ar4 + real, allocatable, dimension(:,:) :: TP_IN, GHT_IN, FICE, TP_OUT + real, allocatable, dimension(:) :: swe_in, depth_in, areasc_in, areasc_out, depth_out + + type(Netcdf4_fileformatter) :: formatter + integer :: i, filetype, n + integer :: status + + ntiles = this%ntiles + + n =count((this%old_catdef .gt. this%old_cdcr1)) + + print*, "Scale tile and pesentile :", n,100*n/ntiles + +! Scale rxexc regardless of CDCR1, CDCR2 differences +! -------------------------------------------------- + this%rzexc = this%old_rzexc * ( this%vgwmax / this%old_vgwmax ) + +! Scale catdef regardless of whether CDCR2 is larger or smaller in the new situation +! ---------------------------------------------------------------------------------- + where (this%old_catdef .gt. this%old_cdcr1) + + this%catdef = this%cdcr1 + & + ( this%old_catdef - this%old_cdcr1 ) / & + ( this%old_cdcr2 - this%old_cdcr1 ) * & + ( this%cdcr2 - this%cdcr1) + end where + +! Scale catdef also for the case where catdef le cdcr1. +! ----------------------------------------------------- + where( this%old_catdef .le. this%old_cdcr1) + this%catdef = this%old_catdef * (this%cdcr1 / this%old_cdcr1) + end where + +! Sanity Check (catch_calc_soil_moist() forces consistency betw. srfexc, rzexc, catdef) +! ------------ + print *, 'Performing Sanity Check ...' + + allocate ( dzsf(ntiles), source = SURFLAY ) + allocate ( ar1( ntiles) ) + allocate ( ar2( ntiles) ) + allocate ( ar4( ntiles) ) + + call catch_calc_soil_moist( ntiles, dzsf, & + this%vgwmax, this%cdcr1, this%cdcr2, & + this%psis, this%bee, this%poros, this%wpwet, & + this%ars1, this%ars2, this%ars3, & + this%ara1, this%ara2, this%ara3, this%ara4, & + this%arw1, this%arw2, this%arw3, this%arw4, & + this%bf1, this%bf2, & + this%srfexc, this%rzexc, this%catdef, & + ar1, ar2, ar4 ) + + + allocate (TP_IN (N_GT, Ntiles)) + allocate (GHT_IN (N_GT, Ntiles)) + allocate (FICE (N_GT, NTILES)) + allocate (TP_OUT (N_GT, Ntiles)) + + GHT_IN (1,:) = this%old_ghtcnt1 + GHT_IN (2,:) = this%old_ghtcnt2 + GHT_IN (3,:) = this%old_ghtcnt3 + GHT_IN (4,:) = this%old_ghtcnt4 + GHT_IN (5,:) = this%old_ghtcnt5 + GHT_IN (6,:) = this%old_ghtcnt6 + + call catch_calc_tp ( NTILES, this%old_poros, GHT_IN, tp_in, FICE) + + do n = 1, ntiles + do i = 1, N_GT + call catch_calc_ght(dzgt(i), this%poros(n), tp_in(i,n), fice(i,n), GHT_IN(i,n)) + end do + end do + + this%ghtcnt1 = GHT_IN (1,:) + this%ghtcnt2 = GHT_IN (2,:) + this%ghtcnt3 = GHT_IN (3,:) + this%ghtcnt4 = GHT_IN (4,:) + this%ghtcnt5 = GHT_IN (5,:) + this%ghtcnt6 = GHT_IN (6,:) + +! Deep soil temp sanity check +! --------------------------- + + call catch_calc_tp ( NTILES, this%poros, GHT_IN, tp_out, FICE) + + print *, 'Percent tiles TP Layer 1 differ : ', 100.* count(ABS(tp_out(1,:) - tp_in(1,:)) > 1.e-5) /float (Ntiles) + print *, 'Percent tiles TP Layer 2 differ : ', 100.* count(ABS(tp_out(2,:) - tp_in(2,:)) > 1.e-5) /float (Ntiles) + print *, 'Percent tiles TP Layer 3 differ : ', 100.* count(ABS(tp_out(3,:) - tp_in(3,:)) > 1.e-5) /float (Ntiles) + print *, 'Percent tiles TP Layer 4 differ : ', 100.* count(ABS(tp_out(4,:) - tp_in(4,:)) > 1.e-5) /float (Ntiles) + print *, 'Percent tiles TP Layer 5 differ : ', 100.* count(ABS(tp_out(5,:) - tp_in(5,:)) > 1.e-5) /float (Ntiles) + print *, 'Percent tiles TP Layer 6 differ : ', 100.* count(ABS(tp_out(6,:) - tp_in(6,:)) > 1.e-5) /float (Ntiles) + +! SNOW scaling +! ------------ + + if(abs(wemin_out-wemin_in) >1.e-6 ) then + + allocate (swe_in (Ntiles)) + allocate (depth_in (Ntiles)) + allocate (depth_out (Ntiles)) + allocate (areasc_in (Ntiles)) + allocate (areasc_out (Ntiles)) + swe_in = this%wesnn1 + this%wesnn2 + this%wesnn3 + depth_in = this%sndzn1 + this%sndzn2 + this%sndzn3 + areasc_in = min(swe_in/wemin_in, 1.) + areasc_out= min(swe_in/wemin_out,1.) + + where (swe_in .gt. 0.) + where (areasc_in .lt. 1. .or. areasc_out .lt. 1.) + ! density_in= swe_in/(areasc_in * depth_in + 1.e-20) + ! depth_out = swe_in/(areasc_out*density_in) + depth_out = areasc_in * depth_in/(areasc_out + 1.e-20) + this%sndzn1 = depth_out/3. + this%sndzn2 = depth_out/3. + this%sndzn3 = depth_out/3. + endwhere + endwhere + + print *, 'Snow scaling summary' + print *, '....................' + print *, 'Percent tiles SNDZ scaled : ', 100.* count (this%sndzn3 .ne. this%old_sndzn3) /float (count (this%sndzn3 > 0.)) + + endif + + ! PEATCLSM - ensure low CATDEF on peat tiles where "old" restart is not also peat + ! ------------------------------------------------------------------------------- + + where ( (this%old_poros < PEATCLSM_POROS_THRESHOLD) .and. (this%poros >= PEATCLSM_POROS_THRESHOLD) ) + this%catdef = 25. + this%rzexc = 0. + this%srfexc = 0. + end where + + end subroutine re_scale + + subroutine set_scale_var(this ) + class(CatchmentRst), intent(inout) :: this + this%old_catdef = this%catdef + this%old_poros = this%poros + this%old_cdcr1 = this%cdcr1 + this%old_cdcr2 = this%cdcr2 + this%old_rzexc = this%rzexc + this%old_vgwmax = this%vgwmax + this%old_sndzn3 = this%sndzn3 + this%old_ghtcnt1 = this%ghtcnt1 + this%old_ghtcnt2 = this%ghtcnt2 + this%old_ghtcnt3 = this%ghtcnt3 + this%old_ghtcnt4 = this%ghtcnt4 + this%old_ghtcnt5 = this%ghtcnt5 + this%old_ghtcnt6 = this%ghtcnt6 + end subroutine set_scale_var + +end module CatchmentRstMod diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/README b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/README index 325e28cee..35eeaaa4f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/README +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/README @@ -1,28 +1,6 @@ This directory contains software to create restarts and bcs for all tile components, not just catch. -The five restarts are created in ./OutData/ by simply -running mk_Restarts (no arguments) from the command line. -This builds the two executables, if necessary. One executable -does all the land grid files (two restarts, an lai-grn file, and -two albedo files). The other executable can be used, with appropriate -arguments (see mk_Restarts) to create the othe three restarts. - -Before running mk_Restarts you must checkout the InData and OutData -subdirectories, and place the input data there. -(see READMEs in these two directories). - -For convenience we have place copies on the MAPL hash software in -this directory, so that the restarts can be created without building -MAPL. - -The directory also contains legacy software by M Kistler that may -still be needed for Eros and for early Fortuna tags. - -On 11/23/2009 Max made changes to sense if the input restart is old (i.e., has the 2 -pairs of prev/next variables or not. By default it produces an new -restart, compatible with Fortuna-2_0 and later. By adding a 4th, optional, -argument of "OutIsOld", it will produce a valid old restart. Note that -any other string as a 4th argument, or no 4th argument will produce -the new restart. +mk_Restart script is called by regrid.pl in GMAO_Shared. Please refer regrid.pl +for the parameters that are passed to mk_Restarts diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_Restarts b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_Restarts index 238396bf1..0dde27468 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_Restarts +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_Restarts @@ -17,8 +17,8 @@ my ($surflay, $rsttime, $grpID, $numtasks, $walltime, $rescale, $qos, $partition my ($mk_catch_j, $mk_catch_log, $weminIN, $weminOUT, $weminDFLT); my ($zoom); -# mk_catchcn job and log file names -#---------------------------------- +# mk_catch job and log file names (also applies to catchcn) +#---------------------------------------------------------- $mk_catch_j = "mk_catch.j"; $mk_catch_log = "mk_catch.log"; @@ -376,11 +376,11 @@ option flags -wemin weminIN minimum snow water equivalent threshold for input catch/cn [$weminDFLT] -wemout weminOUT minimum snow water equivalent threshold for output catch/cn [$weminDFLT] -route create the route internal restart - - -surflay n number of surface layers (catch & catchcn) + -surflay n thickness [mm] of surface soil moisture layer (catch & catchcn) + Ganymed-3 and earlier: SURFLAY=20 + Ganymed-4 and later : SURFLAY=50 -rsttime n10 restart time in format, yyyymmddhh (catchcn) or yyyymm (route) -grpID grpID group ID for batch submittal (catchcn) - -ntasks nt number of tasks to assign to catchcn batch job [112] -walltime wt walltime in format \"hh:mm:ss\" for catchcn batch job [1:00:00] -rescale diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catchANDcnRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catchANDcnRestarts.F90 new file mode 100644 index 000000000..931f97ffa --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catchANDcnRestarts.F90 @@ -0,0 +1,137 @@ +#define I_AM_MAIN +#include "MAPL_Generic.h" + +PROGRAM mk_catchANDcnRestarts + + use mpi + use MAPL + use CatchmentRstMod + use CatchmentCNRstMod + + implicit none + + character(len=:), allocatable :: out_bcsdir, out_dir, in_tilefile, out_tilefile, YYYYMMDDHHMM + character(len=:), allocatable :: model, in_rstfile, out_rstfile + character(len=:), allocatable :: out_File + real :: surflay, wemin_in, wemin_out + integer :: rc, status + integer :: myid, numprocs, mpierr + class (CatchmentRst), allocatable :: catch + + call MPI_INIT(mpierr) + call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) + call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, mpierr ) + + call process_cmd() + + if (index(model, 'catchcn') /=0 ) then + catch = CatchmentCNRst(in_rstfile, model, yyyymmddhhmm, __RC__) + else + catch = CatchmentRst(in_rstfile, yyyymmddhhmm, __RC__) + endif + + call catch%re_tile(in_tilefile, out_bcsdir, out_tilefile, surflay, __RC__) + + if (myid == 0) then + call catch%add_bcs_to_rst(surflay, out_bcsdir, rc) + call catch%re_scale(surflay, wemin_in, wemin_out, __RC__) + call catch%write_nc4(out_file, __RC__) + endif + + call MPI_FINALIZE(mpierr) + + contains + ! process commands + subroutine process_cmd() + integer :: nxt + character(len=256) :: arg + nxt = 1 + call getarg(nxt,arg) + + do while(trim(arg) /= '') + select case (trim(arg)) + case ('-h') + call print_usage() + call exit(0) + case ('-out_bcs') + nxt = nxt + 1 + call getarg(nxt,arg) + out_bcsdir = trim(arg) + case ('-time') + nxt = nxt + 1 + call getarg(nxt,arg) + YYYYMMDDHHMM = trim(arg) + case ('-out_dir') + nxt = nxt + 1 + call getarg(nxt,arg) + out_dir = trim(arg) + case ('-model') + nxt = nxt + 1 + call getarg(nxt,arg) + model = trim(arg) + case ('-surflay') + nxt = nxt + 1 + call getarg(nxt,arg) + read(arg,*) surflay + case ('-in_tilefile') + nxt = nxt + 1 + call getarg(nxt,arg) + in_tilefile = trim(arg) + case ('-out_tilefile') + nxt = nxt + 1 + call getarg(nxt,arg) + out_tilefile = trim(arg) + case ('-in_rst') + nxt = nxt + 1 + call getarg(nxt,arg) + in_rstfile = trim(arg) + case ('-out_rst') + nxt = nxt + 1 + call getarg(nxt,arg) + out_rstfile = trim(arg) + case ('-in_wemin') + nxt = nxt + 1 + call getarg(nxt,arg) + read(arg,*) wemin_in + case ('-out_wemin') + nxt = nxt + 1 + call getarg(nxt,arg) + read(arg,*) wemin_out + case default + print*, trim(arg) + call print_usage() + print*, "wong command line" + call exit(1) + end select + nxt = nxt + 1 + call getarg(nxt,arg) + end do + if (index(model, 'catchcn') /=0 ) then + if((INDEX(out_bcsdir, 'NL') == 0).AND.(INDEX(out_bcsdir, 'OutData') == 0)) then + print *,'Land BCs in : ',trim(out_bcsdir) + print *,'do not support ',trim (model) + stop + endif + endif + out_file = out_dir //'/'//out_rstfile + + end subroutine + + subroutine print_usage() + print *,' ' + print *, 'This program can create catchment or catchmentCN restarts' + print *, 'depending on the command line option "model" ' + print *,' ' + print *,'-out_bcs : BC directory for output restart file' + print *,'-time : time for restart, format (yyyymmddhhmm)' + print *,'-out_dir : directory for output restasrt file' + print *,'-model : model ( catch, catchcnclm40, catchcnclm45)' + print *,'-surflay : surflay value' + print *,'-in_wemin : wemin for input restart' + print *,'-out_wemin : wemin for output restart' + print *,'-in_tilefile : tile_file for input restart' + print *,'-out_tilefile : tile_file for output restart, if none, it will search out_bcs' + print *,'-in_rst : input restart file name WITH path' + print *,'-out_rst : output restart file name WITHOUT path' + end subroutine +end program diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/catchplt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/catchplt similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/catchplt rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/catchplt diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/check_land_restarts.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/check_land_restarts.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/check_land_restarts.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/check_land_restarts.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catch_restart b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_catch_restart similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catch_restart rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_catch_restart diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catch_restart.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_catch_restart.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_catch_restart.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_catch_restart.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_vegdyn_restart b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_vegdyn_restart similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_vegdyn_restart rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_vegdyn_restart diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_vegdyn_restart.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_vegdyn_restart.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_vegdyn_restart.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/mk_vegdyn_restart.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/new_catch.ctl b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/new_catch.ctl similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/new_catch.ctl rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/new_catch.ctl diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/newcatch.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/newcatch.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/newcatch.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/newcatch.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/newvegdyn.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/newvegdyn.f90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/newvegdyn.f90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/newvegdyn.f90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/old_catch.ctl b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/old_catch.ctl similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/old_catch.ctl rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/old_catch.ctl diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/replace_params.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/replace_params.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/replace_params.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/replace_params.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/strip_vegdyn.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/strip_vegdyn.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/strip_vegdyn.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/obsolete/strip_vegdyn.F90