From 780efe613c0bbe82212e8ba5bd6dd9403b520448 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Paulo=20Chambel=20Leit=C3=A3o?= Date: Tue, 23 Mar 2021 10:26:08 +0000 Subject: [PATCH] Time series/discharges & curvilinear ghost corners Consistent Time series/discharges using X,Y location in curvilinear grids with ghost corners (not defined with a real location). --- .../ModuleNetCDFCF_2_HDF5MOHID.F90 | 249 +++++++++++++++++- Software/MOHIDBase1/ModuleGlobalData.F90 | 6 +- Software/MOHIDBase1/ModuleTimeSerie.F90 | 84 +++++- Software/MOHIDBase2/ModuleHorizontalGrid.F90 | 103 +++++++- Software/MOHIDWater/ModuleHydrodynamic.F90 | 59 +++-- Software/MOHIDWater/ModuleSand.F90 | 10 +- .../Convert2netcdf/Convert2netcdf.F90 | 16 ++ 7 files changed, 472 insertions(+), 55 deletions(-) diff --git a/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 b/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 index dbd389aea..9d45828d7 100644 --- a/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 +++ b/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 @@ -128,6 +128,8 @@ Module ModuleNetCDFCF_2_HDF5MOHID real, dimension(:,:), pointer :: LongOut, LatOut, RotationX, RotationY integer :: imax, jmax logical :: Imposed + logical :: ImposedRegular + character(Len=PathLength) :: MohidGrid real :: LongOrig, dLong real :: LatOrig, dLat logical :: Starts180W @@ -160,6 +162,7 @@ Module ModuleNetCDFCF_2_HDF5MOHID logical :: NetCDFNameFaceOff = .false. integer :: RemoveNsurfLayers = 0 real :: Offset = 0. + logical :: ZXYT = .false. end type T_Depth private :: T_Bathym @@ -3318,6 +3321,15 @@ subroutine ReadGridOptions STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ReadGridOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR170' + call GetData(Me%Depth%ZXYT, & + Me%ObjEnterData, iflag, & + SearchType = FromBlockInBlock, & + keyword = 'ZXYT', & + ClientModule = 'ModuleNetCDFCF_2_HDF5MOHID', & + default = .false., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR175' + call GetData(Me%LongLat%Imposed, & Me%ObjEnterData, iflag, & @@ -3330,6 +3342,21 @@ subroutine ReadGridOptions if (Me%LongLat%Imposed) then + call GetData(Me%LongLat%MohidGrid, & + Me%ObjEnterData, iflag, & + SearchType = FromBlockInBlock, & + keyword = 'MOHID_GRID', & + ClientModule = 'ModuleNetCDFCF_2_HDF5MOHID', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR183' + + if (iflag == 0) then + Me%LongLat%ImposedRegular = .true. + else + Me%LongLat%ImposedRegular = .false. + endif + + call GetData(Me%LongLat%LongOrig, & Me%ObjEnterData, iflag, & SearchType = FromBlockInBlock, & @@ -4642,6 +4669,98 @@ subroutine ImposedRegularGrid(ncid) end subroutine ImposedRegularGrid + !------------------------------------------------------------------------ + subroutine ImposedMohidGrid(ncid) + !Arguments------------------------------------------------------------- + integer :: ncid + + !Local----------------------------------------------------------------- + integer :: status, numDims, i, j + integer :: RhVarIdLong + integer, dimension(nf90_max_var_dims) :: rhDimIdsLong + integer :: ObjHorizontalGrid + + !Begin----------------------------------------------------------------- + + + + status=nf90_inq_varid(ncid,trim(Me%LongLat%NetCDFNameLong),RhVarIdLong) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR10' + + status = nf90_inquire_variable(ncid, RhVarIdLong, ndims = numDims) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR20' + + status = nf90_inquire_variable(ncid, RhVarIdLong, dimids = rhDimIdsLong(:numDims)) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR30' + + status=NF90_INQUIRE_DIMENSION(ncid, rhDimIdsLong(1), len = Me%LongLat%jmax) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR30' + + status=NF90_INQUIRE_DIMENSION(ncid, rhDimIdsLong(2), len = Me%LongLat%imax) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR50' + + if (Me%ReadInvertXY) then + + status=NF90_INQUIRE_DIMENSION(ncid, rhDimIdsLong(2), len = Me%LongLat%jmax) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR60' + + status=NF90_INQUIRE_DIMENSION(ncid, rhDimIdsLong(1), len = Me%LongLat%imax) + if (status /= nf90_noerr) stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR70' + + endif + + + !Build HDF5 MOHID Grid + Me%WorkSize%ILB = Me%LongLat%dij + Me%WorkSize%IUB = Me%LongLat%imax - 1 - Me%LongLat%dij + Me%WorkSize%JLB = Me%LongLat%dij + Me%WorkSize%JUB = Me%LongLat%jmax - 1 - Me%LongLat%dij + + + !to warn the user before the model crashes + !cant use a NetCDF with one of the dimension as 2 (or lower) because IUB or JUB would be zero (or lower). + if ((Me%WorkSize%IUB < 1) .or. (Me%WorkSize%JUB < 1)) then + write (*,*) + write (*,*) 'Please use a NETCDF file with more than' + write (*,*) '2x2 points so that the grid can be correctly extracted' + stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR80' + endif + + Me%Size%ILB = Me%WorkSize%ILB - 1 + Me%Size%IUB = Me%WorkSize%IUB + 1 + Me%Size%JLB = Me%WorkSize%JLB - 1 + Me%Size%JUB = Me%WorkSize%JUB + 1 + + allocate(Me%LongLat%LongOut (Me%Size%ILB:Me%Size%IUB,Me%Size%JLB:Me%Size%JUB)) + allocate(Me%LongLat%LatOut (Me%Size%ILB:Me%Size%IUB,Me%Size%JLB:Me%Size%JUB)) + + call ConstructHorizontalGrid(HorizontalGridID = ObjHorizontalGrid, & + DataFile = Me%LongLat%MohidGrid, & + STAT = status) + if (status /= SUCCESS_) then + stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR90' + endif + + call GetCornersCoordinates(HorizontalGridID = ObjHorizontalGrid, & + CoordX = Me%LongLat%LongOut, & + CoordY = Me%LongLat%LatOut, & + STAT = status) + if (status /= SUCCESS_) then + stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR100' + endif + + call KillHorizontalGrid(HorizontalGridID = ObjHorizontalGrid, & + STAT = status) + if (status /= SUCCESS_) then + stop 'ImposedMohidGrid - ModuleNetCDFCF_2_HDF5MOHID - ERR110' + endif + + + end subroutine ImposedMohidGrid + + !------------------------------------------------------------------------ + + !------------------------------------------------------------------------ subroutine ReadWriteGrid2D(ncid) !Arguments------------------------------------------------------------- @@ -4657,7 +4776,11 @@ subroutine ReadWriteGrid2D(ncid) !Begin----------------------------------------------------------------- if (Me%LongLat%Imposed) then + if (Me%LongLat%ImposedRegular) then call ImposedRegularGrid(ncid) + else + call ImposedMohidGrid (ncid) + endif else call ReadGrid2DNetCDF(ncid) endif @@ -6235,7 +6358,8 @@ subroutine ReadTimeNetCDF(ncid) integer :: ncid !Local----------------------------------------------------------------- - character(Len=StringLength) :: ref_date + !character(Len=StringLength) :: ref_date + character(Len=60) :: ref_date real, dimension(6) :: AuxTime real(8) :: Aux, HundredDays, Aux1 integer :: n, status, dimid, i, tmax, jmax @@ -6861,12 +6985,22 @@ subroutine ReadGrid3DNetCDF(ncid) Dim3 = Me%Date%NumberInst) else if (numDims == 4) then + !ZXYT + if (Me%Depth%ZXYT) then + + call AllocateValueIn(Me%Mapping%ValueIn, Dim1 = 1, & + Dim2 = Me%LongLat%jmax, & + Dim3 = Me%LongLat%imax, & + Dim4 = Me%Date%NumberInst) + + else call AllocateValueIn(Me%Mapping%ValueIn, Dim1 = Me%LongLat%jmax, & Dim2 = Me%LongLat%imax, & Dim3 = 1, & Dim4 = Me%Date%NumberInst) endif endif + endif call GetNetCDFMatrix(ncid, mn, Me%Mapping%ValueIn) @@ -6982,8 +7116,13 @@ subroutine ReadGrid3DNetCDF(ncid) else if (Me%Mapping%ValueIn%Dim == 2) then Aux = GetNetCDFValue(Me%Mapping%ValueIn, Dim1 = j+1, Dim2 = i+1) else if (Me%Mapping%ValueIn%Dim == 4) then + !ZXYT + if (Me%Depth%ZXYT) then + Aux = GetNetCDFValue(Me%Mapping%ValueIn, Dim1 = 1, Dim2 = j+1, Dim3 = i+1, Dim4 = Me%Mapping%Instant) + else Aux = GetNetCDFValue(Me%Mapping%ValueIn, Dim1 = j+1, Dim2 = i+1, Dim3 = 1, Dim4 = Me%Mapping%Instant) endif + endif if (Me%Mapping%Limit <= 1) then if (Aux > Me%Mapping%Limit) then @@ -7281,13 +7420,13 @@ subroutine ReadWaterLevelNetCDF(ncid, inst) if (status == nf90_noerr) then status = nf90_inquire_variable(ncid, pn, ndims = numDims) - if (status /= nf90_noerr) stop 'ReadFieldNetCDF - ModuleNetCDFCF_2_HDF5MOHID - ERR50' + if (status /= nf90_noerr) stop 'ReadWaterLevelNetCDF - ModuleNetCDFCF_2_HDF5MOHID - ERR50' Me%Depth%WLValueIn%Dim = numDims else write(*,*) 'possible error causes:' write(*,*) ' Variable not found.=', trim(Me%Depth%NetCDFNameWL) write(*,*) ' The specified netCDF ID does not refer to an open netCDF dataset.' - stop 'ReadFieldNetCDF - ModuleNetCDFCF_2_HDF5MOHID - ERR100' + stop 'ReadWaterLevelNetCDF - ModuleNetCDFCF_2_HDF5MOHID - ERR100' endif @@ -7890,6 +8029,7 @@ subroutine GetNetCDFMatrix(ncid, n, ValueIn, inst, layer, layer_dim) integer :: ILB, IUB, JLB, JUB integer :: Dim1, Dim2, Dim3, Dim4 integer :: iSt1, iSt2, iSt3, iSt4 + integer :: i, j !Begin----------------------------------------------------------------- @@ -8009,18 +8149,47 @@ subroutine GetNetCDFMatrix(ncid, n, ValueIn, inst, layer, layer_dim) else if (Dim==4) then + !ZXYT + if (Me%Depth%ZXYT) then + status = nf90_inquire_dimension(ncid, dimIDs(2), len = xdim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR170' + + status = nf90_inquire_dimension(ncid, dimIDs(3), len = ydim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR172' + + status = nf90_inquire_dimension(ncid, dimIDs(1), len = zdim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR174' + + endif + if (present(inst)) then + + !ZXYT + if (Me%Depth%ZXYT) then + status = nf90_inquire_dimension(ncid, dimIDs(2), len = xdim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR170' + + status = nf90_inquire_dimension(ncid, dimIDs(3), len = ydim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR172' + + status = nf90_inquire_dimension(ncid, dimIDs(1), len = zdim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR174' + + else status = nf90_inquire_dimension(ncid, dimIDs(1), len = xdim) - if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR170' + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR176' status = nf90_inquire_dimension(ncid, dimIDs(2), len = ydim) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR178' + + status = nf90_inquire_dimension(ncid, dimIDs(3), len = zdim) if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR180' - status = nf90_inquire_dimension(ncid, dimIDs(3), len = zdim) - if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR190' + endif status = nf90_inquire_dimension(ncid, dimIDs(4), len = tdim) - if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR191' + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR182' + if (JUB-JLB+1 /= xdim) then stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR192' @@ -8066,6 +8235,25 @@ subroutine GetNetCDFMatrix(ncid, n, ValueIn, inst, layer, layer_dim) else + !ZXYT + if (Me%Depth%ZXYT) then + + allocate(AuxR8D4(1:zdim,1:xdim,1:ydim,1:1)) + + status = nf90_get_var(ncid,n,AuxR8D4, & + start = (/ 1, 1, 1, inst /), & + count = (/ zdim, xdim, ydim, 1 /)) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR204' + + do i =ILB,IUB + do j =JLB, JUB + ValueIn%R84D(j,i,1:zdim,1:1) = AuxR8D4(1:zdim,j,i,1:1) + enddo + enddo + + + else + allocate(AuxR8D4(1:xdim,1:ydim,1:zdim,1:1)) status = nf90_get_var(ncid,n,AuxR8D4, & @@ -8075,6 +8263,9 @@ subroutine GetNetCDFMatrix(ncid, n, ValueIn, inst, layer, layer_dim) ValueIn%R84D(JLB:JUB,ILB:IUB,1:zdim,1:1) = AuxR8D4(1:xdim,1:ydim,1:zdim,1:1) + endif + + deallocate(AuxR8D4) endif @@ -8120,8 +8311,41 @@ subroutine GetNetCDFMatrix(ncid, n, ValueIn, inst, layer, layer_dim) else - status = NF90_GET_VAR(ncid,n,ValueIn%R44D) - if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR216' + !ZXYT + if (Me%Depth%ZXYT) then + + allocate(AuxR4D4(1:zdim,1:xdim,1:ydim,1:1)) + + status = nf90_get_var(ncid,n,AuxR4D4, & + start = (/ 1, 1, 1, 1 /), & + count = (/ zdim, xdim, ydim, 1 /)) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR204' + + do i =ILB,IUB + do j =JLB, JUB + ValueIn%R44D(j,i,1:zdim,1:1) = AuxR4D4(1:zdim,j,i,1:1) + enddo + enddo + + + else + + allocate(AuxR4D4(1:xdim,1:ydim,1:zdim,1:1)) + + status = nf90_get_var(ncid,n,AuxR4D4, & + start = (/ 1, 1, 1, 1 /), & + count = (/ xdim, ydim, zdim, 1 /)) + if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR202' + + ValueIn%R44D(JLB:JUB,ILB:IUB,1:zdim,1:1) = AuxR4D4(1:xdim,1:ydim,1:zdim,1:1) + + endif + + + deallocate(AuxR4D4) + + !status = NF90_GET_VAR(ncid,n,ValueIn%R44D) + !if (status /= nf90_noerr) stop 'GetNetCDFMatrix - ModuleNetCDFCF_2_HDF5MOHID - ERR216' endif @@ -8211,6 +8435,13 @@ function GetNetCDFValue(ValueIn, Dim1, Dim2, Dim3, Dim4) if (Me%ReadMirrorLon .and. Dim >=2) then Dim1_ = Me%LongLat%jmax - Dim1_ + 1 endif + + !!ZXYT + !if (Me%Depth%ZXYT .and. Dim > 3) then + ! Dim2_ = Dim4 + ! Dim3_ = Dim2 + ! Dim4_ = Dim3 + !endif if (Dim==1) then diff --git a/Software/MOHIDBase1/ModuleGlobalData.F90 b/Software/MOHIDBase1/ModuleGlobalData.F90 index db1b942d0..743740c33 100644 --- a/Software/MOHIDBase1/ModuleGlobalData.F90 +++ b/Software/MOHIDBase1/ModuleGlobalData.F90 @@ -76,7 +76,7 @@ Module ModuleGlobalData end interface SetError !Parameter----------------------------------------------------------------- - integer, parameter :: MaxModules = 99 + integer, parameter :: MaxModules = 100 #if defined(_INCREASE_MAXINSTANCES) integer, parameter :: MaxInstances = 2000 @@ -1983,6 +1983,7 @@ Module ModuleGlobalData integer, parameter :: mTwoWay_ = 97 integer, parameter :: mOutputGrid_ = 98 integer, parameter :: mMeshGlue_ = 99 + integer, parameter :: mDelftFM_2_MOHID_ = 100 !Domain decomposition integer, parameter :: WestSouth = 1 @@ -2111,7 +2112,8 @@ Module ModuleGlobalData T_Module(mSediment_ , "Sediment" ), T_Module(mReservoirs_ , "Reservoirs" ), & T_Module(mIrrigation_ , "Irrigation" ), T_Module(mTURBINE_ , "Turbine" ), & T_Module(mLitter_ , "Litter" ), T_Module(mTwoWay_ , "TwoWay" ), & - T_Module(mOutputGrid_ , "OutputGrid" ), T_Module(mMeshGlue_ , "MeshGlue" )/) + T_Module(mOutputGrid_ , "OutputGrid" ), T_Module(mMeshGlue_ , "MeshGlue" ), & + T_Module(mDelftFM_2_MOHID_ , "DelftFM_2_MOHID" )/) !Variables diff --git a/Software/MOHIDBase1/ModuleTimeSerie.F90 b/Software/MOHIDBase1/ModuleTimeSerie.F90 index 1f6999d8a..8371b4682 100644 --- a/Software/MOHIDBase1/ModuleTimeSerie.F90 +++ b/Software/MOHIDBase1/ModuleTimeSerie.F90 @@ -64,6 +64,7 @@ Module ModuleTimeSerie public :: WriteSpecificTimeSerieLine public :: CorrectsCellsTimeSerie public :: TryIgnoreTimeSerie + public :: SetIgnoreTimeSerie public :: ReseatCurrentIndex !Selector @@ -183,6 +184,7 @@ Module ModuleTimeSerie logical :: TimeSerie1D = .false. logical :: ComputeResidual = .true. !logical :: IgnoreON = .false. + logical :: GhostCorners = .false. !TimeSerieInput logical :: TimeCycle = .false. @@ -239,7 +241,8 @@ subroutine StartTimeSerie(TimeSerieID, ObjTime, TimeSerieDataFile, PropertyList, Extension, WaterPoints3D, & WaterPoints2D, WaterPoints1D, ResultFileName, Instance, & ModelName, CoordX, CoordY, UseTabulatedData, & - HavePath, Comment, ModelDomain, ReplacePath, STAT) + HavePath, Comment, ModelDomain, ReplacePath, GhostCorners, & + STAT) !Arguments------------------------------------------------------------- integer :: TimeSerieID @@ -260,6 +263,7 @@ subroutine StartTimeSerie(TimeSerieID, ObjTime, character(len=*), optional, intent(IN ) :: Comment type (T_Polygon), pointer, optional :: ModelDomain character(len=*), optional, intent(IN ) :: ReplacePath + logical, optional, intent(IN ) :: GhostCorners integer, optional, intent(OUT) :: STAT !Local----------------------------------------------------------------- @@ -317,6 +321,12 @@ subroutine StartTimeSerie(TimeSerieID, ObjTime, Me%ModelDomainON = .false. endif + if (present(GhostCorners)) then + Me%GhostCorners = GhostCorners + else + Me%GhostCorners = .false. + endif + !Constructs EnterData call ConstructEnterData(Me%ObjEnterData, TimeSerieDataFile, STAT = STAT_CALL) @@ -366,17 +376,6 @@ subroutine StartTimeSerie(TimeSerieID, ObjTime, endif - !call GetData(Me%IgnoreON, & - ! Me%ObjEnterData, & - ! flag, & - ! SearchType = FromFile, & - ! keyword ='IGNORE_ON', & - ! Default = .false., & - ! ClientModule ='ModuleTimeSerie', & - ! STAT = STAT_CALL) - !if (STAT_CALL .NE. SUCCESS_) & - ! call SetError(FATAL_, KEYWORD_, "Subroutine StartTimeSerie - ModuleTimeSerie. ERR50") - !Stores the number of properties Me%NumberOfProperties = size(PropertyList) @@ -903,10 +902,15 @@ integer function ReadTimeSeriesLocation () TimeSeriesXY%X = Me%TimeSerie(iTimeSerie)%CoordX TimeSeriesXY%Y = Me%TimeSerie(iTimeSerie)%CoordY + Me%TimeSerie(iTimeSerie)%IgnoreON = .false. + + if (Me%GhostCorners) then + Me%TimeSerie(iTimeSerie)%IgnoreON = .false. + else if (.not. IsPointInsidePolygon(TimeSeriesXY, Me%ModelDomain)) then Me%TimeSerie(iTimeSerie)%IgnoreON = .true. endif - + endif deallocate(TimeSeriesXY) endif i12 @@ -3523,6 +3527,60 @@ subroutine TryIgnoreTimeSerie(TimeSerieID, Number, IgnoreOK, STAT) !---------------------------------------------------------------------- end subroutine TryIgnoreTimeSerie + + + !-------------------------------------------------------------------------- + + subroutine SetIgnoreTimeSerie(TimeSerieID, Number, IgnoreOK, STAT) + + !Arguments-------------------------------------------------------------- + integer :: TimeSerieID + integer, intent(IN) :: Number + logical, intent(IN) :: IgnoreOK + integer, intent(OUT), optional :: STAT + + + + !Local----------------------------------------------------------------- + integer :: ready_ + integer :: STAT_ + integer :: STAT_CALL + + !---------------------------------------------------------------------- + + STAT_ = UNKNOWN_ + + call Ready(TimeSerieID, ready_) + +cd1 : if (ready_ .EQ. IDLE_ERR_) then + + Me%TimeSerie(Number)%IgnoreON = IgnoreOK + + if (IgnoreOK) then + call UnitsManager(unit = Me%TimeSerie(Number)%UnitNumber, & + OPENCLOSE = CLOSE_FILE, & + STATUS ='DELETE', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + write(*,*) 'SetIgnoreTimeSerie - ModuleTimeSerie - WRN10' + endif + endif + + STAT_ = SUCCESS_ + else + + STAT_ = ready_ + + end if cd1 + + + if (present(STAT)) & + STAT = STAT_ + + !---------------------------------------------------------------------- + + end subroutine SetIgnoreTimeSerie + !-------------------------------------------------------------------------- !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/Software/MOHIDBase2/ModuleHorizontalGrid.F90 b/Software/MOHIDBase2/ModuleHorizontalGrid.F90 index 60a5237e8..5b99fc7c4 100644 --- a/Software/MOHIDBase2/ModuleHorizontalGrid.F90 +++ b/Software/MOHIDBase2/ModuleHorizontalGrid.F90 @@ -148,6 +148,7 @@ Module ModuleHorizontalGrid public :: GetDDecompMapping2D public :: GetDDecompMPI_ID public :: GetDDecompON + public :: GetGhostCorners public :: GetSonWindow @@ -590,7 +591,7 @@ Module ModuleHorizontalGrid logical :: CornersXYInput = .false. logical :: Distortion = .false. logical :: RegularRotation = .false. - logical :: GostCorners = .false. + logical :: GhostCorners = .false. integer, dimension(:,:), pointer :: DefineCellsMap => null() integer, dimension(:,:), pointer :: DefineFacesUMap => null() @@ -3480,7 +3481,7 @@ subroutine ConstructGlobalVariables endif if (Me%XX_IE(ii, jj) < FillValueReal/2) then - Me%GostCorners = .true. + Me%GhostCorners = .true. Cycle endif @@ -11811,7 +11812,6 @@ integer function GetDDecompMPI_ID(HorizontalGridID, STAT) end function GetDDecompMPI_ID !-------------------------------------------------------------------------- - !-------------------------------------------------------------------------- logical function GetDDecompON(HorizontalGridID, STAT) @@ -11854,6 +11854,46 @@ end function GetDDecompON !-------------------------------------------------------------------------- + + logical function GetGhostCorners(HorizontalGridID, STAT) + + !Arguments------------------------------------------------------------- + integer, intent(IN ) :: HorizontalGridID + integer, optional, intent(OUT) :: STAT + + !External-------------------------------------------------------------- + + integer :: ready_ + + !Local----------------------------------------------------------------- + + integer :: STAT_ !Auxiliar local variable + + !---------------------------------------------------------------------- + + STAT_ = UNKNOWN_ + + GetGhostCorners = .false. + + call Ready(HorizontalGridID, ready_) + +cd1 : if ((ready_ .EQ. IDLE_ERR_ ) .OR. & + (ready_ .EQ. READ_LOCK_ERR_)) then + + GetGhostCorners = Me%GhostCorners + + STAT_ = SUCCESS_ + else + STAT_ = ready_ + end if cd1 + + + if (present(STAT)) STAT = STAT_ + + !---------------------------------------------------------------------- + + end function GetGhostCorners + !-------------------------------------------------------------------------- function GetDDecompOpenBorders(HorizontalGridID, STAT) @@ -18416,6 +18456,50 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, Point%X = XPoint Point%Y = YPoint + f0: if (Me%GhostCorners) then + + dj1: do j = JLB, JUB + di1: do i = ILB, IUB + + if(DefinedPoint(i,j) == 1)then + + Me%AuxPolygon%VerticesF(1)%X = XX2D_Z(i ,j ) + Me%AuxPolygon%VerticesF(1)%Y = YY2D_Z(i ,j ) + + Me%AuxPolygon%VerticesF(2)%X = XX2D_Z(i+1,j ) + Me%AuxPolygon%VerticesF(2)%Y = YY2D_Z(i+1,j ) + + Me%AuxPolygon%VerticesF(3)%X = XX2D_Z(i+1,j+1) + Me%AuxPolygon%VerticesF(3)%Y = YY2D_Z(i+1,j+1) + + Me%AuxPolygon%VerticesF(4)%X = XX2D_Z(i ,j+1) + Me%AuxPolygon%VerticesF(4)%Y = YY2D_Z(i ,j+1) + + !close polygon + Me%AuxPolygon%VerticesF(5)%X = Me%AuxPolygon%VerticesF(1)%X + Me%AuxPolygon%VerticesF(5)%Y = Me%AuxPolygon%VerticesF(1)%Y + + call SetLimits(Me%AuxPolygon) + + if (IsPointInsidePolygon(Point, Me%AuxPolygon)) then + + SearchCell = .false. + CellFound = .true. + + IZ = i + JZ = j + + exit dj1 + + endif + + endif + + enddo di1 + enddo dj1 + + else f0 + !check if f1: if (present(Iold) .and. present(Jold)) then @@ -18485,8 +18569,8 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, IUpper = IUB end if - do j = JLower, JUpper - do i = ILower, IUpper +dj2: do j = JLower, JUpper +di2: do i = ILower, IUpper if(i == iold .and. j == jold) cycle @@ -18518,15 +18602,14 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, IZ = i JZ = j - exit + exit dj2 endif endif - enddo - if(CellFound) exit !exit do J loop - enddo + enddo di2 + enddo dj2 end if f3 @@ -18690,6 +18773,8 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, endif f4 endif f2 + + endif f0 if (present(CellLocated)) CellLocated = CellFound diff --git a/Software/MOHIDWater/ModuleHydrodynamic.F90 b/Software/MOHIDWater/ModuleHydrodynamic.F90 index 411111ca8..e513e5d08 100644 --- a/Software/MOHIDWater/ModuleHydrodynamic.F90 +++ b/Software/MOHIDWater/ModuleHydrodynamic.F90 @@ -126,7 +126,8 @@ Module ModuleHydrodynamic GetNumberOfTimeSeries, TryIgnoreTimeSerie, & GetTimeSerieValue, WriteTimeSerie, & GetTimeSerieName, WriteTimeSerieLine, & - GetTimeSerieNextOutput, KillTimeSerie + GetTimeSerieNextOutput, SetIgnoreTimeSerie, & + KillTimeSerie use ModuleHorizontalMap, only : GetWaterPoints2D, GetBoundaries, GetBoundaryFaces,& GetExteriorBoundaryFaces, UnGetHorizontalMap, & GetLandPoints2D @@ -148,7 +149,7 @@ Module ModuleHydrodynamic GetGridOutBorderPolygon, & GetDDecompWorkSize2D, WriteHorizontalGrid_UV, & GetCellRotation, GetGridCellArea, GetConnections, & - GetCornersCoordinates + GetCornersCoordinates, GetGhostCorners use ModuleTwoWay, only : ConstructTwoWayHydrodynamic, ModifyTwoWay, & AllocateTwoWayAux, PrepTwoWay, & UngetTwoWayExternal_Vars, & @@ -12344,7 +12345,7 @@ subroutine Construct_Time_Serie !Local----------------------------------------------------------------- real :: CoordX, CoordY - logical :: CoordON, IgnoreOK + logical :: CoordON, IgnoreOK, GhostCorners integer :: nProperties integer :: STAT_CALL integer :: iflag, dn, Id, Jd, TimeSerieNumber, i, j @@ -12386,11 +12387,14 @@ subroutine Construct_Time_Serie ClientModule = 'ModuleHydrodynamic', & Default = Me%Files%ConstructData, & STAT = STAT_CALL) - if (STAT_CALL .NE. SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR30' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR30' call GetGridOutBorderPolygon(HorizontalGridID = Me%ObjHorizontalGrid, & Polygon = ModelDomainLimit, STAT = STAT_CALL) - if (STAT_CALL .NE. SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR35' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR40' + + GhostCorners = GetGhostCorners(HorizontalGridID = Me%ObjHorizontalGrid, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR50' !Constructs TimeSerie call StartTimeSerie(Me%ObjTimeSerie, Me%ObjTime, & @@ -12398,19 +12402,22 @@ subroutine Construct_Time_Serie PropertyList, "srh", & WaterPoints3D = Me%External_Var%WaterPoints3D, & ModelName = Me%ModelName, & - ModelDomain = ModelDomainLimit, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR40") + ModelDomain = ModelDomainLimit, & + GhostCorners = GhostCorners, & + STAT = STAT_CALL) + + if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR60") call UngetHorizontalGrid(HorizontalGridID = Me%ObjHorizontalGrid, Polygon= ModelDomainLimit, STAT = STAT_CALL) - if (STAT_CALL .NE. SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR45' + if (STAT_CALL .NE. SUCCESS_) stop 'Construct_Time_Serie - Hydrodynamic - ERR70' !Deallocates PropertyList deallocate(PropertyList, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR50") + if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR80") !Corrects if necessary the cell of the time serie based in the time serie coordinates call GetNumberOfTimeSeries(Me%ObjTimeSerie, TimeSerieNumber, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR60' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR90' do dn = 1, TimeSerieNumber @@ -12421,24 +12428,35 @@ subroutine Construct_Time_Serie call GetTimeSerieLocation(Me%ObjTimeSerie, dn, & CoordX = CoordX, CoordY = CoordY, CoordON = CoordON, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR70' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR110' call GetTimeSerieName(Me%ObjTimeSerie, dn, TimeSerieName, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR80' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR120' i1: if (CoordON) then call GetXYCellZ(Me%ObjHorizontalGrid, CoordX, CoordY, Id, Jd, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_ .and. STAT_CALL /= OUT_OF_BOUNDS_ERR_) & - stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR90' + if (STAT_CALL == OUT_OF_BOUNDS_ERR_) then + if (GhostCorners) then + + call SetIgnoreTimeSerie(Me%ObjTimeSerie, dn, IgnoreOK = .true., STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR130' + + cycle + endif + else + if (STAT_CALL /= SUCCESS_) then + stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR140' + endif + endif call CorrectsCellsTimeSerie(Me%ObjTimeSerie, dn, Id, Jd, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR120' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR150' endif i1 call GetTimeSerieLocation(Me%ObjTimeSerie, dn, LocalizationI = Id, LocalizationJ = Jd, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR130' + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR160' if (Me%External_Var%WaterPoints3D(Id, Jd, Me%WorkSize%KUB) /= WaterPoint) & write(*,*) 'Time Serie in a land cell - ',trim(TimeSerieName),' - ',trim(Me%ModelName) @@ -12447,7 +12465,7 @@ subroutine Construct_Time_Serie call UnGetMap(Me%ObjMap, Me%External_Var%WaterPoints3D, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR140") + if (STAT_CALL /= SUCCESS_) call SetError (FATAL_, OUT_OF_MEM_, "Construct_Time_Serie - Hydrodynamic - ERR170") end subroutine Construct_Time_Serie @@ -52351,6 +52369,7 @@ subroutine OutPut_TimeSeries integer :: STAT_CALL, TimeSerieNumber, dn, id, jd, kd logical :: DepthON, IgnoreOK integer :: IUB,ILB,JUB,JLB,KUB,KLB + character(len=StringLength) :: TimeSerieName !Begin------------------------------------------------------------------- IUB = Me%WorkSize%IUB; JUB = Me%WorkSize%JUB; KUB = Me%WorkSize%KUB ILB = Me%WorkSize%ILB; JLB = Me%WorkSize%JLB; KLB = Me%WorkSize%KLB @@ -52398,7 +52417,11 @@ subroutine OutPut_TimeSeries if (Me%External_Var%WaterPoints3D(Id, JD, kd) /= WaterPoint .and. Me%FirstIteration) then - write(*,*) 'Time serie station I=',Id, 'J=',Jd,'K=',Kd,'is located in land' + call GetTimeSerieName(Me%ObjTimeSerie, dn, TimeSerieName, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'Construct_Time_Serie - ModuleHydrodynamic - ERR120' + + write(*,*) 'Time serie station ', trim(TimeSerieName), ' is located in land' + write(*,*) 'Time serie station I=',Id, 'J=',Jd,'K=',Kd write(*,*) 'OutPut_TimeSeries - ModuleHydrodynamic - WRN100' endif diff --git a/Software/MOHIDWater/ModuleSand.F90 b/Software/MOHIDWater/ModuleSand.F90 index 1e9c40b75..deebe6346 100644 --- a/Software/MOHIDWater/ModuleSand.F90 +++ b/Software/MOHIDWater/ModuleSand.F90 @@ -2098,7 +2098,7 @@ subroutine ConstructGlobalParameters Me%ObjEnterData,iflag, & SearchType = FromFile, & keyword = 'POROSITY', & - default = 0.1, & + default = 0.4, & ClientModule = 'ModuleSand', & STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalParameters - ModuleSand - ERR200' @@ -3707,7 +3707,7 @@ end subroutine ComputeBiHamonicFilter2D subroutine MeyerPeterTransport !Local----------------------------------------------------------------- real :: Tr1,Tr2,Tr3, DeltaTau, Miu, CChezy, Clinha, Depth - real :: DeltaTauNoDim + real :: DeltaTauNoDim, Tau integer :: i, j !---------------------------------------------------------------------- @@ -3729,8 +3729,10 @@ subroutine MeyerPeterTransport Clinha= 18.*LOG10(12.*Depth/Me%D90%Field2D(I,J)) Miu = (CChezy/Clinha)**1.5 + Tau = Me%ExternalVar%TauTotal(I,J) + !N/m2 = [kg * m/s^2 * m] - DeltaTau = min(Miu*Me%ExternalVar%TauTotal(I,J),Me%TauMax)-Me%TauCritic(I,J) + DeltaTau = min(Miu*Tau,Me%TauMax)-Me%TauCritic(I,J) ![ ] = [kg * m/s^2 * m] / [kg] / [m/s^2] / [m] DeltaTauNoDim = DeltaTau / (Me%RhoSl * Gravity * Me%D50%Field2D(I,J)) @@ -3745,7 +3747,7 @@ subroutine MeyerPeterTransport Tr2 = Me%D50%Field2D(I,J)**1.5 ![ ] = [ ] Tr3 = DeltaTauNoDim**1.5 - ![ ] + ![m2/s] Me%TransportCapacity(i, j) = 8.*Tr1*Tr2*Tr3 endif diff --git a/Software/SmallTools/Convert2netcdf/Convert2netcdf.F90 b/Software/SmallTools/Convert2netcdf/Convert2netcdf.F90 index 1d3c6c533..0ae227b29 100644 --- a/Software/SmallTools/Convert2netcdf/Convert2netcdf.F90 +++ b/Software/SmallTools/Convert2netcdf/Convert2netcdf.F90 @@ -2439,18 +2439,34 @@ subroutine BuildAttributes(Name, NCDFName, LongName, StandardName, Units, ! Next 2 keywords are not CF compliant case("mean_wave_direction_X") + NCDFName = "wave_X" + + if (Me%OdysseaProject) then + LongName = "eastward wave direction unit vector" + StandardName = "eastward_wave_direction_unit_vector" + else LongName = "mean wave direction X" StandardName = "mean_wave_direction_X" + endif + Units_ = "-" ValidMin = -1. ValidMax = 1. MissingValue = Me%MissingValue case("mean_wave_direction_Y") + NCDFName = "wave_Y" + + if (Me%OdysseaProject) then + LongName = "northward wave direction unit vector" + StandardName = "northward_wave_direction_unit_vector" + else LongName = "mean wave direction Y" StandardName = "mean_wave_direction_Y" + endif + Units_ = "-" ValidMin = -1. ValidMax = 1.