diff --git a/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 b/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 index f7e907e3b..9d45828d7 100644 --- a/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 +++ b/Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90 @@ -128,11 +128,14 @@ 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 integer :: CorrectJDown, CorrectJUp, BreakJ integer :: dij + logical :: CorrectLongON end type T_LongLat @@ -159,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 @@ -281,6 +285,8 @@ Module ModuleNetCDFCF_2_HDF5MOHID type(T_Bathym) :: Bathym logical :: OutHDF5, OutNetcdf, ReadInvertXY logical :: ReadInvertLat + logical :: ReadInvertLat_V2 + logical :: ReadMirrorLon logical :: Nan_2_Null integer :: OutCountProp = 0 type(T_WindowOut) :: WindowOut @@ -2540,6 +2546,24 @@ subroutine ReadOptions STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ReadOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR130' + call GetData(Me%ReadInvertLat_V2, & + Me%ObjEnterData, iflag, & + SearchType = FromBlock, & + keyword = 'READ_INVERT_LAT_V2', & + default = .false., & + ClientModule = 'ModuleNetCDFCF_2_HDF5MOHID', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR132' + + call GetData(Me%ReadMirrorLon, & + Me%ObjEnterData, iflag, & + SearchType = FromBlock, & + keyword = 'READ_MIRROR_LON', & + default = .false., & + ClientModule = 'ModuleNetCDFCF_2_HDF5MOHID', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR135' + call GetData(Aux4, & Me%ObjEnterData, iflag, & @@ -2949,6 +2973,15 @@ subroutine ReadGridOptions STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ReadGridOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR30' + call GetData(Me%LongLat%CorrectLongON, & + Me%ObjEnterData, iflag, & + SearchType = FromBlockInBlock, & + keyword = 'CORRECT_LONG_ON', & + ClientModule = 'ModuleNetCDFCF_2_HDF5MOHID', & + default = .true., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridOptions - ModuleNetCDFCF_2_HDF5MOHID - ERR35' + Me%LongLat%LatIn%DataType = Real8_ Me%LongLat%LongIn%DataType = Real8_ @@ -3288,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, & @@ -3300,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, & @@ -4612,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------------------------------------------------------------- @@ -4627,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 @@ -4727,6 +4880,7 @@ subroutine ReadWriteGrid3D(ncid) endif if (.not. Me%Bathym%FromMapping) then + call WriteBathymASCII if (Me%OutHDF5 ) call WriteBathymHDF5 if (Me%OutNetCDF) call WriteBathymNetCDF endif @@ -5500,7 +5654,7 @@ subroutine WriteFieldAllInst(iFinal, iP, iT) write(*,*) 'Netcdf value', Me%Field(iP)%Value2DOut(i, j) write(*,*) 'Above Max limit=', Me%Field(iP)%MaxLimit - stop 'WriteFieldAllInst - ModuleNetCDFCF_2_HDF5MOHID - ERR20' + stop 'WriteFieldAllInst - ModuleNetCDFCF_2_HDF5MOHID - ERR40' endif @@ -5511,7 +5665,7 @@ subroutine WriteFieldAllInst(iFinal, iP, iT) write(*,*) 'Netcdf value', Me%Field(iP)%Value2DOut(i, j) write(*,*) 'Below min limit=', Me%Field(iP)%MinLimit - stop 'WriteFieldAllInst - ModuleNetCDFCF_2_HDF5MOHID - ERR30' + stop 'WriteFieldAllInst - ModuleNetCDFCF_2_HDF5MOHID - ERR50' endif @@ -6204,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 @@ -6675,6 +6830,8 @@ subroutine ReadGrid2DNetCDF(ncid) enddo enddo + if (Me%LongLat%CorrectLongON) then + do j=1, Me%LongLat%jmax-1 Aux1 = GetNetCDFValue(Me%LongLat%LongIn, Dim1 = j, Dim2 = 1 ) Aux2 = GetNetCDFValue(Me%LongLat%LongIn, Dim1 = j+1, Dim2 = 1 ) @@ -6699,7 +6856,7 @@ subroutine ReadGrid2DNetCDF(ncid) exit endif enddo - + endif !Build HDF5 MOHID Grid Me%WorkSize%ILB = Me%LongLat%dij @@ -6739,13 +6896,19 @@ subroutine ReadGrid2DNetCDF(ncid) Y3 = GetNetCDFValue(Me%LongLat%LatIn, Dim1 = j, Dim2 = i+1) Y4 = GetNetCDFValue(Me%LongLat%LatIn, Dim1 = j+1, Dim2 = i+1) - Me%LongLat%LongOut(i, j) = (X1 + X2 + X3 + X4) / 4. + if (Me%ReadInvertLat) then Me%LongLat%LatOut (Me%WorkSize%IUB+1+Me%WorkSize%ILB-i, j) = (Y1 + Y2 + Y3 + Y4) / 4. else Me%LongLat%LatOut (i, j) = (Y1 + Y2 + Y3 + Y4) / 4. endif + if (Me%ReadMirrorLon) then + Me%LongLat%LongOut(i, j) = -(X1 + X2 + X3 + X4) / 4. + else + Me%LongLat%LongOut(i, j) = (X1 + X2 + X3 + X4) / 4. + endif + enddo enddo @@ -6822,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) @@ -6943,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 @@ -7242,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 @@ -7851,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----------------------------------------------------------------- @@ -7970,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' @@ -8027,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, & @@ -8036,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 @@ -8081,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 @@ -8165,10 +8428,21 @@ function GetNetCDFValue(ValueIn, Dim1, Dim2, Dim3, Dim4) Dim2_ = Dim1 endif - if (Me%ReadInvertLat .and. Dim >=2) then + if ((Me%ReadInvertLat .or. Me%ReadInvertLat_V2).and. Dim >=2) then Dim2_ = Me%LongLat%imax - Dim2_ + 1 endif + 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 if (DataTypeIn == Real8_ ) then @@ -8292,10 +8566,14 @@ subroutine SetNetCDFValue(ValueIn, SetValue, Dim1, Dim2, Dim3, Dim4) Dim2_ = Dim1 endif - if (Me%ReadInvertLat .and. Dim >=2) then + if ((Me%ReadInvertLat .or. Me%ReadInvertLat_V2) .and. Dim >=2) then Dim2_ = Me%LongLat%imax - Dim2_ + 1 endif + if (Me%ReadMirrorLon .and. Dim >=2) then + Dim1_ = Me%LongLat%jmax - Dim1_ + 1 + endif + if (Dim==1) then if (DataTypeIn == Real8_ ) then diff --git a/Software/MOHIDBase1/ModuleDrawing.F90 b/Software/MOHIDBase1/ModuleDrawing.F90 index 513b9843c..b4ae4d9bc 100644 --- a/Software/MOHIDBase1/ModuleDrawing.F90 +++ b/Software/MOHIDBase1/ModuleDrawing.F90 @@ -969,17 +969,25 @@ subroutine SetLimitsPolygon(Polygon) d1: do i=1, Polygon%Count - if (Polygon%VerticesF(i)%X < Polygon%Limits%Left) & + if (Polygon%VerticesF(i)%X < Polygon%Limits%Left .and. & + Polygon%VerticesF(i)%X > FillValueReal / 2.) then Polygon%Limits%Left = Polygon%VerticesF(i)%X + endif - if (Polygon%VerticesF(i)%X > Polygon%Limits%Right) & + if (Polygon%VerticesF(i)%X > Polygon%Limits%Right .and. & + Polygon%VerticesF(i)%X > FillValueReal / 2.) then Polygon%Limits%Right = Polygon%VerticesF(i)%X + endif - if (Polygon%VerticesF(i)%Y < Polygon%Limits%Bottom) & + if (Polygon%VerticesF(i)%Y < Polygon%Limits%Bottom .and. & + Polygon%VerticesF(i)%Y > FillValueReal / 2.) then Polygon%Limits%Bottom = Polygon%VerticesF(i)%Y + endif - if (Polygon%VerticesF(i)%Y > Polygon%Limits%Top) & + if (Polygon%VerticesF(i)%Y > Polygon%Limits%Top .and. & + Polygon%VerticesF(i)%Y > FillValueReal / 2.) then Polygon%Limits%Top = Polygon%VerticesF(i)%Y + endif enddo d1 @@ -1504,6 +1512,16 @@ logical function IsPointInsidePolygon(Point, Polygon) do1: do CurrentVertix = 1, Polygon%Count - 1 + if (Polygon%VerticesF(CurrentVertix)%X == FillValueReal) then + write (*,*) 'Vertice with invalid X coordinate' + stop 'IsPointInsidePolygon - ModuleDrawing - ERR10' + endif + + if (Polygon%VerticesF(CurrentVertix)%Y == FillValueReal) then + write (*,*) 'Vertice with invalid Y coordinate' + stop 'IsPointInsidePolygon - ModuleDrawing - ERR20' + endif + !construct segment Segment%StartAt%X = Polygon%VerticesF(CurrentVertix)%X Segment%StartAt%Y = Polygon%VerticesF(CurrentVertix)%Y diff --git a/Software/MOHIDBase1/ModuleFunctions.F90 b/Software/MOHIDBase1/ModuleFunctions.F90 index 102ff9f11..6a8b5e29a 100644 --- a/Software/MOHIDBase1/ModuleFunctions.F90 +++ b/Software/MOHIDBase1/ModuleFunctions.F90 @@ -249,6 +249,7 @@ Module ModuleFunctions !Polygon public :: RelativePosition4VertPolygon + public :: RelativePosition4VertPolygonIWD public :: PolygonArea public :: FromGeo2Meters @@ -7928,95 +7929,172 @@ subroutine RelativePosition4VertPolygon(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, !Local------------------------------------------------------------------- real(8) :: DXdc, DYac, DXba, DYbd, DXef, DYeg - real(8) :: MinDx, SumAux + real(8) :: MinDx real(8) :: a1, b1, a2, b2, a3, b3, a4, b4 real(8) :: Seg_ac, Seg_dc,Seg_hc, Seg_ic real(8) :: Xf, Yf, Xg, Yg, Xh, Yh, Xi, Yi, TgX, TgY + real(8) :: SumAux real(8) :: XaR8, YaR8, XbR8, YbR8, XcR8, YcR8, XdR8, YdR8, XeR8, YeR8 !Begin------------------------------------------------------------------- XaR8 = dble(Xa); YaR8 = dble(Ya); XbR8 = dble(Xb); YbR8 = dble(Yb); XcR8 = dble(Xc) YcR8 = dble(Yc); XdR8 = dble(Xd); YdR8 = dble(Yd); XeR8 = dble(Xe); YeR8 = dble(Ye) + + + !the four segments of the cell DXdc = XdR8 - XcR8 DXba = XbR8 - XaR8 DYac = YaR8 - YcR8 DYbd = YbR8 - YdR8 - - + SumAux = abs(DXdc + DXba + DYac + DYbd) - MinDx = SumAux * 1.e-16 - + + if (SumAux > 0) then + MinDx = SumAux * 1.e-16 + else + MinDx = 1.e-16 + endif + if (abs(DXdc) 1 .or. Yey < 0. .or. Yey > 1) then + call RelativePosition4VertPolygonIWD(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey) + endif + end subroutine RelativePosition4VertPolygon !-------------------------------------------------------------------------- + + subroutine RelativePosition4VertPolygonIWD(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey) + + !Arguments--------------------------------------------------------------- + real :: Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey + + !Local------------------------------------------------------------------- + real(8) :: SumAux + real(8) :: XaR8, YaR8, XbR8, YbR8, XcR8, YcR8, XdR8, YdR8, XeR8, YeR8 + real(8) :: da, db, dc, dd, pw + !Begin------------------------------------------------------------------- + + XaR8 = dble(Xa); YaR8 = dble(Ya); XbR8 = dble(Xb); YbR8 = dble(Yb); XcR8 = dble(Xc) + YcR8 = dble(Yc); XdR8 = dble(Xd); YdR8 = dble(Yd); XeR8 = dble(Xe); YeR8 = dble(Ye) + + pw = 2 + + !using a IDW - Inverse distance weighting method + !(Xa,Ya) = (Xex = 0, Yey =1) + !(Xb,Yb) = (Xex = 1, Yey =1) + !(Xc,Yc) = (Xex = 0, Yey =0) + !(Xd,Yd) = (Xex = 1, Yey =0) + ! + ! a_______b + ! | | + ! | | + ! c_______d + + da = sqrt((XeR8-XaR8)**2+(YeR8-YaR8)**2) + db = sqrt((XeR8-XbR8)**2+(YeR8-YbR8)**2) + dc = sqrt((XeR8-XcR8)**2+(YeR8-YcR8)**2) + dd = sqrt((XeR8-XdR8)**2+(YeR8-YdR8)**2) + + if (da == 0) then + Xex = 0 + Yey = 1 + elseif (db == 0) then + Xex = 1 + Yey = 1 + elseif (dc == 0) then + Xex = 0 + Yey = 0 + elseif (dd == 0) then + Xex = 1 + Yey = 0 + else + SumAux = ((1./da)**pw+(1./db)**pw+(1./dc)**pw+(1./dd)**pw) + Xex = ((1./db)**pw+(1./dd)**pw) / SumAux + Yey = ((1./da)**pw+(1./db)**pw) / SumAux + endif + + + if (Xex < 0. .or. Xex > 1) then + stop 'Function - RelativePosition4VertPolygonIWD - ERR10' + endif + + if (Yey < 0. .or. Yey > 1) then + stop 'Function - RelativePosition4VertPolygonIWD - ERR20' + endif + + + end subroutine RelativePosition4VertPolygonIWD + !-------------------------------------------------------------------------- + !! Public-domain function by Darel Rex Finley, 2006. @@ -13746,7 +13824,8 @@ end subroutine DischargeFluxV !>@param[in] Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold subroutine Offline_DischargeFluxU(Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold) !Arguments-------------------------------------------------------------------------- - real, dimension(:, :, :), pointer, intent(IN) :: VelFather, VelSon, AreaU + real, dimension(:, :, :), pointer, intent(IN) :: VelFather, AreaU + real, dimension(:, :, :), allocatable, intent(IN) :: VelSon real(8), dimension(:) , intent(INOUT) :: Flow integer, dimension(:, :) , allocatable, intent(IN) :: DischargeConnection real , dimension(:, :) , pointer, intent(IN) :: DecayTime @@ -13801,7 +13880,8 @@ end subroutine Offline_DischargeFluxU !>@param[in] Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold subroutine Offline_DischargeFluxV(Flow, DischargeConnection, VelFather, VelSon, AreaV, DecayTime, VelDT, CoefCold) !Arguments-------------------------------------------------------------------------- - real, dimension(:, :, :), pointer, intent(IN) :: VelFather, VelSon, AreaV + real, dimension(:, :, :), pointer, intent(IN) :: VelFather, AreaV + real, dimension(:, :, :), allocatable, intent(IN) :: VelSon real(8), dimension(:) , intent(OUT) :: Flow integer, dimension(:, :) , allocatable, intent(IN) :: DischargeConnection real , dimension(:, :) , pointer, intent(IN) :: DecayTime diff --git a/Software/MOHIDBase1/ModuleGlobalData.F90 b/Software/MOHIDBase1/ModuleGlobalData.F90 index 7d9ebb515..743740c33 100644 --- a/Software/MOHIDBase1/ModuleGlobalData.F90 +++ b/Software/MOHIDBase1/ModuleGlobalData.F90 @@ -76,10 +76,12 @@ Module ModuleGlobalData end interface SetError !Parameter----------------------------------------------------------------- - integer, parameter :: MaxModules = 98 + integer, parameter :: MaxModules = 100 -#ifdef _INCREASE_MAXINSTANCES +#if defined(_INCREASE_MAXINSTANCES) integer, parameter :: MaxInstances = 2000 +#elif defined(_INCREASE_MAXINSTANCES_EXTRA) + integer, parameter :: MaxInstances = 20000 #else integer, parameter :: MaxInstances = 500 #endif @@ -1980,6 +1982,8 @@ Module ModuleGlobalData integer, parameter :: mLitter_ = 96 integer, parameter :: mTwoWay_ = 97 integer, parameter :: mOutputGrid_ = 98 + integer, parameter :: mMeshGlue_ = 99 + integer, parameter :: mDelftFM_2_MOHID_ = 100 !Domain decomposition integer, parameter :: WestSouth = 1 @@ -2108,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_ , "OuputGrid" )/) + T_Module(mOutputGrid_ , "OutputGrid" ), T_Module(mMeshGlue_ , "MeshGlue" ), & + T_Module(mDelftFM_2_MOHID_ , "DelftFM_2_MOHID" )/) !Variables diff --git a/Software/MOHIDBase1/ModuleHDF5.F90 b/Software/MOHIDBase1/ModuleHDF5.F90 index 079c4be48..17160c5a7 100644 --- a/Software/MOHIDBase1/ModuleHDF5.F90 +++ b/Software/MOHIDBase1/ModuleHDF5.F90 @@ -7458,7 +7458,7 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn) allocate (offset_in (rank)) allocate (count_in (rank)) - offset_in(:) = 1 + offset_in(:) = 0 count_in (:) = 1 @@ -7506,7 +7506,7 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn) NumType = H5T_NATIVE_REAL - call h5dread_f (dset_id, NumType, ValueOut(1),& + call h5dread_f (dset_id, NumType, ValueOut(1:1),& dims_mem, STAT, memspace_id, space_id) if (STAT /= SUCCESS_) then write(*,*) "FileName =",trim(Me%Filename) @@ -7529,6 +7529,11 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn) if (STAT /= SUCCESS_) then stop 'CheckAllDataSets - ModuleHDF5 - ERR110' endif + + call h5sclose_f (memspace_id, STAT) + if (STAT /= SUCCESS_) then + stop 'CheckAllDataSets - ModuleHDF5 - ERR115' + endif @@ -8119,7 +8124,7 @@ subroutine LocateObjHDF5 (HDF5ID) enddo if (.not. associated(Me)) then - write(*,*) Me%FileName + write(*,*) "HDF5ID =", HDF5ID stop 'ModuleHDF5 - LocateObjHDF5 - ERR01' endif 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/ModuleBoxDif.F90 b/Software/MOHIDBase2/ModuleBoxDif.F90 index 382337ac4..65f371719 100644 --- a/Software/MOHIDBase2/ModuleBoxDif.F90 +++ b/Software/MOHIDBase2/ModuleBoxDif.F90 @@ -1013,7 +1013,11 @@ subroutine ConstructBox2D (NewBox, FirstLine, LastLine) SearchType = FromBlock, & Buffer_Line = Line, & STAT = STAT_CALL) - if(STAT_CALL .ne. SUCCESS_)stop 'ConstructBox2D - ModuleBoxDif - ERR20' + + if(STAT_CALL /= SUCCESS_) then + write(*,*) 'Line=',Line + stop 'ConstructBox2D - ModuleBoxDif - ERR20' + endif select case (Me%CoordinateType) diff --git a/Software/MOHIDBase2/ModuleField4D.F90 b/Software/MOHIDBase2/ModuleField4D.F90 index 8e73569f1..6a2a02f94 100644 --- a/Software/MOHIDBase2/ModuleField4D.F90 +++ b/Software/MOHIDBase2/ModuleField4D.F90 @@ -914,6 +914,12 @@ subroutine ReadGridFromFile() #endif endif + if (LatStag( 1, 1) /= LatStag( 1,Jmax+1) .or. & + LatStag(Imax+1, 1) /= LatStag(Imax+1,Jmax+1)) then + !Grid with some degree of distortion - not able to read a window + Me%ReadWindow = .false. + endif + irw: if (Me%ReadWindow) then @@ -976,21 +982,31 @@ subroutine ReadGridFromFile() ILB = Me%WindowLimitsJI%ILB if (ILB < 1) then + write (*,*) "ILB =", ILB + write (*,*) "Filename=", trim(Me%File%FileName) stop 'ReadGridFromFile - ModuleField4D - ERR80' endif IUB = Me%WindowLimitsJI%IUB if (IUB > imax) then + write (*,*) "IUB =", IUB + write (*,*) "imax=", imax + write (*,*) "Filename=", trim(Me%File%FileName) stop 'ReadGridFromFile - ModuleField4D - ERR90' endif JLB = Me%WindowLimitsJI%JLB if (JLB < 1) then + write (*,*) "JLB =", JLB + write (*,*) "Filename=", trim(Me%File%FileName) stop 'ReadGridFromFile - ModuleField4D - ERR100' endif JUB = Me%WindowLimitsJI%JUB if (JUB > jmax) then + write (*,*) "JUB =", JUB + write (*,*) "jmax=", jmax + write (*,*) "Filename=", trim(Me%File%FileName) stop 'ReadGridFromFile - ModuleField4D - ERR110' endif @@ -1042,7 +1058,7 @@ subroutine ReadGridFromFile() deallocate(Lon ) deallocate(LatStag) deallocate(LonStag) - if (Me%ReadWindowXY .and. Me%WindowWithData) then + if (Me%ReadWindow .and. Me%ReadWindowXY .and. Me%WindowWithData) then deallocate(LatStagW) deallocate(LonStagW) endif @@ -5739,7 +5755,9 @@ subroutine Interpolate2DCloud (PropField, X, Y, Field, NoData) endif call GetXYCellZ(Me%ObjHorizontalGrid, X(nP), Y(nP), i, j, PercI, PercJ, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'Interpolate2DCloud - ModuleValida4D - ERR20' + if (STAT_CALL /= SUCCESS_) then + stop 'Interpolate2DCloud - ModuleValida4D - ERR20' + endif if (PropField%InterpolMethod == NoInterpolation2D_) then @@ -6504,6 +6522,17 @@ subroutine Interpolater3D(PropField, Matrix3D, Depth3D, Mask3D, HorizontalGrid, call GetXYCellZ(HorizontalGrid, X(nP), Y(nP), i, j, PercI, PercJ, STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'Interpolate3DCloud - ModuleValida4D - ERR50' + + if (PercI < 0 .or. PercI >1) then + write(*,*) 'Wrong 01) then + write(*,*) 'Wrong 0 null() integer, dimension(:,:), pointer :: DefineFacesUMap => null() @@ -3089,6 +3092,7 @@ subroutine ConstructGlobalVariables integer :: ClientNumber logical :: BlockFound, ConstantSpacingX, ConstantSpacingY integer :: FirstLine, LastLine, line, i, j, ii, jj, iflag + logical :: CornersOverlap !---------------------------------------------------------------------- @@ -3448,6 +3452,75 @@ subroutine ConstructGlobalVariables end do end do + CornersOverlap = .false. + + do i = Me%GlobalWorkSize%ILB, Me%GlobalWorkSize%IUB + do j = Me%GlobalWorkSize%JLB, Me%GlobalWorkSize%JUB + + + if (Me%DDecomp%MasterOrSlave) then + if (i>= Me%DDecomp%HaloMap%ILB .and. & + i<= Me%DDecomp%HaloMap%IUB+1) then + ii = i - Me%DDecomp%HaloMap%ILB + 1 + else + cycle + endif + else + ii = i + endif + + if (Me%DDecomp%MasterOrSlave) then + if (j>= Me%DDecomp%HaloMap%JLB .and. & + j<= Me%DDecomp%HaloMap%JUB+1) then + jj = j - Me%DDecomp%HaloMap%JLB + 1 + else + cycle + endif + else + jj = j + endif + + if (Me%XX_IE(ii, jj) < FillValueReal/2) then + Me%GhostCorners = .true. + Cycle + endif + + line = FirstLine + (i-1)*(Me%GlobalWorkSize%JUB+1) + j + + if (Me%XX_IE(ii, jj) == Me%XX_IE(ii+1, jj) .and. & + Me%YY_IE(ii, jj) == Me%YY_IE(ii+1, jj)) then + + write(*,*) 'Corners (i,j) overlap',i,j,'vs',i+1,j + write(*,*) 'Line =', Line + CornersOverlap = .true. + + endif + + if (Me%XX_IE(ii, jj) == Me%XX_IE(ii, jj+1) .and. & + Me%YY_IE(ii, jj) == Me%YY_IE(ii, jj+1)) then + + write(*,*) 'Corners (i,j) overlap',i,j,'vs',i,j+1 + write(*,*) 'Line =', Line + CornersOverlap = .true. + + endif + + if (Me%XX_IE(ii, jj) == Me%XX_IE(ii+1, jj+1) .and. & + Me%YY_IE(ii, jj) == Me%YY_IE(ii+1, jj+1)) then + + write(*,*) 'Corners (i,j) overlap',i,j,'vs',i+1,j+1 + write(*,*) 'Line =', Line + CornersOverlap = .true. + + endif + + end do + end do + + if (CornersOverlap) then + stop 'ConstructGlobalVariables - HorizontalGrid - ERR235' + endif + if (Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then @@ -3635,6 +3708,13 @@ subroutine ConstructGlobalVariables Me%YY(ii) = Aux + if (ii >1) then + if (Me%YY(ii) < Me%YY(ii-1)) then + write(*,*) 'Block / need to be written in an ascending way' + stop 'ConstructGlobalVariables - HorizontalGrid - ERR305' + endif + endif + end do end if @@ -3729,6 +3809,14 @@ subroutine ConstructGlobalVariables Me%XX(jj) = Aux + + if (jj >1) then + if (Me%XX(jj) < Me%XX(jj-1)) then + write(*,*) 'Block / need to be written in an ascending way' + stop 'ConstructGlobalVariables - HorizontalGrid - ERR375' + endif + endif + end do end if @@ -4092,7 +4180,6 @@ subroutine ConstructGlobalVariablesV2() Me%Distortion = .false. Me%RegularRotation = .false. - Me%CornersXYInput = .false. Me%CornersXYInput = .true. @@ -4347,7 +4434,9 @@ subroutine Mercator() !real(8) :: EarthRadius, Rad_Lat, CosenLat real :: MaxLon, MinLon, MaxLat, MinLat, MinX, MinY + real(8) :: X1_r8, Y1_r8, X2_r8, Y2_r8 + !Begin----------------------------------------------------------------- !Worksize ILB = Me%WorkSize%ILB @@ -4567,10 +4656,15 @@ subroutine Mercator() ! CosenLat = cos(Rad_Lat) ! XX_IE(i, j) = CosenLat * EarthRadius * (Me%LongitudeConn(i, j) - Me%Longitude) * radians ! YY_IE(i, j) = EarthRadius * (Me%LatitudeConn (i, j) - Me%Latitude ) * radians - call FromSphericalToCart(Lat = Me%LatitudeConn (i, j), & - Lon = Me%LongitudeConn(i, j), & - X = XX_IE(i, j), & - Y = YY_IE(i, j)) + + X1_r8 = Me%LongitudeConn(i, j) + Y1_r8 = Me%LatitudeConn (i, j) + + call FromSphericalToCart(Lat = Y1_r8, Lon = X1_r8, X = X2_r8, Y = Y2_r8) + + XX_IE(i, j) = X2_r8 + YY_IE(i, j) = Y2_r8 + enddo enddo @@ -4849,6 +4943,7 @@ subroutine FromSphericalToCart(Lat, Lon, X, Y) !Local----------------------------------------------------------------- !real(8) :: radians, EarthRadius, Rad_Lat, CosenLat + real(8) :: LatRef, LonRef !Begin----------------------------------------------------------------- @@ -4860,7 +4955,10 @@ subroutine FromSphericalToCart(Lat, Lon, X, Y) !X = CosenLat * EarthRadius * (Lon - Me%Longitude) * radians !Y = EarthRadius * (Lat - Me%Latitude ) * radians - call SphericalToCart(Lat, Lon, X, Y, Me%Longitude, Me%Latitude) + LonRef = Me%Longitude + LatRef = Me%Latitude + + call SphericalToCart(Lat, Lon, X, Y, LonRef, LatRef) end subroutine FromSphericalToCart @@ -10152,7 +10250,8 @@ subroutine GetXYCellZ(HorizontalGridID, XPoint, YPoint, I, J, PercI, PercJ, !stop 'GetXYCellZ - ModuleHorizontalGrid - ERR10' else if (present(PercI) .and. present(PercJ)) then - call RelativePosition4VertPolygon(Xa = XX2D(I+1, J ), Ya = YY2D(I+1, J ), & + call RelativePosition4VertPolygon( & + Xa = XX2D(I+1, J ), Ya = YY2D(I+1, J ), & Xb = XX2D(I+1, J+1), Yb = YY2D(I+1, J+1), & Xc = XX2D(I , J ), Yc = YY2D(I , J ), & Xd = XX2D(I , J+1), Yd = YY2D(I , J+1), & @@ -10323,8 +10422,8 @@ subroutine GetXYCellZ_ThreadSafe(HorizontalGridID, XPoint, YPoint, I, J, PercI, else if (present(PercI) .and. present(PercJ)) then - ! - call RelativePosition4VertPolygon(Xa = XX2D(I+1, J ), Ya = YY2D(I+1, J ), & + call RelativePosition4VertPolygon( & + Xa = XX2D(I+1, J ), Ya = YY2D(I+1, J ), & Xb = XX2D(I+1, J+1), Yb = YY2D(I+1, J+1), & Xc = XX2D(I , J ), Yc = YY2D(I , J ), & Xd = XX2D(I , J+1), Yd = YY2D(I , J+1), & @@ -10449,6 +10548,7 @@ subroutine GetCellZ_XY(HorizontalGridID, I, J, PercI, PercJ, XPoint, YPoint, & integer :: ready_ real :: xac, yac, xbd, ybd integer :: Referential_ + real(8) :: Xin_r8, Yin_r8, XPoint_r8, YPoint_r8 !------------------------------------------------------------------------ STAT_ = UNKNOWN_ @@ -10468,22 +10568,31 @@ subroutine GetCellZ_XY(HorizontalGridID, I, J, PercI, PercJ, XPoint, YPoint, & endif + if (Me%CoordType == SIMPLE_GEOG_ .and. present(Xin) .and. present(Yin)) then + Xin_r8 = Xin + Yin_r8 = Yin + if (Referential_ == GridCoord_) then - call FromCartToSpherical(Lat = YPoint, & - Lon = XPoint, & - X = Xin, & - Y = Yin) + call FromCartToSpherical(Lat = YPoint_r8, & + Lon = XPoint_r8, & + X = Xin_r8, & + Y = Yin_r8) else - call FromSphericalToCart(Lat = Yin, & - Lon = Xin, & - X = XPoint, & - Y = YPoint) + call FromSphericalToCart(Lat = Yin_r8, & + Lon = Xin_r8, & + X = XPoint_r8, & + Y = YPoint_r8) endif + XPoint = XPoint_r8 + YPoint = YPoint_r8 + + + else XX2D => Me%XX_IE YY2D => Me%YY_IE @@ -11703,7 +11812,6 @@ integer function GetDDecompMPI_ID(HorizontalGridID, STAT) end function GetDDecompMPI_ID !-------------------------------------------------------------------------- - !-------------------------------------------------------------------------- logical function GetDDecompON(HorizontalGridID, STAT) @@ -11746,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) @@ -12257,10 +12405,10 @@ real function InterpolXYPoint(HorizontalGridID, Field2DFather, YYLower = YYFather(Ilower) YYUpper = YYFather(Iupper) - XPosition = XInput - YPosition = YInput + XPosition = XInput - Me%Xorig + YPosition = YInput - Me%Yorig - call RODAXY(-Me%Xorig, -Me%Yorig, -Me%Grid_Angle, XPosition, YPosition) + call RODAXY(0., 0., -Me%Grid_Angle, XPosition, YPosition) else @@ -13961,6 +14109,329 @@ end subroutine WriteHorizontalGrid_UV !-------------------------------------------------------------------------- + !---------------------------------------------------------------------------- + + subroutine WriteGridData_ASCII (FileName, & + COMENT1, COMENT2, & + HorizontalGridID, & + FillValue, Overwrite, & + DistortionYes, & + STAT) + + !Arguments--------------------------------------------------------------- + character(LEN = *), intent(IN) :: FileName + character(LEN = *), intent(IN) :: COMENT1 + character(LEN = *), intent(IN) :: COMENT2 + integer :: HorizontalGridID + real :: FillValue + logical, intent(IN) :: Overwrite + logical, optional, intent(IN) :: DistortionYes + integer, optional, intent(OUT) :: STAT + + !Local------------------------------------------------------------------- + real :: Xorig + real :: Yorig + real :: GRID_ANGLE + integer :: Zone + real :: Longitude + real :: Latitude + integer :: CoordType + + real, pointer, dimension(: ) :: XX, YY + real, pointer, dimension(:,:) :: XX_IE, YY_IE + real, pointer, dimension(:,:) :: LatConn, LongConn + + real, pointer, dimension(: ) :: XX_Out, YY_Out + real, pointer, dimension(:,:) :: XX_IE_Out, YY_IE_Out + real, pointer, dimension(:,:) :: LatConn_Out, LongConn_Out + + type (T_Size2D) :: WorkSize, Global + + integer :: Unit + integer :: STAT_CALL + integer :: UTM_, MIL_PORT_ + integer :: PORTUGUESE_UTM_ZONE_, SIMPLE_GEOG_ + logical :: exist, DistortionYes_ + integer :: i, j + + logical :: ON, Master, MasterOrSlave + + !------------------------------------------------------------------------ + + + !Gets Origin + call GetGridOrigin(HorizontalGridID, Xorig, Yorig, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR20' + + !Gets GridAngle + call GetGridAngle(HorizontalGridID, GRID_ANGLE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR30' + + !Gets GridZone + call GetGridZone (HorizontalGridID, Zone, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR40' + + !Gets Lat/lon + call GetLatitudeLongitude (HorizontalGridID, Latitude, Longitude, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR50' + + !Gets Lat/lon + call GetGridCoordType (HorizontalGridID, CoordType, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR60' + + !Gets XX/ YY and XX_IE, YY_IE + call GetHorizontalGrid(HorizontalGridID, XX = XX, YY = YY, XX_IE = XX_IE, YY_IE = YY_IE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR70' + + call GetGridLatitudeLongitude(HorizontalGridID, & + GridLatitudeConn = LatConn, & + GridLongitudeConn = LongConn, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR80' + + if (present(DistortionYes)) then + + DistortionYes_ = DistortionYes + + else + + call GetCheckDistortion(HorizontalGridID, DistortionYes_, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR90' + + endif + + !Gets WorkSize + call GetHorizontalGridSize(HorizontalGridID, WorkSize = WorkSize, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR100' + + call GetDDecompParameters & + (HorizontalGridID = HorizontalGridID, & + ON = ON, & + Master = Master, & + MasterOrSlave = MasterOrSlave, & + Global = Global, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR110' + + + +if3: if (.not. MasterOrSlave .or. Master) then + + call UnitsManager(Unit, OPEN_FILE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR120' + + if(.not. Overwrite)then + !Verifies if file exists + inquire (file = FileName, exist = exist) + if (exist) then + write(*,* )'Cannot write Grid Data. File Exists already!' + write(*, *)trim(adjustl(FileName)) + stop 'WriteGridData_ASCII - ModuleGridData - ERR130' + endif + + !Opens file + open (Unit, file = FileName, status = 'new', IOSTAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR140' + else + !Opens file + open (Unit, file = FileName, status = 'replace', IOSTAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + write(*,*) 'Was not possible to write to file ', FileName + write(*,*) 'OPEN returned IOSTAT = ', STAT_CALL + stop 'WriteGridData_ASCII - ModuleGridData - ERR150' + endif + end if + + endif if3 + + +if1: if (MasterOrSlave) then + + WorkSize = Global + +#ifdef _USE_MPI + + if (DistortionYes_) then + if (CoordType .EQ. SIMPLE_GEOG_) then + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In2D = LongConn, & + Out2D = LongConn_Out, & + dj = 1, & + di = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR210' + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In2D = LatConn, & + Out2D = LatConn_Out, & + dj = 1, & + di = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR220' + else + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In2D = XX_IE, & + Out2D = XX_IE_Out, & + dj = 1, & + di = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR230' + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In2D = YY_IE, & + Out2D = YY_IE_Out, & + dj = 1, & + di = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR240' + end if + else + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In1D = XX, & + Out1D = XX_Out, & + Type1DIJ = Type1DJ, & + dj = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR250' + + call JoinGridData(HorizontalGridID = HorizontalGridID, & + In1D = YY, & + Out1D = YY_Out, & + Type1DIJ = Type1DI, & + di = 1, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR260' + endif + +#endif _USE_MPI + + else if1 + + XX_Out => XX + YY_Out => YY + XX_IE_Out => XX_IE + YY_IE_Out => YY_IE + LatConn_Out => LatConn + LongConn_Out => LongConn + + endif if1 + +if2: if (.not. MasterOrSlave .or. Master) then + + call WriteDataLine(unit, 'COMENT1', COMENT1) + call WriteDataLine(unit, 'COMENT2', COMENT2) + write(unit, *) + write(unit, *) + + write(unit, *) 'ILB_IUB :', WorkSize%ILB, WorkSize%IUB + write(unit, *) 'JLB_JUB :', WorkSize%JLB, WorkSize%JUB + + write(unit, *) 'COORD_TIP :', CoordType + + call GetCoordTypeList(UTM = UTM_, MIL_PORT = MIL_PORT_, & + PORTUGUESE_UTM_ZONE = PORTUGUESE_UTM_ZONE_, & + SIMPLE_GEOG = SIMPLE_GEOG_) + if ((CoordType .EQ. UTM_) .OR. (CoordType .EQ. MIL_PORT_)) then + if (CoordType .EQ. MIL_PORT_) then + write(unit, *) 'ZONE :', PORTUGUESE_UTM_ZONE_ + else + write(unit, *) 'ZONE :', Zone + endif + end if + + write(unit, *) 'ORIGIN :', Xorig, Yorig + write(unit, *) 'GRID_ANGLE :', GRID_ANGLE + write(unit, *) 'LATITUDE :', Latitude + write(unit, *) 'LONGITUDE :', Longitude + write(unit, *) 'FILL_VALUE :', FillValue + write(unit, *) + write(unit, *) + + if (DistortionYes_) then + + write(unit, *) "" + + if (CoordType .EQ. SIMPLE_GEOG_) then + + do i = WorkSize%ILB, WorkSize%IUB+1 + do j = WorkSize%JLB, WorkSize%JUB+1 + write(unit, *) LongConn_Out(i,j), LatConn_Out(i,j) + end do + end do + + else + + do i = WorkSize%ILB, WorkSize%IUB+1 + do j = WorkSize%JLB, WorkSize%JUB+1 + write(unit, *) XX_IE_Out(i,j),YY_IE_Out(i,j) + end do + end do + + end if + + write(unit, *) "<"//backslash//"CornersXY>" + ! write(unit, *) "<\CornersXY>" + + + else + + write(unit, *) "" + do j = WorkSize%JLB, WorkSize%JUB+1 + write(unit, *) XX_Out(j) + end do + write(unit, *) "" + + + write(unit, *) "" + do i = WorkSize%ILB, WorkSize%IUB+1 + write(unit, *) YY_Out(i) + end do + write(unit, *) "" + + endif + + !Closes Files + call UnitsManager(Unit, CLOSE_FILE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR270' + + endif if2 + + !Ungets XX, YY + call UngetHorizontalGrid (HorizontalGridID, XX, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR280' + + call UngetHorizontalGrid (HorizontalGridID, YY, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR290' + + call UngetHorizontalGrid(HorizontalGridID, XX_IE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR300' + + call UngetHorizontalGrid(HorizontalGridID, YY_IE, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR310' + + call UngetHorizontalGrid(HorizontalGridID, LatConn, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR320' + + call UngetHorizontalGrid(HorizontalGridID, LongConn, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteGridData_ASCII - ModuleGridData - ERR330' + + +if31: if (MasterOrSlave) then + + if (associated(XX_Out )) deallocate(XX_Out ) + if (associated(YY_Out )) deallocate(YY_Out ) + if (associated(XX_IE_Out )) deallocate(XX_IE_Out ) + if (associated(YY_IE_Out )) deallocate(YY_IE_Out ) + if (associated(LatConn_Out )) deallocate(LatConn_Out ) + if (associated(LongConn_Out )) deallocate(LongConn_Out ) + + endif if31 + + if (present(STAT)) STAT = SUCCESS_ + + !------------------------------------------------------------------------ + + end subroutine WriteGridData_ASCII + + subroutine ReadHDF5HorizontalGrid () @@ -17985,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 @@ -18054,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 @@ -18087,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 @@ -18119,7 +18633,7 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, IsPointInside = .false. - !Construct 4 polygons NW, NE, SW, SE + !Construct 4 polygons SW, SE, NW, WE do pi = 1, 4 !All polygons are defined anti-Clockwise @@ -18259,6 +18773,8 @@ subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint, endif f4 endif f2 + + endif f0 if (present(CellLocated)) CellLocated = CellFound @@ -18306,6 +18822,9 @@ logical function InsideDomainPolygon(DomainPolygon, XPoint, YPoint) end function InsideDomainPolygon + + + #ifdef _OPENMI_ diff --git a/Software/MOHIDBase2/ModuleLitter.F90 b/Software/MOHIDBase2/ModuleLitter.F90 index 306b573ea..a8d3b5d4f 100644 --- a/Software/MOHIDBase2/ModuleLitter.F90 +++ b/Software/MOHIDBase2/ModuleLitter.F90 @@ -327,7 +327,7 @@ subroutine ConstructLitterLag(ObjLitterID, StartTime, ModelDomain, STAT) call AllocateInstance - Me%Files%Nomfich = "Nomfich.dat" + Me%Files%Nomfich = "nomfich.dat" allocate(VectorX(5), VectorY(5)) diff --git a/Software/MOHIDBase2/ModuleMap.F90 b/Software/MOHIDBase2/ModuleMap.F90 index 959c04f0b..17cccb110 100644 --- a/Software/MOHIDBase2/ModuleMap.F90 +++ b/Software/MOHIDBase2/ModuleMap.F90 @@ -900,7 +900,7 @@ subroutine CorrectIsolatedCells(GridDataID, HorizontalGridID, FirstIsolatedCell) write(*,*)'' if (Me%StopOnBathymetryChange) then - write(*,*)'Modify the file Nomfich.dat and Re-run the model' + write(*,*)'Modify the file nomfich.dat and Re-run the model' stop 'CorrectIsolatedCells - ModuleMap - ERR60' endif diff --git a/Software/MOHIDBase2/ModuleStatistic.F90 b/Software/MOHIDBase2/ModuleStatistic.F90 index 20620ad46..fbde05d9d 100644 --- a/Software/MOHIDBase2/ModuleStatistic.F90 +++ b/Software/MOHIDBase2/ModuleStatistic.F90 @@ -3607,6 +3607,7 @@ subroutine WriteValuesToFileHDF5 (WriteGlobal, WriteDaily, & real, dimension(:, :, :), pointer :: AuxMatrix3D real, dimension(:, : ), pointer :: AuxMatrix2D real :: Aux, d1, d2 + real :: AuxFreq real(8), dimension(:), allocatable :: P,C real(8) :: Px,Cx, sumFreq integer :: nc, n @@ -4351,7 +4352,9 @@ subroutine WriteValuesToFileHDF5 (WriteGlobal, WriteDaily, & do iClass = 1, nc - sumFreq = sumFreq + Me%Classification%Frequency(i, j, k, iClass) + AuxFreq = Me%Classification%Frequency(i, j, k, iClass) + + sumFreq = sumFreq + AuxFreq enddo @@ -4359,7 +4362,9 @@ subroutine WriteValuesToFileHDF5 (WriteGlobal, WriteDaily, & do iClass = 1, nc - P(iClass + 1) = P(iClass) + Me%Classification%Frequency(i, j, k, iClass) / sumFreq * 100. + AuxFreq = Me%Classification%Frequency(i, j, k, iClass) + + P(iClass + 1) = P(iClass) + AuxFreq / sumFreq * 100. enddo @@ -4367,7 +4372,9 @@ subroutine WriteValuesToFileHDF5 (WriteGlobal, WriteDaily, & do iClass = 1, nc - P(iClass + 1) = P(iClass) + Me%Classification%Frequency(i, j, k, iClass) * 100. + AuxFreq = Me%Classification%Frequency(i, j, k, iClass) + + P(iClass + 1) = P(iClass) + AuxFreq * 100. enddo @@ -4377,7 +4384,9 @@ subroutine WriteValuesToFileHDF5 (WriteGlobal, WriteDaily, & doClass1: do iClass = 1, nc - P(iClass + 1) = P(iClass) + Me%Classification%Frequency(i, j, k, iClass) * 100. + AuxFreq = Me%Classification%Frequency(i, j, k, iClass) + + P(iClass + 1) = P(iClass) + AuxFreq * 100. enddo doClass1 diff --git a/Software/MOHIDBase2/ModuleTwoWay.F90 b/Software/MOHIDBase2/ModuleTwoWay.F90 index 991b38408..87fc82e39 100644 --- a/Software/MOHIDBase2/ModuleTwoWay.F90 +++ b/Software/MOHIDBase2/ModuleTwoWay.F90 @@ -1520,7 +1520,8 @@ subroutine Offline_Upscaling_Discharge (TwoWayID, Flow, VelFather, VelSon, Decay !Arguments-------------------------------------------------------------- integer, intent(IN ) :: TwoWayID real, dimension(:,:,:), pointer, intent(INOUT) :: Flow - real, dimension(:,:,:), pointer, intent(IN ) :: VelFather, VelSon + real, dimension(:,:,:), pointer, intent(IN ) :: VelFather + real, dimension(:,:,:), allocatable, intent(IN ) :: VelSon real, dimension(:,: ), pointer, intent(IN ) :: DecayTime real , intent(IN ) :: CoefCold, VelDT integer , intent(IN ) :: VelID diff --git a/Software/MOHIDWater/ModuleHydrodynamic.F90 b/Software/MOHIDWater/ModuleHydrodynamic.F90 index b8c19c404..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, & @@ -1641,6 +1642,7 @@ Module ModuleHydrodynamic logical :: AssimilaOneField = .false. logical :: MatrixesOutputOpt = .false. + logical :: OperationalDefault = .false. end type T_HydroOptions @@ -1794,6 +1796,8 @@ Module ModuleHydrodynamic character(len=PathLength) :: FileName = null_str !Instance of ModuleGridData integer :: ObjGridData = 0 + !Wave generate + logical :: WaveGenerate = .false. ! FOCAL DEPTH, MEASURED FROM MEAN EARTH SURFACE TO THE TOP EDGE OF FAULT PLANE real :: HH = null_real ! LENGTH OF THE FAULT PLANE @@ -1821,7 +1825,8 @@ Module ModuleHydrodynamic private :: T_Tsunami type T_Tsunami logical :: ON = .true. - type (T_Fault) :: Fault + type (T_Fault), dimension(:), pointer :: Fault + integer :: Fault_Number = null_int end type T_Tsunami private :: T_Relaxation @@ -2370,8 +2375,8 @@ Subroutine Construct_Hydrodynamic(DischargesID, AssimilationID) !Construct the Time Serie Obj if (Me%OutPut%TimeSerieON) call Construct_Time_Serie - if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & - Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON.or. & + if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & + Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON.or. & Me%OutW%OutPutWindowsON) then call ConstructMatrixesOutput endif @@ -4329,6 +4334,8 @@ subroutine ConstructTsunami !Local----------------------------------------------------------------- integer :: STAT_CALL, iflag + integer :: ifl, ClientNumber + logical :: AtLeastOneBlock, BlockFound !---------------------------------------------------------------------- @@ -4343,226 +4350,302 @@ subroutine ConstructTsunami i1: if (Me%Tsunami%ON) then + call GetData(Me%Tsunami%Fault_Number, & + Me%ObjEnterData, iflag, & + Keyword = 'FAULT_NUMBER', & + Default = 1, & + SearchType = FromFile, & + ClientModule = 'ModuleHydrodynamic', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR20' + + allocate(Me%Tsunami%Fault(Me%Tsunami%Fault_Number)) + + AtLeastOneBlock = .false. + +do1 : do ifl =1,Me%Tsunami%Fault_Number + call ExtractBlockFromBuffer(Me%ObjEnterData, & + ClientNumber = ClientNumber, & + block_begin = "", & + block_end = "", & + BlockFound = BlockFound, & + STAT = STAT_CALL) + cd1 : if (STAT_CALL .EQ. SUCCESS_ ) then + cd2 : if (BlockFound) then + + AtLeastOneBlock = .true. + + call TsunamiKeywords(FromBlock, ifl) + + if (ifl==Me%Tsunami%Fault_Number) then + exit + endif + + else cd2 + + if (ifl==1) then + call TsunamiKeywords(FromFile, ifl) + exit + else + stop 'ConstructTsunami - ModuleHydrodynamic - ERR030' + endif + + endif cd2 + + else cd1 + stop 'ConstructTsunami - ModuleHydrodynamic - ERR040' + endif cd1 + + enddo do1 + + if (AtLeastOneBlock) then + call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR050' + endif + + + endif i1 + + +i2: if (Me%Tsunami%ON) then + + do ifl = 1, Me%Tsunami%Fault_Number + + if (Me%Tsunami%Fault(ifl)%T0 <= Me%BeginTime .and. & + .not. Me%Tsunami%Fault(ifl)%WaveGenerate) then + + call GetOpenPoints3D(Me%ObjMap, Me%External_Var%OpenPoints3D, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR060' + + call GetGridLatitudeLongitude(Me%ObjHorizontalGrid, & + GridLatitude = Me%External_Var%LatitudeZ,& + GridLongitude = Me%External_Var%LongitudeZ,& + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) & + stop 'ConstructTsunami - ModuleHydrodynamic - ERR070' + + call WaterLevel_Tsunami(ifl) + + Me%Tsunami%Fault(ifl)%WaveGenerate = .true. + + call UnGetMap(Me%ObjMap, Me%External_Var%OpenPoints3D, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR080' + + call UnGetHorizontalGrid(Me%ObjHorizontalGrid, & + Me%External_Var%LatitudeZ, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) & + stop 'ConstructTsunami - ModuleHydrodynamic - ERR090' + + call UnGetHorizontalGrid(Me%ObjHorizontalGrid, & + Me%External_Var%LongitudeZ, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) & + stop 'ConstructTsunami - ModuleHydrodynamic - ERR100' + + endif + + enddo + endif i2 + + + end subroutine ConstructTsunami + + !-------------------------------------------------------------------------- + + subroutine TsunamiKeywords(FromWhere, ifl) + + !Arguments------------------------------------------------------------- + integer :: FromWhere, ifl + + !Local----------------------------------------------------------------- + integer :: STAT_CALL, iflag + + !---------------------------------------------------------------------- + + ! TIME WHEN THE RUTPURE STARTS [YYYY MM DD HH MM SS] - call GetData(Me%Tsunami%Fault%T0, & + call GetData(Me%Tsunami%Fault(ifl)%T0, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_RUTPURE_START_TIME', & default = Me%BeginTime, & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR20' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR20' ! Fault input method (0 - Okada 1985 Model, 1 - input grid data file) - call GetData(Me%Tsunami%Fault%InputMethod, & + call GetData(Me%Tsunami%Fault(ifl)%InputMethod, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_INPUT_METHOD', & default = FaultOkada1985_, & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR30' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR30' - if (Me%Tsunami%Fault%InputMethod /= FaultFile_ .and. & - Me%Tsunami%Fault%InputMethod /= FaultOkada1985_) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR40' + if (Me%Tsunami%Fault(ifl)%InputMethod /= FaultFile_ .and. & + Me%Tsunami%Fault(ifl)%InputMethod /= FaultOkada1985_) then + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR40' endif ! AMPLIFICATION OF THE TSUNAMI - call GetData(Me%Tsunami%Fault%Amplification, & + call GetData(Me%Tsunami%Fault(ifl)%Amplification, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_AMPLIFICATION', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & default = 1., & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR245' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR245' -i3: if (Me%Tsunami%Fault%InputMethod == FaultFile_ ) then +i3: if (Me%Tsunami%Fault(ifl)%InputMethod == FaultFile_ ) then ! Fault filemane of the input grid data file - call GetData(Me%Tsunami%Fault%FileName, & + call GetData(Me%Tsunami%Fault(ifl)%FileName, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_FILENAME', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR50' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR50' if (iflag == 0) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR60' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR60' endif - call ConstructGridData(GridDataID = Me%Tsunami%Fault%ObjGridData, & + call ConstructGridData(GridDataID = Me%Tsunami%Fault(ifl)%ObjGridData,& HorizontalGridID = Me%ObjHorizontalGrid, & - FileName = Me%Tsunami%Fault%FileName, & + FileName = Me%Tsunami%Fault(ifl)%FileName, & DefaultValue = 0., & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR70' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR70' - elseif (Me%Tsunami%Fault%InputMethod == FaultOkada1985_) then i3 + elseif (Me%Tsunami%Fault(ifl)%InputMethod == FaultOkada1985_) then i3 ! FOCAL DEPTH, MEASURED FROM MEAN EARTH SURFACE TO THE TOP EDGE OF FAULT PLANE [m] - call GetData(Me%Tsunami%Fault%HH, & + call GetData(Me%Tsunami%Fault(ifl)%HH, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_FOCAL_DEPTH', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR80' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR80' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR90' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR90' endif ! LENGTH OF THE FAULT PLANE [m] - call GetData(Me%Tsunami%Fault%L, & + call GetData(Me%Tsunami%Fault(ifl)%L, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_LENGTH', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR100' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR100' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR110' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR110' endif ! WIDTH OF THE FAULT PLANE [m] - call GetData(Me%Tsunami%Fault%W, & + call GetData(Me%Tsunami%Fault(ifl)%W, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_WIDTH', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR120' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR120' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR130' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR130' endif ! DISLOCATION [m] - call GetData(Me%Tsunami%Fault%D, & + call GetData(Me%Tsunami%Fault(ifl)%D, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_DISLOCATION', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR50' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR50' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR140' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR140' endif ! (=THETA) STRIKE DIRECTION [º] - call GetData(Me%Tsunami%Fault%TH, & + call GetData(Me%Tsunami%Fault(ifl)%TH, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_STRIKE_DIRECTION', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR150' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR150' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR160' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR160' endif ! (=DELTA) DIP ANGLE [º] - call GetData(Me%Tsunami%Fault%DL, & + call GetData(Me%Tsunami%Fault(ifl)%DL, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_DIP_ANGLE', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR170' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR170' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR180' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR180' endif ! (=LAMDA) SLIP ANGLE [º] - call GetData(Me%Tsunami%Fault%RD, & + call GetData(Me%Tsunami%Fault(ifl)%RD, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_SLIP_ANGLE', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR190' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR190' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR200' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR200' endif ! EPICENTER (LATITUDE)[º] - call GetData(Me%Tsunami%Fault%Y0, & + call GetData(Me%Tsunami%Fault(ifl)%Y0, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_EPICENTER_Y', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR210' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR210' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR220' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR220' endif ! EPICENTER (LONGITUDE)[º] - call GetData(Me%Tsunami%Fault%X0, & + call GetData(Me%Tsunami%Fault(ifl)%X0, & Me%ObjEnterData, iflag, & Keyword = 'FAULT_EPICENTER_X', & - SearchType = FromFile, & + SearchType = FromWhere, & ClientModule = 'ModuleHydrodynamic', & STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR230' + if (STAT_CALL /= SUCCESS_) stop 'TsunamiKeywords - ModuleHydrodynamic - ERR230' if (iflag /= 1) then - stop 'ConstructTsunami - ModuleHydrodynamic - ERR240' + stop 'TsunamiKeywords - ModuleHydrodynamic - ERR240' endif endif i3 - endif i1 - - -i2: if (Me%Tsunami%ON) then - if (Me%Tsunami%Fault%T0 <= Me%BeginTime) then - - call GetOpenPoints3D(Me%ObjMap, Me%External_Var%OpenPoints3D, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR250' - call GetGridLatitudeLongitude(Me%ObjHorizontalGrid, & - GridLatitude = Me%External_Var%LatitudeZ,& - GridLongitude = Me%External_Var%LongitudeZ,& - STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) & - stop 'Subroutine ReadLock_ModuleHorizontalGrid - ModuleHydrodynamic. ERR260.' - - call WaterLevel_Tsunami - - Me%Tsunami%ON = .false. - - call UnGetMap(Me%ObjMap, Me%External_Var%OpenPoints3D, STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructTsunami - ModuleHydrodynamic - ERR270' - - call UnGetHorizontalGrid(Me%ObjHorizontalGrid, & - Me%External_Var%LatitudeZ, & - STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) & - stop 'Subroutine ConstructTsunami - ModuleHydrodynamic. ERR280.' - - call UnGetHorizontalGrid(Me%ObjHorizontalGrid, & - Me%External_Var%LongitudeZ, & - STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) & - stop 'Subroutine ConstructTsunami - ModuleHydrodynamic. ERR290.' - - endif - endif i2 - - - end subroutine ConstructTsunami + end subroutine TsunamiKeywords !-------------------------------------------------------------------------- @@ -8734,13 +8817,12 @@ subroutine OperationalModelDefaultOptions integer :: iflag, FromFile integer :: STAT_CALL integer :: i, j - logical :: OperationalModel !---------------------------------------------------------------------- call GetExtractType(FromFile = FromFile) - call GetData(OperationalModel, & + call GetData(Me%ComputeOptions%OperationalDefault, & Me%ObjEnterData, iflag, & keyword = 'OPERATIONAL_MODEL_DEFAULT', & default = .false., & @@ -8754,7 +8836,7 @@ subroutine OperationalModelDefaultOptions Me%ComputeOptions%AssimilaOneField = .false. - if (OperationalModel) then + if (Me%ComputeOptions%OperationalDefault) then call GetData(Me%ComputeOptions%AssimilaOneField, & Me%ObjEnterData, iflag, & @@ -8804,8 +8886,14 @@ subroutine OperationalModelDefaultOptions !RADIATION : 2 Me%ComputeOptions%BarotropicRadia = FlatherLocalSolution_ + + if (Me%ComputeOptions%InvertBarometer) then + !LOCAL_SOLUTION : 8 + Me%ComputeOptions%LocalSolution = AssimilaGauge_ + else !LOCAL_SOLUTION : 3 Me%ComputeOptions%LocalSolution = AssimilationField_ + endif Me%ComputeOptions%MinLeavingVelocity = 1e-6 Me%ComputeOptions%MinLeavingComponent = 1e-3 @@ -10311,6 +10399,8 @@ subroutine Construct_Sub_Modules(DischargesID, AssimilationID) if (STAT_CALL /= SUCCESS_) stop 'Construct_Sub_Modules - ModuleHydrodynamic - ERR70' if (IgnoreOK) then + write(*,*) 'I , J =', Id, Jd + write(*,*) 'STAT_CALL =', STAT_CALL write(*,*) 'Discharge outside the domain - ',trim(DischargeName),' - ',trim(Me%ModelName) cycle else @@ -12255,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 @@ -12297,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, & @@ -12309,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 @@ -12332,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) @@ -12358,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 @@ -12945,6 +13052,7 @@ Subroutine ConstructRelaxation integer :: ILB, IUB, JLB, JUB, KLB, KUB integer :: ILBWork, IUBWork, JLBWork, JUBWork, KLBWork, KUBWork integer :: STATUS, iflag,fromfile + logical :: LogDefault !Begin -------------------------------------------------------------------- IUB = Me%Size%IUB; JUB = Me%Size%JUB; KUB = Me%Size%KUB @@ -13142,10 +13250,17 @@ Subroutine ConstructRelaxation Me%Relaxation%ReferenceVelocity /= BaroclVel_) & call SetError(FATAL_, KEYWORD_, 'ConstructRelaxation - Hydrodynamic - ERR31.') + + if (Me%ComputeOptions%OperationalDefault) then + LogDefault = Me%Relaxation%Force + else + LogDefault = .false. + endif + call GetData(Me%Relaxation%Force, & Me%ObjEnterData, iflag, & keyword = 'BRFORCE', & - default = .false., & + default = LogDefault, & SearchType = FromFile, & ClientModule ='ModuleHydrodynamic', & STAT = status) @@ -24671,6 +24786,7 @@ Subroutine Compute_WaterLevel integer :: IUB, ILB, JUB, JLB integer :: IJmin, IJmax, JImin, JImax integer :: di, dj, DirectionXY + integer :: ifl !Begin---------------------------------------------------------------- @@ -24760,10 +24876,15 @@ Subroutine Compute_WaterLevel endif if (Me%Tsunami%ON) then - if (Me%Tsunami%Fault%T0 <= Me%CurrentTime) then - call WaterLevel_Tsunami - Me%Tsunami%ON = .false. + do ifl = 1, Me%Tsunami%Fault_Number + + if (Me%Tsunami%Fault(ifl)%T0 <= Me%CurrentTime .and. & + .not. Me%Tsunami%Fault(ifl)%WaveGenerate) then + call WaterLevel_Tsunami(ifl) +! Me%Tsunami%ON = .false. + Me%Tsunami%Fault(ifl)%WaveGenerate = .true. endif + enddo endif if (Me%Relaxation%WaterLevel) & @@ -24968,7 +25089,10 @@ end subroutine STEREO_PROJECTION ! #. UPDATED ON FEB 16 2009 (XIAOMING WANG, GNS) ! 1. ADD AN OPTION TO SELECT THE FOCUS LOCATION !---------------------------------------------------------------------- - subroutine FaultOkada1985 + subroutine FaultOkada1985(ifl) + + !Arguments---------------------------------------------------------------------- + integer :: ifl !Local-------------------------------------------------------------------------- @@ -24994,22 +25118,22 @@ subroutine FaultOkada1985 KUB = Me%WorkSize%KUB - ANG_L = RAD_DEG*Me%Tsunami%Fault%DL - ANG_R = RAD_DEG*Me%Tsunami%Fault%RD - ANG_T = RAD_DEG*Me%Tsunami%Fault%TH - HALFL = 0.5*Me%Tsunami%Fault%L + ANG_L = RAD_DEG*Me%Tsunami%Fault(ifl)%DL + ANG_R = RAD_DEG*Me%Tsunami%Fault(ifl)%RD + ANG_T = RAD_DEG*Me%Tsunami%Fault(ifl)%TH + HALFL = 0.5*Me%Tsunami%Fault(ifl)%L !.....CALCULATE FOCAL DEPTH USED FOR OKADA'S MODEL - HH = Me%Tsunami%Fault%HH+0.5*Me%Tsunami%Fault%W*SIN(ANG_L) + HH = Me%Tsunami%Fault(ifl)%HH+0.5*Me%Tsunami%Fault(ifl)%W*SIN(ANG_L) !.....DISPLACEMENT DUE TO DIFFERENT EPICENTER DEFINITION ! EPICENTER IS DEFINED AT THE CENTER OF FAULT PLANE - DEL_X = 0.5*Me%Tsunami%Fault%W*COS(ANG_L)*COS(ANG_T) - DEL_Y = 0.5*Me%Tsunami%Fault%W*COS(ANG_L)*SIN(ANG_T) + DEL_X = 0.5*Me%Tsunami%Fault(ifl)%W*COS(ANG_L)*COS(ANG_T) + DEL_Y = 0.5*Me%Tsunami%Fault(ifl)%W*COS(ANG_L)*SIN(ANG_T) H1 = HH/SIN(ANG_L) - H2 = HH/SIN(ANG_L)+Me%Tsunami%Fault%W + H2 = HH/SIN(ANG_L)+Me%Tsunami%Fault(ifl)%W - DS = Me%Tsunami%Fault%D*COS(ANG_R) - DD = Me%Tsunami%Fault%D*SIN(ANG_R) + DS = Me%Tsunami%Fault(ifl)%D*COS(ANG_R) + DD = Me%Tsunami%Fault(ifl)%D*SIN(ANG_R) SN = SIN(ANG_L) CS = COS(ANG_L) @@ -25017,8 +25141,8 @@ subroutine FaultOkada1985 X_SHIFT = 0.0 Y_SHIFT = 0.0 - LON0 = Me%Tsunami%Fault%X0 - LAT0 = Me%Tsunami%Fault%Y0 + LON0 = Me%Tsunami%Fault(ifl)%X0 + LAT0 = Me%Tsunami%Fault(ifl)%Y0 do j = JLB, JUB do i = ILB, IUB @@ -25040,14 +25164,14 @@ subroutine FaultOkada1985 P = X2*CS+HH*SN CALL Strike_Slip (X2,X1+HALFL,P ,ANG_L,HH,F1) - CALL Strike_Slip (X2,X1+HALFL,P-Me%Tsunami%Fault%W,ANG_L,HH,F2) + CALL Strike_Slip (X2,X1+HALFL,P-Me%Tsunami%Fault(ifl)%W,ANG_L,HH,F2) CALL Strike_Slip (X2,X1-HALFL,P ,ANG_L,HH,F3) - CALL Strike_Slip (X2,X1-HALFL,P-Me%Tsunami%Fault%W,ANG_L,HH,F4) + CALL Strike_Slip (X2,X1-HALFL,P-Me%Tsunami%Fault(ifl)%W,ANG_L,HH,F4) CALL Dip_Slip (X2,X1+HALFL,P ,ANG_L,HH,G1) - CALL Dip_Slip (X2,X1+HALFL,P-Me%Tsunami%Fault%W,ANG_L,HH,G2) + CALL Dip_Slip (X2,X1+HALFL,P-Me%Tsunami%Fault(ifl)%W,ANG_L,HH,G2) CALL Dip_Slip (X2,X1-HALFL,P ,ANG_L,HH,G3) - CALL Dip_Slip (X2,X1-HALFL,P-Me%Tsunami%Fault%W,ANG_L,HH,G4) + CALL Dip_Slip (X2,X1-HALFL,P-Me%Tsunami%Fault(ifl)%W,ANG_L,HH,G4) US = (F1-F2-F3+F4)*DS UD = (G1-G2-G3+G4)*DD @@ -25066,7 +25190,7 @@ subroutine FaultOkada1985 if (abs(UD) < 1e-5) UD = 0. !Water level actualization - Me%WaterLevel%New(i, j) = Me%WaterLevel%New(i, j) + (US + UD) * Me%Tsunami%Fault%Amplification + Me%WaterLevel%New(i, j) = Me%WaterLevel%New(i, j) + (US + UD) * Me%Tsunami%Fault(ifl)%Amplification endif enddo @@ -25077,7 +25201,10 @@ subroutine FaultOkada1985 end subroutine FaultOkada1985 !---------------------------------------------------------------------- - subroutine FaultInputFile + subroutine FaultInputFile (ifl) + + !Arguments---------------------------------------------------------------------- + integer :: ifl !Local-------------------------------------------------------------------------- real, dimension(:,:), pointer :: Tsunami2D @@ -25091,7 +25218,7 @@ subroutine FaultInputFile JUB = Me%WorkSize%JUB KUB = Me%WorkSize%KUB - call GetGridData(GridDataID = Me%Tsunami%Fault%ObjGridData, & + call GetGridData(GridDataID = Me%Tsunami%Fault(ifl)%ObjGridData, & GridData2D = Tsunami2D, & STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'FaultInputFile - ModuleHydrodynamic - ERR10' @@ -25100,17 +25227,17 @@ subroutine FaultInputFile do i = ILB, IUB if (Me%External_Var%OpenPoints3D(i, j, KUB) == OpenPoint) then !Water level actualization - Me%WaterLevel%New(i, j) = Me%WaterLevel%New(i, j) + Tsunami2D(i, j) * Me%Tsunami%Fault%Amplification + Me%WaterLevel%New(i, j) = Me%WaterLevel%New(i, j) + Tsunami2D(i, j) * Me%Tsunami%Fault(ifl)%Amplification endif enddo enddo - call UnGetGridData(GridDataID = Me%Tsunami%Fault%ObjGridData, & + call UnGetGridData(GridDataID = Me%Tsunami%Fault(ifl)%ObjGridData, & Array = Tsunami2D, & STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'FaultInputFile - ModuleHydrodynamic - ERR20' - call KillGridData(GridDataID = Me%Tsunami%Fault%ObjGridData, & + call KillGridData(GridDataID = Me%Tsunami%Fault(ifl)%ObjGridData, & STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'FaultInputFile - ModuleHydrodynamic - ERR30' @@ -25118,12 +25245,18 @@ end subroutine FaultInputFile !----------------------------------------------------------------------------------- - subroutine WaterLevel_Tsunami + subroutine WaterLevel_Tsunami(ifl) - if (Me%Tsunami%Fault%InputMethod == FaultFile_ ) then - call FaultInputFile - elseif (Me%Tsunami%Fault%InputMethod == FaultOkada1985_) then - call FaultOkada1985 + !Arguments---------------------------------------------------------------------- + integer :: ifl + + !Begin-------------------------------------------------------------------------- + + + if (Me%Tsunami%Fault(ifl)%InputMethod == FaultFile_ ) then + call FaultInputFile(ifl) + elseif (Me%Tsunami%Fault(ifl)%InputMethod == FaultOkada1985_) then + call FaultOkada1985(ifl) endif end subroutine WaterLevel_Tsunami @@ -49174,8 +49307,8 @@ subroutine Hydrodynamic_OutPut if (Me%CurrentTime >= NextProfileOutput) ProfileFileOK = .true. endif - if (ProfileFileOK .or. OutPutFileOK .or. TimeSeriesFileOK .or. OutPutWindowFileOK .or. & - OutPutSurfaceFileOK) then + if (ProfileFileOK .or. OutPutFileOK .or. TimeSeriesFileOK .or. & + OutPutWindowFileOK .or. OutPutSurfaceFileOK) then call ModifyMatrixesOutput @@ -49187,9 +49320,9 @@ subroutine Hydrodynamic_OutPut else - if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & - Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON .or. & - Me%OutW%OutPutWindowsON)then + if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & + Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON .or. & + Me%OutW%OutPutWindowsON)then call ModifyMatrixesOutput @@ -49206,6 +49339,11 @@ subroutine Hydrodynamic_OutPut !! $OMP SECTION if (Me%OutPut%hdf5ON) then + !Statistics output + call Statistics_OutPut(Me%OutPut%CenterU, Me%OutPut%CenterV, & + Me%OutPut%CenterW, Me%OutPut%ModulusH, & + Me%WaterLevel%New) + NextOutPut = Me%OutPut%NextOutPut if (NextOutPut <= Me%OutPut%Number) then @@ -49320,7 +49458,7 @@ subroutine Hydrodynamic_OutPut if (Me%OutPut%TurbineON) & call OutPut_Turbine(Me%ObjTurbine) - + !! $OMP SECTION if(Me%OutPut%WriteRestartFile .and. .not. Me%OutPut%Run_End)then @@ -49365,6 +49503,9 @@ subroutine Hydrodynamic_OutPut enddo endif + + + if (MonitorPerformance) call StopWatch ("ModuleHydrodynamic", "Hydrodynamic_OutPut") @@ -50853,11 +50994,11 @@ subroutine Write_HDF5_Format(iW) endif - !Do statistics analysis - if (.not. present(iW)) then - call Statistics_OutPut(Me%OutPut%CenterU, Me%OutPut%CenterV, Me%OutPut%CenterW, & - Me%OutPut%ModulusH, Me%WaterLevel%New) - endif + !!Do statistics analysis + !if (.not. present(iW)) then + ! call Statistics_OutPut(Me%OutPut%CenterU, Me%OutPut%CenterV, Me%OutPut%CenterW, & + ! Me%OutPut%ModulusH, Me%WaterLevel%New) + !endif CHUNK = CHUNK_J(WorkJLB, WorkJUB) @@ -52228,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 @@ -52275,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 @@ -54765,6 +54911,13 @@ subroutine Kill_Sub_Modules if (STAT_CALL /= SUCCESS_) stop 'Subroutine Kill_Sub_Modules - ModuleHydrodynamic. ERR10.' endif cd1 + +i111: if (Me%Tsunami%ON) then + + deallocate(Me%Tsunami%Fault) + + endif i111 + !Disposes the rest of the energy buffer if (Me%ComputeOptions%Energy) call KillEnergy @@ -55688,8 +55841,8 @@ Subroutine DeallocateVariables end if - if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & - Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON.or. & + if (Me%OutPut%TimeSerieON .or. Me%OutPut%hdf5ON .or. & + Me%OutPut%ProfileON .or. Me%OutPut%HDF5_Surface_ON.or. & Me%OutW%OutPutWindowsON) then call KillMatrixesOutput endif diff --git a/Software/MOHIDWater/ModuleInterfaceSedimentWater.F90 b/Software/MOHIDWater/ModuleInterfaceSedimentWater.F90 index 580eb4457..ca3820d29 100644 --- a/Software/MOHIDWater/ModuleInterfaceSedimentWater.F90 +++ b/Software/MOHIDWater/ModuleInterfaceSedimentWater.F90 @@ -1297,6 +1297,8 @@ subroutine ConstructShearStressResidual !Gets File Access Code call GetHDF5FileAccess (HDF5_READ = HDF5_READ) + ObjHDF5 = 0 + !Opens HDF5 File call ConstructHDF5 (ObjHDF5, trim(Me%Files%Initial), HDF5_READ, STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ConstructShearStressResidual - ModuleInterfaceSedimentWater - ERR50' diff --git a/Software/MOHIDWater/ModuleLagrangianGlobal.F90 b/Software/MOHIDWater/ModuleLagrangianGlobal.F90 index 26180b79e..00d129d13 100644 --- a/Software/MOHIDWater/ModuleLagrangianGlobal.F90 +++ b/Software/MOHIDWater/ModuleLagrangianGlobal.F90 @@ -1257,7 +1257,7 @@ Module ModuleLagrangianGlobal logical :: NetCDF = .false. integer :: ObjNETCDF = null_int integer :: NetCDF_DimID = null_int - + logical :: AmbientConc = OFF end type T_OutPut @@ -4185,6 +4185,16 @@ subroutine ConstructOrigins call GetComputeTimeStep(Me%ExternalVar%ObjTime, DT, STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ConstructOrigins - ModuleLagrangianGlobal - ERR70' + + call GetData(Me%OutPut%AmbientConc, & + Me%ObjEnterData, & + flag, & + SearchType = FromFile, & + keyword ='OUTPUT_AMBIENT_CONC', & + ClientModule ='ModuleLagrangianGlobal', & + Default = .false., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ConstructOrigins - ModuleLagrangianGlobal - ERR75' !Get time step for particle from file DT_PARTIC = null_real @@ -25784,7 +25794,9 @@ subroutine ParticleOutput () enddo if (CurrentOrigin%State%VariableGeom) then - call HDF5WriteParticAmbientConc(CurrentOrigin, em, OutPutNumber,Matrix1D) + if(Me%Output%AmbientConc)then + call HDF5WriteParticAmbientConc(CurrentOrigin, em, OutPutNumber,Matrix1D) + endif endif deallocate (Matrix1D) @@ -26845,11 +26857,12 @@ subroutine ParticleOutput () enddo - if (Me%State%VariableGeom) then - call HDF5WriteAllGroupParticAmbientConc(GroupName, ig, em, OutPutNumber,Matrix1D) + if (Me%State%VariableGeom) then + if(Me%Output%AmbientConc)then + call HDF5WriteAllGroupParticAmbientConc(GroupName, ig, em, OutPutNumber,Matrix1D) + endif endif - deallocate (Matrix1D) enddo dig 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. diff --git a/Software/SmallTools/ConvertToTimeSerie/ConvertToTimeSerie.f90 b/Software/SmallTools/ConvertToTimeSerie/ConvertToTimeSerie.f90 index d1b24fe35..049db62fd 100644 --- a/Software/SmallTools/ConvertToTimeSerie/ConvertToTimeSerie.f90 +++ b/Software/SmallTools/ConvertToTimeSerie/ConvertToTimeSerie.f90 @@ -14,7 +14,7 @@ ! !------------------------------------------------------------------------------ -!Nomfich.dat: +!nomfich.dat: ! ! IN_MODEL : char [-] !Path to the input file with ! !user's instructions (DataFile) diff --git a/Software/SmallTools/Delft3DGrid2Mohid/MainDelft3DGrid2Mohid.F90 b/Software/SmallTools/Delft3DGrid2Mohid/MainDelft3DGrid2Mohid.F90 new file mode 100644 index 000000000..c39067220 --- /dev/null +++ b/Software/SmallTools/Delft3DGrid2Mohid/MainDelft3DGrid2Mohid.F90 @@ -0,0 +1,259 @@ +!------------------------------------------------------------------------------ +! IST/MARETEC, Water Modelling Group, Mohid modelling system +!------------------------------------------------------------------------------ +! +! TITLE : Mohid Model +! PROJECT : Delft3DGrid2Mohid +! PROGRAM : MainDelft3DGrid2Mohid +! URL : http://www.mohid.com +! AFFILIATION : IST/MARETEC, Marine Modelling Group +! DATE : May 2003 +! REVISION : Frank Braunschweig /Luis Fernandes - v4.0 +! DESCRIPTION : Delft3DGrid2Mohid to create main program to use MOHID modules +! +!------------------------------------------------------------------------------ + +program MohidDelft3DGrid2Mohid + + use ModuleGlobalData + use ModuleTime + use ModuleEnterData + use ModuleFunctions + use ModuleHorizontalGrid + + implicit none + + type (T_Time) :: InitialSystemTime, FinalSystemTime + real :: TotalCPUTime, ElapsedSeconds + integer, dimension(8) :: F95Time + + integer :: STAT_CALL + character(len=PathLength) :: Input_Delft3D_Grid, Output_Mohid_Grid + real :: Delft3D_Grid_FillValue + + integer, parameter :: FileOpen = 1, FileClose = 0 + + + call ConstructMohidDelft3DGrid2Mohid + call ModifyMohidDelft3DGrid2Mohid + call KillMohidDelft3DGrid2Mohid + + contains + + !-------------------------------------------------------------------------- + + subroutine ConstructMohidDelft3DGrid2Mohid + + call StartUpMohid("MohidDelft3DGrid2Mohid") + + call StartCPUTime + + call ReadKeywords + + + end subroutine ConstructMohidDelft3DGrid2Mohid + + !-------------------------------------------------------------------------- + + subroutine ModifyMohidDelft3DGrid2Mohid + + !Local----------------------------------------------------------------- + + !Begin----------------------------------------------------------------- + + call ReadWriteGrids + + end subroutine ModifyMohidDelft3DGrid2Mohid + + !-------------------------------------------------------------------------- + + subroutine KillMohidDelft3DGrid2Mohid + + call StopCPUTime + + call ShutdownMohid ("MohidDelft3DGrid2Mohid", ElapsedSeconds, TotalCPUTime) + + end subroutine KillMohidDelft3DGrid2Mohid + + !-------------------------------------------------------------------------- + + subroutine StartCPUTime + + call date_and_time(Values = F95Time) + + call SetDate (InitialSystemTime, float(F95Time(1)), float(F95Time(2)), & + float(F95Time(3)), float(F95Time(5)), & + float(F95Time(6)), float(F95Time(7))+ & + F95Time(8)/1000.) + + end subroutine StartCPUTime + + !-------------------------------------------------------------------------- + + subroutine StopCPUTime + + call date_and_time(Values = F95Time) + + call SetDate (FinalSystemTime, float(F95Time(1)), float(F95Time(2)), & + float(F95Time(3)), float(F95Time(5)), & + float(F95Time(6)), float(F95Time(7))+ & + F95Time(8)/1000.) + + call cpu_time(TotalCPUTime) + + ElapsedSeconds = FinalSystemTime - InitialSystemTime + + end subroutine StopCPUTime + + !-------------------------------------------------------------------------- + + subroutine ReadKeywords + + !Local----------------------------------------------------------------- + character(PathLength) :: DataFile + integer :: STAT_CALL + integer :: ObjEnterData = 0 + integer :: FromFile, iflag + + !Begin----------------------------------------------------------------- + + call ReadFileName('IN_MODEL', DataFile, "MohidDelft3DGrid2Mohid", STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR10' + + call ConstructEnterData (ObjEnterData, DataFile, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR20' + + call GetExtractType (FromFile = FromFile) + + call GetData(Input_Delft3D_Grid, & + ObjEnterData, iflag, & + keyword = 'INPUT_DELFT3D_GRID', & + SearchType = FromFile, & + ClientModule ='MohidDelft3DGrid2Mohid', & + STAT = STAT_CALL) + + if (STAT_CALL /= SUCCESS_ .or. iflag == 0) then + call SetError(FATAL_, INTERNAL_, 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR30') + endif + + call GetData(Delft3D_Grid_FillValue, & + ObjEnterData, iflag, & + keyword = 'DELFT3D_GRID_FILLVALUE', & + SearchType = FromFile, & + ClientModule ='MohidDelft3DGrid2Mohid', & + STAT = STAT_CALL) + + if (STAT_CALL /= SUCCESS_ .or. iflag == 0) then + call SetError(FATAL_, INTERNAL_, 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR40') + endif + + call GetData(Output_Mohid_Grid, & + ObjEnterData, iflag, & + keyword = 'OUTPUT_MOHID_GRID', & + SearchType = FromFile, & + ClientModule ='MohidDelft3DGrid2Mohid', & + STAT = STAT_CALL) + + if (STAT_CALL /= SUCCESS_ .or. iflag == 0) then + call SetError(FATAL_, INTERNAL_, 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR50') + endif + + + call KillEnterData (ObjEnterData, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidDelft3DGrid2Mohid - ERR60' + + end subroutine ReadKeywords + + !-------------------------------------------------------------------------- + + subroutine ReadWriteGrids + + !Local----------------------------------------------------------------- + real, dimension(:,:), pointer :: XX, YY + integer :: Imax, Jmax, InputID, OutputID, i, j + integer :: STAT_CALL, DummyInt + character(13) :: DummyChar + + !Begin----------------------------------------------------------------- + + call UnitsManager(InputID, FileOpen, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + call SetError(FATAL_, INTERNAL_, 'ReadWriteGrids - MohidDelft3DGrid2Mohid - ERR10') + endif + + open(unit = InputID, file = Input_Delft3D_Grid, form = 'FORMATTED', status = 'UNKNOWN') + + do i =1, 5 + read(InputID,*) + enddo + + read(InputID,*) jmax, imax + read(InputID,*) + + allocate(XX(0:imax,0:jmax),YY(0:imax,0:jmax)) + + do i=1,imax + read(InputID,*) DummyChar,DummyInt,(XX(i,j),j=1,jmax) + enddo + + do i=1,imax + read(InputID,*) DummyChar,DummyInt,(YY(i,j),j=1,jmax) + enddo + + do i=1,imax + do j=1,jmax + if (XX(i,j) == Delft3D_Grid_FillValue) then + XX(i,j) = FillValueReal + YY(i,j) = FillValueReal + endif + enddo + enddo + + call UnitsManager(InputID, FileClose, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + call SetError(FATAL_, INTERNAL_, 'ReadWriteGrids - MohidDelft3DGrid2Mohid - ERR20') + endif + + + call UnitsManager(OutputID, FileOpen, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + call SetError(FATAL_, INTERNAL_, 'ReadWriteGrids - MohidDelft3DGrid2Mohid - ERR30') + endif + + open(unit = OutputID, file = Output_Mohid_Grid, form = 'FORMATTED', status = 'UNKNOWN') + + + write(OutputID,*) "COMENT1 : From Delft3D Grid" + write(OutputID,*) "COMENT2 : to MOHID Grid" + write(OutputID,*) "ILB_IUB : 1", Imax-1 + write(OutputID,*) "JLB_JUB : 1", Jmax-1 + write(OutputID,*) "COORD_TIP : 5" + write(OutputID,*) "ORIGIN : 0 0" + write(OutputID,*) "GRID_ANGLE : 0" + write(OutputID,*) "LATITUDE : 42" + write(OutputID,*) "LONGITUDE : -9" + + write(OutputID,*) "FILL_VALUE : -99.0" + + + write(OutputID,*) "" + + do i=1,imax + do j=1,jmax + write(OutputID,*) XX(i,j), YY(i,j) + enddo + enddo + + + write(OutputID,*) "<\CornersXY>" + + call UnitsManager(OutputID, FileClose, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) then + call SetError(FATAL_, INTERNAL_, 'ReadWriteGrids - MohidDelft3DGrid2Mohid - ERR50') + endif + + + end subroutine ReadWriteGrids + + +end program MohidDelft3DGrid2Mohid diff --git a/Software/SmallTools/HDF5_2_EsriGridData/HDF5_2_EsriGridData.F90 b/Software/SmallTools/HDF5_2_EsriGridData/HDF5_2_EsriGridData.F90 index bec251cea..83c9033f1 100644 --- a/Software/SmallTools/HDF5_2_EsriGridData/HDF5_2_EsriGridData.F90 +++ b/Software/SmallTools/HDF5_2_EsriGridData/HDF5_2_EsriGridData.F90 @@ -82,6 +82,8 @@ Program HDF5_2_EsriGridData logical :: TransferToOutGrid logical :: ExportXYZ logical :: ExportXY_Vector + logical :: VectorON + character(len=StringLength) :: Vector_X, Vector_Y end type T_HDF5_2_EsriGridData type(T_HDF5_2_EsriGridData), pointer :: Me @@ -158,10 +160,10 @@ subroutine ConstructHDF5_2_EsriGridData call GetHorizontalGridSize(Me%ObjHorizontalGridOut, Me%SizeOut, Me%WorkSizeOut, STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_) stop 'ConstructHDF5_2_EsriGridData - HDF5_2_EsriGridData - ERR40' - call ConstructFatherGridLocation(Me%ObjHorizontalGridOut, Me%ObjHorizontalGrid, & - OkCross = .false., OkZ = .true., & - OkU = .false., OkV = .false., STAT = STAT_CALL) - if (STAT_CALL /= SUCCESS_) stop 'ConstructHDF5_2_EsriGridData - HDF5_2_EsriGridData - ERR50' + !call ConstructFatherGridLocation(Me%ObjHorizontalGridOut, Me%ObjHorizontalGrid, & + ! OkCross = .false., OkZ = .true., & + ! OkU = .false., OkV = .false., STAT = STAT_CALL) + !if (STAT_CALL /= SUCCESS_) stop 'ConstructHDF5_2_EsriGridData - HDF5_2_EsriGridData - ERR50' endif @@ -373,6 +375,37 @@ subroutine ReadGlobalOptions if (STAT_CALL /= SUCCESS_) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR120' + call GetData(Me%VectorON, & + Me%ObjEnterData, iflag, & + SearchType = FromFile, & + keyword = 'VECTOR_ON', & + ClientModule = 'HDF5ToASCIIandBIN', & + default = .false., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR130' + + if (Me%VectorON) then + + call GetData(Me%Vector_X, & + Me%ObjEnterData, iflag, & + SearchType = FromFile, & + keyword = 'VECTOR_X', & + ClientModule = 'HDF5ToASCIIandBIN', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR140' + if (iflag == 0 ) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR150' + + call GetData(Me%Vector_Y, & + Me%ObjEnterData, iflag, & + SearchType = FromFile, & + keyword = 'VECTOR_Y', & + ClientModule = 'HDF5ToASCIIandBIN', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR160' + if (iflag == 0 ) stop 'ReadGlobalOptions - HDF5_2_EsriGridData - ERR170' + + endif + end subroutine ReadGlobalOptions !-------------------------------------------------------------------------- @@ -597,6 +630,7 @@ subroutine ModifyHDF5_2_EsriGridData logical :: Found2Blanks real, dimension(:,:), pointer :: CoordX, CoordY real, dimension(:,:), pointer :: CoordXout, CoordYout + integer, dimension(:,:), pointer :: InPutMap2D logical :: Exist type (T_PointF), pointer :: Point @@ -639,6 +673,7 @@ subroutine ModifyHDF5_2_EsriGridData allocate(Aux2DOut(Me%SizeOut%ILB:Me%SizeOut%IUB, Me%SizeOut%JLB:Me%SizeOut%JUB)) allocate(Me%InPutMap (Me%Size%ILB :Me%Size%IUB , Me%Size%JLB :Me%Size%JUB,1: Me%KUB)) + allocate( InPutMap2D (Me%Size%ILB :Me%Size%IUB , Me%Size%JLB :Me%Size%JUB )) allocate(Me%OutPutMap(Me%SizeOut%ILB:Me%SizeOut%IUB, Me%SizeOut%JLB:Me%SizeOut%JUB)) @@ -750,18 +785,45 @@ subroutine ModifyHDF5_2_EsriGridData Aux3DOut(:,:,:) = Me%FillValue Aux2D (:,: ) = Aux3D(:,:,k) + if (Me%LandMaskON) then + InPutMap2D (:,:) = Me%InPutMap(:,:,k) + else + InPutMap2D (:,:) = 1 + endif if (Me%InterpolOut) then - call InterpolRegularGrid(HorizontalGridSonID = Me%ObjHorizontalGridOut,& - HorizontalGridFatherID = Me%ObjHorizontalGrid, & + !if (Me%Regular) then + ! + ! call InterpolRegularGrid(HorizontalGridSonID = Me%ObjHorizontalGridOut,& + ! HorizontalGridFatherID = Me%ObjHorizontalGrid, & + ! Field2DFather = Aux2D, & + ! Field2DSon = Aux2DOut, & + ! ComputeFather = Me%InPutMap, & + ! KUBFather = k, & + ! STAT = STAT_CALL) + ! if (STAT_CALL /= SUCCESS_)stop 'ModifyHDF5_2_EsriGridData - HDF5_2_EsriGridData - ERR80' + ! + !else + + + + do j = Me%WorkSizeOut%JLB, Me%WorkSizeOut%JUB + do i = Me%WorkSizeOut%ILB, Me%WorkSizeOut%IUB + + Aux2DOut(i,j) = InterpolXYPoint(HorizontalGridID = Me%ObjHorizontalGrid, & Field2DFather = Aux2D, & - Field2DSon = Aux2DOut, & - ComputeFather = Me%InPutMap, & - KUBFather = k, & + ComputeFather = InPutMap2D, & + XInput = CoordXout(i, j), & + YInput = CoordYout(i, j), & STAT = STAT_CALL) if (STAT_CALL /= SUCCESS_)stop 'ModifyHDF5_2_EsriGridData - HDF5_2_EsriGridData - ERR80' + enddo + enddo + + !endif + Aux3DOut(:,:,k) = Aux2DOut(:,:) do j = Me%WorkSizeOut%JLB, Me%WorkSizeOut%JUB @@ -796,7 +858,7 @@ subroutine ModifyHDF5_2_EsriGridData if (Me%ExportXYZ) then - call Export_To_XYZ(Aux3DOut, k, ILB, IUB, JLB, JUB, Me%OutputESRI(l)) + call Export_To_XYZ(Aux3DOut, k, ILBout, IUBout, JLBout, JUBout, Me%OutputESRI(l)) endif @@ -832,6 +894,7 @@ subroutine ModifyHDF5_2_EsriGridData if (Me%TransferToOutGrid) then deallocate(Aux3DOut) deallocate(Aux2DOut) + deallocate(InPutMap2D) endif @@ -1178,9 +1241,9 @@ subroutine Export_To_XYZ(Aux3D, k, ILB, IUB, JLB, JUB, Filename) do i = ILB, IUB do j = JLB, JUB - if (Aux3D(i,j,k) == FillValueReal) then - write(Unit,*)CoordX(i, j), CoordY(i, j), Me%FillValue - else + if (Aux3D(i,j,k) > FillValueReal .and. Aux3D(i,j,k) /= Me%FillValue) then + ! write(Unit,*)CoordX(i, j), CoordY(i, j), Me%FillValue + !else write(Unit,*)CoordX(i, j), CoordY(i, j), Aux3D(i,j,k) endif enddo diff --git a/Software/SmallTools/MeshGlue/MainMeshGlue.F90 b/Software/SmallTools/MeshGlue/MainMeshGlue.F90 new file mode 100644 index 000000000..73c3913f7 --- /dev/null +++ b/Software/SmallTools/MeshGlue/MainMeshGlue.F90 @@ -0,0 +1,134 @@ +!------------------------------------------------------------------------------ +! IST/MARETEC, Water Modelling Group, Mohid modelling system +!------------------------------------------------------------------------------ +! +! TITLE : Mohid Model +! PROJECT : MeshGlue +! PROGRAM : MainMeshGlue +! URL : http://www.mohid.com +! AFFILIATION : Hidromod +! DATE : December 2020 +! REVISION : Paulo Leitao - v1.0 +! DESCRIPTION : MeshGlue to create main program to use MOHID modules +! +!------------------------------------------------------------------------------ + +program MohidMeshGlue + + use ModuleGlobalData + use ModuleEnterData + use ModuleFunctions + use ModuleTime + use ModuleMeshGlue + + implicit none + + type (T_Time) :: InitialSystemTime, FinalSystemTime + real :: TotalCPUTime, ElapsedSeconds + integer, dimension(8) :: F95Time + + integer :: STAT_CALL + integer :: MeshGlueID = 0 + + + call ConstructMohidMeshGlue + call ModifyMohidMeshGlue + call KillMohidMeshGlue + + contains + + !-------------------------------------------------------------------------- + + subroutine ConstructMohidMeshGlue + + call StartUpMohid("MohidMeshGlue") + + call StartCPUTime + + !call ReadKeywords + + call ConstructMeshGlue(ObjMeshGlueID = MeshGlueID, STAT = STAT_CALL) + + + end subroutine ConstructMohidMeshGlue + + !-------------------------------------------------------------------------- + + subroutine ModifyMohidMeshGlue + + !Local----------------------------------------------------------------- + logical :: Running + + + call ModifyMeshGlue(ObjMeshGlueID = MeshGlueID, STAT = STAT_CALL) + + + + end subroutine ModifyMohidMeshGlue + + !-------------------------------------------------------------------------- + + subroutine KillMohidMeshGlue + + call KillMeshGlue(ObjMeshGlueID = MeshGlueID, STAT = STAT_CALL) + + call StopCPUTime + + call ShutdownMohid ("MohidMeshGlue", ElapsedSeconds, TotalCPUTime) + + end subroutine KillMohidMeshGlue + + !-------------------------------------------------------------------------- + + subroutine StartCPUTime + + call date_and_time(Values = F95Time) + + call SetDate (InitialSystemTime, float(F95Time(1)), float(F95Time(2)), & + float(F95Time(3)), float(F95Time(5)), & + float(F95Time(6)), float(F95Time(7))+ & + F95Time(8)/1000.) + + end subroutine StartCPUTime + + !-------------------------------------------------------------------------- + + subroutine StopCPUTime + + call date_and_time(Values = F95Time) + + call SetDate (FinalSystemTime, float(F95Time(1)), float(F95Time(2)), & + float(F95Time(3)), float(F95Time(5)), & + float(F95Time(6)), float(F95Time(7))+ & + F95Time(8)/1000.) + + call cpu_time(TotalCPUTime) + + ElapsedSeconds = FinalSystemTime - InitialSystemTime + + end subroutine StopCPUTime + + !-------------------------------------------------------------------------- + + subroutine ReadKeywords + + !Local----------------------------------------------------------------- + character(PathLength) :: DataFile + integer :: STAT_CALL + integer :: ObjEnterData = 0 + integer :: FromFile + + call ReadFileName('IN_MODEL', DataFile, "MohidMeshGlue", STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidMeshGlue - ERR01' + + call ConstructEnterData (ObjEnterData, DataFile, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidMeshGlue - ERR02' + + call GetExtractType (FromFile = FromFile) + + call KillEnterData (ObjEnterData, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadKeywords - MohidMeshGlue - ERR03' + + end subroutine ReadKeywords + +end program MohidMeshGlue diff --git a/Software/SmallTools/MeshGlue/ModuleMeshGlue.F90 b/Software/SmallTools/MeshGlue/ModuleMeshGlue.F90 new file mode 100644 index 000000000..91b4ef23d --- /dev/null +++ b/Software/SmallTools/MeshGlue/ModuleMeshGlue.F90 @@ -0,0 +1,836 @@ +!------------------------------------------------------------------------------ +! IST/MARETEC, Water Modelling Group, Mohid modelling system +!------------------------------------------------------------------------------ +! +! TITLE : Mohid Model +! PROJECT : Mohid Base 1 +! MODULE : MeshGlue +! URL : http://www.mohid.com +! AFFILIATION : Hidromod +! DATE : December 2020 +! REVISION : Paulo Leitao - v1.0 +! DESCRIPTION : Module to serve as MeshGlue to create new modules +! +!------------------------------------------------------------------------------ + + +Module ModuleMeshGlue + + use ModuleGlobalData + use ModuleEnterData + use ModuleHorizontalGrid + + implicit none + + private + + !Subroutines--------------------------------------------------------------- + + !Constructor + public :: ConstructMeshGlue + private :: AllocateInstance + + !Selector + + + !Modifier + public :: ModifyMeshGlue + + !Destructor + public :: KillMeshGlue + private :: DeAllocateInstance + + !Management + private :: Ready + private :: LocateObjMeshGlue + + !Interfaces---------------------------------------------------------------- + + !Types--------------------------------------------------------------------- + private :: T_Grid + type T_Grid + ! Filename of the input grid data file + character(len=PathLength) :: FileName = null_str + !Instance of ModuleHorizontalGrid + integer :: ObjHorizontalGrid= 0 + integer :: IO = null_int + integer :: JO = null_int + logical :: Flip_IJ = .false. + logical :: RowFlip = .false. + logical :: ColumnFlip = .false. + integer :: ILB = null_int + integer :: IUB = null_int + integer :: JLB = null_int + integer :: JUB = null_int + end type T_Grid + + + private :: T_MeshGlue + type T_MeshGlue + type (T_Grid), dimension(:), pointer :: GridIn + integer :: GridIn_Number = null_int + ! Filename of the input grid data file + character(len=PathLength) :: FileNameGrid_Out_Start = null_str + character(len=PathLength) :: FileNameGrid_Out = null_str + !Instance of ModuleHorizontalGrid + integer :: ObjHorizontalGridOut1 = 0 + integer :: ObjHorizontalGridOut2 = 0 + + logical :: ResetGrid = .false. + + real(8), pointer, dimension(:,:) :: XX_IE => null() + real(8), pointer, dimension(:,:) :: YY_IE => null() + + type(T_Size2D) :: Size, WorkSize + + integer :: InstanceID + !Instance of ObjEnterData + integer :: ObjEnterData = 0 + + type(T_MeshGlue), pointer :: Next + + + end type T_MeshGlue + + !Global Module Variables + type (T_MeshGlue), pointer :: FirstObjMeshGlue + type (T_MeshGlue), pointer :: Me + + + !-------------------------------------------------------------------------- + + contains + + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONS + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine ConstructMeshGlue(ObjMeshGlueID, STAT) + + !Arguments--------------------------------------------------------------- + integer :: ObjMeshGlueID + integer, optional, intent(OUT) :: STAT + + !External---------------------------------------------------------------- + integer :: ready_ + + !Local------------------------------------------------------------------- + integer :: STAT_ + + !------------------------------------------------------------------------ + + STAT_ = UNKNOWN_ + + !Assures nullification of the global variable + if (.not. ModuleIsRegistered(mMeshGlue_)) then + nullify (FirstObjMeshGlue) + call RegisterModule (mMeshGlue_) + endif + + call Ready(ObjMeshGlueID, ready_) + +cd0 : if (ready_ .EQ. OFF_ERR_) then + + call AllocateInstance + + !Returns ID + ObjMeshGlueID = Me%InstanceID + + call ReadGridIn + + call ReadGridOut + + STAT_ = SUCCESS_ + + else cd0 + + stop 'ModuleMeshGlue - ConstructMeshGlue - ERR01' + + end if cd0 + + + if (present(STAT)) STAT = STAT_ + + !---------------------------------------------------------------------- + + end subroutine ConstructMeshGlue + + !-------------------------------------------------------------------------- + + subroutine AllocateInstance + + !Arguments------------------------------------------------------------- + + !Local----------------------------------------------------------------- + type (T_MeshGlue), pointer :: NewObjMeshGlue + type (T_MeshGlue), pointer :: PreviousObjMeshGlue + + + !Allocates new instance + allocate (NewObjMeshGlue) + nullify (NewObjMeshGlue%Next) + + !Insert New Instance into list and makes Current point to it + if (.not. associated(FirstObjMeshGlue)) then + FirstObjMeshGlue => NewObjMeshGlue + Me => NewObjMeshGlue + else + PreviousObjMeshGlue => FirstObjMeshGlue + Me => FirstObjMeshGlue%Next + do while (associated(Me)) + PreviousObjMeshGlue => Me + Me => Me%Next + enddo + Me => NewObjMeshGlue + PreviousObjMeshGlue%Next => NewObjMeshGlue + endif + + Me%InstanceID = RegisterNewInstance (mMeshGlue_) + + + end subroutine AllocateInstance + + + !-------------------------------------------------------------------------- + + subroutine ReadGridIn + + !Local----------------------------------------------------------------- + integer :: STAT_CALL, iflag + integer :: ifl, ClientNumber + logical :: BlockFound + + !---------------------------------------------------------------------- + + !Construct enter data + call ConstructEnterData(Me%ObjEnterData, "MeshGlue.dat", STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridIn - ModuleMeshGlue - ERR10' + + call GetData(Me%GridIn_Number, & + Me%ObjEnterData, iflag, & + Keyword = 'GRID_IN_NUMBER', & + Default = 1, & + SearchType = FromFile, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridIn - ModuleMeshGlue - ERR20' + + allocate(Me%GridIn(Me%GridIn_Number)) + +do1 : do ifl =1, Me%GridIn_Number + + call ExtractBlockFromBuffer(Me%ObjEnterData, & + ClientNumber = ClientNumber, & + block_begin = "", & + block_end = "", & + BlockFound = BlockFound, & + STAT = STAT_CALL) +cd1 : if (STAT_CALL == SUCCESS_ ) then +cd2 : if (BlockFound) then + + call ReadGridKeywords(FromBlock, Me%GridIn(ifl)) + + else cd2 + + stop 'ReadGridIn - ModuleMeshGlue - ERR030' + + endif cd2 + + else cd1 + stop 'ReadGridIn - ModuleMeshGlue - ERR040' + endif cd1 + + enddo do1 + + call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridIn - ModuleMeshGlue - ERR050' + + + end subroutine ReadGridIn + + !-------------------------------------------------------------------------- + + subroutine ReadGridKeywords(FromWhere, Grid) + + !Arguments------------------------------------------------------------- + integer :: FromWhere + type (T_Grid) :: Grid + + !Local----------------------------------------------------------------- + integer :: STAT_CALL, iflag + type(T_Size2D) :: WorkSize + + !---------------------------------------------------------------------- + + !Read grid filemane of the input grid data file + call GetData(Grid%FileName, & + Me%ObjEnterData, iflag, & + Keyword = 'FILENAME', & + SearchType = FromWhere, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR10' + + if (iflag == 0) then + stop 'ReadGridKeywords - ModuleMeshGlue - ERR20' + endif + + call ConstructHorizontalGrid(HorizontalGridID = Grid%ObjHorizontalGrid, & + DataFile = Grid%FileName, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR30' + + + call GetHorizontalGridSize(HorizontalGridID = Grid%ObjHorizontalGrid, & + WorkSize = WorkSize, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR40' + + call GetData(Grid%IO, & + Me%ObjEnterData, iflag, & + Keyword = 'IO', & + SearchType = FromWhere, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR50' + + if (iflag /= 1) then + stop 'ReadGridKeywords - ModuleMeshGlue - ERR60' + endif + + call GetData(Grid%JO, & + Me%ObjEnterData, iflag, & + Keyword = 'JO', & + SearchType = FromWhere, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR70' + + if (iflag /= 1) then + stop 'ReadGridKeywords - ModuleMeshGlue - ERR80' + endif + + call GetData(Grid%Flip_IJ, & + Me%ObjEnterData, iflag, & + Keyword = 'FLIP_IJ', & + SearchType = FromWhere, & + default = .false., & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR90' + + call GetData(Grid%RowFlip, & + Me%ObjEnterData, iflag, & + Keyword = 'ROW_FLIP', & + SearchType = FromWhere, & + default = .false., & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR100' + + call GetData(Grid%ColumnFlip, & + Me%ObjEnterData, iflag, & + Keyword = 'COLUMN_FLIP', & + SearchType = FromWhere, & + default = .false., & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR110' + + call GetData(Grid%ILB, & + Me%ObjEnterData, iflag, & + Keyword = 'ILB', & + SearchType = FromWhere, & + default = WorkSize%ILB, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR120' + + call GetData(Grid%IUB, & + Me%ObjEnterData, iflag, & + Keyword = 'IUB', & + SearchType = FromWhere, & + default = WorkSize%IUB, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR130' + + call GetData(Grid%JLB, & + Me%ObjEnterData, iflag, & + Keyword = 'JLB', & + SearchType = FromWhere, & + default = WorkSize%JLB, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR140' + + call GetData(Grid%JUB, & + Me%ObjEnterData, iflag, & + Keyword = 'JUB', & + SearchType = FromWhere, & + default = WorkSize%JUB, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR150' + + + end subroutine ReadGridKeywords + + !-------------------------------------------------------------------------- + + subroutine ReadGridOut + + !Local----------------------------------------------------------------- + integer :: STAT_CALL, iflag + integer :: ifl, ClientNumber + logical :: BlockFound + + !---------------------------------------------------------------------- + + + !Read grid filemane of the first iteration of the output grid + call GetData(Me%FileNameGrid_Out_Start, & + Me%ObjEnterData, iflag, & + Keyword = 'FILENAME_OUT_START', & + SearchType = FromFile, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR10' + + if (iflag == 0) then + stop 'ReadGridKeywords - ModuleMeshGlue - ERR20' + endif + + call ConstructHorizontalGrid(HorizontalGridID = Me%ObjHorizontalGridOut1, & + DataFile = Me%FileNameGrid_Out_Start, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR30' + + !Read grid filemane of the first iteration of the output grid + call GetData(Me%FileNameGrid_Out, & + Me%ObjEnterData, iflag, & + Keyword = 'FILENAME_OUT', & + SearchType = FromFile, & + ClientModule = 'ModuleMeshGlue', & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR40' + + if (iflag == 0) then + stop 'ReadGridKeywords - ModuleMeshGlue - ERR50' + endif + + call GetData(Me%ResetGrid, & + Me%ObjEnterData, iflag, & + Keyword = 'RESET_GRID', & + SearchType = FromFile, & + ClientModule = 'ModuleMeshGlue', & + default = .false., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR50' + + + Me%ObjHorizontalGridOut2 = 0 + + + !Kill enter data + call KillEnterData(Me%ObjEnterData, STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'ReadGridKeywords - ModuleMeshGlue - ERR60' + + end subroutine ReadGridOut + !-------------------------------------------------------------------------- + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SE + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !-------------------------------------------------------------------------- + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODI + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine ModifyMeshGlue(ObjMeshGlueID, STAT) + + !Arguments------------------------------------------------------------- + integer :: ObjMeshGlueID + integer, intent(OUT), optional :: STAT + + !Local----------------------------------------------------------------- + integer :: STAT_, ready_ + + !---------------------------------------------------------------------- + + STAT_ = UNKNOWN_ + + call Ready(ObjMeshGlueID, ready_) + + if (ready_ .EQ. IDLE_ERR_) then + + call WriteNewGrid + + STAT_ = SUCCESS_ + else + STAT_ = ready_ + end if + + if (present(STAT)) STAT = STAT_ + + end subroutine ModifyMeshGlue + + !-------------------------------------------------------------------------- + + subroutine WriteNewGrid + + !Local----------------------------------------------------------------- + real, pointer, dimension(:,:) :: XX_IE, YY_IE + real, pointer, dimension(:) :: Dummy + real :: Latitude, Longitude + integer :: STAT_CALL + + !Begin----------------------------------------------------------------- + + + call GetHorizontalGridSize(HorizontalGridID = Me%ObjHorizontalGridOut1, & + Size = Me%Size, & + WorkSize = Me%WorkSize, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR10' + + call GetLatitudeLongitude(HorizontalGridID = Me%ObjHorizontalGridOut1, & + Latitude = Latitude, & + Longitude = Longitude, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR20' + + allocate(Me%XX_IE(Me%Size%ILB:Me%Size%IUB,Me%Size%JLB:Me%Size%JUB)) + allocate(Me%YY_IE(Me%Size%ILB:Me%Size%IUB,Me%Size%JLB:Me%Size%JUB)) + + call GetGridLatitudeLongitude(HorizontalGridID = Me%ObjHorizontalGridOut1, & + GridLatitudeConn = YY_IE, & + GridLongitudeConn = XX_IE, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR30' + + Me%XX_IE(:,:) = XX_IE(:,:) + Me%YY_IE(:,:) = YY_IE(:,:) + + if (Me%ResetGrid) then + Me%XX_IE(:,:) = FillValueReal + Me%YY_IE(:,:) = FillValueReal + endif + + call UpdateGrid + + call ConstructHorizontalGrid(HorizontalGridID = Me%ObjHorizontalGridOut2, & + LatitudeConn = Me%YY_IE, & + LongitudeConn = Me%XX_IE, & + XX = Dummy, & + YY = Dummy, & + Latitude = Latitude, & + Longitude = Longitude, & + ILB = Me%WorkSize%ILB, & + IUB = Me%WorkSize%IUB, & + JLB = Me%WorkSize%JLB, & + JUB = Me%WorkSize%JUB, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR40' + + call WriteGridData_ASCII (FileName = Me%FileNameGrid_Out, & + COMENT1 = "Mesh Glue", & + COMENT2 = "Line 2", & + HorizontalGridID = Me%ObjHorizontalGridOut2, & + FillValue = FillValueReal, & + Overwrite = .true., & + DistortionYes = .true., & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR50' + + + end subroutine WriteNewGrid + !-------------------------------------------------------------------------- + + subroutine UpdateGrid() + + !Local----------------------------------------------------------------- + real, pointer, dimension(:,:) :: XX_IE, YY_IE + real, pointer, dimension(:,:) :: XX_Aux, YY_Aux, Aux2D + real, pointer, dimension(:) :: Dummy + real :: Latitude, Longitude + integer :: STAT_CALL, ifl, io, jo, di, dj + integer :: ILB, IUB, JLB, JUB, i, j + + !Begin----------------------------------------------------------------- + + + do ifl =1, Me%GridIn_Number + + + call GetGridLatitudeLongitude(HorizontalGridID = Me%GridIn(ifl)%ObjHorizontalGrid, & + GridLatitudeConn = YY_IE, & + GridLongitudeConn = XX_IE, & + STAT = STAT_CALL) + if (STAT_CALL /= SUCCESS_) stop 'WriteNewGrid - ModuleMeshGlue - ERR30' + + + io = Me%GridIn(ifl)%IO + jo = Me%GridIn(ifl)%JO + + if (Me%GridIn(ifl)%Flip_IJ) then + + ILB = Me%GridIn(ifl)%JLB + IUB = Me%GridIn(ifl)%JUB + JLB = Me%GridIn(ifl)%ILB + JUB = Me%GridIn(ifl)%IUB + + allocate(XX_Aux(ILB:IUB+1,JLB:JUB+1)) + allocate(YY_Aux(ILB:IUB+1,JLB:JUB+1)) + + do i = ILB, IUB + 1 + do j = JLB, JUB + 1 + XX_Aux(i, j) = XX_IE(j, i) + YY_Aux(i, j) = YY_IE(j, i) + enddo + enddo + + else + + ILB = Me%GridIn(ifl)%ILB + IUB = Me%GridIn(ifl)%IUB + JLB = Me%GridIn(ifl)%JLB + JUB = Me%GridIn(ifl)%JUB + + allocate(XX_Aux(ILB:IUB+1,JLB:JUB+1)) + allocate(YY_Aux(ILB:IUB+1,JLB:JUB+1)) + + XX_Aux(ILB:IUB+1,JLB:JUB+1) = XX_IE(ILB:IUB+1,JLB:JUB+1) + YY_Aux(ILB:IUB+1,JLB:JUB+1) = YY_IE(ILB:IUB+1,JLB:JUB+1) + + endif + + if (Me%GridIn(ifl)%RowFlip) then + + allocate(Aux2D(ILB:IUB+1,JLB:JUB+1)) + + Aux2D(:,:) = XX_Aux(:,:) + + do i = ILB, IUB + 1 + do j = JLB, JUB + 1 + XX_Aux(i, j) = Aux2D(IUB + 2 - i, j) + enddo + enddo + + Aux2D(:,:) = YY_Aux(:,:) + + do i = ILB, IUB + 1 + do j = JLB, JUB + 1 + YY_Aux(i, j) = Aux2D(IUB + 2 - i, j) + enddo + enddo + + endif + + if (Me%GridIn(ifl)%ColumnFlip) then + + allocate(Aux2D(ILB:IUB+1,JLB:JUB+1)) + + Aux2D(:,:) = XX_Aux(:,:) + + do i = ILB, IUB + 1 + do j = JLB, JUB + 1 + XX_Aux(i, j) = Aux2D(i, JUB + 2 - j) + enddo + enddo + + Aux2D(:,:) = YY_Aux(:,:) + + do i = ILB, IUB + 1 + do j = JLB, JUB + 1 + YY_Aux(i, j) = Aux2D(i, JUB + 2 - j) + enddo + enddo + + endif + + dj = JUB + 1 - JLB + di = IUB + 1 - ILB + + Me%XX_IE(io:io+di,jo:jo+dj) = XX_Aux(ILB:IUB+1, JLB:JUB+1) + + Me%YY_IE(io:io+di,jo:jo+dj) = YY_Aux(ILB:IUB+1, JLB:JUB+1) + + deallocate(XX_Aux, YY_Aux) + + if (associated(Aux2D)) deallocate(Aux2D) + + enddo + + + + end subroutine UpdateGrid + + !-------------------------------------------------------------------------- + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + + subroutine KillMeshGlue(ObjMeshGlueID, STAT) + + !Arguments--------------------------------------------------------------- + integer :: ObjMeshGlueID + integer, optional, intent(OUT) :: STAT + + !External---------------------------------------------------------------- + integer :: ready_ + + !Local------------------------------------------------------------------- + integer :: STAT_, nUsers + + !------------------------------------------------------------------------ + + STAT_ = UNKNOWN_ + + call Ready(ObjMeshGlueID, ready_) + +cd1 : if (ready_ .NE. OFF_ERR_) then + + nUsers = DeassociateInstance(mMeshGlue_, Me%InstanceID) + + if (nUsers == 0) then + + !Deallocates Instance + call DeallocateInstance () + + ObjMeshGlueID = 0 + STAT_ = SUCCESS_ + + end if + else + STAT_ = ready_ + end if cd1 + + if (present(STAT)) STAT = STAT_ + + + !------------------------------------------------------------------------ + + end subroutine KillMeshGlue + + + !------------------------------------------------------------------------ + + + subroutine DeallocateInstance () + + !Arguments------------------------------------------------------------- + + !Local----------------------------------------------------------------- + type (T_MeshGlue), pointer :: AuxObjMeshGlue + type (T_MeshGlue), pointer :: PreviousObjMeshGlue + + !Updates pointers + if (Me%InstanceID == FirstObjMeshGlue%InstanceID) then + FirstObjMeshGlue => FirstObjMeshGlue%Next + else + PreviousObjMeshGlue => FirstObjMeshGlue + AuxObjMeshGlue => FirstObjMeshGlue%Next + do while (AuxObjMeshGlue%InstanceID /= Me%InstanceID) + PreviousObjMeshGlue => AuxObjMeshGlue + AuxObjMeshGlue => AuxObjMeshGlue%Next + enddo + + !Now update linked list + PreviousObjMeshGlue%Next => AuxObjMeshGlue%Next + + endif + + !Deallocates instance + deallocate (Me) + nullify (Me) + + + end subroutine DeallocateInstance + + !-------------------------------------------------------------------------- + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + !MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEME + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + !-------------------------------------------------------------------------- + + subroutine Ready (ObjMeshGlue_ID, ready_) + + !Arguments------------------------------------------------------------- + integer :: ObjMeshGlue_ID + integer :: ready_ + + !---------------------------------------------------------------------- + + nullify (Me) + +cd1: if (ObjMeshGlue_ID > 0) then + call LocateObjMeshGlue (ObjMeshGlue_ID) + ready_ = VerifyReadLock (mMeshGlue_, Me%InstanceID) + else + ready_ = OFF_ERR_ + end if cd1 + + !---------------------------------------------------------------------- + + end subroutine Ready + + !-------------------------------------------------------------------------- + + subroutine LocateObjMeshGlue (ObjMeshGlueID) + + !Arguments------------------------------------------------------------- + integer :: ObjMeshGlueID + + !Local----------------------------------------------------------------- + + Me => FirstObjMeshGlue + do while (associated (Me)) + if (Me%InstanceID == ObjMeshGlueID) exit + Me => Me%Next + enddo + + if (.not. associated(Me)) stop 'ModuleMeshGlue - LocateObjMeshGlue - ERR01' + + end subroutine LocateObjMeshGlue + + !-------------------------------------------------------------------------- + +end module ModuleMeshGlue + + + + + + + + + diff --git a/Software/SmallTools/TimeSeriesAnalyser/ModuleTimeSeriesAnalyser.f90 b/Software/SmallTools/TimeSeriesAnalyser/ModuleTimeSeriesAnalyser.f90 index 1c0be94e8..82f03ee4d 100644 --- a/Software/SmallTools/TimeSeriesAnalyser/ModuleTimeSeriesAnalyser.f90 +++ b/Software/SmallTools/TimeSeriesAnalyser/ModuleTimeSeriesAnalyser.f90 @@ -1493,7 +1493,7 @@ subroutine ModifyTimeSeriesSmooth jaux = int(real(iaux)/2.) + 1 - write(Me%iSmooth,*) Me%TimeTS(i)-Me%BeginTime, SortArray(jaux) + write(Me%iSmooth,*) Me%TimeTS(i)-Me%BeginTime, SortArray(jaux), Me%DataMatrix(i,Me%DataColumn) !write(Me%iSmooth,*) Me%TimeTS(i)-Me%BeginTime, sum(SortArray(1:iaux))/iaux deallocate(SortArray) diff --git a/Solutions/VisualStudio2017_IntelFortran19/MOHIDNumerics/MOHIDWater/MOHIDWater.vfproj b/Solutions/VisualStudio2017_IntelFortran19/MOHIDNumerics/MOHIDWater/MOHIDWater.vfproj index b36bff3a1..b0a9332f3 100644 --- a/Solutions/VisualStudio2017_IntelFortran19/MOHIDNumerics/MOHIDWater/MOHIDWater.vfproj +++ b/Solutions/VisualStudio2017_IntelFortran19/MOHIDNumerics/MOHIDWater/MOHIDWater.vfproj @@ -35,8 +35,8 @@ - - + + @@ -96,7 +96,7 @@ - + @@ -143,14 +143,14 @@ - + - +